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

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

0 引言

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

1 理論方法

1.1 等效渦流理論

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

圖1 導(dǎo)電薄板內(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

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

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

(1)

(2)

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

1.2 異常場提取

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

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

(3)

(4)

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

(5)

(6)

(7)

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

(8)

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

(9)

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

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

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

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

圖3 GUI功能及程序運(yùn)行流程簡圖Fig.3 GUI function and program flow diagram

2.1 人機(jī)交互界面設(shè)計(jì)和功能實(shí)現(xiàn)

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

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

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

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

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

圖4 子GUI運(yùn)行效果示意圖Fig.4 Sketch map of running effect of sub GUI

圖5 主GUI控件設(shè)計(jì)示意圖Fig.5 Schematic diagram of the main GUI control

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

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

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

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

2.4 異常場擬合

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

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

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

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

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

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

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

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

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

圖6 GUI主程序運(yùn)行界面Fig.6 GUI main program running interface

圖7 副GUI程序運(yùn)行界面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ù)值的給定。點(diǎn)擊filamentNumber菜單，在彈出對(duì)話框中輸入?yún)⑴c異常場擬合的電流環(huán)個(gè)數(shù)，點(diǎn)擊后載入到下拉菜單filament choose中，然后在其正下方表征電流環(huán)七參數(shù)的七個(gè)可編輯輸入文本框中給出每個(gè)電流環(huán)的參數(shù)，注意每完成一個(gè)電流環(huán)參數(shù)的輸入，需要在下拉菜單中選中對(duì)應(yīng)的電流環(huán)名稱已完成參數(shù)的載入，完成后點(diǎn)擊Filament plot按鈕將電流環(huán)繪制于axes5中完成顯示。

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

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

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

3 理論算例

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

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

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

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

4 結(jié)論

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

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

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

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