基于重力異常特征點的等軸狀場源中心埋深自動化估計算法

張 攀，杜勁松,2，邱 峰

(1.中國地質(zhì)大學 地球物理與空間信息學院 地球內(nèi)部多尺度成像湖北省重點實驗室，湖北 武漢 430074；2.中國地質(zhì)大學 地質(zhì)過程與礦產(chǎn)資源國家重點實驗室，湖北 武漢 430074)

1 引 言

在重力勘探中，特征點法是根據(jù)重力異常曲線上的一些特征點(如極大值點、1/n極值點、零點、拐點等)(其中n>1)的重力異常值以及它的相對坐標位置求取場源的幾何參數(shù)和物性參數(shù)的一種方法[1,2]。

當?shù)刭|(zhì)體的最大與最小延伸尺度之差ΔL與最小埋藏深度D滿足ΔL<0.2ΔD時，可將場源近似為球體。在實際應用中，鹽丘、礦囊、巖體等近似等軸狀的地質(zhì)體均可以等效為球體進行處理。本文根據(jù)特征點法的特點，設計自動化算法，當輸入是一個二維的重力數(shù)據(jù)時，首先自動找出異常極值點，其次判斷其是否為等軸狀，然后對反演可用的1/n極值點范圍進行搜索，挑取多個1/n極值點進行特征點法深度反演，最后對反演結(jié)果進行分析，確定最終結(jié)果。隨后，將該算法分別應用于模型數(shù)據(jù)中，得到了比較精確的結(jié)果，在加噪聲的模型數(shù)據(jù)和實測數(shù)據(jù)中也得到了良好的結(jié)果[3-8]。

2 自動化反演方法

2.1 特征點反演方法

對于一個各向同性的均勻球體來言，其引起的重力異常與同質(zhì)量的點源引起的重力異常是完全相同的[10-14]。將球心在地面的投影處記為O點，由于球體具有對稱性，只需研究通過O點的任意一個方向的剖面重力異常即可，離O點的距離為x的一點處的重力異常Δg，單位 m/s2，計算表達式為

(1)

式中，D為球體中心的埋藏深度，單位 m；ΔM為球體的剩余質(zhì)量，單位t；x為測點離球心在地面的投影O點之間的距離，單位 m；G為萬有引力常數(shù)(6.672×10-11N·m2/kg2)。

在異常曲線上任意取出兩個點，設它們的坐標為x1和x2，相應的異常值為Δg1和Δg2，則根據(jù)式(1)可得

(2-1)

(2-2)

將上式求比值以消去G和ΔM，將其比值記為w，則

(3)

進而，可以推導得到中心埋深D的顯示表達式為

(4)

對式(4)進行簡化，令x1=0，在此處異常具有極大值，即Δgx1=Δgmax，設x2為1/n極值，則有

(5)

(6)

式(5)和式(6)即為球體中心埋藏深度的反演計算通式。

在得到球體中心埋藏深度之后，若已知剩余密度Δσ，單位 g/cm3；進而可以得到球體的剩余質(zhì)量ΔM、半徑R、上頂埋深hu、下底埋深hl分別為

2.2 自動化算法

2.2.1 極值點搜索

平面重力數(shù)據(jù)為一個二維的數(shù)組，在尋找極大值時，先預分配兩個大小與原始數(shù)組相同的零矩陣[15-18]，記其為A1和A2。首先沿著數(shù)組的每一行進行極大值的搜索，如果當前位置的值為極大值，則將A1中相同的位置標記為1，之后沿著數(shù)組的每一列進行極大值的搜索，如果當前位置的值為極大值，則將A2中相應的位置標記為1，最后將A1和A2進行點乘，其結(jié)果記為A。在A矩陣中，1出現(xiàn)的地方即為平面數(shù)據(jù)中極大值點的位置。

但是，在數(shù)據(jù)含有一定的噪聲時，這種標定方法會將所有噪聲造成的小突或凹點均標定出來。為了解決這個問題，本算法在標定極值的過程中引入了三種判定準則。

1)當標定的極值的絕對值小于一定的閾值之后，則舍去這個標定出的極值；

2)以極值為中心取一個滑動窗口，如果該極值的絕對值不是當前窗口中的最大值時，則舍棄不采用，滑動窗口的大小視具體情況而定；

3)重力異常主體如果一部分在邊界之外，則不適合進行特征點法反演，則對離邊界很近的一些極值點不進行標定。

通過這三種準則的約束，即使重力異常數(shù)據(jù)存在一定的噪音，也能對數(shù)據(jù)中的極大值點進行準確的標定。

2.2.2 判斷異常是否為(似)等軸狀

以標定出的極大值點為中心，沿N個方向獲取N條剖面，在這N個方位通過線性插值尋找1/n極值點，在正常情況下，將會搜索到2N個1/n極值點，在本算法中，引入兩種評價標準用于判斷異常是否能近似為等軸狀。

1)對這2N個1/n極值點的位置進行橢圓擬合[3]，求出它們組成的橢圓的離心率，若離心率過大，則證明該異常不能近似為等軸狀；

2)根據(jù)2N個1/n極值點距離極值點的距離，求出距離的極差，由于重力異常有大有小，極差不能進行比較，則將極差與距離的最大值進行比值計算，將其定義為相對極差，若相對極差較大，則證明該異常不能近似為等軸狀。

如果異常比較特殊，沿某一個方向沒有搜索到1/n極值點，則直接將異常認定為非等軸狀異常。這一步可以選擇一個1/n極值點或者多個1/n極值點進行綜合評價。經(jīng)過這一步，數(shù)據(jù)中的非等軸狀異常則會被排除，而不參與后面的反演計算。

2.2.3 確定反演的1/n極值點

在這個環(huán)節(jié)，需要挑選多個1/n極值點共同反演，最后的結(jié)果挑選比較可靠的一部分進行均值計算。設起始搜索的1/n極值為1/n_start 極值，截止搜索的1/n極值為1/n_end 極值。

n_start的選取準則為：極值與1/n_start 極值之間的距離不能小于測點平均點距。令n從n_start開始，步長為dn，不斷地增大n，并且求出當前搜索到的1/n極值擬合橢圓的離心率和相對極差。本算法引入四種判斷準則來確定n_end，如下：

1)如果沒有搜索到當前的1/n極值，則在前一個步長停止；

2)如果當前搜索到的1/n極值擬合橢圓的離心率過大，則證明異常當前的形態(tài)已經(jīng)不符合等軸狀的特征，則在前一個步長停止；

3)如果當前搜索到的1/n極值的相對極差過大，則證明異常當前的形態(tài)已經(jīng)不符合等軸狀的特征，則在前一個步長停止；

4)如果當前搜索到的1/n極值小于3倍的觀測數(shù)據(jù)誤差標準差，則在前一個步長停止。

2.2.4 在n_start 到n_end范圍內(nèi)進行特征點法深度反演

將n_start到n_end范圍中每條剖面取出M個1/n極值點進行特征點法深度反演，得到N×M個深度點，同時求出當前1/n極值搜索到的1/n極值點的擬合橢圓的離心率和距離的相對極差作為結(jié)果是否準確的評判標準。

沿每條測線得出的深度值可能不同，若直接將結(jié)果進行平均計算，則一些由噪聲引起的“奇異值”會影響結(jié)果的準確性。為了解決這個問題，本算法在這一步引入了離散系數(shù)Vs作為評判準則，離散系數(shù)的定義為：

(8)

3 模型試驗

3.1 簡單模型試驗

設計點源質(zhì)量為3×1011kg，中心點埋深為1 000 m。數(shù)據(jù)共201條測線，每條測線上面具有201個采樣點，數(shù)據(jù)的點距和線距均為100 m，正演得到地面上的重力異常，如圖1所示。

將模擬數(shù)據(jù)輸入算法，先后經(jīng)過極值點的搜索、異常是否為等軸狀的判斷、反演范圍的搜索以及反演。在計算過程中，總共沿四個方向均勻的拉出四條剖面，將n_start到n_end范圍中均勻取出30個1/n極值點進行特征點法深度反演，n_end的截止最大值為6，計算結(jié)果如圖2所示。

圖2 反演結(jié)果Fig.2 Inversion results

將所有的反演結(jié)果直接進行平均，得到的深度反演結(jié)果為1 004.4 m，相對誤差為0.44 %。根據(jù)式(5)可見，若n值較小，可能引起埋深偏大。因此，根據(jù)離散系數(shù)，將第一個反演深度數(shù)據(jù)舍去之后再求平均，得到的深度反演結(jié)果為1 001.2 m，相對誤差降為0.12 %，結(jié)果也更為準確。

3.1.1 噪聲的影響分析

對于實際數(shù)據(jù)，其不可避免地會含有噪聲干擾。為了更好地符合實際情況，將上述重力數(shù)據(jù)中加入一定量的噪聲再進行計算，用于檢驗算法的穩(wěn)定性。由表1所示，由于每次加入的隨機噪聲幅度不同，計算結(jié)果也存在差異?？梢园l(fā)現(xiàn)，隨著噪聲的增大，結(jié)果的相對誤差總體上呈現(xiàn)出增大的趨勢，而且計算結(jié)果的離散系數(shù)也呈現(xiàn)出增大的趨勢。雖然算法顯示具有比較穩(wěn)定的性能，但是在實際應用之前，還是建議預先對實測數(shù)據(jù)進行適當?shù)娜ピ胩幚怼?/p>

表1 添加不同幅度噪聲情況之下的反演結(jié)果

3.1.2 背景場的影響分析

特征點法具有其應用條件，該方法是基于單個異常體的重力異常表達式推導得出的，而實測異常是多個異常體疊加而成的，一般也含有背景場，致使極值和1/n極值點之間的相對關系就會發(fā)生變化，故運用特征點法之前必須去除背景場。但是，在實際應用中，目前的位場分離技術(shù)并不能完美地將局部場和背景場分離，或多或少還是會存在一定的剩余。因此，為了更好地符合實際情況，將上述重力數(shù)據(jù)中加入背景場進行計算，用于分析背景場對特征點算法的影響。

在前述模擬的重力異常數(shù)據(jù)中加入不同的常值作為背景場，計算結(jié)果如表2所示。如果加入背景場的數(shù)據(jù)為正值，則隨著1/n極值中n的增大，反演的深度會越來越大；如果加入的背景場的數(shù)值為負值，則隨著1/n極值中n的增大，反演的深度將會越來越小。當背景場幅度為局部異常幅值的10 %時，反演埋深的相對誤差可達18 %左右。因此，背景場對于反演深度結(jié)果的影響非常大，建議在實際應用之前，采用不同位場分離技術(shù)，進行綜合判斷。

表2 不同背景場下的反演結(jié)果

表3 不同非等軸狀場源的反演結(jié)果

3.1.3 場源非等軸狀的影響

在實測重力異常數(shù)據(jù)中，很少能夠找到完美的等軸狀異常。為了更好地符合實際情況，將非等軸狀場源正演得到的重力異常數(shù)據(jù)輸入算法進行測試，用于分析場源非等軸狀的影響。

采用寬與高分別為1 000 m與1 000 m，剩余密度為0.4 g/cm3，埋深為2 000 m的直立棱柱體作為場源，不斷改變棱柱體的長，最后得到的反演深度結(jié)果如表3所示?？梢钥闯?，場源非等軸狀對反演深度造成的影響很大，當長寬比超過2∶1之后，反演深度的相對誤差會超過7 %。因此，在實際數(shù)據(jù)反演時，必須結(jié)合離心率和相對極差評判結(jié)果的可靠性。

3.2 組合模型實驗

對于實測重力異常數(shù)據(jù)，往往包含多個異常。因此，設計組合模型實驗，采用兩個點源模型模擬等軸狀異常，采用一個棱柱模型模擬非等軸狀異常。第一個點源質(zhì)量為3×1011kg，中心埋深為2 000 m；第二個點源質(zhì)量為1.2×1011kg，中心埋深為700 m；采用的棱柱模型的長寬高分別為5 000 m、200 m、40 m，中心埋深為700 m；三個模型的剩余密度均為0.4 g/cm3。模擬得到的重力異常如圖3所示。

圖3 組合模型模擬的重力異常Fig.3 Synthetic gravity anomaly by combined density models

首先進行重力異常極大值點的標定，接著在極大值附近搜索1/5極值的位置用于橢圓擬合和相對極差的計算。經(jīng)過計算，左上角非等軸狀異常的擬合離心率為0.894 7，右上角等軸狀異常的擬合離心率為0.191 4，下方等軸狀異常的擬合離心率為0.127 8；左上角非等軸狀異常的相對極差為0.566 9，右上角等軸狀異常的相對極差為0.041 8，下方等軸狀異常的相對極差為0.016 3?？梢姡笊辖堑闹亓Ξ惓５碾x心率和相對極差均較大，因此舍去該異常不對其進行后續(xù)反演。

之后，對標定出的異常進行反演范圍的搜索和深度反演，得到下方相對寬緩的等軸狀異常的深度為2 011.0 m，相對誤差為0.55 %，根據(jù)離散系數(shù)將一些不穩(wěn)定結(jié)果剔除之后，得到的反演深度為2 007.0 m，相對誤差為0.35 %；右上角相對尖銳的異常的深度為704.2 m，相對誤差為0.60 %，同樣根據(jù)離散系數(shù)將一些不穩(wěn)定結(jié)果進行剔除，最后得到的反演深度為702.8 m，相對誤差降為0.40 %。

組合模型試驗表明，對于擁有多個等軸狀和非等軸狀異常的數(shù)據(jù)，算法可以挑選出等軸狀的異常并且進行反演，反演結(jié)果也較為準確，從而證明了本文算法的可行性。

4 實際應用

將上述算法應用于位于美國路易斯安那州的Vinton鹽丘，測區(qū)位于路易斯安那州與德克薩斯州交界處，測區(qū)以南即為著名地石油產(chǎn)地墨西哥灣。測區(qū)多為沉積巖，底部發(fā)育鹽丘構(gòu)造，鹽丘區(qū)經(jīng)常形成油氣圈閉，因此對鹽丘進行研究具有實際意義[4，5]。

圖4 Vinton鹽丘重力異常Fig.4 Gravity anomaly of Vinton salt dome

文中使用的數(shù)據(jù)為航空重力測量數(shù)據(jù)，總測線長度為1 087.5 km，測區(qū)面積為192.6 km2，數(shù)據(jù)的點距和線距均為30 m，地形改正密度值為1.8 g/cm3。本文截取面積為61 69 m×6 169 m的具有目標重力異常的區(qū)域(圖4)，使用本文算法進行埋深估計，進而與前人計算結(jié)果進行對比。

由模型試驗可知，噪聲以及背景場對反演計算均會產(chǎn)生一定的影響。因此，首先采用遞歸濾波方法，截斷波長設置為294.19 m，補償系數(shù)為1，對Vinton鹽丘重力異常數(shù)據(jù)進行了去噪處理，結(jié)果如圖5所示；然后，采用匹配濾波法，分頻點波長設置為5 589.68 m，對去噪之后的重力異常數(shù)據(jù)進行局部場和趨勢場分離，結(jié)果如圖6所示。

圖5 Vinton鹽丘重力異常去噪處理Fig.5 Noise reduction for gravity anomaly of Vinton salt dome

圖6 Vinton鹽丘重力異常位場分離結(jié)果Fig.6 Potential field separation for gravity anomaly of Vinton salt dome

圖7 反演結(jié)果Fig.7 Inversion results

表4 本文與前人反演的Vinton鹽丘幾何參數(shù)

計算得到的結(jié)果如圖7所示。設鹽丘的密度[6,10]為2.7 g/cm3，選取離散系數(shù)小于0.1的計算結(jié)果。表4對比了本文與前人反演的Vinton鹽丘的幾何參數(shù)，對比顯示本文所得中心與上頂埋深與前人計算結(jié)果一致，但是下底埋深大于前人結(jié)果。邱峰等(2020)[17]認為該鹽丘的構(gòu)造指數(shù)為1.198，其形態(tài)比較接近于柱體狀而非球狀(球體的構(gòu)造指數(shù)為3)，重力異常由鹽丘高密度、扁平的蓋層引起而非鹽丘本身。場源的非等軸狀導致本文計算所得的下底埋深偏大。

5 結(jié) 論

本文將重力勘探中的特征點反演方法進行了計算機自動化處理，使其能夠識別重力異常平面數(shù)據(jù)中的(似)等軸狀異常，進而進行場源中心埋深反演與誤差統(tǒng)計，并且輸出擬合離心率、相對極差和離散系數(shù)作為評價結(jié)果準確性的參數(shù)。

理論模型試驗表明，數(shù)據(jù)包含的噪聲對反演結(jié)果存在一定的影響，但是影響較弱，而背景場對反演結(jié)果影響較大。因此，建議在實際反演之前，對實測數(shù)據(jù)進行適當?shù)娜ピ肱c分場處理。此外，理論模型試驗顯示，本文算法可以準確地找出數(shù)據(jù)中的(似)等軸狀異常，并且比較準確地計算得出異常體的中心埋深。將本文算法應用于美國路易斯安那州的Vinton鹽丘的航空重力數(shù)據(jù)，得到了較好的結(jié)果，與前人計算埋深一致。理論模型試驗和實際應用均證明了本文算法的可行性與實用性。

致謝

感謝美國Bell Geospace公司提供Vinton dome地區(qū)的重力異常數(shù)據(jù)！