重力向下延拓的迭代法對比分析研究

吳姿穎， 鄧居智， 王彥國， 羅 瀟， 王 娥

(1. 東華理工大學(xué) 放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南昌 330013； 2. 廣東省核工業(yè)地質(zhì)局 二九三大隊(duì)，廣州 510800)

0 引言

向下延拓是位場數(shù)據(jù)處理及反演解釋的重要手段，它能夠突出淺部地質(zhì)異常，劃分水平疊加異常。但向下延拓是典型的不適定問題，直接使用會導(dǎo)致計(jì)算結(jié)果出現(xiàn)明顯的不穩(wěn)定現(xiàn)象，通常采用添加低通濾波器的方法來提高向下延拓穩(wěn)定性。在改善向下延拓穩(wěn)定性方法中，正則化法[1-2]、Wiener濾波法[3]、匹配組合濾波法[4]、泰勒級數(shù)展開法[5]及廣義逆算法[6]等在一定程度上提高了計(jì)算穩(wěn)定性，但向下延拓深度有限。徐世浙[7]提出的積分迭代法實(shí)現(xiàn)了大跨度向下延拓，之后各種不同模式的迭代法應(yīng)運(yùn)而生，如泰勒級數(shù)迭代法[8]、導(dǎo)數(shù)迭代法[9]、相關(guān)系數(shù)法[10]、補(bǔ)償延拓法[11]、迭代維納濾波法[12]。學(xué)者們還針對空間域積分迭代法計(jì)算效率低的問題，將積分迭代法引入到了波數(shù)域中，實(shí)現(xiàn)了快速計(jì)算[13]；針對積分迭代法壓制高頻干擾不足問題還提出了一系列改進(jìn)措施[14-15]；針對泰勒級數(shù)迭代法迭代次數(shù)少的優(yōu)點(diǎn)和高頻干擾壓制能力不足問題，提出了正則-積分迭代法[16]。

雖然迭代法在解決向下延拓上的良好應(yīng)用效果已被普遍認(rèn)同，但不同迭代模式之間關(guān)聯(lián)性及差異性，尤其不同迭代模式的應(yīng)用效果，需要進(jìn)行深入研究與分析，以便選擇更為適合的迭代模式進(jìn)行有效的向下延拓?；诖耍P者對常用的積分迭代法、泰勒級數(shù)迭代法和補(bǔ)償向下延拓法這三種方法進(jìn)行了理論模擬和應(yīng)用實(shí)例的對比研究，總結(jié)典型迭代法的適用性。

1 基本原理

在波數(shù)域中，向下延拓可以表示為式(1)。

UB=UAψ

(1)

其中:UA、UB分別為觀測面和延拓面的位場波譜；ψ=exp(wz)為理論向下延拓算子；z為向下延拓深度；w為圓波數(shù)。由于向下延拓是不適定問題，直接使用理論延拓因子會導(dǎo)致計(jì)算不穩(wěn)定性。在解決這一問題上，目前普遍認(rèn)為迭代法具有較好的應(yīng)用效果。這里選取了常用的積分迭代法、泰勒級數(shù)迭代法和補(bǔ)償向下延拓法三種方法進(jìn)行對比研究分析，以評判不同迭代模式的應(yīng)用效果。

1.1 積分迭代法

積分迭代法是由徐世浙院士[7]在空間域中提出的，劉東甲等[13]將其引入到了波數(shù)域中，并嚴(yán)格證明了迭代法的收斂性。該方法是將原平面位場值作為延拓面初值，再將原平面位場值與延拓面上初值向上延拓到原平面結(jié)果之差作為剩余值對延拓面位場值進(jìn)行修正，采用逐次循環(huán)迭代實(shí)現(xiàn)的，該方法的波數(shù)域迭代遞推公式為：

(2)

其中:φ=exp(-wz)為向上延拓算子，利用數(shù)學(xué)歸納法可得到積分迭代法的通式為式(3)。

(1-φ)+1]=

UAψ[1-(1-φ)n+1]

(3)

1.2 泰勒級數(shù)迭代法

泰勒級數(shù)迭代法是將波數(shù)域中的向下延拓算子進(jìn)行泰勒級數(shù)展開，取前N項(xiàng)和代替向下延拓理論算子進(jìn)行近似向下延拓，并采用迭代法進(jìn)行逐步逼近得到的，其迭代公式為式(4)。

(4)

(5)

1.3 補(bǔ)償向下延拓法

補(bǔ)償法是基于侯重初[17]的補(bǔ)償圓滑濾波和陳生昌等[6]提出的廣義逆算法推導(dǎo)出來的，其補(bǔ)償向下延拓的迭代公式為式(6)。

(6)

(7)

從式(3)、式(5)、式(7)可以看出，迭代事實(shí)上是在理論延拓算子上添加了一個與迭代次數(shù)有關(guān)的低通濾波器，只是不同迭代法出發(fā)點(diǎn)不同。另外，積分迭代法、泰勒級數(shù)迭代法和補(bǔ)償向下延拓法的迭代通式具有很強(qiáng)的相似性，不同之處在于通式中向上延拓因子φ之前所使用的初始濾波算子不同，積分迭代法是“1”，為全通濾波器；泰勒級數(shù)迭代法的是φ，為高通濾波器；而補(bǔ)償法的是HD，為中低通濾波器。

1.4 濾波特性

圖1是積分迭代法、泰勒級數(shù)迭代法和補(bǔ)償向下延拓法這三種方法的濾波特性曲線，可以看出，理論下延算子隨著波數(shù)的增加呈現(xiàn)指數(shù)型放大，而三種迭代法的濾波因子對高頻成分均有一定的壓制，但積分迭代法和泰勒級數(shù)迭代法的響應(yīng)算子仍是高通濾波，對高頻干擾仍有放大作用。補(bǔ)償法的濾波因子則在中低頻段更能較好地逼近理論延拓因子，在高頻段則更有效地壓制高頻成分，即補(bǔ)償法在低頻段能更精確地逼近理論下延因子而對高頻成分則更有效地壓制。因此，補(bǔ)償法在理論上具有更高的計(jì)算精度和計(jì)算穩(wěn)定性，筆者認(rèn)為根本原因在于迭代法使用的初始濾波器不同。

圖1 不同迭代法的濾波特性Fig.1 Filter characteristics of different iterative methods

2 模型試算

為驗(yàn)證積分迭代法、泰勒級數(shù)迭代法和補(bǔ)償向下延拓法的有效性，構(gòu)建了由四個不同參數(shù)球體組成的疊加模型，各個模型體相對位置見圖2，模型體A、B、C、D中心坐標(biāo)分別為(10，10，1.5)、(20，10，6)、(10，20，6)和(20，20，4)，半徑分別為0.5 km、2.5 km、2 km、1 km，四個模型體剩余密度均為0.5 g/cm3，并在組合模型產(chǎn)生的重力數(shù)據(jù)上添加了0.5%隨機(jī)噪聲(圖2)，數(shù)據(jù)網(wǎng)格點(diǎn)距為0.1 km。從圖2可以看出，受疊加異常影響，模型體A、B上方的重力極值與模型體中心位置不對應(yīng)，模型體D上方的重力異常則是以等值線扭曲為標(biāo)志，難以從原始異常圖中精確識別出地質(zhì)體的位置。為了提高異常分辨率，對原始含噪異常采用三種迭代法進(jìn)行了向下延拓10倍、20倍點(diǎn)距計(jì)算。

圖2 含0.5%噪聲的疊加重力異常Fig.2 Combined gravity anomaly including 0.5% random noise

圖3是向下延拓10倍時三種迭代法均方誤差與迭代次數(shù)的關(guān)系曲線，可以看出，泰勒級數(shù)迭代法迭代誤差隨著迭代次數(shù)的增加而快速遞增，最小誤差0.42 mGal在n=1時出現(xiàn)；積分迭代法和補(bǔ)償法迭代誤差則是隨著迭代次數(shù)的增加先減少后增加，存在一個誤差極小值，分別為0.10 mGal和0.07 mGal，均顯著小于泰勒級數(shù)迭代法的。圖4是三種迭代法在迭代誤差最小時對應(yīng)的向下延拓10倍的計(jì)算結(jié)果，可以看出，補(bǔ)償法無論是在計(jì)算穩(wěn)定性上還是在計(jì)算精度上均高于積分迭代法和泰勒級數(shù)迭代法，而泰勒級數(shù)迭代法的穩(wěn)定性最差，噪聲干擾幾乎淹沒了地質(zhì)體A、B、D上方的有效異常。

圖3 不同迭代法向下延拓10倍點(diǎn)距誤差 與迭代次數(shù)關(guān)系曲線Fig.3 Relationship curves of iterative errors and iterative numbers using different iterative methods for downward continuation of 1.0 km

圖5是積分迭代法、泰勒級數(shù)迭代法和補(bǔ)償向下延拓法選取迭代次數(shù)分別為15、1、2 000時的向下延拓20倍點(diǎn)距的計(jì)算結(jié)果，由于向下延拓深度已超過地質(zhì)體A的中心埋深，因此理論延拓值無法與實(shí)際延拓值進(jìn)行對比，即無法獲得不同迭代方法的誤差曲線，此時迭代次數(shù)選取是以計(jì)算結(jié)果圖中開始出現(xiàn)不穩(wěn)定現(xiàn)象為依據(jù)。從圖5中同樣可以看出，補(bǔ)償法同樣具有更高的計(jì)算精度和計(jì)算穩(wěn)定性，另外從圖5(c)中還可以看出，當(dāng)延拓深度大于地質(zhì)體的實(shí)際深度時，地質(zhì)體上方的極大值圈閉外圍存在明顯的極小負(fù)值環(huán)狀圈閉。為清晰說明這種現(xiàn)象，圖6給出了過地質(zhì)體A、C中心位置剖面的重力及向下延拓重力異常，可以看出向下延拓較好地提高了對地質(zhì)體A的異常識別能力，向下延拓20倍點(diǎn)距的重力異常在地質(zhì)體A上方的重力正異常外圍出現(xiàn)了邊瓣負(fù)異常，而邊瓣負(fù)異常外圍則又出現(xiàn)了邊瓣正異常(原始異常和向下延拓10倍點(diǎn)距的則無此異常特征)，因此在實(shí)際資料處理中可以此異常特征來評價地質(zhì)體的大致埋深。

3 應(yīng)用實(shí)例

為驗(yàn)證積分迭代法、泰勒級數(shù)迭代法和補(bǔ)償法這三種方法的實(shí)用性，選取江西相山鈾礦田重力數(shù)據(jù)進(jìn)行方法對比試驗(yàn)。相山鈾多金屬礦田是我國最大的火山巖型鈾礦田，位于贛-杭火山巖型鈾多金屬成礦帶與大王山-于山花崗巖型鈾多金屬成礦帶交匯區(qū)，其經(jīng)歷了晉寧期、加里東期、海西-印支期造山作用，燕山期時處于NEE向贛杭構(gòu)造火山巖帶西南端與近NS向贛中南花崗巖帶的交接地帶，發(fā)生了強(qiáng)烈的構(gòu)造-巖漿-成礦活動。相山盆地?cái)嗔褬?gòu)造(圖7)較為復(fù)雜，主構(gòu)造方向?yàn)镹E向，而NS向、近EW向、NW向及NEE向斷裂同樣較為發(fā)育。相山盆地中心部位地表巖性比較單一，主要以低密度的碎斑熔巖為主；盆地基底為高密度的青白口紀(jì)變質(zhì)巖，盆地外圍四周均有出露；高密度的晚侏羅世砂巖砂礫巖呈現(xiàn)出明顯的環(huán)形條帶狀分布。

圖4 不同迭代法向下延拓10倍點(diǎn)距的計(jì)算結(jié)果Fig.4 Results of downward continuation of 1.0 km using different iterative methods(a)理論向下延拓；(b)積分迭代法；(c)泰勒級數(shù)迭代法；(d)補(bǔ)償向下延拓

圖5 不同迭代法向下延拓20倍點(diǎn)距的計(jì)算結(jié)果Fig.5 Results of downward continuation of 2.0 km using different iterative methods(a)積分迭代法；(b)泰勒級數(shù)迭代法；(c)補(bǔ)償向下延拓

圖6 過地質(zhì)體A、C中心位置剖面上的重力及向下延拓重力異常Fig.6 Original gravity and downward continuation gravity anomalies of the profile over the center position of geological body A and C

圖7 相山礦田地質(zhì)構(gòu)造綱要圖[18]Fig.7 Geological structure map of Xiangshan orefield

圖8(a)是網(wǎng)格點(diǎn)距為250 m的相山布格重力異常，由圖8(a)可以看出，布格重力異常整體上呈現(xiàn)為中心重力低、外圍重力高，與低密度的碎斑熔巖及高密度的變質(zhì)巖表現(xiàn)出了良好的對應(yīng)關(guān)系，但重力異常中局部細(xì)節(jié)特征受區(qū)域場影響較大而難以清晰識別。為了提高異常分辨能力，這里采用了積分迭代法、泰勒級數(shù)迭代法及補(bǔ)償向下延拓法分別迭代100次、5次和200次向下延拓5倍點(diǎn)距的計(jì)算結(jié)果分別見圖8(b)～圖8(d)，可以看出，積分迭代法和泰勒級數(shù)迭代法存在著明顯的不穩(wěn)定現(xiàn)象，僅能識別部分有效異常，而補(bǔ)償法獲得的計(jì)算結(jié)果顯然穩(wěn)定性較好，并未出現(xiàn)明顯的振蕩現(xiàn)象，同時，異常分辨率得到了明顯提高，局部異常間的特征關(guān)系得到了清晰地展示。異常整體上呈現(xiàn)出NE走向，局部異常伴隨著NS走向、EW向及NW向，這也恰與該地區(qū)的構(gòu)造特征相一致。另外，向下延拓1.25 km的重力異常在研究區(qū)大部分區(qū)域均表現(xiàn)為明顯的正、負(fù)相間出現(xiàn)，表明大多淺部地質(zhì)體中心埋深則是小于這個延拓深度的。

圖8 相山礦田重力數(shù)據(jù)迭代法向下延拓5倍點(diǎn)距的計(jì)算結(jié)果Fig.8 Downward continuation 1.25 km height of gravity data in Xiangshan orefield using different iterative methods(a)相山重力異常；(b)積分迭代法；(c)泰勒級數(shù)迭代法；(d)補(bǔ)償向下延拓

4 結(jié)論

1)積分迭代法和泰勒級數(shù)迭代法仍是高通濾波，在實(shí)際資料處理中，需要對原始數(shù)據(jù)進(jìn)行濾噪后再使用。補(bǔ)償向下延拓法為中低通濾波，對高頻成分具有較強(qiáng)的壓制能力，在實(shí)際資料處理中直接使用也可以獲得穩(wěn)定的計(jì)算結(jié)果。

2)向下延拓可以直接用于估計(jì)地質(zhì)體的大致埋深，可以清晰展示出異常間的相互關(guān)系和異常等值線圈閉走向，這些異常特征可為地質(zhì)體平面賦存狀態(tài)、邊界位置及深度信息等方面提供了借鑒。

3)構(gòu)造壓制高頻干擾能力強(qiáng)于廣義逆算子的低通濾波器，結(jié)合迭代法則在理論上可以獲得更為理想的向下延拓結(jié)果。

參考文獻(xiàn)：

[1] 欒文貴. 位場解析延拓的穩(wěn)定化算法[J]. 地球物理學(xué)報，1983, 26(3): 263-274.

LUAN W G. The stabilized algorithm of the analytic continuation for the potential filed[J]. Chinese Journal of Geophysics, 1983, 26(3): 263-274.(In Chinese)

[2] 梁錦文. 位場向下延拓的正則化方法[J]. 地球物理學(xué)報，1989, 32(5): 600-608.

LIANG J W. Downward continuation of regularization methods for potential fields[J].Chinese Journal of Geophysics, 1989, 32(5): 600-608. (In Chinese)

[3] PAWLOWSKI R S.Preferential continuation for potential-field anomaly enhancement[J]. Geophysics,1995,60(2):390-398.

[4] 王延忠，熊光楚. 位場向下延拓組合濾波器的設(shè)計(jì)和應(yīng)用[J]. 地球物理學(xué)報，1985, 28(5): 537-543.

WANG Y Z, XIONG G C. A combined filter used for the downward continuation of the potential filed[J].Chinese Journal of Geophysics, 1985, 28(5): 537-543. (In Chinese)

[5] FEDI M, FLORIO G. A stable downward continuation by using the ISVD method[J].Geophysical Journal International, 2002, 151(1): 146-156.

[6] 陳生昌，肖鵬飛. 位場向下延拓的波數(shù)域廣義逆算法[J].地球物理學(xué)報，2007, 50(6): 1816-1822.

CHEN S C, XIAO P F. Wavenumber domain generalized inverse algorithm for potential field downward continuation[J].Chinese Journal of Geophysics, 2007, 50(6): 1816-1822. (In Chinese)

[7] 徐世浙. 位場延拓的積分-迭代法[J]. 地球物理學(xué)報，2006, 49(4): 1176-1182.

XU S Z. The integral-iteration method for continuation of potential fields[J]. Chinese Journal of Geophysics, 2006, 49(4): 1176-1182. (In Chinese)

[8] 王彥國，張鳳旭，王祝文，等.位場向下延拓的泰勒級數(shù)迭代法[J]. 石油地球物理勘探，2011, 46(4): 657-662.

WANG Y G, ZHANG F X, WANG Z W, et al. Taylor series iteration for downward continuation of potential field [J]. Oil Geophysical prospecting, 2011, 46(4): 657-662.(In Chinese)

[9] 王彥國，王祝文，張鳳旭，等.位場向下延拓的導(dǎo)數(shù)迭代法[J].吉林大學(xué)學(xué)報(地球科學(xué)版), 2012, 42(1): 240-245.

WANG Y G, WANG Z W, ZHANG F X, et al. Derivative-iteration method for downward continuation of potential fields[J]. Journal of Jilin University (Earth science Edition), 2012, 42(1): 240-245. (In Chinese)

[10] 張志厚，吳樂園. 位場向下延拓的相關(guān)系數(shù)法 [J].吉林大學(xué)學(xué)報(地球科學(xué)版), 2012, 42(6): 1912-1919.

ZHANG Z H, WU L Y. Correlation coefficient method for downward continuation of potential field [J]. Journal of Jilin University (Earth Science Edition), 2012, 42(6): 1912-1919. (In Chinese)

[11] 高玉文，駱遙，文武. 補(bǔ)償向下延拓方法研究及應(yīng)用 [J]. 地球物理學(xué)報，2012, 55(8): 2747-2756.

GAO Y W, LUO Y, WEN W. The compensation method for downward continuation of potential field from horizontal plane and its application [J]. Chinese Journal of Geophysics, 2012, 55(8): 2747-2756. (In Chinese)

[12] 曾小牛，劉代志，李夕海，等.位場向下延拓的改進(jìn)迭代維納濾波法 [J]. 地球物理學(xué)報，2014, 57(6): 1958-1967.

ZENG X N, LIU D Z, LI X H, et al.An improved iterative Wiener filter for downward continuation of potential fields[J]. Chinese Journal of Geophysics, 2014, 57(6): 1958-1967. (In Chinese)

[13] 劉東甲，洪天求，賈志海，等. 位場向下延拓的波數(shù)域迭代法及其收斂性[J]. 地球物理學(xué)報，2009, 52(6): 1599-1605.

LIU D J, HONG T Q, JIA Z H, et al. Wave number domain iteration method for downward continuation of potential fields and its convergence [J]. Chinese Journal of Geophysics, 2009, 52(6): 1599-1605. (In Chinese)

[14] 張輝，陳龍偉，任治新，等. 位場向下延拓迭代法收斂性分析及穩(wěn)健向下延拓方法研究 [J]. 地球物理學(xué)報, 2009, 52(4): 1107-1113.

ZHANG H, CHEN L W, REN Z X, et al. Analysis on convergence of iteration method for potential fields downward continuation and research on robust downward continuation method [J]. Chinese Journal of Geophysics, 2009, 52(4): 1107-1113. (In Chinese)

[15] 于波，翟國君，劉雁春，等. 噪聲對磁場向下延拓迭代法的計(jì)算誤差影響分析[J]. 地球物理學(xué)報，2009, 52(8): 2182-2188.

YU B, ZHAI G J, LIU Y C, et al. Analysis of noise effect on the calculation error of downward continuation with iteration method [J]. Chinese Journal of Geophysics, 2009, 52(8): 2182-2188. (In Chinese)

[16] 曾小牛，李夕海，牛超，等. 位場向下延拓的波數(shù)域正則-積分迭代法[J]. 石油地球物理學(xué)報，2013, 48(4): 643-650.

ZENG X N, LI X H, NIU C, et al. Regularization-integral iteration in wavenumber domain for downward continuation of potential fields [J]. Oil Geophysical Prospecting, 2013, 48(4): 643-650. (In Chinese)

[17] 侯重初. 補(bǔ)償圓滑濾波方法[J]. 石油物探，1981，20(2): 22-29.

HOU Z C. Filtering of smooth compensation [J]. Geophysical prospecting for Petroleum, 1981, 20(2): 22-29. (In Chinese)

[18] 陳越. 相山鈾礦田地球物理特征及深部地質(zhì)結(jié)構(gòu)研究[D]：北京：核工業(yè)北京地質(zhì)研究所, 2014.

CHEN Y. Study on geophysical characteristics and deep geological structures of the Xiangshan uranium orefield. [D].Beijing: thesis of Beijing Research insititute of Uranium Geology，2014. (In Chinese)