曹創(chuàng)華， 鄧 專， 康方平,， 柳建新， 童孝忠

(1. 湖南省地質調查院, 長沙 410116;2. 中南大學 地球科學與信息物理學院, 長沙 410083;3. 中南大學 有色資源與地質災害探查湖南省重點實驗室, 長沙 410083)

0 引言

激發(fā)極化法在有色金屬礦床勘探領域發(fā)揮了重要的作用，自使用以來發(fā)現(xiàn)了大批硫化物礦床[1-3]。特別是在大于萬分之一圖幅礦產調查研究階段，因效率高、可以旁側等特點，激發(fā)極化法中間梯度法裝置在面積性工作時發(fā)揮了很大的作用[4-9]。但在實際應用之中，目前中間梯度法測線布設方位問題鮮有文獻資料進行敘述參考文獻[10]也僅僅直接給出結論：測線方位需垂直研究區(qū)主要目標體方向，但在實際工作中不可能完全垂直地質異常體。究其原因，筆者在對承擔的項目實際探測資料分析以及查閱其他技術人員的實測資料時發(fā)現(xiàn)：礦區(qū)地質目標體在地層中走向往往不一，或者由于地形地貌原因，使得在布設測線時不同的技術人員測線方位不盡一致，隨之其在地表探測出的數據不盡相同，最后導致地質解釋不同而存在了不同程度的偏差。從這一現(xiàn)象出發(fā)，利用有限單元法進行了激發(fā)極化法中間梯度裝置條件下的三維正演，接著分析了數值計算的有效性，最后根據實際的常見地質條件和施工參數設計了模型進行了激電中間梯度法剖面和異常走向不同夾角的模擬計算，得出了具有實際意義的結論，為勘查地球物理技術工作者提供了借鑒。

如圖1所示，激電法勘探在實際布設電極時的供電源A極(或者B極)，在數據采集區(qū)域的中間2/3區(qū)間可近似均勻電流場的區(qū)域進行測量[10]，若中間區(qū)域地層中有異常體存在，則均勻場發(fā)生變化，在野外通過測量M極和N極之間電性參數來研究地質體電阻率變化。理論上的研究方法具體有：①物理模擬；②數值模擬。物理模擬是通過實驗室物理實驗模擬真實物理過程的方法，將實際地形物理的縮小模型置于實驗體(如風洞、水槽等)內，在滿足基本相似條件(包括幾何、運動、熱力、動力和邊界條件相似)的基礎上，模擬真實過程的主要特征[11]；而數值模擬是依靠電子計算機，結合有限元、有限差分和有限容積的概念，通過數值計算和圖像顯示的方法，達到對工程問題和物理問題乃至自然界各類問題的研究[12]。近年來，隨著電子計算機技術的發(fā)展，因室內數值模擬成本低，速度快而受推崇，數值模擬目前利用最多的是有限單元法(FEM)和有限差分法(FDM)法，有限差分法以計算方法簡單速度快為特點，但對異常體的邊界模擬精度有限，而有限單元法對異常體和測量剖面的邊界擬合精度較高，其中王家林教授[13]課題組利用積分方程法模擬了低電阻率差情況下偶極-偶極剖面三維激電法正演，加快了計算速度；強建科博士[21]利用有限單元法數值模擬了單個剖面上偶極-偶極、單極-偶極、對稱四極等裝置條件下地質異常體響應規(guī)律；劉海飛等[14]實現(xiàn)了起伏地形電導率連續(xù)變化條件下利用單極-偶極裝置的響應三維模擬，較好地反映了異常體的規(guī)模和空間形態(tài)。隨著計算機性能的提高，譚捍東教授[15-16]課題組以CUDA為基礎，利用MPI并行技術對電法勘探進行了三維數值模擬，大大提高了計算速度；戴前偉等[17]利用PSO-BP 非線性反演方法實現(xiàn)了二維Wenner-Schlumberger裝置的規(guī)則異常體成像效果，證實了其方法垂直方向上的反演分辨率更高。在三維響應和成像規(guī)律的研究方面，崔益安等[18]利用二極裝置研究了電法勘探三維正反演技術，把地球物理技術應用范圍拓展到了污染監(jiān)測領域，發(fā)現(xiàn)其具有良好的實用性和時效性。

圖1 激電法二維主剖面測量示意圖Fig.1 The schematic diagram of two-dimensional principal profile measurement by IP method

我們發(fā)現(xiàn)目前在有色金屬找礦技術領域，最為常用的中間梯度法裝置掃面形式的數據研究甚有學者問津，加之筆者對以往工作中遇到的迷惑，因此利用有限單元法，針對供電電源為三維場源、計算了激電中間梯度法地表面積性測量條件下進行了三維數值，以期從本質上對激電掃面中間梯度法測線布設有一個客觀的認知。首先參考黃俊革、強建科等[19-22]已有成果推導了三維場源的變分公式，然后設計層狀模型計算了解析解和數值結果進行對比；最后設計了符合實際工作地電參數的模型，計算了激電中間梯度裝置條件下二度地質異常體的走向和測量剖面方向之間不同夾角面上響應，分析并討論了理論計算結果。

1 正演響應計算原理

正演響應是利用設計的地電模型，根據電場的特征和邊界條件進行室內數值計算，得出其在特定裝置下在地表的響應(圖2)。

采用有限單元法，經過單元剖分、插值，單元積分、單元集成等來計算電場值，根據裝置參數等地電觀測數據計算激電中間梯度法條件下的面積性視電阻率和視極化率平面分布。

圖2 二度體正演模型三維模型示意圖Fig.2 The model diagram two-dimensional body forward model in three-dimensional

1.1 變分方程

為了更好地對一定角度的異常體模擬，且網格能夠均勻的分布在測量區(qū)域，利用三角剖分進行插值，通過構造變分進行推導，得到三維電場的邊值問題與下列變分問題等價[21]：

δF(u)=0

(1)

式中:I為電流強度；j表示電流密度矢量；U為總電位；σ表示介質的電導率；Ω為計算內部區(qū)域閉合空間；Γ計算區(qū)域內部的電阻率界面；r為點電源和計算點之間的距離；A為地表計算的點位置。

1.2 有限單元法

利用有限單元法對式(1)進行計算，對圖3(a)中的三維地質體進行剖分，在測量區(qū)域剖分為三角柱狀體。

(2)

將kij(其中i,j=1,2,…,6)集合呈矩陣形式，可得式(3)。

(3)

圖3 有限元剖分模型示意圖Fig.3 Typical finite-element mesh(a)三維地質體剖分示意圖;(b)單個單元剖分示意圖

依據上面系數組成單元剛度矩陣如下：

(4)

則式(1)中的第二項可變換為式(5)。

(5)

式中積分只與電源點有關。第三項積分與邊界有關，若單元e的一個面1254落在無窮遠邊界上時，其系數矩陣為式(6)～式(7)。

(6)

(7)

同理，當單元e的面1463落在無窮遠邊界時，系數矩陣為式(8)。

當單元e的面123落在向下無窮遠邊界時，系數矩陣為式(9)。

(8)

(9)

F(u)= ∑Fe(u)=

(10)

KU=P

(11)

式中:K為系數矩陣；U為電位向量；P為與點源有關的向量，解方程組得各節(jié)點的電位值U。然后通過計算響應場值，然后按照式(12)～式(14)計算激電中間梯度法視電阻率和視極化率值[23]。

(12)

(13)

即可計算MN測量點之間的視電阻率，其中AM、AN、BM、BN、MN為極距之間的距離，ΔUMN為M和N極之間的電位差(圖1)。為了衡量礦化物的直接激發(fā)極化大小特征，根據異常體響應條件下的二次場和總電壓比值來計算其大小[24]。

(14)

式中:ΔU2為MN之間的供電電源斷開后的二次場電位差。

2 正演算法正確性驗證

設計了三層H型地層結構的一維介質，具體參數為：第一層電阻率50 Ω·m，厚度5 m；第二層電阻率10 Ω·m，厚度5 m；第三層電阻率200 Ω·m。正演響應的解析解用一維層狀模型遞推公式進行迭代計算，數值解按照模型設計利用本文中的計算方法計算，其計算結果見表1。由表1可以看出，在18個極距中，最大誤差為3.26%，說明精度相對較高，數值計算結果較為可信。

表1 層狀模型數值解和解析解計算結果

3 三維正演模擬

按照高斯坐標系度量，沿x軸、y軸和z軸方向分布的地質體為1 200 m×800 m×400 m，設計的二度體模型及相關參數見圖2。模擬空間坐標為x=[0, 1200]m，y=[-200, 600]m，z=[0, -400]m。異常體以[600,200,-60] m為中心，其大小為200 m×30 m×30 m的低電阻率和高極化率異常體，電阻率圍巖為4 000 Ω·m、異常體為10 Ω·m；極化率圍巖1%、異常體10%。

圖4 中間梯度法測線方位和異常體夾角變化下掃面正演結果Fig.4 Forward result under the change of the angle between the azimuth line and the anomaly body in the middle gradient method(a)夾角為0°視電阻率掃面模擬結果;(b)夾角為0°視極化率掃面模擬結果;(c)夾角為30°視電阻率掃面模擬結果;(d)夾角為30°視極化率掃面模擬結果;(e)夾角為60°視電阻率掃面模擬結果;(f)夾角為60°視極化率掃面模擬結果;(g)夾角為90°視電阻率掃面模擬結果;(h)夾角為90°視極化率掃面模擬結果

研究了測線MN和異常體走向在地表夾角0°(測線方位和異常體走向一致)、夾角30°、夾角60°和相互垂直時的激電掃面結果(圖4)，其中點距為5 m，線距為10 m，保證跨越異常體的為主至少有三個響應測點。圖4中，隨著異常體和測線夾角的增加，掃面結果和異常體的地表投影逐漸趨于一致，且視電阻率對測線方位不同的測量結果敏感程度大于視極化率。在夾角為0°時，視電阻率被分為兩個對稱的相對低阻體；當夾角為30°時，兩個對稱的低阻異常峰值減小且與兩個峰值之間的電阻率差異越來越??；當夾角增加到60°時，兩個對稱的低阻異常峰值幾乎未能發(fā)現(xiàn)，但異常體的邊界與地表下的異常體投影有一定的偏差；當夾角增加到90°時，異常體的地面投影和掃面結果完全對應，說明夾角過小時測量的異常是不可信的。雖然視極化率沒有視電阻率明顯，但是隨著角度的增加異常的幅值逐漸增大，異常得到強化，說明在實際野外測線進行布置時需要盡量的垂直主要構造或者礦(化)體。

為了量化測線和異常地質體走向對對異常體測量的影響，特抽取了主剖面數據(圖4中的藍色實線)，研究了隨著異常的走向的變化得到的結果(圖5)。當異常體和測線平行時，其面上視電阻率不僅出現(xiàn)了兩個峰值且位置與中心測點不對稱，其視極化率也出現(xiàn)了負值，隨著兩者夾角的增加，異常逐漸得到正確的顯示。明顯的當兩者夾角大于60°時對異常的定性影響較小。

圖5 中間梯度法測線主剖面不同角度測量數據正演結果Fig.5 Forward grades of survey data at different angles in main profile of intermediate gradient method(a)測量主軸視電阻率掃面模擬結果;(b)測量主軸視極化率掃面模擬結果

4 結論

通過以上分析，我們利用有限單元法正確模擬了激電法三維模擬，特別是從實際勘查角度出發(fā)，研究了激電中間梯度法測量時其剖面和異常體地質走向呈不同夾角時的響應規(guī)律，主要結論如下：

1)實現(xiàn)了激電中間梯度法條件下的三維數值模擬，設計了模型計算了地表掃面數據的響應。

2)通過計算發(fā)現(xiàn)，不同正演響應數據其對測線布設方位敏感不一，通過對視電阻率和視極化率的正演響應結果對比來看，隨著兩者夾角的減小視參數數值差異越小，對后期的正確解釋困難也隨之增加。

3)激電中間梯度法的測線方向應盡可能的垂直異常體走向，如果研究區(qū)域異常沿其走向存在連續(xù)性和對稱性，可考慮變換方位進行測量以研究地質異常體是否為一個異常?？梢酝ㄟ^改變測線方向，用異常相交的方法定位異常源的頂板中心位置。

4)一個研究區(qū)一般情況要求布置相同的測線方向，不同時期的物探資料如果測線方位不同，即不能簡單的進行拼接和比較，以免在綜合研究分析時產生誤導。