位場(chǎng)垂向高階導(dǎo)數(shù)的Tikhonov正則化迭代法

杜 威，許家姝，吳燕岡，郝夢(mèng)成

吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

0 引言

在位場(chǎng)數(shù)據(jù)的處理中，垂向?qū)?shù)可以壓制區(qū)域場(chǎng)、圈定局部場(chǎng)[1]、分離水平疊加異常。并且在位場(chǎng)其他處理方法中，垂向?qū)?shù)的求取也至關(guān)重要，例如歐拉反褶積、向下延拓[2]、重力歸一化梯度等常規(guī)位場(chǎng)數(shù)據(jù)處理方法中都需要先進(jìn)行垂向?qū)?shù)的計(jì)算。可見(jiàn)導(dǎo)數(shù)換算作為處理過(guò)程的重要組成部分，其結(jié)果的好壞直接影響到這些位場(chǎng)數(shù)據(jù)處理方法的計(jì)算精度。位場(chǎng)垂向?qū)?shù)的換算方法可以分為：空間域法和波數(shù)域法。空間域法包括有限元法[3]、樣條函數(shù)法[4-6]等；波數(shù)域法包括常規(guī)FFT算子、維納濾波法[7]、補(bǔ)償圓滑濾波法[8]、正則化方法[9]、新正則化方法[10]及波數(shù)域迭代法[11]。其他方法包括Fedi等[12]提出的ISDV算法，該方法在空間域與波數(shù)域相結(jié)合換算出垂向?qū)?shù)。目前，ISVD算法已經(jīng)成為穩(wěn)定求解位場(chǎng)各階垂向?qū)?shù)的重要方法，并廣泛應(yīng)用到位場(chǎng)的數(shù)據(jù)處理中[13-15]。

由于空間域的算法普遍存在計(jì)算速度慢等缺點(diǎn)，因此目前常在頻率域?qū)ξ粓?chǎng)的高階垂向?qū)?shù)進(jìn)行計(jì)算；但是頻率域的導(dǎo)數(shù)響應(yīng)因子具有高頻放大作用，當(dāng)所求導(dǎo)數(shù)的階數(shù)越高，數(shù)據(jù)對(duì)噪聲的敏感程度就越大，導(dǎo)致所求結(jié)果不穩(wěn)定。為了壓制高階導(dǎo)數(shù)計(jì)算結(jié)果的不穩(wěn)定性，本文提出一種基于Tikhonov正則化迭代法計(jì)算位場(chǎng)各階垂向?qū)?shù)的方法，并通過(guò)模型試驗(yàn)及實(shí)際數(shù)據(jù)的應(yīng)用，對(duì)該方法計(jì)算垂向?qū)?shù)的穩(wěn)定性進(jìn)行了驗(yàn)證。

1 方法原理

位場(chǎng)T(x,y)與其各階垂向?qū)?shù)Dm(x,y)(m代表階數(shù))在波數(shù)域的關(guān)系式為

(1)

針對(duì)式(1)的不穩(wěn)定問(wèn)題，Tikhonov正則化是一種廣泛應(yīng)用的方法，即求解一個(gè)極小化的正則化泛函：

(2)

其中：‖‖2表示L2范數(shù)；α為正則化參數(shù)，用于平衡不穩(wěn)定性及光滑性。上述極小化問(wèn)題的解為

(3)

(4)

(5)

(6)

(7)

由式(1)可知：

(8)

代入式(7)得：

(9)

由數(shù)學(xué)歸納法很容易得到迭代通式：

(10)

(11)

當(dāng)?shù)鸁o(wú)窮大時(shí)，

(12)

上式表明，本文提出的Tikhonov正則化迭代法進(jìn)行高階導(dǎo)數(shù)換算是可以收斂到理論解的。其中正則化參數(shù)α的常用取值范圍為(0.001，1) ，迭代次數(shù)在10以內(nèi)即可。

2 濾波特性分析

圖1和圖2分別是垂向一階導(dǎo)數(shù)和垂向二階導(dǎo)數(shù)的Tikhonov正則化迭代算子濾波特性。其中，圖1a和圖2a分別是不同迭代次數(shù)的垂向一階導(dǎo)數(shù)和垂向二階導(dǎo)數(shù)Tikhonov正則化迭代算子的濾波響應(yīng)特征曲線。其中，橫坐標(biāo)ω為徑向圓頻率，數(shù)據(jù)的點(diǎn)距均為1 km，垂向一階導(dǎo)數(shù)的正則化參數(shù)α的取值為1.0，垂向二階導(dǎo)數(shù)的正則化參數(shù)α的取值為0.1。對(duì)理論垂向?qū)?shù)算子φm、Tikhonov正則化迭代法垂向?qū)?shù)換算算子的濾波特性進(jìn)行對(duì)比分析可得：理論垂向?qū)?shù)因子φm隨著頻率的增加急劇上升，且導(dǎo)數(shù)階次越高，高頻成分被放大得越強(qiáng)；Tikhonov正則化迭代算子在低頻段逼近理論垂向?qū)?shù)算子，保證了低頻有用信號(hào)處導(dǎo)數(shù)的換算，在高頻階段卻能對(duì)其進(jìn)行壓制；隨著迭代次數(shù)的增加，Tikhonov正則化迭代算子逐漸趨近于理論垂向?qū)?shù)算子。為了分析正則化參數(shù)α的取值對(duì)Tikhonov正則化迭代算子濾波特性的影響，圖1b和圖2b分別為迭代3次時(shí)，垂向一階和二階導(dǎo)數(shù)Tikhonov正則化迭代算子的濾波響應(yīng)特征曲線。可以看出：當(dāng)?shù)螖?shù)不變時(shí)，隨著導(dǎo)數(shù)階次的增加，Tikhonov正則化迭代算子對(duì)高頻成分的壓制作用逐漸增強(qiáng)，保證了高階導(dǎo)數(shù)換算的穩(wěn)定性；正則化參數(shù)α取值越小，Tikhonov正則化迭代算子在高頻段的壓制作用越弱，當(dāng)其值為0時(shí)，則結(jié)果為理論導(dǎo)數(shù)算子；當(dāng)正則化參數(shù)α的取值過(guò)大時(shí)，對(duì)低頻段也產(chǎn)生一定的壓制作用，則會(huì)影響濾波效果。

a.α=1.0，不同迭代次數(shù)；b. 迭代3次，不同α取值。ω為徑向圓頻率。圖1 垂向一階導(dǎo)數(shù)的Tikhonov正則化迭代算子濾波特性Fig.1 Filtering properties of vertical derivative with Tikhonov regularization iteration factor

a.α=0.1，不同迭代次數(shù)；b. 迭代3次，不同α取值。圖2 垂向二階導(dǎo)數(shù)的Tikhonov正則化迭代算子濾波特性Fig.2 Filtering properties of second vertical derivative with Tikhonov regularization iteration factor

可見(jiàn)，本文提出的Tikhonov正則化迭代算子在低頻段能精確地逼近理論導(dǎo)數(shù)算子，對(duì)高頻成分則能有效壓制，因而該方法在理論上具有保幅性和穩(wěn)定性的優(yōu)點(diǎn)。

3 模型試驗(yàn)

文中采用2個(gè)二度直立柱體的疊加異常進(jìn)行數(shù)值檢驗(yàn)。模型體1和2長(zhǎng)均為2.0 km，上頂埋深h為2.0 km，下底埋深h為2.5 km;其中心位置x分別為10、15 km，點(diǎn)距為0.1 km，2個(gè)模型體的剩余密度均為1.0 g/cm3。組合模型體相對(duì)位置及其產(chǎn)生的重力異常Δg如圖3a所示，圖3b和圖3c分別是無(wú)噪聲情況下理論的垂向二階導(dǎo)數(shù)和垂向三階導(dǎo)數(shù)。添加重力異常幅值0.1%的隨機(jī)噪聲，圖4和圖5分別是含噪聲情況下常規(guī)FFT法和Tikhonov正則化迭代法所計(jì)算出的垂向二階以及垂向三階導(dǎo)數(shù)。其中，二階垂向?qū)?shù)的正則化參數(shù)為0.1，迭代8次；垂向三階導(dǎo)數(shù)的正則化參數(shù)為0.5，迭代8次。

圖3 模型重力異常及模型位置(a)、垂向二階導(dǎo)數(shù)(b) 和垂向三階導(dǎo)數(shù)(c)Fig.3 Gravity anomaly and the position of the model (a), the second vertical derivative (b) and the third vertical derivative (c)

圖4 FFT垂向二階導(dǎo)數(shù)(a)和Tikhonov正則化迭代法垂向二階導(dǎo)數(shù)(b)Fig.4 Second vertical derivative using FFT method (a) and the second vertical derivative using Tikhonov regularization iteration method (b)

圖5 常規(guī)FFT垂向三階導(dǎo)數(shù)(a)和Tikhonov正則化迭代法垂向三階導(dǎo)數(shù)(b)Fig.5 Third vertical derivative using FFT method (a) and the third vertical derivative using Tikhonov regularization iteration method (b)

圖4和圖5分別是常規(guī)FFT法和Tikhonov正則化迭代法計(jì)算垂向二階和三階導(dǎo)數(shù)的對(duì)比圖。從圖中可以看出，由于噪聲的干擾，用常規(guī)FFT變換計(jì)算的垂向?qū)?shù)結(jié)果存在著明顯的波動(dòng)性，且導(dǎo)數(shù)的階次越高，穩(wěn)定性越差，無(wú)法識(shí)別有效異常，這是由于理論導(dǎo)數(shù)算子對(duì)高頻成分的放大作用。Tikhonov正則化迭代法計(jì)算垂向二階及三階導(dǎo)數(shù)結(jié)果穩(wěn)定，對(duì)噪聲的壓制能力強(qiáng)，這是因?yàn)檎齽t化迭代法的濾波特性，在低頻段較好地逼近理論因子，而在高頻段可以較好地壓制噪聲。對(duì)比圖3b和3c中不含噪聲計(jì)算得到的垂向二階導(dǎo)數(shù)和三階導(dǎo)數(shù)可以看出，隨著求導(dǎo)次數(shù)的增加，高頻部分的影響越來(lái)越大，即圖3c中曲線的震蕩幅度大于圖3b中。從圖4b和圖5b可以看出，Tikhonov正則化迭代法計(jì)算得到的導(dǎo)數(shù)，在數(shù)量級(jí)上和換算得到的導(dǎo)數(shù)(圖3b和3c)數(shù)量級(jí)是相同的，垂向三階導(dǎo)數(shù)的結(jié)果圖中并沒(méi)有出現(xiàn)明顯的震蕩，再次說(shuō)明了Tikhonov正則化迭代法的準(zhǔn)確性和穩(wěn)定性。圖4b模型體1和2中，零值點(diǎn)對(duì)應(yīng)的橫坐標(biāo)分別為8.5和11.3、13.7和16.5，圖5b模型體1和2中，零值點(diǎn)對(duì)應(yīng)的橫坐標(biāo)分別為8.6和11.2、13.8和16.4，可以看出，隨著所計(jì)算導(dǎo)數(shù)的階次增加，曲線的零值點(diǎn)位置與地質(zhì)體邊界的對(duì)應(yīng)效果越來(lái)越好。

4 實(shí)際應(yīng)用

為了試驗(yàn)Tikhonov正則化迭代法求垂向?qū)?shù)對(duì)實(shí)際資料的處理效果，我們對(duì)老撾萬(wàn)象盆地的布格重力異常Δg數(shù)據(jù)進(jìn)行垂向二階導(dǎo)數(shù)的計(jì)算，其中，實(shí)際數(shù)據(jù)的點(diǎn)距為0.2 km。該研究區(qū)的地質(zhì)情況復(fù)雜，如圖6。圖7a為該研究區(qū)的布格重力異常，比例尺為1∶50 000，局部細(xì)化為1∶25 000。

從圖7b中可以看出，由于研究區(qū)數(shù)據(jù)中噪聲的干擾，常規(guī)FFT法得到的垂向二階導(dǎo)數(shù)結(jié)果穩(wěn)定性差，單點(diǎn)異常多且雜亂，幾乎不能從中看出研究區(qū)的淺部異常特征；圖7c為Tikhonov正則化迭代法對(duì)該研究區(qū)的垂向二階導(dǎo)數(shù)計(jì)算結(jié)果，其中，正則化參數(shù)α為0.001，迭代5次。從圖7c中可以看出Tikhonov正則化迭代法計(jì)算結(jié)果穩(wěn)定，噪聲得到了有效壓制，研究區(qū)的淺部異常清晰明顯、連續(xù)性高，這有利于對(duì)研究區(qū)斷裂的劃分以及成礦遠(yuǎn)景區(qū)的預(yù)測(cè)。

為了驗(yàn)證Tikhonov正則化迭代法計(jì)算的垂向二階導(dǎo)數(shù)結(jié)果，我們對(duì)圖7c進(jìn)行二次積分，得到了通過(guò)重力垂向二階導(dǎo)數(shù)計(jì)算出的重力異常圖7d。通過(guò)對(duì)原始重力異常圖7a以及積分計(jì)算出的重力異常圖7d對(duì)比可以看出，計(jì)算所得重力異常和原始重力異常在形態(tài)與數(shù)值上吻合得非常好，說(shuō)明該方法對(duì)垂向二階導(dǎo)數(shù)的計(jì)算精度高。

據(jù)文獻(xiàn)[16]修編。圖6 研究區(qū)地質(zhì)圖 Fig.6 Geological map of the study area

圖7 原始重力異常(a)、常規(guī)FFT法垂向二階導(dǎo)數(shù)(b)、Tikhonov正則化迭代法得到的垂向二階導(dǎo)數(shù)(c)及垂向二次積分得到的重力異常(d)Fig.7 Original gravity anomaly (a), the second vertical derivative using FFT (b), the second vertical derivative using Tikhonov regularization iteration method (c) and the gravity anomaly result by the secondary vertical integration (d)

5 結(jié)論

1)本文將Tikhonov正則化與迭代法相結(jié)合，提出一種計(jì)算位場(chǎng)高階導(dǎo)數(shù)的新方法，通過(guò)分析該算子的濾波特性，表明該方法具有一定的穩(wěn)定性及保幅性。

2)通過(guò)模型試驗(yàn)也可以看出，該方法計(jì)算出的高階導(dǎo)數(shù)比常規(guī)FFT求導(dǎo)法更穩(wěn)定，因此對(duì)異常體邊界的識(shí)別效果更好。Tikhonov正則化迭代法對(duì)實(shí)際數(shù)據(jù)的處理結(jié)果說(shuō)明了該方法較常規(guī)FFT求導(dǎo)法有更高的穩(wěn)定性和實(shí)用價(jià)值。

3)該方法的提出可以為實(shí)際數(shù)據(jù)的處理提供一種參考，提高后續(xù)異常解釋的準(zhǔn)確性。

[1] 肖鋒，吳燕岡，孟令順. 重力異常圖中的邊界增強(qiáng)和提取技術(shù)[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2011，41(4)：1197-1203.

Xiao Feng, Wu Yangang, Meng Lingshun.Edge Enhancement and Detection Technology in Gravity Anomaly Map [J]. Journal of Jilin University (Earth Science Edition), 2011，41(4)：1197-1203.

[2] 張沖，黃大年，秦朋波，等．重力場(chǎng)向下延拓的三階Adams-Bashforth公式法[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2017，47(5)：1533-1542.

Zhang Chong, Huang Danian, Qin Pengbo, et al. Third-Order Adams-Bashforth Formula Method for Downward Continuation of Gravity Field [J]. Journal of Jilin University (Earth Science Edition), 2017，47 (5)：1533-1542.

[3] 徐世浙. 用邊界單元法計(jì)算任意形體的重力異常及其導(dǎo)數(shù)[J]. 石油物探, 1984，23(2):22-37.

Xu Shizhe. The Calculation of the Gravitational Anomaly and Its Derivative of the Geologic Body with Arbitrary Configuration by Boundary-Elements Method [J]. Geophysical Prospecting for Petroleum, 1984, 23(2):22-37.

[4] 汪炳柱. 用樣條函數(shù)法求重力異常二階垂向?qū)?shù)和向上延拓計(jì)算[J]. 石油地球物理勘探, 1996, 31(3):415-422.

Wang Bingzhu. Computing the Vertical Second Derivative and Upward Continuation of Gravity Anomaly by Spline Function Method[J]. Oil Geophysical Prospecting, 1996, 31(3):415-422.

[5] 姚長(zhǎng)利，黃衛(wèi)寧，管志寧. 綜合利用位場(chǎng)及其垂直梯度的快速樣條曲化平方法[J]. 石油地球物理勘探, 1997, 32(2):229-236.

Yao Changli, Huang Weining,Guan Zhining. Fast Splines Conversion of Curvedsurface Potential Field and Vertical Gradient Data into Horizontal-Plane Data [J]. Oil Geophysical Prospecting, 1997, 32(2):229-236.

[6] Wang B, Krebes E S, Ravat D. High-Precision Poten-tial-Field and Gradient-Component Transformations and Derivative Computations Using Cubic B-Splines [J]. Geophysics, 2008, 73(5):135-142.

[7] Clarke G K C. Optimum Second-Derivative and Down-ward-Continuation Filters [J]. Geophysics, 1969, 34(3): 424-437.

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

Hou Zhongchu. Filtering of Smooth Compensation [J]. Geophysical Prospecting for Petroleum, 1981, 20(2): 22-29.

[9] Pa?teka R, Richter F P, Karcol R, et al. Regularized Derivatives of Potential Fields and Their Role in Semi-Automated Interpretation Methods [J]. Geophysical Prospecting, 2009, 57(4): 507-516.

[10] 曾小牛, 李夕海, 賈維敏,等. 位場(chǎng)各階垂向?qū)?shù)換算的新正則化方法[J]. 地球物理學(xué)報(bào), 2015, 58(4):1400-1410.

Zeng Xiaoniu, Li Xihai, Jia Weimin, et al. A New Regularization Method for Calculating the Vertical Derivatives of the Potential Field [J].Chinese Journal of Geophysics, 2015, 58(4): 1400-1410.

[11] 王彥國(guó), 張瑾. 位場(chǎng)高階導(dǎo)數(shù)的波數(shù)域迭代法[J]. 物探與化探, 2016, 40(1):143-147.

Wang Yanguo, Zhang Jin. The Iterative Method for Higher-Order Derivative Calculation of Potential Fields in Wave Number Domain [J]. Geophysical and Geochemical Exploration, 2016, 40(1):143-147.

[12] Fedi M, Florio G. Detection of Potential Fields Sou-rce Boundaries by Enhanced Horizontal Derivative Method [J]. Geophysical Prospecting, 2001, 49(1):40-58.

[13] Fedi M, Florio G. A Stable Downward Continuation by Using the ISVD Method [J]. Geophysical Journal International, 2002, 151(1):146-156.

[14] Cooper G R J, Cowan D R. Filtering Using Variable Order Vertical Derivatives [J]. Computers & Geosciences, 2004, 30(5):455-459.

[15] 卞光浪, 翟國(guó)君, 高金耀,等. 總強(qiáng)度磁場(chǎng)穩(wěn)健向下延拓的改進(jìn)泰勒級(jí)數(shù)法[J]. 測(cè)繪學(xué)報(bào), 2014,43(1):30-36.

Bian Guanglang, Zhai Guojun, Gao Jinyao, et al. Improved Taylor Series Approach for Stable Downward Continuation of Total Strength of Geomagnetic Field [J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(1):30-36.

[16] Xu M L, Yang C B, Wu Y G, et al. Edge Detection in the Potential Field Using the Correlation Coefficients of Multidirectional Standard Deviations [J]. Applied Geophysics, 2015, 12: 23-34.