重力異常曲化平的改進(jìn)插值—迭代法

楊婧,郭良輝

(中國(guó)地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院，北京 100083)

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

0 引言

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

當(dāng)前重力數(shù)據(jù)的曲化平方法較多，Pilkington等[6]將其分為基于源和基于場(chǎng)的兩大類。基于源的方法主要有等效源法[7-12]，它先對(duì)起伏觀測(cè)面異常反演獲得地下等效場(chǎng)源，然后再正演計(jì)算等效源在任意水平面引起的異常，該類方法通常計(jì)算量較大，適用于中、小規(guī)模數(shù)據(jù)的曲化平處理。基于場(chǎng)的方法有有限調(diào)和級(jí)數(shù)法[13]、泰勒級(jí)數(shù)法[14-15]、有限元法[16-17]、積分法[18-19]、逐次逼近法[20]等，該類方法無(wú)需反演場(chǎng)源，通常需要向上延拓或向下延拓到某中間層面進(jìn)而獲取最終平面的異常。向上延拓是穩(wěn)定的，但向下延拓通常是不穩(wěn)定的，這是由于它具有放大作用，對(duì)高頻噪聲干擾比較敏感，延拓跨度越大，高頻干擾越嚴(yán)重。因此，基于向下延拓的曲化平方法一般迭代收斂速度慢、精度有限，僅適用于延拓跨度小的曲化平。

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

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

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

1 方法原理

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

(1)

(2)

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

(3)

其中，T為起伏觀測(cè)面高程，Tmin和Tmax分別是起伏觀測(cè)面高程的最小值和最大值?？芍^測(cè)面起伏越大，修正因子S將越大，曲化平收斂速度將越快，反之，觀測(cè)面起伏越小，修正因子S將越小，曲化平收斂速度將越慢。一般，n取非負(fù)數(shù)，若n取值過(guò)大，則修正因子S將越大；反之，n取值越小，修正因子S將越?。划?dāng)n=0時(shí)，修正因子S=1，等效于徐世浙等[21]的算法。

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

2 理論模型數(shù)據(jù)試驗(yàn)

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

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

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

a、b—常規(guī)插值-迭代法曲化平結(jié)果及與海拔0 m平面理論重力異常殘差；c、d—本文改進(jìn)的插值-迭代法曲化平結(jié)果及與海拔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 不同曲化平方法的結(jié)果及與海拔0 m平面理論重力異常對(duì)比Fig.4 Comparison between the results of different interpolation-iteration methods and the gravity forward values of plane surface

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

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

改進(jìn)的插值—迭代法有效增強(qiáng)了信號(hào)，曲化平后其最小值為0.23 mGal，最大值為12.33 mGal，結(jié)果與海拔0 m平面理論重力異常差別較小。異常殘差最小值為-0.86 mGal，最大值為1.03 mGal，標(biāo)準(zhǔn)差為0.17 mGal，我們分析這一微小差異是由高程差距引起收斂速度不一致產(chǎn)生的。因此，與常規(guī)插值—迭代法相比，本文改進(jìn)的插值—迭代法曲化平結(jié)果更接近海拔0 m平面理論重力異常值，效果更優(yōu)，適用于大規(guī)模數(shù)據(jù)、復(fù)雜起伏觀測(cè)面、大延拓跨度的曲化平。

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

3 實(shí)際數(shù)據(jù)試驗(yàn)

川滇地區(qū)構(gòu)造上位于印支塊體、青藏高原塊體和揚(yáng)子塊體的交匯處，區(qū)域構(gòu)造動(dòng)力作用復(fù)雜，在緬甸弧東向俯沖、印度板塊藏東構(gòu)造結(jié)NE向楔體擠出和俯沖、青藏高原隆升和SE向擠出等共同作用下，該地區(qū)整體向E—SE方向擠出，形成了大陸內(nèi)部典型的剪切、拉張和推覆構(gòu)造，表現(xiàn)出地殼變形速率高、斷層運(yùn)動(dòng)劇烈、強(qiáng)震原地復(fù)發(fā)周期短等區(qū)域特征。因此，川滇地區(qū)是研究青藏高原隆升過(guò)程與大陸強(qiáng)震孕育機(jī)理的天然實(shí)驗(yàn)場(chǎng)，也是我國(guó)區(qū)域防震減災(zāi)的重要地區(qū)。高精度重力數(shù)據(jù)是川滇地區(qū)深部結(jié)構(gòu)與構(gòu)造研究的基礎(chǔ)數(shù)據(jù)。由于該地區(qū)地形上呈現(xiàn)青藏高原、云貴高原、四川盆地等復(fù)雜、多樣化的起伏形態(tài)，因此對(duì)該地區(qū)起伏地形面的重力異常數(shù)據(jù)作常規(guī)頻率域處理和解釋?xiě)?yīng)用之前，值得優(yōu)先作曲化平處理。

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

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

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

本文以常規(guī)頻率域垂直導(dǎo)數(shù)換算為例作說(shuō)明曲化平的意義。導(dǎo)數(shù)換算是重力異常分析解釋的一個(gè)常用步驟，用于突出淺部場(chǎng)源、分析構(gòu)造邊界等。圖8a和圖8b顯示了曲化平前后的布格重力異常的垂直導(dǎo)數(shù)，曲化平前的垂直導(dǎo)數(shù)數(shù)值范圍為-0.002 9～0.002 3 mGal/m，曲化平后的數(shù)值范圍轉(zhuǎn)為-0.004 4～0.003 7 mGal/m?？梢?jiàn)，曲化平后的異常幅值得到明顯增強(qiáng)，而且刻畫(huà)出更多、更清晰的異常變化細(xì)節(jié)，從而反映出更多的淺部場(chǎng)源細(xì)節(jié)和構(gòu)造邊界信息。

a—川滇地區(qū)地形高程；b—去噪后的布格重力異常；c—本文改進(jìn)插值—迭代法曲化平后布格重力異常；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 川滇地區(qū)曲化平前后重力異常對(duì)比Fig.6 The Bouguer gravity comparison before and after continuation from undulating surface to plane of Sichuan-Yunnan region

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

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

圖8中標(biāo)注了2008年汶川8級(jí)地震、2013年蘆山7級(jí)地震、2014年魯?shù)?.5級(jí)地震、2021年漾濞6.4級(jí)地震的震中位置，易見(jiàn)經(jīng)過(guò)曲化平這些地震震中周緣地區(qū)的重力異常變化特征更為鮮明，這有助于后續(xù)的深部結(jié)構(gòu)與構(gòu)造解釋研究。

4 結(jié)論與討論

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