潘北斗

(1.中國冶金地質(zhì)總局 第一地質(zhì)勘查院,燕郊 065201；2.中國冶金地質(zhì)總局 礦產(chǎn)資源研究院，北京 101300)

0 引言

內(nèi)蒙古東烏旗地區(qū)主要有斑巖型、侵染型礦床的鉬礦，熱液脈狀、熱液蝕變型礦床的銀鉛鋅礦，矽卡巖型礦床的鐵鋅礦等[1-2]。找礦標(biāo)志幾乎都是金屬硫化物，而金屬硫化物較之于圍巖有較強(qiáng)的激電效應(yīng)[3-4]。2001年至今第一地質(zhì)勘查院在東烏旗地區(qū)開展了大量的激電工作。在阿爾哈達(dá)、迪彥欽阿木和花腦特等礦區(qū)取得了較好的找礦成果,而在其他礦區(qū)沒有突破性找礦成果，其主要原因之一是物探工作的多樣性、反演解釋的多解性和反演精度的不穩(wěn)定性造成的[5-7]。依據(jù)本次研究工作目標(biāo)和工作內(nèi)容，綜合研究區(qū)的具體情況，確定本次研究工作技術(shù)路線。研究利用相關(guān)資料、通過軟件進(jìn)行數(shù)值分析，推導(dǎo)出激電異常相對(duì)于目標(biāo)體真實(shí)空間位置的偏移量(Δa)與目標(biāo)體埋深(h0)、測(cè)量極距(AN)之間的關(guān)系，對(duì)模型加以修改，并找出影響模型的相關(guān)因素定量分析，通過實(shí)際運(yùn)用對(duì)其可行性和實(shí)用效果加以客觀評(píng)價(jià)，結(jié)合本區(qū)情況提出下一步找礦方向。項(xiàng)目研究的總體技術(shù)路線見圖1。

1 激發(fā)極化法反演問題

地球物理反演是通過地球物理異常的分布特征確定地質(zhì)體的賦存狀態(tài)，即在數(shù)學(xué)基礎(chǔ)上，通過實(shí)際觀測(cè)值使其與模型理論值達(dá)到最大限度的擬合，以此來求解地球物理模型。

圖1 技術(shù)路線圖Fig.1 Technology roadmap

1.1 廣義線性反演

通常把地球物理反演問題用一個(gè)泛涵方程組表示，其離散化后取得一些非線性方程組，但是地球物理數(shù)據(jù)是有限的，地球物理的解通常具有多解性，此時(shí)只能使用某種可接受的解估計(jì)，該方法稱為廣義反演方法，利用泰勒級(jí)數(shù)展開的形式形成的方程組反演方法都稱為廣義線性反演[8]。

1.2 時(shí)間域激發(fā)極化原理

時(shí)間域激電法作為電法勘探的重要分支，廣泛應(yīng)用在資源勘查中，通過人工場(chǎng)向地下以脈沖直流供電，并以一種裝置的形式對(duì)地下電流的分布規(guī)律和激發(fā)效應(yīng)進(jìn)行觀測(cè)研究，通過反演得出異常體的埋深、產(chǎn)狀和位置等，從而查明礦產(chǎn)資源賦存狀態(tài)和地下地質(zhì)體情況[9-10]。

1.3 激發(fā)極化法的數(shù)值模擬

隨著電子計(jì)算技術(shù)的飛速發(fā)展，為電場(chǎng)的數(shù)值模擬提供了技術(shù)支撐，使之得以實(shí)現(xiàn)。目前，國內(nèi)、外已經(jīng)有一些學(xué)者和單位，廣泛地將數(shù)值模擬技術(shù)應(yīng)用于研究激電異常，用來對(duì)穩(wěn)定電流場(chǎng)作數(shù)值模擬的方法主要采用有限單元法和有限差分法。其中有限差分法一般適用于小規(guī)模模擬計(jì)算，而有限單元法通常用于大型模型分析。

在激電法模擬中，利用有限單元法求解地電斷面為真電阻率和有效電阻率分布時(shí)的穩(wěn)定電流場(chǎng)，再依據(jù)等效電阻率法計(jì)算激電異常[11]。

1.4 最小二乘約束反演方法

地球物理反演是非線性的，病態(tài)的，反演解釋存在多解性以及非唯一性。目前比較成功的做法是將光滑約束、先驗(yàn)信息等加入反演算法中，建立基于某種約束的最小二乘反演算法。

非線性問題線性化，并引入光滑約束和已知先驗(yàn)信息，構(gòu)造出最小二乘反演目標(biāo)函數(shù)為式(1)[12-13]。

ψ=‖Wd(δd-Aδm)‖2+

‖Wm(m-m0+δm)‖2

(1)

式中：‖Wd(δd-Aδm)‖2為常規(guī)的最小二乘法;‖Wm(m-m0+δm)‖2為先驗(yàn)信息項(xiàng)。其中，δd為數(shù)據(jù)殘差矢量，其值為實(shí)測(cè)數(shù)據(jù)對(duì)數(shù)值與模型正演計(jì)算數(shù)據(jù)對(duì)數(shù)值之差；m為預(yù)測(cè)模型向量，其值為模型參數(shù)的對(duì)數(shù)值；m0為基本模型向量，其值為模型參數(shù)的對(duì)數(shù)值；A為偏導(dǎo)數(shù)矩陣；Wd為觀測(cè)數(shù)據(jù)加權(quán)矩陣。Wm為稱模型加權(quán)矩陣。目標(biāo)函數(shù)式對(duì)δm求導(dǎo)，并令其等于零,解線性方程組便可得模型參數(shù)帶約束條件時(shí)的最小二乘反演結(jié)果。

2 精細(xì)化模型反演研究

由于激電測(cè)量受體積效應(yīng)的影響，激電異常往往是各種異常的疊加反應(yīng)，隨著極距的增大，探測(cè)成果中疊加的上部異常信息越多，反演解釋難度隨之增大。精細(xì)化反演處理的目的是通過反演將深部的疊加異常信息分離開來，從而獲得異常源的有效信息。精細(xì)化反演的特征是力求同時(shí)使觀測(cè)數(shù)據(jù)與計(jì)算數(shù)據(jù)間的差值、背景與反演模擬參數(shù)間的差值以及反演模擬粗糙度最小化，通常利用均方根誤差準(zhǔn)則量度數(shù)值的擬合性，偏離先前模擬的距離以及模擬粗糙度。設(shè)計(jì)如圖2和圖3所示的簡(jiǎn)單模型來進(jìn)行激電測(cè)深精細(xì)化反演模擬。

圖2 電阻率模型Fig.2 Resistivity model

圖3 充電率模型Fig.3 Charging rate model

2.1 有限元網(wǎng)格的選擇

2.1.1 初始模型層數(shù)

一般來講，當(dāng)保持網(wǎng)格內(nèi)單元格大小不變情況下，網(wǎng)格越大越好；反之保持網(wǎng)格不變情況下，單元格越小越好，這兩種情況下都會(huì)增加大量的計(jì)算工作量。為了解決精度和工作量之間的矛盾，往往采用非均勻網(wǎng)格，即網(wǎng)格中心單元小、節(jié)點(diǎn)密，邊界單元大、節(jié)點(diǎn)稀，由中心到邊緣單元逐漸放大。優(yōu)點(diǎn)是在保證了網(wǎng)格有足夠的大小前提下，又保證地電斷面的復(fù)雜部位位于網(wǎng)格中心，從而滿足求解要求。反演實(shí)踐中通常會(huì)給定一組初始模型作為擬合初值，假定為層狀地電結(jié)構(gòu)，從表層向深部層厚依次增加，增加比例在1.01～1.5之間，首層厚度一般設(shè)置為0.3電極距，保證模擬深度不小于實(shí)際探測(cè)深度。

反演采用了10層(圖4)、20層(圖5)、30層(圖6)模型進(jìn)行擬合。從反演結(jié)果看，10層模型反演結(jié)果比較粗糙，尤其是下部控制稀疏，已經(jīng)形成拖拽現(xiàn)象，造成異常放大，顯然不足以擬合觀測(cè)數(shù)據(jù)，20層與30層模型反演結(jié)果區(qū)別不大，由此可以推定，20層模型已經(jīng)完全可以擬合觀測(cè)數(shù)據(jù)。

圖4 10層模型的反演結(jié)果Fig.4 Inversion results of 10-layer model(a)電阻率模型;(b)充電率模型

圖5 20層模型的反演結(jié)果Fig.5 Inversion results of 20-layer model(a)電阻率模型;(b)充電率模型

圖6 30層模型的反演結(jié)果Fig.6 Inversion results of 30-layer model(a)電阻率模型;(b)充電率模型

2.1.2 網(wǎng)格密度

網(wǎng)格密度增大會(huì)提高計(jì)算精度，但會(huì)增加計(jì)算規(guī)模，所以在工作中要綜合這兩個(gè)因素考慮。網(wǎng)格較少時(shí)，在計(jì)算時(shí)間增加不大的情況下，通過增加網(wǎng)格數(shù)量提高計(jì)算精度提高。當(dāng)網(wǎng)格密度增加到一定程度后，再繼續(xù)增加時(shí)精度提高很小，網(wǎng)格密度應(yīng)增加到計(jì)算結(jié)果在誤差允許范圍以內(nèi)便可。在網(wǎng)格密度上分別采用了等倍電極距(510個(gè)模塊)圖7，0.5倍電極距(1 020個(gè)模塊)圖8和0.25倍電極距(2 040個(gè)模塊)圖9進(jìn)行反演，從反演結(jié)果看，等倍電極距的橫向分辨率不夠，造成了異常橫向拉長(zhǎng)，低充電率區(qū)域被掩蓋，0.5倍電極距和0.25倍電極距在異常顯示上并沒有較大形變，但是0.25倍電極距斷面圖引入了許多虛假細(xì)節(jié)，反演結(jié)果不夠光滑，這是過擬合地反映。

圖7 等電極距劃分510個(gè)模塊的反演Fig.7 Inversion results of 510 modules divided by equidistant electrodes(a)電阻率模型;(b)充電率模型

圖8 0.5電極距劃分1 020個(gè)模塊的反演Fig.8 Inversion results of 1 020 modules divided by 0.5 electrode spacing(a)電阻率模型;(b)充電率模型

圖9 0.25電極距劃分2 040個(gè)模塊的反演結(jié)果Fig.9 Inversion results of 2 040 modules divided by 0.25 electrode spacing(a)電阻率模型;(b)充電率模型

2.2 平滑度約束

平滑度約束是反演模型的重要參數(shù)，擬合觀測(cè)數(shù)據(jù)和維持平滑度變化是反演過程當(dāng)中同時(shí)需要考慮的。如果一個(gè)新模型非常平滑，但其計(jì)算數(shù)據(jù)不能與野外數(shù)據(jù)很好擬合，則這個(gè)反演是失敗的，然而一味追求數(shù)據(jù)擬合度則會(huì)使得模型斷面失去趨勢(shì)性，過多的虛假異常往往會(huì)掩蓋了目標(biāo)異常。這時(shí)候可試用一定范圍的模擬約束權(quán)值進(jìn)行反演，改變數(shù)據(jù)擬合與維持平滑模擬間的折中，直到用較低的模擬約束權(quán)產(chǎn)生對(duì)觀測(cè)數(shù)據(jù)的完美擬合，要避免過高約束形成的粗糙模型和過低約束形成的過擬合現(xiàn)象。

反演采用了0.1(圖10)、1(圖11)和10(圖12)作為平滑約束因子。從反演結(jié)果來看:0.1的反演結(jié)果幾乎沒有平滑約束，數(shù)據(jù)有很大的靈活性，結(jié)果形成類一維地電結(jié)構(gòu)，二維約束缺失；平滑因子為10的結(jié)果過多地模糊了細(xì)節(jié)異常，兩個(gè)水平的異常體不能分辨，顯然對(duì)于精細(xì)化反演來講是背道而馳的，平滑因子為1的結(jié)果則顯示了較好的平滑折中，達(dá)到了精細(xì)反演的目的。

圖10 無平滑約束的反演結(jié)果(平滑因子0.1)Fig.10 Inversion results without smoothing constraints (smoothing factor 0.1)(a)電阻率模型;(b)充電率模型

圖11 適度平滑約束的反演結(jié)果(平滑因子1)Fig.11 Inversion results of moderate smoothing constraints (smoothing factor 1)(a)電阻率模型;(b)充電率模型

圖12 過度平滑約束的反演結(jié)果(平滑因子10)Fig.12 Inversion results of excessive smoothing constraints (smoothing factor 10)(a)電阻率模型;(b)充電率模型

圖13 5次迭代的反演結(jié)果Fig.13 The inversion results of 5 iterations(a)電阻率模型;(b)充電率模型

圖14 20次迭代的反演結(jié)果Fig.14 Inversion results of 20 iterations(a)電阻率模型;(b)充電率模型

2.3 收斂標(biāo)準(zhǔn)

收斂標(biāo)準(zhǔn)是影響反演結(jié)果的主要因素，主要指標(biāo)是最小殘差，通常以百分比來衡量。如果殘差大于10%即可認(rèn)為反演的可信度較低，百分比通常是越小越好。但是這將需要較長(zhǎng)的反演時(shí)間，最好地反演是能在最小殘差和觀測(cè)數(shù)據(jù)擬合約束間找到平衡。這需要在反演迭代次數(shù)，最小步長(zhǎng)和最小殘差上的多次試驗(yàn)，以找到合理的擬合模型，從而改善反演結(jié)果。

圖15 梯度約束控制圖Fig.15 Gradient constrained control chart

反演實(shí)踐中迭代3次以上即可以得到擬合模型的大致輪廓，隨著迭代次數(shù)的增加，反演精度有明顯提高。圖13、圖14中使用1%殘差和1.5%步長(zhǎng)分別迭代5次與20次的結(jié)果，圖13所示5次迭代的模型雖然顯示了主要異常，但是分辨率很低，模型粗糙，迭代20次的結(jié)果較完美的還原了原始模型，可見多次迭代可以有效改善迭代效果。

2.4 梯度約束

在某些已知地質(zhì)條件的情況下，初步反演之后，在明顯高對(duì)比度接觸帶的模型斷面區(qū)域，引入梯度約束來切斷平滑度約束是有益的。尤其在取得了鉆孔資料后，有必要在鉆孔部位約束反演，這可以大大減少反演的多解性，在一些導(dǎo)電覆蓋區(qū)則必須進(jìn)行梯度約束來抵制低阻屏蔽的影響。初步反演能給出淺部信息的良好指示，但是平滑度約束與電法體積效應(yīng)的作用可能模糊深部特征。因此在數(shù)據(jù)靈敏度降低的深部切斷平滑度約束，然后重新反演可改善深部特征的分辨力，這可使得基巖或者接觸帶類型的激電異常變得更加明顯。

鑒于低阻的拖拽效應(yīng)，筆者在低阻下方和高低阻接觸帶部位編輯梯度約束。從圖16可以看出，低阻和高低阻接觸部位得到了很好的約束，準(zhǔn)確地反映了原始模型的結(jié)構(gòu)，反演效果良好。

圖16 帶梯度約束的反演結(jié)果Fig.16 Inversion results with gradient constraints(a)電阻率模型;(b)充電率模型

圖17 69線反演結(jié)果Fig.17 Inversion results of 69 lin(a)電阻率;(b)充電率

3 精細(xì)化反演

3.1 迪彥欽阿木激電測(cè)深精細(xì)化反演

根據(jù)已知資料，選取較有代表性的69線激電測(cè)深工作來研究。由鉆孔資料得知，ZK6904和ZK6903鉆孔見礦位置圍巖主要為凝灰?guī)r和安山巖等。由圖17可以看出，69線地電結(jié)構(gòu)比較簡(jiǎn)單，淺部激電等值線比較平緩，視電阻率偏低，與第四系埋深較厚有關(guān)；深部特征比較明顯，其中視電阻率在剖面兩側(cè)較高，中間相對(duì)低。充電率在形態(tài)上呈現(xiàn)“桃”狀形態(tài)，測(cè)線大號(hào)方向衰減較快，異常體中間稍有凹陷，與已控制的礦體比較，反演充電率和礦(化)體吻合很好，礦體細(xì)節(jié)形態(tài)也有顯現(xiàn)，精細(xì)化反演效果良好。

從圖18可以看出， 21線地電結(jié)構(gòu)比較復(fù)雜，第四系埋深較淺，深部橫向上高低阻相間，形態(tài)復(fù)雜，充電率向大小號(hào)方向放射，底部有凹陷和已控制礦體對(duì)比，總體形態(tài)吻合，但是已控制礦體形態(tài)異常復(fù)雜，礦體多呈細(xì)脈狀，反演斷面對(duì)細(xì)節(jié)的顯示有限，說明精細(xì)化反演對(duì)于復(fù)雜斷面的反演仍然有待提高。

圖18 21線反演結(jié)果Fig.18 Inversion results of 21 line(a)電阻率;(b)充電率

3.2 哈巴特蓋激電測(cè)深精細(xì)化反演

2014年哈巴特蓋工區(qū)共完成激電測(cè)深剖面6條，其中在19線進(jìn)行了工程驗(yàn)證。在19線點(diǎn)號(hào)920處施工鉆孔，孔深為400.70 m，全孔共分4層，地層上部主要為安格爾音烏拉組角巖化砂巖，9.20 m～54.85 m為氧化帶，巖心幾處破碎較強(qiáng)，具褐鐵礦化，54.85 m～302.90 m為砂巖原生帶，多處云英巖脈穿插，厚度為0.3 m～8 m不等。與原巖接觸面中軸夾角為30°～50°。下部302.90 m～400.70 m為花崗巖體，中粗粒至細(xì)粒不等。

圖19 19線反演結(jié)果Fig.19 Inversion results of 19 line

在已有鉆孔資料的基礎(chǔ)上，對(duì)原有數(shù)據(jù)進(jìn)行了精細(xì)化重新反演，從圖20可以看出，中部的低阻凹陷有了較大的變化，且充電率形態(tài)有了細(xì)節(jié)反映，新斷面也與已知鉆孔吻合，可見精細(xì)化反演對(duì)本條測(cè)線有較大地改善。

圖20 19線精細(xì)化反演結(jié)果Fig.20 Refined inversion results of 19 line

4 結(jié)果與討論

4.1 精細(xì)化反演選擇

1)網(wǎng)格選擇。在劃分網(wǎng)格時(shí)一般采用非均勻網(wǎng)格，網(wǎng)格的中心部分單元小，節(jié)點(diǎn)密，邊界單元大，節(jié)點(diǎn)稀，由中心到邊緣單元逐漸放大，初始模型層數(shù)建議一般不少于16層，過少的層數(shù)會(huì)使擬合模型非常粗糙，然而過多的層數(shù)對(duì)精細(xì)化反演改善并不大，且反演消耗將會(huì)指數(shù)增加，甚至產(chǎn)生假異常，橫向上最少要保持半極距分辨力以上的網(wǎng)格數(shù)量才能保證擬合的準(zhǔn)確度。

2)平滑度選擇。平滑過高會(huì)導(dǎo)致計(jì)算數(shù)據(jù)不能與野外數(shù)據(jù)很好擬合，最終形成簡(jiǎn)單粗糙的模型，較低的平滑度會(huì)提高數(shù)據(jù)擬合度，但同時(shí)引入噪聲與假異常，通常需要在這兩者之間選擇折中。

3)收斂標(biāo)準(zhǔn)。收斂標(biāo)準(zhǔn)以數(shù)據(jù)擬合誤差作為評(píng)價(jià)，在平滑度保證的前提下擬合誤差越小越好，但是誤差設(shè)置過小可能會(huì)造成模型不收斂，在設(shè)置誤差門檻的過程中迭代次數(shù)和搜索步長(zhǎng)是兩個(gè)關(guān)鍵參數(shù)，實(shí)踐證明搜索步長(zhǎng)維持在1%，迭代次數(shù)在5次以上，大多數(shù)情況下模型擬合殘差會(huì)小于5%，而且迭代次數(shù)在20次以上會(huì)對(duì)反演結(jié)果有一定改善，更多的迭代會(huì)消耗大量資源，且改善效果不明顯。

4)梯度約束。在已知地質(zhì)條件的情況下，或者斷面內(nèi)有鉆孔資料后，引入梯度約束的效果是明顯的。另外一種情況是在一些導(dǎo)電覆蓋區(qū)進(jìn)行梯度約束來抵制低阻屏蔽的影響。在初步反演之后，在明顯梯度帶部位編輯梯度約束重新反演，能更好地修正反演結(jié)果，可使接觸帶類型的激電異常變得更加明顯。

4.2 關(guān)于過擬合

1)進(jìn)行精細(xì)化反演之前，物探數(shù)據(jù)的預(yù)處理顯得尤為重要，如果反演數(shù)據(jù)存在較多噪音干擾，這將使精細(xì)化工作變得異常困難，更有可能反演出假異常，數(shù)據(jù)的前期濾波就變得重要了。

2)目前的反演軟件通常以不等邊四邊形作為網(wǎng)格劃分，在進(jìn)行反演時(shí)充分考慮數(shù)據(jù)靈敏度和分辨率寬度，可以提供合理的網(wǎng)格劃分思路，一味追求細(xì)分網(wǎng)格不僅會(huì)使反演開銷成倍增長(zhǎng)，而且會(huì)使得噪聲和假異常泛濫，這種過擬合現(xiàn)象對(duì)反演弊大于利。

3)平滑度對(duì)反演的控制異常重要，追求數(shù)據(jù)擬合度往往形成過擬合反演模型，從而淹沒主要異常，使反演結(jié)果變得難以分辨。

5 結(jié)論

筆者分別從有限元網(wǎng)格、平滑度約束、收斂標(biāo)準(zhǔn)、梯度約束等幾方面研究，對(duì)精細(xì)化反演做了深入探索，主要得出以下結(jié)論。

1)網(wǎng)格劃分時(shí)一般采用非均勻網(wǎng)格，即網(wǎng)格的中心部分單元小，節(jié)點(diǎn)密，邊界單元大，節(jié)點(diǎn)稀，由中心到邊緣單元逐漸放大；平滑度不宜過高也不宜太低，通常需要在兩者之間選擇折中；在平滑度保證的前提下擬合誤差越小越好，但誤差設(shè)置不能過小，以免造成模型不收斂。

2)在已知地質(zhì)條件的情況下，引入梯度約束能夠明顯提高反演精度。

3)復(fù)雜的地電結(jié)構(gòu)仍然是困擾反演精度的主要因素，后續(xù)研究將針對(duì)這些問題進(jìn)行改進(jìn)。