杜 威，許家姝，吳燕岡，郝夢成

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

0 引言

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

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

1 方法原理

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

(1)

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

(2)

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

(3)

(4)

(5)

(6)

(7)

由式(1)可知：

(8)

代入式(7)得：

(9)

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

(10)

(11)

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

(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ù)的點距均為1 km，垂向一階導(dǎo)數(shù)的正則化參數(shù)α的取值為1.0，垂向二階導(dǎo)數(shù)的正則化參數(shù)α的取值為0.1。對理論垂向?qū)?shù)算子φm、Tikhonov正則化迭代法垂向?qū)?shù)換算算子的濾波特性進(jìn)行對比分析可得：理論垂向?qū)?shù)因子φm隨著頻率的增加急劇上升，且導(dǎo)數(shù)階次越高，高頻成分被放大得越強；Tikhonov正則化迭代算子在低頻段逼近理論垂向?qū)?shù)算子，保證了低頻有用信號處導(dǎo)數(shù)的換算，在高頻階段卻能對其進(jìn)行壓制；隨著迭代次數(shù)的增加，Tikhonov正則化迭代算子逐漸趨近于理論垂向?qū)?shù)算子。為了分析正則化參數(shù)α的取值對Tikhonov正則化迭代算子濾波特性的影響，圖1b和圖2b分別為迭代3次時，垂向一階和二階導(dǎo)數(shù)Tikhonov正則化迭代算子的濾波響應(yīng)特征曲線?？梢钥闯觯寒?dāng)?shù)螖?shù)不變時，隨著導(dǎo)數(shù)階次的增加，Tikhonov正則化迭代算子對高頻成分的壓制作用逐漸增強，保證了高階導(dǎo)數(shù)換算的穩(wěn)定性；正則化參數(shù)α取值越小，Tikhonov正則化迭代算子在高頻段的壓制作用越弱，當(dāng)其值為0時，則結(jié)果為理論導(dǎo)數(shù)算子；當(dāng)正則化參數(shù)α的取值過大時，對低頻段也產(chǎn)生一定的壓制作用，則會影響濾波效果。

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

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

3 模型試驗

文中采用2個二度直立柱體的疊加異常進(jìn)行數(shù)值檢驗。模型體1和2長均為2.0 km，上頂埋深h為2.0 km，下底埋深h為2.5 km;其中心位置x分別為10、15 km，點距為0.1 km，2個模型體的剩余密度均為1.0 g/cm3。組合模型體相對位置及其產(chǎn)生的重力異常Δg如圖3a所示，圖3b和圖3c分別是無噪聲情況下理論的垂向二階導(dǎo)數(shù)和垂向三階導(dǎo)數(shù)。添加重力異常幅值0.1%的隨機噪聲，圖4和圖5分別是含噪聲情況下常規(guī)FFT法和Tikhonov正則化迭代法所計算出的垂向二階以及垂向三階導(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正則化迭代法計算垂向二階和三階導(dǎo)數(shù)的對比圖。從圖中可以看出，由于噪聲的干擾，用常規(guī)FFT變換計算的垂向?qū)?shù)結(jié)果存在著明顯的波動性，且導(dǎo)數(shù)的階次越高，穩(wěn)定性越差，無法識別有效異常，這是由于理論導(dǎo)數(shù)算子對高頻成分的放大作用。Tikhonov正則化迭代法計算垂向二階及三階導(dǎo)數(shù)結(jié)果穩(wěn)定，對噪聲的壓制能力強，這是因為正則化迭代法的濾波特性，在低頻段較好地逼近理論因子，而在高頻段可以較好地壓制噪聲。對比圖3b和3c中不含噪聲計算得到的垂向二階導(dǎo)數(shù)和三階導(dǎo)數(shù)可以看出，隨著求導(dǎo)次數(shù)的增加，高頻部分的影響越來越大，即圖3c中曲線的震蕩幅度大于圖3b中。從圖4b和圖5b可以看出，Tikhonov正則化迭代法計算得到的導(dǎo)數(shù)，在數(shù)量級上和換算得到的導(dǎo)數(shù)(圖3b和3c)數(shù)量級是相同的，垂向三階導(dǎo)數(shù)的結(jié)果圖中并沒有出現(xiàn)明顯的震蕩，再次說明了Tikhonov正則化迭代法的準(zhǔn)確性和穩(wěn)定性。圖4b模型體1和2中，零值點對應(yīng)的橫坐標(biāo)分別為8.5和11.3、13.7和16.5，圖5b模型體1和2中，零值點對應(yīng)的橫坐標(biāo)分別為8.6和11.2、13.8和16.4，可以看出，隨著所計算導(dǎo)數(shù)的階次增加，曲線的零值點位置與地質(zhì)體邊界的對應(yīng)效果越來越好。

4 實際應(yīng)用

為了試驗Tikhonov正則化迭代法求垂向?qū)?shù)對實際資料的處理效果，我們對老撾萬象盆地的布格重力異常Δg數(shù)據(jù)進(jìn)行垂向二階導(dǎo)數(shù)的計算，其中，實際數(shù)據(jù)的點距為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)定性差，單點異常多且雜亂，幾乎不能從中看出研究區(qū)的淺部異常特征；圖7c為Tikhonov正則化迭代法對該研究區(qū)的垂向二階導(dǎo)數(shù)計算結(jié)果，其中，正則化參數(shù)α為0.001，迭代5次。從圖7c中可以看出Tikhonov正則化迭代法計算結(jié)果穩(wěn)定，噪聲得到了有效壓制，研究區(qū)的淺部異常清晰明顯、連續(xù)性高，這有利于對研究區(qū)斷裂的劃分以及成礦遠(yuǎn)景區(qū)的預(yù)測。

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

據(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é)合，提出一種計算位場高階導(dǎo)數(shù)的新方法，通過分析該算子的濾波特性，表明該方法具有一定的穩(wěn)定性及保幅性。

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

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

[1] 肖鋒，吳燕岡，孟令順. 重力異常圖中的邊界增強和提取技術(shù)[J]. 吉林大學(xué)學(xué)報(地球科學(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] 張沖，黃大年，秦朋波，等．重力場向下延拓的三階Adams-Bashforth公式法[J]. 吉林大學(xué)學(xué)報(地球科學(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] 徐世浙. 用邊界單元法計算任意形體的重力異常及其導(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]. 石油地球物理勘探, 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] 姚長利，黃衛(wèi)寧，管志寧. 綜合利用位場及其垂直梯度的快速樣條曲化平方法[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] 侯重初. 補償圓滑濾波方法[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] 曾小牛, 李夕海, 賈維敏,等. 位場各階垂向?qū)?shù)換算的新正則化方法[J]. 地球物理學(xué)報, 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] 王彥國, 張瑾. 位場高階導(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] 卞光浪, 翟國君, 高金耀,等. 總強度磁場穩(wěn)健向下延拓的改進(jìn)泰勒級數(shù)法[J]. 測繪學(xué)報, 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.