潘北斗
(1.中國冶金地質(zhì)總局 第一地質(zhì)勘查院,燕郊 065201;2.中國冶金地質(zhì)總局 礦產(chǎn)資源研究院,北京 101300)
內(nèi)蒙古東烏旗地區(qū)主要有斑巖型、侵染型礦床的鉬礦,熱液脈狀、熱液蝕變型礦床的銀鉛鋅礦,矽卡巖型礦床的鐵鋅礦等[1-2]。找礦標志幾乎都是金屬硫化物,而金屬硫化物較之于圍巖有較強的激電效應(yīng)[3-4]。2001年至今第一地質(zhì)勘查院在東烏旗地區(qū)開展了大量的激電工作。在阿爾哈達、迪彥欽阿木和花腦特等礦區(qū)取得了較好的找礦成果,而在其他礦區(qū)沒有突破性找礦成果,其主要原因之一是物探工作的多樣性、反演解釋的多解性和反演精度的不穩(wěn)定性造成的[5-7]。依據(jù)本次研究工作目標和工作內(nèi)容,綜合研究區(qū)的具體情況,確定本次研究工作技術(shù)路線。研究利用相關(guān)資料、通過軟件進行數(shù)值分析,推導(dǎo)出激電異常相對于目標體真實空間位置的偏移量(Δa)與目標體埋深(h0)、測量極距(AN)之間的關(guān)系,對模型加以修改,并找出影響模型的相關(guān)因素定量分析,通過實際運用對其可行性和實用效果加以客觀評價,結(jié)合本區(qū)情況提出下一步找礦方向。項目研究的總體技術(shù)路線見圖1。
地球物理反演是通過地球物理異常的分布特征確定地質(zhì)體的賦存狀態(tài),即在數(shù)學(xué)基礎(chǔ)上,通過實際觀測值使其與模型理論值達到最大限度的擬合,以此來求解地球物理模型。
圖1 技術(shù)路線圖Fig.1 Technology roadmap
通常把地球物理反演問題用一個泛涵方程組表示,其離散化后取得一些非線性方程組,但是地球物理數(shù)據(jù)是有限的,地球物理的解通常具有多解性,此時只能使用某種可接受的解估計,該方法稱為廣義反演方法,利用泰勒級數(shù)展開的形式形成的方程組反演方法都稱為廣義線性反演[8]。
時間域激電法作為電法勘探的重要分支,廣泛應(yīng)用在資源勘查中,通過人工場向地下以脈沖直流供電,并以一種裝置的形式對地下電流的分布規(guī)律和激發(fā)效應(yīng)進行觀測研究,通過反演得出異常體的埋深、產(chǎn)狀和位置等,從而查明礦產(chǎn)資源賦存狀態(tài)和地下地質(zhì)體情況[9-10]。
隨著電子計算技術(shù)的飛速發(fā)展,為電場的數(shù)值模擬提供了技術(shù)支撐,使之得以實現(xiàn)。目前,國內(nèi)、外已經(jīng)有一些學(xué)者和單位,廣泛地將數(shù)值模擬技術(shù)應(yīng)用于研究激電異常,用來對穩(wěn)定電流場作數(shù)值模擬的方法主要采用有限單元法和有限差分法。其中有限差分法一般適用于小規(guī)模模擬計算,而有限單元法通常用于大型模型分析。
在激電法模擬中,利用有限單元法求解地電斷面為真電阻率和有效電阻率分布時的穩(wěn)定電流場,再依據(jù)等效電阻率法計算激電異常[11]。
地球物理反演是非線性的,病態(tài)的,反演解釋存在多解性以及非唯一性。目前比較成功的做法是將光滑約束、先驗信息等加入反演算法中,建立基于某種約束的最小二乘反演算法。
非線性問題線性化,并引入光滑約束和已知先驗信息,構(gòu)造出最小二乘反演目標函數(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為先驗信息項。其中,δd為數(shù)據(jù)殘差矢量,其值為實測數(shù)據(jù)對數(shù)值與模型正演計算數(shù)據(jù)對數(shù)值之差;m為預(yù)測模型向量,其值為模型參數(shù)的對數(shù)值;m0為基本模型向量,其值為模型參數(shù)的對數(shù)值;A為偏導(dǎo)數(shù)矩陣;Wd為觀測數(shù)據(jù)加權(quán)矩陣。Wm為稱模型加權(quán)矩陣。目標函數(shù)式對δm求導(dǎo),并令其等于零,解線性方程組便可得模型參數(shù)帶約束條件時的最小二乘反演結(jié)果。
由于激電測量受體積效應(yīng)的影響,激電異常往往是各種異常的疊加反應(yīng),隨著極距的增大,探測成果中疊加的上部異常信息越多,反演解釋難度隨之增大。精細化反演處理的目的是通過反演將深部的疊加異常信息分離開來,從而獲得異常源的有效信息。精細化反演的特征是力求同時使觀測數(shù)據(jù)與計算數(shù)據(jù)間的差值、背景與反演模擬參數(shù)間的差值以及反演模擬粗糙度最小化,通常利用均方根誤差準則量度數(shù)值的擬合性,偏離先前模擬的距離以及模擬粗糙度。設(shè)計如圖2和圖3所示的簡單模型來進行激電測深精細化反演模擬。
圖2 電阻率模型Fig.2 Resistivity model
圖3 充電率模型Fig.3 Charging rate model
2.1.1 初始模型層數(shù)
一般來講,當保持網(wǎng)格內(nèi)單元格大小不變情況下,網(wǎng)格越大越好;反之保持網(wǎng)格不變情況下,單元格越小越好,這兩種情況下都會增加大量的計算工作量。為了解決精度和工作量之間的矛盾,往往采用非均勻網(wǎng)格,即網(wǎng)格中心單元小、節(jié)點密,邊界單元大、節(jié)點稀,由中心到邊緣單元逐漸放大。優(yōu)點是在保證了網(wǎng)格有足夠的大小前提下,又保證地電斷面的復(fù)雜部位位于網(wǎng)格中心,從而滿足求解要求。反演實踐中通常會給定一組初始模型作為擬合初值,假定為層狀地電結(jié)構(gòu),從表層向深部層厚依次增加,增加比例在1.01~1.5之間,首層厚度一般設(shè)置為0.3電極距,保證模擬深度不小于實際探測深度。
反演采用了10層(圖4)、20層(圖5)、30層(圖6)模型進行擬合。從反演結(jié)果看,10層模型反演結(jié)果比較粗糙,尤其是下部控制稀疏,已經(jīng)形成拖拽現(xiàn)象,造成異常放大,顯然不足以擬合觀測數(shù)據(jù),20層與30層模型反演結(jié)果區(qū)別不大,由此可以推定,20層模型已經(jīng)完全可以擬合觀測數(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)格密度增大會提高計算精度,但會增加計算規(guī)模,所以在工作中要綜合這兩個因素考慮。網(wǎng)格較少時,在計算時間增加不大的情況下,通過增加網(wǎng)格數(shù)量提高計算精度提高。當網(wǎng)格密度增加到一定程度后,再繼續(xù)增加時精度提高很小,網(wǎng)格密度應(yīng)增加到計算結(jié)果在誤差允許范圍以內(nèi)便可。在網(wǎng)格密度上分別采用了等倍電極距(510個模塊)圖7,0.5倍電極距(1 020個模塊)圖8和0.25倍電極距(2 040個模塊)圖9進行反演,從反演結(jié)果看,等倍電極距的橫向分辨率不夠,造成了異常橫向拉長,低充電率區(qū)域被掩蓋,0.5倍電極距和0.25倍電極距在異常顯示上并沒有較大形變,但是0.25倍電極距斷面圖引入了許多虛假細節(jié),反演結(jié)果不夠光滑,這是過擬合地反映。
圖7 等電極距劃分510個模塊的反演Fig.7 Inversion results of 510 modules divided by equidistant electrodes(a)電阻率模型;(b)充電率模型
圖8 0.5電極距劃分1 020個模塊的反演Fig.8 Inversion results of 1 020 modules divided by 0.5 electrode spacing(a)電阻率模型;(b)充電率模型
圖9 0.25電極距劃分2 040個模塊的反演結(jié)果Fig.9 Inversion results of 2 040 modules divided by 0.25 electrode spacing(a)電阻率模型;(b)充電率模型
平滑度約束是反演模型的重要參數(shù),擬合觀測數(shù)據(jù)和維持平滑度變化是反演過程當中同時需要考慮的。如果一個新模型非常平滑,但其計算數(shù)據(jù)不能與野外數(shù)據(jù)很好擬合,則這個反演是失敗的,然而一味追求數(shù)據(jù)擬合度則會使得模型斷面失去趨勢性,過多的虛假異常往往會掩蓋了目標異常。這時候可試用一定范圍的模擬約束權(quán)值進行反演,改變數(shù)據(jù)擬合與維持平滑模擬間的折中,直到用較低的模擬約束權(quán)產(chǎn)生對觀測數(shù)據(jù)的完美擬合,要避免過高約束形成的粗糙模型和過低約束形成的過擬合現(xiàn)象。
反演采用了0.1(圖10)、1(圖11)和10(圖12)作為平滑約束因子。從反演結(jié)果來看:0.1的反演結(jié)果幾乎沒有平滑約束,數(shù)據(jù)有很大的靈活性,結(jié)果形成類一維地電結(jié)構(gòu),二維約束缺失;平滑因子為10的結(jié)果過多地模糊了細節(jié)異常,兩個水平的異常體不能分辨,顯然對于精細化反演來講是背道而馳的,平滑因子為1的結(jié)果則顯示了較好的平滑折中,達到了精細反演的目的。
圖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)充電率模型
收斂標準是影響反演結(jié)果的主要因素,主要指標是最小殘差,通常以百分比來衡量。如果殘差大于10%即可認為反演的可信度較低,百分比通常是越小越好。但是這將需要較長的反演時間,最好地反演是能在最小殘差和觀測數(shù)據(jù)擬合約束間找到平衡。這需要在反演迭代次數(shù),最小步長和最小殘差上的多次試驗,以找到合理的擬合模型,從而改善反演結(jié)果。
圖15 梯度約束控制圖Fig.15 Gradient constrained control chart
反演實踐中迭代3次以上即可以得到擬合模型的大致輪廓,隨著迭代次數(shù)的增加,反演精度有明顯提高。圖13、圖14中使用1%殘差和1.5%步長分別迭代5次與20次的結(jié)果,圖13所示5次迭代的模型雖然顯示了主要異常,但是分辨率很低,模型粗糙,迭代20次的結(jié)果較完美的還原了原始模型,可見多次迭代可以有效改善迭代效果。
在某些已知地質(zhì)條件的情況下,初步反演之后,在明顯高對比度接觸帶的模型斷面區(qū)域,引入梯度約束來切斷平滑度約束是有益的。尤其在取得了鉆孔資料后,有必要在鉆孔部位約束反演,這可以大大減少反演的多解性,在一些導(dǎo)電覆蓋區(qū)則必須進行梯度約束來抵制低阻屏蔽的影響。初步反演能給出淺部信息的良好指示,但是平滑度約束與電法體積效應(yīng)的作用可能模糊深部特征。因此在數(shù)據(jù)靈敏度降低的深部切斷平滑度約束,然后重新反演可改善深部特征的分辨力,這可使得基巖或者接觸帶類型的激電異常變得更加明顯。
鑒于低阻的拖拽效應(yīng),筆者在低阻下方和高低阻接觸帶部位編輯梯度約束。從圖16可以看出,低阻和高低阻接觸部位得到了很好的約束,準確地反映了原始模型的結(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)充電率
根據(jù)已知資料,選取較有代表性的69線激電測深工作來研究。由鉆孔資料得知,ZK6904和ZK6903鉆孔見礦位置圍巖主要為凝灰?guī)r和安山巖等。由圖17可以看出,69線地電結(jié)構(gòu)比較簡單,淺部激電等值線比較平緩,視電阻率偏低,與第四系埋深較厚有關(guān);深部特征比較明顯,其中視電阻率在剖面兩側(cè)較高,中間相對低。充電率在形態(tài)上呈現(xiàn)“桃”狀形態(tài),測線大號方向衰減較快,異常體中間稍有凹陷,與已控制的礦體比較,反演充電率和礦(化)體吻合很好,礦體細節(jié)形態(tài)也有顯現(xiàn),精細化反演效果良好。
從圖18可以看出, 21線地電結(jié)構(gòu)比較復(fù)雜,第四系埋深較淺,深部橫向上高低阻相間,形態(tài)復(fù)雜,充電率向大小號方向放射,底部有凹陷和已控制礦體對比,總體形態(tài)吻合,但是已控制礦體形態(tài)異常復(fù)雜,礦體多呈細脈狀,反演斷面對細節(jié)的顯示有限,說明精細化反演對于復(fù)雜斷面的反演仍然有待提高。
圖18 21線反演結(jié)果Fig.18 Inversion results of 21 line(a)電阻率;(b)充電率
2014年哈巴特蓋工區(qū)共完成激電測深剖面6條,其中在19線進行了工程驗證。在19線點號920處施工鉆孔,孔深為400.70 m,全孔共分4層,地層上部主要為安格爾音烏拉組角巖化砂巖,9.20 m~54.85 m為氧化帶,巖心幾處破碎較強,具褐鐵礦化,54.85 m~302.90 m為砂巖原生帶,多處云英巖脈穿插,厚度為0.3 m~8 m不等。與原巖接觸面中軸夾角為30°~50°。下部302.90 m~400.70 m為花崗巖體,中粗粒至細粒不等。
圖19 19線反演結(jié)果Fig.19 Inversion results of 19 line
在已有鉆孔資料的基礎(chǔ)上,對原有數(shù)據(jù)進行了精細化重新反演,從圖20可以看出,中部的低阻凹陷有了較大的變化,且充電率形態(tài)有了細節(jié)反映,新斷面也與已知鉆孔吻合,可見精細化反演對本條測線有較大地改善。
圖20 19線精細化反演結(jié)果Fig.20 Refined inversion results of 19 line
1)網(wǎng)格選擇。在劃分網(wǎng)格時一般采用非均勻網(wǎng)格,網(wǎng)格的中心部分單元小,節(jié)點密,邊界單元大,節(jié)點稀,由中心到邊緣單元逐漸放大,初始模型層數(shù)建議一般不少于16層,過少的層數(shù)會使擬合模型非常粗糙,然而過多的層數(shù)對精細化反演改善并不大,且反演消耗將會指數(shù)增加,甚至產(chǎn)生假異常,橫向上最少要保持半極距分辨力以上的網(wǎng)格數(shù)量才能保證擬合的準確度。
2)平滑度選擇。平滑過高會導(dǎo)致計算數(shù)據(jù)不能與野外數(shù)據(jù)很好擬合,最終形成簡單粗糙的模型,較低的平滑度會提高數(shù)據(jù)擬合度,但同時引入噪聲與假異常,通常需要在這兩者之間選擇折中。
3)收斂標準。收斂標準以數(shù)據(jù)擬合誤差作為評價,在平滑度保證的前提下擬合誤差越小越好,但是誤差設(shè)置過小可能會造成模型不收斂,在設(shè)置誤差門檻的過程中迭代次數(shù)和搜索步長是兩個關(guān)鍵參數(shù),實踐證明搜索步長維持在1%,迭代次數(shù)在5次以上,大多數(shù)情況下模型擬合殘差會小于5%,而且迭代次數(shù)在20次以上會對反演結(jié)果有一定改善,更多的迭代會消耗大量資源,且改善效果不明顯。
4)梯度約束。在已知地質(zhì)條件的情況下,或者斷面內(nèi)有鉆孔資料后,引入梯度約束的效果是明顯的。另外一種情況是在一些導(dǎo)電覆蓋區(qū)進行梯度約束來抵制低阻屏蔽的影響。在初步反演之后,在明顯梯度帶部位編輯梯度約束重新反演,能更好地修正反演結(jié)果,可使接觸帶類型的激電異常變得更加明顯。
1)進行精細化反演之前,物探數(shù)據(jù)的預(yù)處理顯得尤為重要,如果反演數(shù)據(jù)存在較多噪音干擾,這將使精細化工作變得異常困難,更有可能反演出假異常,數(shù)據(jù)的前期濾波就變得重要了。
2)目前的反演軟件通常以不等邊四邊形作為網(wǎng)格劃分,在進行反演時充分考慮數(shù)據(jù)靈敏度和分辨率寬度,可以提供合理的網(wǎng)格劃分思路,一味追求細分網(wǎng)格不僅會使反演開銷成倍增長,而且會使得噪聲和假異常泛濫,這種過擬合現(xiàn)象對反演弊大于利。
3)平滑度對反演的控制異常重要,追求數(shù)據(jù)擬合度往往形成過擬合反演模型,從而淹沒主要異常,使反演結(jié)果變得難以分辨。
筆者分別從有限元網(wǎng)格、平滑度約束、收斂標準、梯度約束等幾方面研究,對精細化反演做了深入探索,主要得出以下結(jié)論。
1)網(wǎng)格劃分時一般采用非均勻網(wǎng)格,即網(wǎng)格的中心部分單元小,節(jié)點密,邊界單元大,節(jié)點稀,由中心到邊緣單元逐漸放大;平滑度不宜過高也不宜太低,通常需要在兩者之間選擇折中;在平滑度保證的前提下擬合誤差越小越好,但誤差設(shè)置不能過小,以免造成模型不收斂。
2)在已知地質(zhì)條件的情況下,引入梯度約束能夠明顯提高反演精度。
3)復(fù)雜的地電結(jié)構(gòu)仍然是困擾反演精度的主要因素,后續(xù)研究將針對這些問題進行改進。