重力異常曲化平的改進插值—迭代法

楊婧,郭良輝

(中國地質大學(北京) 地球物理與信息技術學院，北京 100083)

通訊作者： 郭良輝(1980-)，男，教授，博導，主要從事地球物理數據處理反演新方法及應用研究工作。Email: guo_lianghui@163.com

0 引言

高精度重力測量可以為區域地質構造研究、資源與能源勘查等提供重要的基礎數據。野外起伏觀測面采集的重力數據，經過各項校正得到的重力異常仍然是原起伏觀測面上的異常，反映地下密度不均勻體在起伏觀測面上的重力效應[1]。頻率域算法是地球物理領域的一種常用快速算法，已廣泛應用于重力數據處理與反演解釋，比如頻率域延拓[1-2]、頻率域優化濾波[3]、頻率域界面反演[4]、頻率域三維成像[5]等。然而，這些頻率域算法通常要求重力異常所對應的觀測面為平面。在觀測面起伏較大情況下，若對起伏面重力數據直接應用常規頻率域算法做處理與反演，則將帶來不可忽略的處理偏差。因此，在起伏觀測面條件下，重力異常數據從起伏觀測面化到平面上，即曲化平，是開展頻率域精細處理解釋的必要步驟。

當前重力數據的曲化平方法較多，Pilkington等[6]將其分為基于源和基于場的兩大類?；谠吹姆椒ㄖ饕械刃г捶╗7-12]，它先對起伏觀測面異常反演獲得地下等效場源，然后再正演計算等效源在任意水平面引起的異常，該類方法通常計算量較大，適用于中、小規模數據的曲化平處理?；趫龅姆椒ㄓ杏邢拚{和級數法[13]、泰勒級數法[14-15]、有限元法[16-17]、積分法[18-19]、逐次逼近法[20]等，該類方法無需反演場源，通常需要向上延拓或向下延拓到某中間層面進而獲取最終平面的異常。向上延拓是穩定的，但向下延拓通常是不穩定的，這是由于它具有放大作用，對高頻噪聲干擾比較敏感，延拓跨度越大，高頻干擾越嚴重。因此，基于向下延拓的曲化平方法一般迭代收斂速度慢、精度有限，僅適用于延拓跨度小的曲化平。

徐世浙等[21]提出了位場(重、磁場)曲化平的插值—迭代法，屬于基于場的方法類，它計算簡單快速，可以實現延拓跨度達10～20倍點距的大規模數據的穩定曲化平。該方法基本思路是，將起伏觀測面(A)的頂部水平面(Bm)與底部水平面(B)之間的空間剖分為若干等間距的中間層水平面(Bi)，見圖1所示；接著，把起伏觀測面(A)的位場值放到底部水平面(B)并作為該水平面的位場初始值；然后，通過頻率域向上延拓計算出各中間水平面(Bi)的位場值，進而從各中間水平面(Bi)集合的位場值插值出起伏觀測面(A)的位場近似值；之后，計算起伏觀測面(A)的位場實測值與近似值的殘差，并用于對底部水平面(B)的位場值進行修正；如此反復迭代，直至起伏觀測面(A)的位場實測值與近似值的殘差達到誤差容限。

圖1 插值—迭代法曲化平示意Fig.1 The diagram of interpolation-iteration method for gravity anomaly continuation from undulating surface to plane

然而，在實際應用中，本文發現上述插值—迭代方法[21]在觀測面起伏較大、延拓跨度較大的復雜條件下的曲化平效果還是有限的(見后面理論模型數據試驗的圖4所示)。為此，本文對上述插值—迭代方法做進一步的改進，在異常迭代修正過程中引入起伏觀測面修正因子，加快曲化平迭代收斂，促進曲化平效果提升。本文首先介紹改進的插值—迭代法的方法原理和算法流程，然后利用理論模型和川滇地區實際數據試驗驗證方法的有效性。

1 方法原理

本文曲化平方法是在徐世浙等[21]的插值—迭代方法基礎上改進得到， 如圖1所示，假定起伏觀測面A和底部水平面B，曲化平的目標是由起伏觀測面A的位場數據pA通過曲化平處理獲得底部水平面B的位場數據pB。根據起伏觀測面的高程最大值給定頂部水平面Bm，進而將底部水平面和頂部水平面之間的空間按照等間距剖分為m層，各中間層水平面為Bi(i=1,2,…,m)。

(1)

(2)

在徐世浙等[21]的插值—迭代方法中，修正因子S取常數，且通常取為1，即觀測面上各測點的修正因子相同，導致在相同迭代次數下，觀測面起伏變化大的地方曲化平迭代收斂速度慢、效果差些。因此，本文對修正因子S做了改進，采用與起伏觀測面相關的修正因子，從而加快曲化平迭代收斂速度，提升曲化平效果。本文的修正因子公式如下：

(3)

其中，T為起伏觀測面高程，Tmin和Tmax分別是起伏觀測面高程的最小值和最大值?？芍?，觀測面起伏越大，修正因子S將越大，曲化平收斂速度將越快，反之，觀測面起伏越小，修正因子S將越小，曲化平收斂速度將越慢。一般，n取非負數，若n取值過大，則修正因子S將越大；反之，n取值越小，修正因子S將越??；當n=0時，修正因子S=1，等效于徐世浙等[21]的算法。

圖2 本文改進的插值—迭代法算法流程Fig.2 The workflow of the modified interpolation-iteration method

2 理論模型數據試驗

理論模型為一個直立長方體，長、寬、高分別為 2 000、4 000、1 000 m，中心埋深1 500 m，剩余密度為1.0 g/cm3。假定測網E向和N向范圍均為0～10 km、點距均為0.1 km，則理論模型在海拔0 km的水平面上產生的理論重力異常如圖3a所示，異常等值線為近橢圓狀，最小值為0.23 mGal，最大值為12.65 mGal，平均值為2.41 mGal。假定起伏觀測面高程如圖3b所示，高程最小值為1.46 m，最大值為2 023.55 m，觀測面高差是測網點距的20多倍。理論模型在起伏觀測面上產生的理論重力異常如圖3c所示，異常等值線為不規則形狀，且與起伏觀測面形狀有一定相關，異常最小值為0.23 mGal，最大值為7.52 mGal，平均值為1.88 mGal。因此，受觀測面起伏影響，理論重力異常變為復雜，難以直接處理解釋，值得優先進行曲化平處理。

圖4a和圖4b顯示了常規插值—迭代方法[21]迭代20次的曲化平結果及其與海拔0 m平面理論重力異常(圖3a)的殘差，修正因子中n取0，即S=1。常規插值—迭代法曲化平增強了有效信號，最小值為0.23 mGal，最大值為13.78 mGal，但結果與海拔0 m平面理論重力異常有明顯的差別。異常殘差特征與起伏觀測面比較相似，在觀測面平緩區殘差較小，起伏區殘差較大，殘差最小值為-2.15 mGal，最大值為3.73 mGal，標準差為0.39 mGal。因此，在觀測面起伏較大、延拓跨度較大的情況下，常規插值—迭代方法的曲化平效果是有限的。

a—海拔0 m平面的理論重力異常；b—起伏觀測面高程；c—起伏觀測面的理論重力異常a—theoretical gravity anomaly at 0 m altitude；b—elevation of undulating observation surface；c—theoretical gravity anomaly at undulating observation surface圖3 理論模型起伏觀測面與重力異常Fig.3 Undulating observation surface of theoretical model and its gravity anomaly

a、b—常規插值-迭代法曲化平結果及與海拔0 m平面理論重力異常殘差；c、d—本文改進的插值-迭代法曲化平結果及與海拔0 m平面理論重力異常殘差a、b—the result of the routine interpolation-iteration method for continuation from undulating surface to plane and its deviation from the values of plane surface；c、d—the result of the modified interpolation-iteration method and its deviation from the values of plane surface圖4 不同曲化平方法的結果及與海拔0 m平面理論重力異常對比Fig.4 Comparison between the results of different interpolation-iteration methods and the gravity forward values of plane surface

圖4c和圖4d顯示了本文改進的插值—迭代方法迭代20次的曲化平結果及其與海拔0 m平面理論重力異常(圖3a)的殘差，其中修正因子采用式(3)且經過測試，選擇效果較好的n=1.5。

圖5顯示了剖面A不同曲化平方法結果與起伏觀測面、海拔0 m平面的理論重力異常、起伏觀測面高程的對比。觀察發現，水平面的理論重力異常出現一個異常峰值，對應地下密度異常體，起伏觀測面下的重力異常出現兩個異常峰值，其位置與起伏觀測面相對應，這說明起伏觀測面扭曲了地下密度體的異常信息，觀測面起伏越高，扭曲作用越大。n的取值與實際研究區觀測面起伏程度相關，n=0時，S為常數1，等效于徐世浙等[21]的算法，其均方根誤差為0.1551 mGal；n=1時，其均方根誤差0.044 5 mGal；n=1.5時，其均方根誤差最小，為0.028 7 mGal；n=2時，其均方根誤差0.079 4 mGal，故本文選n=1.5。

改進的插值—迭代法有效增強了信號，曲化平后其最小值為0.23 mGal，最大值為12.33 mGal，結果與海拔0 m平面理論重力異常差別較小。異常殘差最小值為-0.86 mGal，最大值為1.03 mGal，標準差為0.17 mGal，我們分析這一微小差異是由高程差距引起收斂速度不一致產生的。因此，與常規插值—迭代法相比，本文改進的插值—迭代法曲化平結果更接近海拔0 m平面理論重力異常值，效果更優，適用于大規模數據、復雜起伏觀測面、大延拓跨度的曲化平。

圖5 剖面A不同修正因子S曲化平的對比Fig.5 Different S comparison along profile A

3 實際數據試驗

川滇地區構造上位于印支塊體、青藏高原塊體和揚子塊體的交匯處，區域構造動力作用復雜，在緬甸弧東向俯沖、印度板塊藏東構造結NE向楔體擠出和俯沖、青藏高原隆升和SE向擠出等共同作用下，該地區整體向E—SE方向擠出，形成了大陸內部典型的剪切、拉張和推覆構造，表現出地殼變形速率高、斷層運動劇烈、強震原地復發周期短等區域特征。因此，川滇地區是研究青藏高原隆升過程與大陸強震孕育機理的天然實驗場，也是我國區域防震減災的重要地區。高精度重力數據是川滇地區深部結構與構造研究的基礎數據。由于該地區地形上呈現青藏高原、云貴高原、四川盆地等復雜、多樣化的起伏形態，因此對該地區起伏地形面的重力異常數據作常規頻率域處理和解釋應用之前，值得優先作曲化平處理。

本文從ICGEM官網(http://icgem.gfz-potsdam.de/calcgrid)下載了川滇地區地形高程網格數據和布格重力異常網格數據，東經99°～108°，北緯24°～32°，網格大小均為0.1°×0.1°。其中，布格重力異常數據是由地球重力場模型 (earth gravitational model 2008)的自由空氣重力異常數據經過標準的地形校正和中間層校正得到的。本文利用窗口大小為3的均值濾波對研究區地形高程數據做了平滑處理，結果見圖6a所示，研究區地形起伏較大，高程最大值為4 882 m，最小值為269.2 m，高差約為 4 612.8 m。由于原始布格重力異常數據存在明顯的高頻噪聲干擾，這主要是由于自由空氣重力異常數據和地形數據的分辨率與精度不一致造成的[22]，因此，本文對原始布格重力異常數據作了截止波長為100 km的低通濾波，得到去噪后的川滇地區布格重力異常，如圖6b所示，異常最大值為-64.51 mGal，最小值為-495.5 mGal，青藏高原地區異常較低，四川盆地等海拔低地區異常較高，云貴高原介于兩者之間。

本文應用改進的插值—迭代法對去噪后的布格重力異常作曲化平處理，修正因子采用式(3)，n=1.5，經過50次迭代，將其化到海拔0 m的位置，結果如圖6c所示。曲化平后的布格重力異常最大值為-64.47 mGal，最小值為-508.9 mGal，山區異常信號有效增強，尤其是在川西地區增強較為明顯、異常細節更為突出。圖6d顯示了曲化平前后的布格重力異常殘差，最大值為13.36 mGal，最小值為-8.45 mGal，殘差變化主要集中在青藏高原和云貴高原等山區，與地形有較強對應關系。

為進一步探討地形對重力異常值的影響及曲化平的效果，本文選取一條斜向穿過研究區的剖面B(位置如圖6a所示)，該剖面地勢復雜，地形起伏變化較大。該剖面曲化平前后的重力異常對比(圖7)顯示，地形起伏較大地區的曲化平效果顯著，地形起伏較為平緩地區的曲化平效果較弱，這一結果與前文的分析相符。

本文以常規頻率域垂直導數換算為例作說明曲化平的意義。導數換算是重力異常分析解釋的一個常用步驟，用于突出淺部場源、分析構造邊界等。圖8a和圖8b顯示了曲化平前后的布格重力異常的垂直導數，曲化平前的垂直導數數值范圍為-0.002 9～0.002 3 mGal/m，曲化平后的數值范圍轉為-0.004 4～0.003 7 mGal/m?？梢?，曲化平后的異常幅值得到明顯增強，而且刻畫出更多、更清晰的異常變化細節，從而反映出更多的淺部場源細節和構造邊界信息。

a—川滇地區地形高程；b—去噪后的布格重力異常；c—本文改進插值—迭代法曲化平后布格重力異常；d—曲化平前后的異常殘差a—the elevation of Sichuan-Yunnan region；b—the denoised Bouguer gravity anomaly；c—the Bouguer gravity anomaly after continuation from undulating surface to plane by using the modified interpolation-iteration method；d—the anomaly deviations before and after continuation from undulating surface to plane圖6 川滇地區曲化平前后重力異常對比Fig.6 The Bouguer gravity comparison before and after continuation from undulating surface to plane of Sichuan-Yunnan region

圖7 剖面B的對比Fig.7 The gravity anomaly comparison and the elevation along profile B

a—布格重力異常的垂向導數；b—曲化平后布格重力異常的垂向導數a—the vertical derivative of Bouguer gravity anomaly；b—the vertical derivative of Bouguer gravity圖8 川滇地區曲化平前后垂向導數對比Fig.8 The vertical derivative of Bouguer gravity comparison before and after continuation from undulating surface to plane of Sichuan-Yunnan region

圖8中標注了2008年汶川8級地震、2013年蘆山7級地震、2014年魯甸6.5級地震、2021年漾濞6.4級地震的震中位置，易見經過曲化平這些地震震中周緣地區的重力異常變化特征更為鮮明，這有助于后續的深部結構與構造解釋研究。

4 結論與討論

本文在重力異常曲化平的常規插值—迭代法基礎上給出了一種改進的插值—迭代法，即在異常迭代修正過程中引入起伏觀測面修正因子，加快曲化平迭代收斂，促進曲化平效果提升，實現適用于觀測面起伏較大、延拓跨度較大的復雜條件下重力異常曲化平。理論模型和川滇地區實際數據試驗驗證了本文方法的有效性，效果優于常規插值—迭代法。本文改進的插值—迭代方法不僅適用于重力異常數據，也適用于磁異常數據以及梯度、三分量等數據。本文曲化平的結果可進一步通過常規平化平、平化曲處理實現任意平面或曲面的延拓。