于德武，龔勝平

中國地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所,河北 廊坊 065000

?

對迭代法位場向下延拓方法的剖析

于德武，龔勝平

中國地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所,河北 廊坊 065000

位場向下延拓迭代法的實質(zhì)內(nèi)容是“向上延拓而不是向下延拓”和以“迭代的結(jié)果趨近于觀測值”為標(biāo)準(zhǔn)的操作過程。根據(jù)數(shù)據(jù)操作流程剖析了位場向下延拓迭代法的運行機制，得到了迭代數(shù)據(jù)在空間域的變化規(guī)律，即用迭代法將觀測高度的位場向下延拓一個深度h。這實際上是通過不同高度的向上延拓來實現(xiàn)的。也就是說，迭代次數(shù)增加一次，涉及的上延平面就增大一個h的高度。一般地，迭代次數(shù)n與上延高度h的關(guān)系為n～(n+1)h。在空間域中，初值、上延結(jié)果、差以及每一次校正后的結(jié)果都能用滿足萊布尼茲定理的交錯級數(shù)表示，從而得出了迭代法能夠收斂的結(jié)論；或者，以“觀測高度上的實測值與計算值的差值小到可以忽略”為標(biāo)準(zhǔn)，從數(shù)學(xué)上也能證明迭代法能夠收斂。數(shù)學(xué)推論和模型試驗結(jié)果說明了迭代的位場初值可以任意給定。在實際操作中，迭代誤差標(biāo)準(zhǔn)的影響和由于迭代誤差標(biāo)準(zhǔn)不恰當(dāng)可能出現(xiàn)不能達(dá)到迭代標(biāo)準(zhǔn)的情況，需引起注意，也值得進(jìn)一步研究。

位場；向下延拓；收斂性；空間域；迭代法；頻率域；快速傅里葉變換

0 引言

通過快速傅里葉變換(FFT)實現(xiàn)位場向下延拓是在重、磁資料處理中常用的方法，但向下延拓的不穩(wěn)定性十分明顯。頻率域向下延拓一般只能向下延拓二三個資料點距。為了改善向下延拓的穩(wěn)定性，提高延拓距離，許多地球物理學(xué)家提出了種種方法，包括改進(jìn)的方法[1-10]。已有文獻(xiàn)[11-13]對此做過詳細(xì)介紹，在此不贅述。其中，徐世浙院士[7-8]對重磁位場向下延拓這一課題進(jìn)行了系統(tǒng)的研究，提出了全新的位場向下延拓方法。其主要步驟是：將水平觀測面上的位場實測值垂直投影至下部的延拓水平面上，作為該水平面上的位場初值；再根據(jù)該水平面上的初值，通過向上延拓計算觀測面上的位場值；最后用觀測面上的實測值與計算值的差值對延拓面上的位場值進(jìn)行校正。如此反復(fù)迭代，直至觀測面上的實測值與計算值的差值小到可以忽略。針對大數(shù)據(jù)量的網(wǎng)格數(shù)據(jù)，可將向下延拓深度提高到10～20倍點距。此后，不斷有人對該方法進(jìn)行補充完善，主要進(jìn)展有：

劉東甲等[11]把迭代下延法的空間域迭代過程變換到波數(shù)域(頻率域)，得到了位場向下延拓的波數(shù)域迭代法，并從波數(shù)域迭代公式出發(fā)證明了其收斂性，即當(dāng)?shù)螖?shù)n趨于無窮時，第n次迭代的位場波譜趨于向下延拓場的波譜。

駱遙[12]提出了一個下延異常逐次分離的過程，闡述了該過程即為迭代法下延的物理實質(zhì)。通過對可下延異常逐次進(jìn)行分離，將位場分離成可下延的區(qū)域場(信息)和剩余異常(噪聲)兩部分，并對分離出的區(qū)域場進(jìn)行向下延拓，其累計向下延拓結(jié)果的頻譜逼近于頻率域下延譜。剩余異常的頻譜經(jīng)不斷分離而被壓制，說明迭代過程本身是不斷收斂的。研究者認(rèn)為如何減少剩余異常的噪聲水平是迭代法下延中有待深入探討和解決的問題。

姚長利等[13]認(rèn)為迭代法是近來在研究中受到普遍重視的方法技術(shù)，雖然也取得了較好的成果，但也存在對迭代法研究還不夠深入、對其存在的缺點認(rèn)識不夠充分和客觀等問題。其推導(dǎo)了迭代法的通式，并分析了對迭代法收斂性影響的各種因素后指出：映射函數(shù)決定著迭代法是否收斂；增加迭代次數(shù)雖然能夠使計算收斂到FFT直接計算理論結(jié)果，但如果該理論結(jié)果本身就是不穩(wěn)定的，則迭代法計算如果收斂，也是收斂到一個不穩(wěn)定的結(jié)果。所以針對位場處理轉(zhuǎn)換中一些不穩(wěn)定計算采用迭代法，并沒有從根本上解決計算的不穩(wěn)定性問題。研究者[13]突出了“向下延拓迭代法主要過程是向上延拓而不是向下延拓”這一概念，筆者認(rèn)為這一概念代表了問題的本質(zhì)。

通過研究我們認(rèn)為，為了深刻理解方法的本質(zhì)，應(yīng)該對“向上延拓而不是向下延拓”做詳細(xì)的剖析，同時還應(yīng)該突出“迭代的結(jié)果趨近于觀測值”這一實際操作過程，并證明其收斂性。因此，研究迭代數(shù)據(jù)在空間域的變化規(guī)律以及迭代數(shù)據(jù)與觀測值之間的關(guān)系,就成了必要的事情。為了簡便起見，并不失一般性，在本文中，筆者以剖面磁異常向下延拓為例，按照向下延拓迭代法的步驟，通過追蹤迭代過程中數(shù)據(jù)的變化(研究迭代過程與觀測值之間的動態(tài)關(guān)系)，得到了迭代數(shù)據(jù)在空間域的變化規(guī)律;并從不同途徑證明空間域迭代法的收斂性，用模型試驗結(jié)果說明迭代的位場初值可以任意給定;最后,對實際操作中可能出現(xiàn)的、不能達(dá)到誤差標(biāo)準(zhǔn)的情況進(jìn)行了探討。以上工作將使得位場向下延拓迭代法理論更加完善。

1 迭代法位場向下延拓原理簡述

迭代法位場向下延拓原理和步驟[7]簡述如下。

高度ΓA與ΓB之間是無源空間(圖1)。用u(x，0)代表高度ΓA(z=0)的位。用

(1)

表示u(x，0)的傅里葉變換。其中，k是波數(shù)。則高度ΓB(z=h)上的位是

(2)

其中，F(xiàn)-1是傅里葉反變換。式(2)是頻率域的延拓公式。當(dāng)h>0 時，為向上延拓，延拓很穩(wěn)定；當(dāng)h<0時，為向下延拓，延拓很不穩(wěn)定。

據(jù)文獻(xiàn)[7]。圖1 觀測高度ΓB 與延拓高度ΓA的位置Fig.1 Location of observation height ΓB and continuation height ΓA

設(shè)ΓB高度的位uB是已知的，其下ΓA高度的位uA是待求的，ΓA，ΓB之間沒有場源。迭代法向下延拓的步驟如下：

1)將ΓB上的位uB垂直放置在ΓA的對應(yīng)點上，作為uA的初值。

2)從uA的初值，用公式(2)計算ΓB上的位，也可以用其他空間域的方法計算。

3)根據(jù)ΓB上原始值與計算值的差值，對ΓA上的uA進(jìn)行校正。

4)如此反復(fù)迭代，直至觀測高度上的實測值與計算值的差值小到可以忽略。一般迭代20～50次即可。

這種方法的位場向下延拓比簡單的FFT穩(wěn)定得多，可以向下延拓較大的距離。

2 向下延拓迭代法的數(shù)據(jù)流分析及收斂性

2.1 迭代法向下延拓的具體操作過程

為了分析向下延拓迭代法的收斂性，先分析其具體操作過程。

將ΓB上的位uB(或?qū)懽鱱B0)作為uA的初值，用公式(2)計算ΓB上的位。于是迭代法向下延拓在空間域的具體執(zhí)行 (數(shù)據(jù)追蹤) 情況如表1所示。表1中各符號的意義如圖2所示。

圖2 迭代過程中涉及的各個高度Fig.2 Altitude involved in iteration

一般地，設(shè)uA的初值為P，上延結(jié)果為Q，原始值與計算值的差為D，用原始值與計算值的差對計算值進(jìn)行校正后的值為C，則對于每一次迭代有

(3)

(4)

(5)

(6)

式中，n為迭代次數(shù)。

從上面的分析可知： 1)用迭代法將觀測高度的位場向下延拓一個深度h，實際上是通過不同高度的向上延拓來實現(xiàn)的，迭代次數(shù)增加一次，涉及的上延平面就增大一個h高度的平面。一般地，迭代次數(shù)n與上延高度h的關(guān)系為n～(n+1)h。 2)如果觀測高度上的實測值與計算值的差值小到可以忽略，則將本次迭代使用的初值作為向下延拓的結(jié)果(根據(jù)向下延拓的步驟 4))。

2.2 向下延拓迭代法的收斂性

根據(jù)公式(1)和公式(2)，有

(7)

也可以從另一方面論述迭代法的收斂性。因為迭代法是以“觀測高度上的實測值與計算值的差值小到可以忽略”為標(biāo)準(zhǔn)，從表1(或式(4))我們期望有

(8)

兩邊做傅里葉變換得

即期望有

(9)

而式(9)中公比為e-kh的交錯級數(shù)的無窮等比數(shù)列的和為

(10)

其中，k的取值范圍為〔0，∞〕。故也能推出上延結(jié)果的和Q≤uB的結(jié)論。

因此，當(dāng)?shù)螖?shù)趨近于無窮時, 原始值與計算值的差為

(11)

可以推論，對于有限的任意初值,每一次的上延結(jié)果是收斂的，且其和Q≤uB。

3 關(guān)于向下延拓迭代法的誤差標(biāo)準(zhǔn)和初值給定

根據(jù)迭代法步驟 4)，向下延拓迭代法的均方誤差可以寫成(用于模型試驗，主要與模型和計算精度有關(guān))

(12)

式中：uB,cont和uB,theo分別是ΓB高度上的位場延拓值和理論值；m是計算所用的點數(shù)?；蛘?，均方誤差也可以寫成(用于實測值，主要與實測值精度以及異常的復(fù)雜程度有關(guān))

(13)

式中：uB,surv是ΓB高度上的位場實測值。

對于式(12)，根據(jù)需要也可以在絕對誤差

(14)

的約束下用模型試驗來檢驗迭代法的收斂速度、結(jié)果的可信程度(例如最大值之比)等。顯然，式(14)較式(12)要嚴(yán)格得多。

圖3即是一個板狀磁性體模型利用不同的誤差標(biāo)準(zhǔn)得到的算例。板體的傾角為60°，磁化傾角為45°。板體半寬度為20 m，上頂埋深為200 m，點距為10 m，剖面長度為1 280 m。把磁性體上方200 m高度的磁異常向下延拓到磁性體上方50 m高度。試驗中所用的不同誤差標(biāo)準(zhǔn)和誤差公式對應(yīng)的迭代次數(shù)、計算異常最大值與理論異常最大值之百分比列于表2，迭代初值為觀測高度上的理論異常值(也試驗了迭代初值為常數(shù)或0值的情況，結(jié)果相同)。從表2看出，給定的誤差標(biāo)準(zhǔn)越嚴(yán)格，向下延拓結(jié)果與理論曲線越逼近，但迭代次數(shù)增多。

1.理論磁異常(h=200 m)；2.理論磁異常(h=50 m)；3.向下延拓磁異常(h=50 m，ε=1.0 nT)；4.向下延拓磁異常(h=50 m，ε=0.01 nT)。圖3 板體模型磁異常向下延拓算例Fig.3 Example of magnetic anomaly downward continuation by a dike model

圖4a為一個球狀磁性體模型按照式(12)給定的誤差標(biāo)準(zhǔn)為0.1 nT、迭代初值數(shù)據(jù)組為0得到的算例。球體半徑為100 m，中心深度為500 m，磁化傾角為50°，磁化偏角為0°，線距和點距均為50 m，計算了邊長為5 000 m 的方形面積的磁異常。圖4b為從平面圖中心部位截取的剖面示意圖(圖中縱橫距離單位未按比例)。把磁性體中心上方500 m高度平面的磁異常向下延拓到磁性體上方300 m高度的平面(即下延200 m)。圖4a中計算值與理論異常值符合得相當(dāng)好，因此曲線沒有用不同線條來表示。圖a邊部的0值等值線，是計算產(chǎn)生的噪聲。

表2 不同誤差標(biāo)準(zhǔn)、誤差公式對應(yīng)的迭代次數(shù)和最大值之百分比

Table 2 Iterate period and data difference with distinct error bounds and error formula

urms/nT按照式(12)nE/%按照式(14)nE/%1．0408089850．122491503940．01126297286698

注：E為計算異常最大值與理論異常最大值之百分比。

對于該球狀磁性體模型，利用不同的初值以及按照式(12)給定不同的誤差標(biāo)準(zhǔn)進(jìn)行試驗， 試驗結(jié)果列于表3。表3也說明，給定的誤差標(biāo)準(zhǔn)越嚴(yán)格，向下延拓結(jié)果與理論曲線越逼近，但迭代次數(shù)增多。

a.平面等值線圖 ； b.剖面示意圖。1.理論磁異常(h=500 m)；2.理論磁異常(h=300 m)；3.向下延拓磁異常(h=300 m，Urms=0.1 nT)。圖4 球體模型磁異常向下延拓算例Fig.4 Example of magnetic anomaly downward continuation by a spherical model

板狀體二維模型和球體三維模型算例還說明，迭代初值不一定必須用地面觀測值uB，可使用任意有限數(shù)據(jù)，包括常數(shù)(包括0值)，即從不同的初值出發(fā)，把磁性體上方的磁異常向下延拓到磁性體上方某個高度，都能得到收斂的結(jié)果：迭代法不受初值的影響。這與第3節(jié)的數(shù)學(xué)推論是一致的。

在用式(12)、式(13)或者式(14)給定誤差標(biāo)準(zhǔn)時，要考慮既不會因過于嚴(yán)格而造成迭代誤差標(biāo)準(zhǔn)得不到滿足，又不會因過于寬松而“欠迭代”從而造成向下延拓“不到位”的結(jié)果。顯然，在迭代法實際應(yīng)用中，誤差標(biāo)準(zhǔn)、實測值精度以及異常的復(fù)雜程度，都是影響迭代次數(shù)和向下延拓結(jié)果的重要因素。對于式(13)，定性地講：受偶然誤差影響的實測值可以用較小的誤差標(biāo)準(zhǔn)，受系統(tǒng)誤差影響的實測值就要用較大的誤差標(biāo)準(zhǔn)；對于簡單異?？梢杂幂^小的

表3 不同初值、不同誤差標(biāo)準(zhǔn)對應(yīng)的迭代次數(shù)和最大值之百分比

誤差標(biāo)準(zhǔn)，對于復(fù)雜異常就要用較大的誤差標(biāo)準(zhǔn)。

4 結(jié)論

1)迭代法向下延拓是通過不同高度的向上延拓來實現(xiàn)的，隨著迭代次數(shù)的增加，涉及的上延高度增大。一般地，迭代次數(shù)n與上延高度h的關(guān)系為n～(n+1)h。

2)迭代過程形成一系列交錯級數(shù)，每一次上延結(jié)果都是收斂的，且其和Q≤uB。模型試驗說明，隨著迭代次數(shù)的增加，迭代結(jié)果收斂于觀測高度上的實測值。

3)數(shù)學(xué)推論和模型試驗還說明迭代初值可使用任意有限數(shù)據(jù)，迭代法的收斂性與初值無關(guān)。

4)誤差標(biāo)準(zhǔn)關(guān)聯(lián)著迭代次數(shù)。誤差標(biāo)準(zhǔn)越嚴(yán)格，下延結(jié)果會更“到位”，而迭加次數(shù)也會增加。值得注意的是，實際中存在著達(dá)不到誤差標(biāo)準(zhǔn)的情況。

5)在迭代法實際應(yīng)用中，誤差標(biāo)準(zhǔn)、實測值精度以及異常的復(fù)雜程度，都是影響迭代次數(shù)和向下延拓結(jié)果的重要因素，這在實際應(yīng)用中需引起注意，也值得今后進(jìn)一步研究。

[1] Fedi M，F(xiàn)lorio G．A Stable Downward Continuation by Using the ISVD Method[J]．Geophysical Journal International，2002，151(1)：146-156.

[2] Cooper G．The Stable Downward Continuation of Potential Field Data[J]．Exploration Geophysics，2004，35(2)：260-265.

[3] 劉保華，焦湘恒，王重德．位場向下延拓的邊界單元法[J]．石油地球物理勘探，1990，25(4)：495-499. Liu Baohua，Jiao Xiangheng，Wang Chongde．Boundary Element Method for Continuation of Potential Field[J]. Oil Geophysical Prospecting, 1990，25(4)：495-499.

[4] 毛小平，吳蓉元，曲贊．頻率域位場下延的振蕩機制及消除方法[J]．石油地球物理勘探，1998，33(2)：230-237. Mao Xiaoping,Wu Rongyuan,Qu Zan. Oscillation Mechanism of Potential Field Downward Continuation in Frequency Domain and Its Elimination[J]. Oil Geophysical Prospecting, 1998，33(2)：230-237.

[5] 徐世浙，沈曉華，鄒樂君，等．將航磁異常從飛行高度向下延拓至地形線[J]．地球物理學(xué)報，2004，47(6)：1127-1130. Xu Shizhe，Shen Xiaohua，Zou Lejun，et al．Downward Continuation of Aeromagnetic Anomaly from Flying Altitude to Terrain[J]．Chinese Journal of Geophysics，2004，47(6)：1127-1130.

[6] 徐世浙．位場延拓的積分-迭代法[J]．地球物理學(xué)報，2006，49(4)：1176-1182. Xu Shizhe．The Integral-Iteration Method for Continuation of Potential Fields[J]．Chinese Journal of Geophysics，2006，49(4)：1176-1182.

[7] 徐世浙．迭代法與FFT法位場向下延拓效果的比較[J]．地球物理學(xué)報，2007，50(1)：285-289. Xu Shizhe．A Comparison of Effects Between the Iteration Method and FFT for Downward Continuation of Potential Fields[J]．Chinese Journal of Geophysics，2007，50(1)：285-289.

[8] 徐世浙，余海龍．位場曲化平的插值-迭代法[J]．地球物理學(xué)報，2007，50(6)：1811-1815. Xu Shizhe, Yu Hailong．The Interpolation-Iteration Method for Potential Field Continuation from Undulating Surface to Plane[J]．Chinese Journal of Geophysics，2007，50(6)：1811-1815.

[9] 王彥國，王祝文，張鳳旭，等. 位場向下延拓的導(dǎo)數(shù)迭代法[J].吉林大學(xué)學(xué)報：地球科學(xué)版，2012,42(1): 240-245. Wang Yanguo,Wang Zhuwen,Zhang Fengxu,et al. Derivative-Iteration Method for Downward Continuation of Potential Fields[J]. Journal of Jilin University：Earth Science Edition, 2012,42(1): 240-245.

[10] 張志厚，吳樂園. 位場向下延拓的相關(guān)系數(shù)法法[J].吉林大學(xué)學(xué)報：地球科學(xué)版，2012,42(6): 1912-1919. Zhang Zhihou,Wu Leyuan. Corelation Coefficient Method for Downward Continuation of Potential Fields[J]. Journal of Jilin University：Earth Science Edition, 2012,42(6): 1912-1919.

[11] 劉東甲，洪天求，賈志海，等．位場向下延拓的波數(shù)域迭代法及其收斂性[J]．地球物理學(xué)報，2009，52(6)：1599-1605. Liu Dongjia，Hong Tianqiu，Jia Zhihai，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.

[12] 駱遙．位場迭代法向下延拓的地球物理含義：以可下延異常逐次分離過程為例[J]．地球物理學(xué)進(jìn)展，2011，26(4)：1197-1200. Luo Yao．The Integral-Iteration Method for Continuation of Potential Fields: A Case by Computing Regional Anomaly Can be Continued Downward[J]．Progress in Geophysics，2011，26(4)：1197-1200.

[13] 姚長利，李宏偉，鄭元滿，等．重磁位場轉(zhuǎn)換計算中迭代法的綜合分析與研究[J]．地球物理學(xué)報，2012，55(6)：2062-2078. Yao Changli, Li Hongwei, Zheng Yuanman, et al. Research on Iteration Method Using in Potential Field Transformations[J]. Chinese Journal of Geophysics，2012，55(6)：2062-2078.

Analysis of the Potential Field Downward Continuation Iteration Method

Yu Dewu, Gong Shengping

InstituteofGeophysical&GeochemicalExploration,ChineseAcademyofGeologicalSciences,Langfang065000,Hebei,China

The substance of the potential field downward continuation iteration method is an operating process of “upward continuation rather than downward continuation” based on the criterion of “iteration result tend to be the observation data”. We have analyzed the operational mechanism of the method by chasing its operation routine, and obtained the dynamic rule of the data flow in spatial domain. It turned out that the iteration method for a downward continuation leads to an upward continuation of the potential field on the individual levels. The relationship between the iteration periodnand the heighthof the upward continuation isn～(n+1)h.Shown in an alternating serie of satisfying Leibniz theorem to each iteration result in spatial domain, the convergence of this method is obtained. The convergence can also be verified by using the criterion of minimizing the difference between measured value and calculated data. The mathematical reasoning and numerical computation results show that an initial value of the iteration method can be arbitrary. Also, we have studied the effect of the iteration tolerance criterion, and pointed out the possibility that a tolerance might never be satisfied in practice if the criterion is not fitted.

potential field；downward continuation；convergence；spatial domain；iteration method；frequency domain；fast Fourier transformation (FFT)

10.13278/j.cnki.jjuese.201503302.

2014-06-25

國家重大科學(xué)儀器設(shè)備開發(fā)專項(2011YQ050060，2011YQ05006011)

于德武(1955--)，男，教授級高級工程師，主要從事重磁勘探數(shù)據(jù)處理和解釋研究，E-mail:yu_dewu@aliyun.com。

10.13278/j.cnki.jjuese.201503302

P631.2

A

于德武，龔勝平. 對迭代法位場向下延拓方法的剖析.吉林大學(xué)學(xué)報:地球科學(xué)版,2015,45(3):934-940.

Yu Dewu, Gong Shengping. Analysis of the Potential Field Downward Continuation Iteration Method.Journal of Jilin University:Earth Science Edition,2015,45(3):934-940.doi:10.13278/j.cnki.jjuese.201503302.