基于Matlab平臺的地-井TEM快速成像

楊 毅， 張 杰， 霍美凈， 鄧曉紅，王興春， 智慶全， 武軍杰

(1.中國地質(zhì)科學院 地球物理地球化學勘查研究所，廊坊 065000;2．國土資源部 地球物理電磁法探測技術重點實驗室，廊坊 065000;3．成都理工大學，成都 610059;4.廊坊市水利勘察規(guī)劃設計院，廊坊 065000)

0 引言

地-井TEM法是瞬變電磁法的一個重要分支方法，興起于上世紀70年代，最早由蘇聯(lián)、加拿大等國開展相關研究。Crone于上世紀70年代研制了第一套井中瞬變電磁測量系統(tǒng)，標志著這項方法技術應用于實際勘察工作。多年來經(jīng)過國內(nèi)外學者的努力，不論是儀器系統(tǒng)[1-3]還是方法技術[4-12]都得到了極大發(fā)展，使得地-井TEM法廣泛應用于國內(nèi)外的硫化物型銅鎳礦的勘察[13-23]，并取得巨大成功。

1 理論方法

1.1 等效渦流理論

地-井TEM法一般采用地面發(fā)射，井中接收的工作方式，能以“一孔之見”對井旁和井底可能存在的異常做出準確預判，為進一步的地質(zhì)鉆孔設計提供重要參考。井中瞬變電磁法最早是通過研究自由空間中局部異常體(薄板體、球體等)的特征來開展的[24-25]，其場的特征研究均以感應渦流為基礎，其中又以導電薄板(圖1)最具代表性，導電薄板內(nèi)感應渦流建立與消失的物理過程，為等效電流環(huán)反演奠定理論基礎。

圖1 導電薄板內(nèi)渦旋電流分布示意圖Fig.1 Schematic diagram of eddy current distribution in conductive thin plate

圖2 圓形載流環(huán)參數(shù)模型Fig.2 Parametric model of circular current

位于一次場中的導電薄板，在發(fā)射突然關斷瞬間，根據(jù)法拉第定律，為了維持關斷之前一次場的作用，在導電薄板體內(nèi)會產(chǎn)生感應渦流以維持關斷之前一次場的作用，感應渦流在板內(nèi)產(chǎn)生與薄板體形態(tài)相近的一系列電流強度和大小均不同的電流環(huán)，對應于瞬變過程早期的電流環(huán)主要集中在板體的邊緣，尺寸相對大，這些電流環(huán)隨著時間推移逐漸向?qū)w中心擴散，直至以熱損耗形式消亡。在這過程中導電薄板體內(nèi)的任意時刻的感應電流分布可以用一個等效電流環(huán)表示，如圖2所示。這樣通過等效電流環(huán)即實現(xiàn)了對局部導電薄板的描述，位于自由空間的電流環(huán)可由七個參數(shù)來表示，即電流環(huán)中心坐標(x，y，z)、電流環(huán)半徑、電流強度、電流環(huán)傾角、電流環(huán)方位角。

等效載流環(huán)包含七個參數(shù)(圖2)，其在自由空間中任意點磁場的數(shù)學表達式為：

(1)

(2)

依據(jù)式 ( 1 ) 和式 ( 2 ) 即可計算自由空間中任意電流強度、大小、展布的電流環(huán)在任一點產(chǎn)生的場。

1.2 異常場提取

在異常場成像之初需要對實測數(shù)據(jù)進行異常提取，采用層狀大地響應來近似等效背景圍巖的響應，從實測數(shù)據(jù)中減去這個背景場，從而得到異常場響應。計算上采用Kerry Key(2009)[26]的方法，從電偶極子出發(fā)推導其在層狀大地空間中任意點的電磁響應，再通過有限個長度的偶極子疊加獲得回線源層狀大地瞬變電磁響應(李建平，2005)[27]。

水平電偶極子位于層狀大地時，除z方向的矢量位之外，還出現(xiàn)x方向的分量，其柱坐標下矢量勢分量分別寫為式(3)與式(4)。

(3)

(4)

其相應的磁場各分量表達式可通過式(5)～式(7)求得。

(5)

(6)

(7)

通過以上步驟求得的磁場響應為包含貝塞爾函數(shù)的復雜積分，求解上采用Hankel數(shù)字濾波計算，其一般表達式可寫為式(8)。

(8)

完成頻率域場的計算再采用余弦變換將場值換算到時間域，以Hz為例其數(shù)值離散近似表達式寫為式(9)。

(9)

依據(jù)以上公式，使用Fortran編寫層狀大地正演計算代碼，為人機交互功能控件提供回調(diào)函數(shù)。

2 人機交互快速成像實現(xiàn)

Matlab圖形用戶界面 ( GUI ) 是指由窗口、菜單、圖標、光標、按鍵、對話框和文本等各種圖形對象組成的用戶界面。它讓用戶定制用戶與Matlab的交互方式。用戶通過鼠標或鍵盤選擇、激活這些圖形對象，使計算機產(chǎn)生某種動作或變化?；緢D形對象分為控件對象和用戶界面菜單對象，簡稱控件和菜單。

Matlab提供了一套可視化的創(chuàng)建圖形窗口的工具，使用圖形用戶界面開發(fā)環(huán)境可方便地創(chuàng)建GUI應用程序，它可以根據(jù)用戶設計的GUI布局，自動生成M文件的框架，用戶使用這一框架編制自己的應用程序。Matlab提供了一套可視化的創(chuàng)建圖形用戶接口 ( GUI ) 的工具，層次架構如圖3。

圖3 GUI功能及程序運行流程簡圖Fig.3 GUI function and program flow diagram

2.1 人機交互界面設計和功能實現(xiàn)

交互界面的布局應該在保證基本功能的前提下盡量簡潔美觀，秉承這一理念，根據(jù)需要將快速成像程序分成三個功能模塊來實現(xiàn)①數(shù)據(jù)載入編輯；②人機交互背景場擬合；③人機交互異常場擬合。

其功能簡圖及其主代碼名稱如下圖3所示，運行主程序PureAnomalyFitting.m文件，點擊菜單Data Raw Process調(diào)用inversion.m程序，進入界面，在此界面完成三分量地-井TEM數(shù)據(jù)的讀入和編輯工作，完成后關閉子GUI界面將數(shù)據(jù)傳遞給主GUI。主GUI界面左側分區(qū)可進行地層參數(shù)、發(fā)射框角點坐標的輸入、發(fā)射電流輸入、下降沿輸入以及場值計算類型的選擇，輸入完畢后，繪制地層厚度和電阻率折線圖，程序中的交互設計允許在此圖上對層厚和電阻率值進行拖動，快速實現(xiàn)層厚和電阻率參數(shù)值修改，再進行層狀大地響應正演擬合，最后保存背景場和異常場。主GUI界面右側分區(qū)可以進行電流環(huán)個數(shù)以及七參數(shù)的快速載入，點擊電流環(huán)正演，進行異常場擬合，同樣在另一個圖形中顯示電流環(huán)空間位置展布圖，在此設計了交互編輯功能，允許用戶在不同視角下對電流環(huán)的參數(shù)進行快速編輯，然后根據(jù)異常特征進行電流環(huán)交互擬合，人工試錯不斷修正，直至達到最佳擬合，最后將異常場不同道擬合最好的電流環(huán)在三維坐標里成圖即可得到異常的三維空間展布。

2.2 數(shù)據(jù)載入編輯

程序開始部分首先應該對實測數(shù)據(jù)進行載入并以剖面曲線的形式展現(xiàn)出來，為了美觀和數(shù)據(jù)操作編輯方便，將其作為一個子GUI單獨呈現(xiàn)，其作用為為初始成像做好基本數(shù)據(jù)處理，包含數(shù)據(jù)載入、數(shù)據(jù)編輯、曲線圓滑等三個大的功能。代碼主M文件為inversion.m，對應的GUI文件名為inversion.fig。

GUI運行后界面如圖4所示，設計了三個axes控件，分別用以顯示地-井TEM測量的x、y、z三個分量數(shù)據(jù)，在Data Raw Process下拉菜單PLOT選項可對需要顯示的道數(shù)進行選擇。菜單的安排上，設置了Data Load、Data Plot、Data Smooth、Data Delete四個功能，分別用以完成數(shù)據(jù)載入、成圖、數(shù)據(jù)圓滑、數(shù)據(jù)點刪除等功能，其中在plot菜單下調(diào)用selectchannel子GUI用于選擇需要操作的時間道。

圖4 子GUI運行效果示意圖Fig.4 Sketch map of running effect of sub GUI

圖5 主GUI控件設計示意圖Fig.5 Schematic diagram of the main GUI control

2.3 背景場擬合和異常場提取

M代碼主文件名PureAnomalyFitting.m文件，對應的GUI文件文PureAnomalyFitting.fig，其布設如圖5所示，按功能分區(qū)，左側為背景場擬合部分，使用的控件有一個Button Group，內(nèi)設兩個互斥的Radio Button，分別代表背景場計算場值為電動勢還是磁場；Layer earth H&R負責層狀大地層厚和層電阻率的錄入，具體方式為在界面上手動輸入；TXI ( A )和RAMP(ms)分別對應發(fā)射電流和下降沿；Loop corner coordinates負責發(fā)射框角點坐標的錄入，可在界面上手動輸入；Axes6負責顯示錄入地層的厚度和電阻率(折線圖)；Axes1、Axes2、Axes3分別用于擬合曲線顯示，axes4用于擬合差的顯示，這幾個axes控件為公共控件，在異常場擬合部分也用于擬合曲線和擬合差的顯示。pure anomaly按鈕用于擬合曲線和擬合差曲線的繪制。

對層狀大地層厚和電阻率的交互編輯功能為本段程序的重點實現(xiàn)部分，在初始給定層數(shù)、層厚度和電阻率之后，程序允許在Axes6界面中對頂?shù)捉缑?、電阻率折線圖進行拖動，實現(xiàn)層狀大地層厚和電阻率參數(shù)的快速修改。

在菜單中設置integral菜單，其功能為在完成層狀大地正演后，從實測數(shù)據(jù)中減去層狀大地響應，獲取異常場響應，并對異常場積分，獲取磁場響應，為異常場的電流環(huán)成像做好準備。

2.4 異常場擬合

在完成背景場擬合和異常場提取之后即進行此步驟，在GUI布局設計上，如圖5所示，從右至左，filament choose為下拉菜單控件，用于電流環(huán)參數(shù)設置后各參數(shù)的載入；下拉菜單之下為七組兩兩橫排的控件，用于電流環(huán)七參數(shù)的輸入，左側為edit text可編輯文本框控件，右側為slide滑塊控件，滑塊和文本框輸入設置為聯(lián)動，即拖動滑塊或輸入文本框，另一項值也隨著同步變化，可以快速實現(xiàn)參數(shù)的修改；Edit Panel Select 控件是一個互斥的radio button組，用于交互操作時，交互區(qū)域的選擇，分別對應背景場擬合時地層層厚、電阻率折線圖和電流環(huán)兩個分區(qū)的選擇和鎖定。Axes5用于發(fā)射框、井、一次場以及給定電流環(huán)的空間位置三維顯示；filament plot按鈕控件用于初次電流環(huán)的繪制以及參數(shù)重置；refresh按鈕用于交互成像電流環(huán)參數(shù)改變后對電流環(huán)的刷新顯示。

對電流環(huán)的交互編輯功能為本段程序的重點實現(xiàn)部分，在初始給定電流環(huán)個數(shù)和基本位置之后，程序允許在x-y，x-z，y-x三個視角界面下對電流環(huán)位置進行拖動，實現(xiàn)參數(shù)的快速修改。

2.5 GUI主程序運行流程及其功能實現(xiàn)

本程序的編制基于Matlab的GUI平臺，根據(jù)需要設置了控件，并逐一完善控件功能。下面將簡單敘述筆者編寫的人機交互程序的工作流程，實現(xiàn)功能的控件、菜單以及回調(diào)函數(shù)。

打開Matlab程序，將路徑切換到目標文件夾，雙擊PureAnomalyFitting.m函數(shù)，點擊運行按鈕主程序開始運行，彈出如下界面(圖6)。

1 ) 完成數(shù)據(jù)的讀入與編輯。單擊Data Raw Process菜單，彈出inversion.fig界面(圖7)，在此界面下載入數(shù)據(jù)，并可對數(shù)據(jù)進行簡單編輯，其調(diào)用的函數(shù)文件為loaddata1.m，主要完成實測數(shù)據(jù)文件的整理和讀入；Data Raw Process點擊后調(diào)用inversion.fig負責數(shù)據(jù)的初步編輯處理，其調(diào)用的函數(shù)文件為loaddata1.m，其主要功能為實測數(shù)據(jù)文件的整理和讀入。

2 ) 完成發(fā)射參數(shù)和層狀大地參數(shù)給定。從Button Group中選擇層狀大地正演計算(背景場)所需場類型，df為感應電動勢，ff為感應磁場；在TXI和RAMP兩個可編輯輸入文本控件中分別輸入發(fā)射電流和下降沿；在layer earth H&R輸入文本框中輸入正演電阻率模型參數(shù)；在Loop corner coordinates可編輯文本框里輸入發(fā)射框角點坐標；完成之后點擊Layer earth model plot菜單，即可在GUI左下角的axes6中繪制出電阻率-層厚折線圖；點擊Primary field plot可在axes5中繪制一次場矢量圖，發(fā)射框和井的位置圖，主要調(diào)用ZY.m完成一次矢量場的空間分布計算(圖8)。

3 )通過人機交互方式完成背景場計算。根據(jù)實測曲線特征，已經(jīng)完成了層狀大地正演所需參數(shù)準備，點擊Layer earth forward按鈕完成層狀大地背景場正演計算，此處主要調(diào)用temFwdLayerBoreholeFast.exe執(zhí)行文件；在完成正演后，按下pure anomaly按鈕繪制實測數(shù)據(jù)與層狀大地正演擬合曲線(圖9)；如果擬合不滿意，則可以點GUI擊右下角Edit panel select下的layer earth edit單選按鈕，啟動GUI左下角axes6上的電阻率-層厚折線圖編輯模式，在axes6范圍內(nèi)左鍵選中線段(圖9)，線段被選中后即可進行快速拖動以改變電阻率或?qū)雍穸戎?水平方向線段可上下拖動，改變地層厚度，垂直方向線段可左右拖動，改變層電阻率值)，然后再點擊 Layer earth forward按鈕進行正演，直到擬合滿意為止。

4 ) 完成異常場轉換。點擊integral菜單，用“差值法”完成異常場提取并將電動勢積分轉換為磁場響應，調(diào)用函數(shù)jifen11.m。

圖6 GUI主程序運行界面Fig.6 GUI main program running interface

圖7 副GUI程序運行界面Fig.7 Vice GUI program running interface

圖8 主GUI程序分區(qū)參數(shù)輸入示意圖Fig.8 Schematic diagram of main GUI program partition parameter input

圖9 主GUI程序背景擬合示及地層參數(shù)編輯意圖Fig.9 Background fitting of main GUI program and editing intention of formation parameters

圖10 異常場擬合與電流環(huán)編輯示意圖Fig.10 Abnormal field fitting and current loop edit diagram

5 ) 完成電流環(huán)初始參數(shù)值的給定。點擊filamentNumber菜單，在彈出對話框中輸入?yún)⑴c異常場擬合的電流環(huán)個數(shù)，點擊后載入到下拉菜單filament choose中，然后在其正下方表征電流環(huán)七參數(shù)的七個可編輯輸入文本框中給出每個電流環(huán)的參數(shù)，注意每完成一個電流環(huán)參數(shù)的輸入，需要在下拉菜單中選中對應的電流環(huán)名稱已完成參數(shù)的載入，完成后點擊Filament plot按鈕將電流環(huán)繪制于axes5中完成顯示。

6 ) 通過人機交互方式完成電流環(huán)擬合。點擊GUI右下角的Filament Edit radio控件選中電流環(huán)編輯模式，則可在x-y，x-z，y-z三個平面視角內(nèi)對電流環(huán)進行拖動操作，完成其空間位置坐標的迅速修改，三個平面視角的實現(xiàn)由三個activity 控件控制，對電流環(huán)大小、強度以及兩個空間角度的修改則可以在圖中點擊選中(圖10)相應的電流環(huán)，然后在電流環(huán)輸入文本里編輯或拖動其后的滑塊實現(xiàn)，完成后點擊GUI下方的refresh按鈕完成數(shù)據(jù)更新以及電流環(huán)的重繪。

完成電流環(huán)設置后，可點擊pure plot按鈕，在axes1、axes2、axes3、axes4中分別顯示各分量異常場與電流環(huán)響應擬合情況(圖10)，如不滿意則可返回上一步采用人機交互方式完成電流環(huán)參數(shù)的快速時時修正。

7 ) 數(shù)據(jù)的輸出與顯示。在完成各道異常場曲線與對應電流環(huán)的擬合后，點擊spatial按鈕，將不同道擬合得到的電流環(huán)會繪制于一個圖上，最后用800 dpi分辨率將圖從四個視角以*.jpeg格式輸出。至此，完成了整個交互成像程序的編制。

3 理論算例

使用Maxwell軟件中設計的理論模型計算結果作為假設實測數(shù)據(jù)進行交互成像。

模型均勻半空間中放置一個薄板體，其中均勻半空間電阻率為300 Ω·m，中心框發(fā)射，框邊長為200 m，采樣道為36道，間隔采用Crone 50 ms間隔，下降沿為0.5 ms，發(fā)射電流10 A，計算單位為pT/s；導電板體大小設置為50 m×50 m，縱向電導為500 S，導電薄板中心坐標為(-50,-50，-100)，傾角為25°，旋轉角為30°。

圖11 異常提取與電流環(huán)成像擬合曲線Fig.11 Fitting curve of anomaly extraction and current loop imaging

采用上一節(jié)所述方法對Maxwell設計的半空間含導電薄板模型計算的結果進行交互成像，其中背景場擬合(圖11(a))在早期相對差，在中晚期擬合較好，自編程序計算的半空間響應能夠完全反映背景場信息，通過“差值法”從實測數(shù)據(jù)中減去自編程序模擬計算的背景場即得到異常場，并采用電流環(huán)響應與之擬合(圖11(b))，對于三個分量擬合結果在早期道相對差，自編程序在給定與Maxwell背景值一致條件下計算得到的場值相對強，原因可能為未考慮背景與導電薄板的耦合所致，中晚期擬合好(圖11(b)、圖11(c)、圖11(d))，所得電流環(huán)與薄板大小相近，形態(tài)一致(圖12)。表明了程序的正確性，為下一步實測數(shù)據(jù)的交互成像奠定了基礎。

4 結論

筆者實現(xiàn)了基于等效電流環(huán)理論的Matlab 人機交互快速成像，編寫了人機交互程序，主要取得如下結論：

1)編寫了基于Matlab GUI平臺的等效電流環(huán)快速成像界面，利用GUI自帶的動作控件函數(shù)實現(xiàn)了層狀大地模型和電流環(huán)模型的快速修改編輯，為異常場提取和成像奠定了基礎。

2)完成了理論數(shù)據(jù)的快速成像，所得電流環(huán)與理論導電薄板模型一致，驗證了程序的正確性，為地-井TEM實測數(shù)據(jù)的定量解釋開辟了新思路。

圖12 等效電流環(huán)空間分布圖Fig.12 Spatial distribution diagram of equivalent current