楊 磊

(中鐵第四勘察設(shè)計(jì)院集團(tuán)有限公司 地質(zhì)路基設(shè)計(jì)研究院，湖北 武漢 430063)

1 引 言

重力勘探資料是研究地質(zhì)構(gòu)造、圈定地下空洞異常區(qū)域和地下成礦帶的重要數(shù)據(jù)集。隨著空間探測(cè)技術(shù)的發(fā)展，衛(wèi)星重力、航空重力和微重力測(cè)量得到廣泛應(yīng)用[1,2]。野外實(shí)測(cè)重力異常通常是由地表到地下深部所有密度不均勻體引起的重力效應(yīng)疊加場(chǎng)，包含著不同深度異常體的豐富信息。重力資料解釋通常把實(shí)測(cè)重力異?？醋饔蓞^(qū)域異常和局部異常組成，區(qū)域異常主要由分布范圍較廣的、相對(duì)深的地質(zhì)因素引起，其異常幅值大、范圍寬、梯度小、體征出現(xiàn)低頻信息；局部異常通常由比區(qū)域地質(zhì)因素范圍小的研究對(duì)象引起，其異常幅度小、范圍小、梯度大、體征出現(xiàn)高頻信息[3]。重力異常包含了地下異常構(gòu)造和不同深度密度不均勻體的疊加信號(hào)，因此，位場(chǎng)分離研究已經(jīng)成為重力勘探解譯工作的重要環(huán)節(jié)[4-6]。

目前，重力異常分離的主要方法有解析延拓法、趨勢(shì)分析法、匹配濾波法、小波分析法等[4-6]。但這些方法在復(fù)雜的重力異常信號(hào)源中提取目標(biāo)異常信號(hào)時(shí)，濾波或信號(hào)分離層次的選取一直是較為棘手的問(wèn)題，處理效果較大地依賴于實(shí)際地質(zhì)資料。近些年來(lái)逐漸發(fā)展起來(lái)的反向傳播(Back-Propagation)神經(jīng)網(wǎng)絡(luò)算法，它具有數(shù)據(jù)映射高度非線性、數(shù)據(jù)處理運(yùn)行速度快和批量處理的并行性等特點(diǎn)[7-9]。其較好的學(xué)習(xí)和訓(xùn)練能力，在面對(duì)大量數(shù)據(jù)集合時(shí)，通過(guò)特定的學(xué)習(xí)訓(xùn)練能力計(jì)算出數(shù)據(jù)的特征分布，在實(shí)際工程勘察研究中得到了系列發(fā)展[10,11]。朱志大等[12]針對(duì)BP神經(jīng)網(wǎng)絡(luò)算法擬合重力異常問(wèn)題開(kāi)展了理論探究；聶建亮等[13]利用BP神經(jīng)網(wǎng)絡(luò)擬合區(qū)域似大地水準(zhǔn)面，結(jié)果表明其能較好地逼近實(shí)際似大地水準(zhǔn)面；郭文斌等[14]利用BP神經(jīng)網(wǎng)絡(luò)實(shí)現(xiàn)了重力的三維物性反演，證明其具有較好的容錯(cuò)和泛化能力；鄧才林等[15]根據(jù)實(shí)測(cè)GPS和水準(zhǔn)數(shù)據(jù)，構(gòu)建BP神經(jīng)網(wǎng)絡(luò)模型對(duì)縣域GPS高程異常進(jìn)行了擬合，證明了其擬合精度優(yōu)于曲面擬合法；王智偉等[16]利用BP神經(jīng)網(wǎng)絡(luò)數(shù)據(jù)融合算法對(duì)多源異構(gòu)監(jiān)測(cè)數(shù)據(jù)的滑坡變形成功進(jìn)行了預(yù)測(cè)，并有效提高了滑坡變形預(yù)測(cè)結(jié)果的精度。

常規(guī)的實(shí)測(cè)重力異常探測(cè)數(shù)據(jù)，區(qū)域重力場(chǎng)信號(hào)幅值相較于局部重力信號(hào)要大得多，并且局部異常的范圍分布較小。因而，利用BP神經(jīng)網(wǎng)絡(luò)進(jìn)行區(qū)域重力異常場(chǎng)擬合，進(jìn)而實(shí)現(xiàn)位場(chǎng)分離，理論合理，實(shí)施也較為便捷。本文的主要思路是通過(guò)比較BP神經(jīng)網(wǎng)絡(luò)方法和純集合的曲面擬合法的精度優(yōu)劣，將BP神經(jīng)網(wǎng)絡(luò)算法引用到實(shí)際重力數(shù)據(jù)資料處理中，為其在地球物理重力勘探資料解譯中位場(chǎng)分離工作提供實(shí)用性的參考。

2 方法原理

2.1 BP神經(jīng)網(wǎng)絡(luò)模型

反向傳播(Back Propagation)神經(jīng)網(wǎng)絡(luò)是一種利用網(wǎng)絡(luò)自適應(yīng)映射能力來(lái)進(jìn)行反向傳播的多層前饋網(wǎng)絡(luò)，可以實(shí)現(xiàn)從輸入到輸出的任意非線性運(yùn)算[17-20]。從網(wǎng)絡(luò)結(jié)構(gòu)分層上分析，網(wǎng)絡(luò)結(jié)構(gòu)一般由輸入層、輸出層和隱含層組成，隱含層可為多層復(fù)合層，各層內(nèi)只有相鄰的神經(jīng)元之間才產(chǎn)生反饋連接。BP神經(jīng)網(wǎng)絡(luò)的的本質(zhì)是求解目標(biāo)結(jié)果最優(yōu)化時(shí)所對(duì)應(yīng)的網(wǎng)絡(luò)權(quán)值，通過(guò)比較網(wǎng)絡(luò)輸出和期望輸出之間誤差，構(gòu)建相應(yīng)的誤差函數(shù)，利用最速下降法算法求解；不斷修正調(diào)整相應(yīng)的網(wǎng)絡(luò)權(quán)值，并重新返回到輸入層進(jìn)行相應(yīng)計(jì)算，直至網(wǎng)絡(luò)的整體誤差值減小到誤差限內(nèi)輸出結(jié)果。利用BP神經(jīng)網(wǎng)絡(luò)算法擬合區(qū)域重力信號(hào)的關(guān)鍵是通過(guò)選取合理比例尺的重力數(shù)據(jù)集作為訓(xùn)練集，建立起訓(xùn)練數(shù)據(jù)集中每一個(gè)的測(cè)點(diǎn)坐標(biāo)和該點(diǎn)重力異常數(shù)據(jù)的映射關(guān)系，從而可以訓(xùn)練輸出得到區(qū)域重力場(chǎng)數(shù)據(jù)集。由于反向傳播網(wǎng)絡(luò)的非線性映射和無(wú)模型估計(jì)的特點(diǎn)，對(duì)于擬合得到測(cè)網(wǎng)低頻區(qū)域異常信號(hào)具有一定的高精度。

BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)和實(shí)際擬合流程如圖1所示。BP神經(jīng)網(wǎng)絡(luò)可以對(duì)任意復(fù)雜函數(shù)進(jìn)行逼近，信號(hào)按照輸入層到輸出層的方向傳遞，而權(quán)值和偏置值的不斷修正是輸出層至輸入層的逆方向。

圖1 BP 神經(jīng)網(wǎng)絡(luò)示意圖

公式(1)表明了輸入層到輸出層之間的信號(hào)傳遞關(guān)系，其中，xi表示輸入層第i個(gè)節(jié)點(diǎn)的輸入值，隱含層第i個(gè)節(jié)點(diǎn)到到輸入層第j個(gè)節(jié)點(diǎn)之間的權(quán)值ωij，隱含層第i個(gè)節(jié)點(diǎn)的偏置值bi，輸入層神經(jīng)元個(gè)數(shù)n1，隱含層神經(jīng)元個(gè)數(shù)n2，隱含層的激勵(lì)函數(shù)f(x)，輸出層第k個(gè)節(jié)點(diǎn)到到輸入層第i個(gè)節(jié)點(diǎn)之間的權(quán)值ωki，yk則表示為第k個(gè)節(jié)點(diǎn)的輸出。

(1)

此外，所謂的誤差反向傳播，是通過(guò)輸出層來(lái)逐層計(jì)算各層神經(jīng)元的輸出誤差，并根據(jù)誤差梯度下降法來(lái)調(diào)節(jié)各層的權(quán)值和偏置值，在訓(xùn)練過(guò)程中不斷修改參數(shù)降低誤差。其總誤差目標(biāo)函數(shù)可表示為：

(2)

式(2)中，N和P分別表示為樣本個(gè)數(shù)和輸出層節(jié)點(diǎn)數(shù)，di為期望輸出，E為誤差目標(biāo)函數(shù)。

2.2 趨勢(shì)分析法基本原理

在重力勘探數(shù)據(jù)分析中，測(cè)點(diǎn)P(x,y)的區(qū)域異?？梢杂蓽y(cè)點(diǎn)坐標(biāo)x和y的n階多項(xiàng)式擬合。趨勢(shì)分析方法通常選用一個(gè)由測(cè)點(diǎn)坐標(biāo)P(x,y)的n階多項(xiàng)式表示的曲面來(lái)表示測(cè)網(wǎng)的區(qū)域重力異常信息。其中n階多項(xiàng)式的一般形式為：

g(x,y)=a0+a1x+a2y+a3x2+a4xy+a5y2+…+aN-1yn

(3)

(4)

g(x,y)即為趨勢(shì)分析方法得到的低頻區(qū)域信息,A為坐標(biāo)參數(shù)矩陣，ζ為真實(shí)測(cè)網(wǎng)重力信號(hào)。其方法是利用n階多項(xiàng)式曲線對(duì)實(shí)測(cè)重力異常數(shù)據(jù)進(jìn)行擬合，求算參數(shù)aN的最小二乘解，在一定誤差范圍內(nèi)對(duì)異常進(jìn)行分離。

3 模型驗(yàn)證

本文設(shè)計(jì)建立了由多個(gè)棱柱體的組合異常體來(lái)模擬二層地質(zhì)模型，淺層用4個(gè)大小相等的正方體來(lái)模擬地表中淺部的地下穿插小構(gòu)造，深部平鋪一個(gè)比較大的長(zhǎng)方體來(lái)模擬區(qū)域莫霍面的巖石圈重力異常。區(qū)域面積為100 km×100 km，采樣間距為0.5 km，模型示意圖和具體參數(shù)由圖2和表1列出。

表1 重力模型正演參數(shù)值

圖2 組合模型平鋪示意圖

對(duì)于BP神經(jīng)網(wǎng)絡(luò)隱含層神經(jīng)元采用對(duì)數(shù)S型傳遞函數(shù)logsig作為傳遞函數(shù)，輸出層神經(jīng)元?jiǎng)t采用purelin線性函數(shù)，網(wǎng)絡(luò)訓(xùn)練函數(shù)為traingdx，網(wǎng)絡(luò)迭代次數(shù)最多為1 000，期望誤差為10-5。局部異常一般是比區(qū)域地質(zhì)因素范圍小的研究對(duì)象引起，其異常范圍明顯小于整體區(qū)域范圍。為了研究BP神經(jīng)網(wǎng)絡(luò)擬合的效果，將重力正演網(wǎng)格偶數(shù)隔點(diǎn)抽稀的部分正演值作為學(xué)習(xí)集，一定程度上減弱了高頻局部信息的干擾。學(xué)習(xí)集用來(lái)對(duì)網(wǎng)絡(luò)進(jìn)行學(xué)習(xí)和訓(xùn)練，以建立平均重力異常與各點(diǎn)位坐標(biāo)的映射關(guān)系。將重力測(cè)網(wǎng)的坐標(biāo)值作為神經(jīng)網(wǎng)絡(luò)的輸入，測(cè)試集的重力異常值作為輸出，神經(jīng)網(wǎng)的學(xué)習(xí)訓(xùn)練相當(dāng)于對(duì)原始數(shù)據(jù)集進(jìn)行低通濾波處理。該低通濾波器作用是為了減弱高頻淺部異常體的干擾，突出深部大范圍區(qū)域異常特征。此外，測(cè)試集的其他加密網(wǎng)格值作為該測(cè)點(diǎn)的區(qū)域異常擬合值，實(shí)現(xiàn)了大區(qū)域范圍的低頻區(qū)域重力異常擬合，從而可以實(shí)現(xiàn)重力場(chǎng)的位場(chǎng)分離。由于輸入的是點(diǎn)位的坐標(biāo)即所建立坐標(biāo)系x,y的值，故輸入層神經(jīng)元取2個(gè)，輸出的是重力異常，輸出層神經(jīng)元取1個(gè)。隱含層神經(jīng)元的選取一般通過(guò)經(jīng)驗(yàn)性選取，通過(guò)測(cè)試集的多次試驗(yàn)，當(dāng)選取雙隱含層神經(jīng)元5個(gè)和1個(gè)時(shí)，擬合的結(jié)果相對(duì)最佳，因此選擇2×5×1×1的三層神經(jīng)網(wǎng)絡(luò)模型較為合適。

組合模型的重力異常圖及模型正演的局部和區(qū)域異常見(jiàn)圖3和圖4，四個(gè)小棱柱體異常和大尺度的區(qū)域異常很好地顯示出來(lái)，滿足趨勢(shì)分析法分析的條件。接下來(lái)利用趨勢(shì)分析法對(duì)上述理論重力模型進(jìn)行分析，比較分析得到二階多項(xiàng)式擬合區(qū)域異常的相似度最大，分離出的局部異常誤差最小，結(jié)果如圖5所示。局部異常往往包含著由地質(zhì)體巖性或密度不同引起的異常信息，備受地質(zhì)研究者的關(guān)注[21-24]。因此，接下來(lái)將這兩種方法進(jìn)行對(duì)比分析，對(duì)這兩種方法分離得到的重力局部異常信息與理論局部異常信息進(jìn)行誤差分析。

圖3 模型疊加異常

利用corr2函數(shù)可以將計(jì)算結(jié)果與理論模型數(shù)據(jù)(圖4a)來(lái)進(jìn)行相似度比較，對(duì)比其分離得出的局部異常與正演局部異常值的的近似程度。該函數(shù)返回值介于-1和1之間，返回值的模為0表示完全不同，模接近1則表示完全相似。經(jīng)過(guò)用趨勢(shì)分析法和BP神經(jīng)網(wǎng)絡(luò)法兩種方法分離重力異常，趨勢(shì)分析法的相似度約為0.61，而B(niǎo)P神經(jīng)網(wǎng)絡(luò)法的相似度約為0.91，對(duì)于重力異常的分離，可以比較出BP神經(jīng)網(wǎng)絡(luò)法分析的結(jié)果(圖6b)極大地優(yōu)于趨勢(shì)分析法的結(jié)果(圖5b)。

圖4 模型局部和區(qū)域重力異常正演

圖5 趨勢(shì)分析法分離異常

圖6 BP神經(jīng)網(wǎng)絡(luò)法分離異常

corr2(A,B)=

(5)

一般而言，局部重力異常包含著高頻短波長(zhǎng)異常信息源，是礦產(chǎn)勘查、斷層結(jié)構(gòu)的重要信息源。不同地質(zhì)體的邊界、斷裂構(gòu)造線等線性構(gòu)造兩側(cè)往往具有明顯的密度差異，重力異常正負(fù)差異帶可以凸顯橫向線性構(gòu)造。Theta map方法[25]利用總梯度振幅歸一化的總水平梯度模，這種歸一化有效地提升了增益控制，利用該方法對(duì)分離得到的局部異常信息分析，能有效刻畫出局部異常體的輪廓邊界。

(6)

對(duì)BP神經(jīng)網(wǎng)絡(luò)分離得到的局部異常體信息的Theta map成像如圖7所示，可以清晰地發(fā)現(xiàn)中淺部四個(gè)棱柱體的水平邊界，勾畫出了局部異常源的橫向特征分布，與重力正演模型參數(shù)表1中的場(chǎng)源參數(shù)比較十分吻合，證實(shí)了理論模型分析的可靠性。

圖7 BP算法分離的局部異常的邊界信息

4 實(shí)測(cè)數(shù)據(jù)分析

接下來(lái)選取了Gooper于2008年于SEG官網(wǎng)上提供的南非Witwatersrand Basin某區(qū)域的航空重力數(shù)據(jù)集進(jìn)行測(cè)試。約翰內(nèi)斯堡區(qū)域地層廣泛外露太古宙的花崗質(zhì)巖石，盆地中沉積大量巨厚的元古宙火山-系沉積巖，區(qū)域黃金礦產(chǎn)資源豐富[26]。此航空測(cè)網(wǎng)大小為200 km×200 km，采樣間隔為1 km，該測(cè)網(wǎng)區(qū)域位于盆地西南方，有一個(gè)類似“火山口”特征的天然隕石坑，周圍地質(zhì)表層構(gòu)造信息豐富。實(shí)際資料處理中選擇區(qū)域異常的樣本往往是需要進(jìn)行多類測(cè)試對(duì)比分析的。此處參考前文正演測(cè)試方法，采用偶數(shù)抽稀隔點(diǎn)選取的數(shù)據(jù)點(diǎn)作為訓(xùn)練學(xué)習(xí)集，擬合的原始分辨率網(wǎng)格數(shù)據(jù)點(diǎn)作為低頻區(qū)域異常輸出，通過(guò)分離結(jié)果的驗(yàn)證，成功建立了平均重力異常與各點(diǎn)位坐標(biāo)的映射關(guān)系。測(cè)網(wǎng)布格重力異常值及BP神經(jīng)網(wǎng)絡(luò)分離得到的重力異常如圖(8)和圖(9)所示。

從測(cè)區(qū)的布格重力異常分布(圖8)可以發(fā)現(xiàn)測(cè)區(qū)的高低密度異常體的大致分布，經(jīng)BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練得到的區(qū)域異常場(chǎng)(圖9a)，其東北部區(qū)域有著明顯的相對(duì)高密度異常，其東南側(cè)角有著明顯的相對(duì)低密度異常。分離得出的局部異常信息(圖9b)在測(cè)區(qū)右上角區(qū)域有著一定范圍的異常源，并廣泛分布著部分相對(duì)高密度異常信息。經(jīng)由Theta map成像(圖10)，測(cè)區(qū)東北部的隕石坑邊界輪廓突出明顯，測(cè)網(wǎng)區(qū)域較為明顯的條帶狀信息推測(cè)為豐富的斷層、褶皺信息。結(jié)合測(cè)區(qū)詳細(xì)的地質(zhì)資料，發(fā)現(xiàn)測(cè)區(qū)的隕石坑邊界信息，與 BP 神經(jīng)網(wǎng)絡(luò)分離得出的火山口位置信息十分吻合，相對(duì)高密度異常區(qū)域邊界可判斷出明顯的斷層、褶皺信息，推測(cè)可能為深部集礦區(qū)域。中淺部的局部異常體信息很好地被提取出來(lái)，驗(yàn)證了該方法在實(shí)際重力資料數(shù)據(jù)處理中的適用性。

圖8 測(cè)區(qū)布格重力異常

圖9 BP神經(jīng)網(wǎng)絡(luò)法分離實(shí)際重力異常

圖10 BP算法分離的局部異常的邊界信息

5 結(jié) 論

本文通過(guò)重力正演構(gòu)造典型地質(zhì)模型，對(duì)于區(qū)域地質(zhì)因素范圍影響較小的局部異常源，通過(guò)BP神經(jīng)網(wǎng)絡(luò)模型和趨勢(shì)分析模型兩類位場(chǎng)分離測(cè)試。BP神經(jīng)網(wǎng)絡(luò)具有較好的非線性映射和無(wú)模型估計(jì)能力，其分離結(jié)果明顯優(yōu)于趨勢(shì)分析模型，提取出的局部異常theta map成像精度更高。在實(shí)測(cè)重力測(cè)網(wǎng)數(shù)據(jù)分析中亦發(fā)現(xiàn)，其在實(shí)際地球物理的數(shù)據(jù)資料解譯也具有一定的適用性。但值得一提的是，BP神經(jīng)網(wǎng)絡(luò)方法在復(fù)雜地質(zhì)模型情況下，較難找出能夠代表全區(qū)范圍內(nèi)的重力場(chǎng)來(lái)作為訓(xùn)練樣本，訓(xùn)練過(guò)程容易陷入局部極值，需要選取不同區(qū)域范圍進(jìn)行多類測(cè)試。

致 謝

感謝Gooper于2008年在SEG官網(wǎng)上提供的南非Witwatersrand Basin區(qū)域內(nèi)的航空重力實(shí)測(cè)數(shù)據(jù)。