激電中梯裝置測量長方異常體的三維正演模擬及啟示

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

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

0 引言

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

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

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

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

1 正演響應(yīng)計算原理

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

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

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

1.1 變分方程

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

δF(u)=0

(1)

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

1.2 有限單元法

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

(2)

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

(3)

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

依據(jù)上面系數(shù)組成單元剛度矩陣如下：

(4)

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

(5)

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

(6)

(7)

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

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

(8)

(9)

F(u)= ∑Fe(u)=

(10)

KU=P

(11)

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

(12)

(13)

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

(14)

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

2 正演算法正確性驗證

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

表1 層狀模型數(shù)值解和解析解計算結(jié)果

3 三維正演模擬

按照高斯坐標(biāo)系度量，沿x軸、y軸和z軸方向分布的地質(zhì)體為1 200 m×800 m×400 m，設(shè)計的二度體模型及相關(guān)參數(shù)見圖2。模擬空間坐標(biāo)為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 中間梯度法測線方位和異常體夾角變化下掃面正演結(jié)果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°視電阻率掃面模擬結(jié)果;(b)夾角為0°視極化率掃面模擬結(jié)果;(c)夾角為30°視電阻率掃面模擬結(jié)果;(d)夾角為30°視極化率掃面模擬結(jié)果;(e)夾角為60°視電阻率掃面模擬結(jié)果;(f)夾角為60°視極化率掃面模擬結(jié)果;(g)夾角為90°視電阻率掃面模擬結(jié)果;(h)夾角為90°視極化率掃面模擬結(jié)果

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

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

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

4 結(jié)論

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

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

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

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

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