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

(1. 東華理工大學 放射性地質與勘探技術國防重點學科實驗室，南昌 330013； 2. 廣東省核工業(yè)地質局 二九三大隊，廣州 510800)

0 引言

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

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

1 基本原理

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

UB=UAψ

(1)

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

1.1 積分迭代法

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

(2)

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

(1-φ)+1]=

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

(3)

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

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

(4)

(5)

1.3 補償向下延拓法

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

(6)

(7)

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

1.4 濾波特性

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

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

2 模型試算

為驗證積分迭代法、泰勒級數(shù)迭代法和補償向下延拓法的有效性，構建了由四個不同參數(shù)球體組成的疊加模型，各個模型體相對位置見圖2，模型體A、B、C、D中心坐標分別為(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，并在組合模型產生的重力數(shù)據(jù)上添加了0.5%隨機噪聲(圖2)，數(shù)據(jù)網(wǎng)格點距為0.1 km。從圖2可以看出，受疊加異常影響，模型體A、B上方的重力極值與模型體中心位置不對應，模型體D上方的重力異常則是以等值線扭曲為標志，難以從原始異常圖中精確識別出地質體的位置。為了提高異常分辨率，對原始含噪異常采用三種迭代法進行了向下延拓10倍、20倍點距計算。

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

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

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

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

3 應用實例

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

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

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

圖6 過地質體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 相山礦田地質構造綱要圖[18]Fig.7 Geological structure map of Xiangshan orefield

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

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

4 結論

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

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

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

參考文獻：

[1] 欒文貴. 位場解析延拓的穩(wěn)定化算法[J]. 地球物理學報，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]. 地球物理學報，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] 王延忠，熊光楚. 位場向下延拓組合濾波器的設計和應用[J]. 地球物理學報，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].地球物理學報，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]. 地球物理學報，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] 王彥國，王祝文，張鳳旭，等.位場向下延拓的導數(shù)迭代法[J].吉林大學學報(地球科學版), 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] 張志厚，吳樂園. 位場向下延拓的相關系數(shù)法 [J].吉林大學學報(地球科學版), 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] 高玉文，駱遙，文武. 補償向下延拓方法研究及應用 [J]. 地球物理學報，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]. 地球物理學報，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]. 地球物理學報，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]. 地球物理學報, 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]. 地球物理學報，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]. 石油地球物理學報，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] 侯重初. 補償圓滑濾波方法[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] 陳越. 相山鈾礦田地球物理特征及深部地質結構研究[D]：北京：核工業(yè)北京地質研究所, 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)