位場垂向高階導數的Tikhonov正則化迭代法

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

吉林大學地球探測科學與技術學院，長春 130026

0 引言

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

由于空間域的算法普遍存在計算速度慢等缺點，因此目前常在頻率域對位場的高階垂向導數進行計算；但是頻率域的導數響應因子具有高頻放大作用，當所求導數的階數越高，數據對噪聲的敏感程度就越大，導致所求結果不穩定。為了壓制高階導數計算結果的不穩定性，本文提出一種基于Tikhonov正則化迭代法計算位場各階垂向導數的方法，并通過模型試驗及實際數據的應用，對該方法計算垂向導數的穩定性進行了驗證。

1 方法原理

位場T(x,y)與其各階垂向導數Dm(x,y)(m代表階數)在波數域的關系式為

(1)

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

(2)

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

(3)

(4)

(5)

(6)

(7)

由式(1)可知：

(8)

代入式(7)得：

(9)

由數學歸納法很容易得到迭代通式：

(10)

(11)

當迭代無窮大時，

(12)

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

2 濾波特性分析

圖1和圖2分別是垂向一階導數和垂向二階導數的Tikhonov正則化迭代算子濾波特性。其中，圖1a和圖2a分別是不同迭代次數的垂向一階導數和垂向二階導數Tikhonov正則化迭代算子的濾波響應特征曲線。其中，橫坐標ω為徑向圓頻率，數據的點距均為1 km，垂向一階導數的正則化參數α的取值為1.0，垂向二階導數的正則化參數α的取值為0.1。對理論垂向導數算子φm、Tikhonov正則化迭代法垂向導數換算算子的濾波特性進行對比分析可得：理論垂向導數因子φm隨著頻率的增加急劇上升，且導數階次越高，高頻成分被放大得越強；Tikhonov正則化迭代算子在低頻段逼近理論垂向導數算子，保證了低頻有用信號處導數的換算，在高頻階段卻能對其進行壓制；隨著迭代次數的增加，Tikhonov正則化迭代算子逐漸趨近于理論垂向導數算子。為了分析正則化參數α的取值對Tikhonov正則化迭代算子濾波特性的影響，圖1b和圖2b分別為迭代3次時，垂向一階和二階導數Tikhonov正則化迭代算子的濾波響應特征曲線?？梢钥闯觯寒數螖挡蛔儠r，隨著導數階次的增加，Tikhonov正則化迭代算子對高頻成分的壓制作用逐漸增強，保證了高階導數換算的穩定性；正則化參數α取值越小，Tikhonov正則化迭代算子在高頻段的壓制作用越弱，當其值為0時，則結果為理論導數算子；當正則化參數α的取值過大時，對低頻段也產生一定的壓制作用，則會影響濾波效果。

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

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

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

3 模型試驗

文中采用2個二度直立柱體的疊加異常進行數值檢驗。模型體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。組合模型體相對位置及其產生的重力異常Δg如圖3a所示，圖3b和圖3c分別是無噪聲情況下理論的垂向二階導數和垂向三階導數。添加重力異常幅值0.1%的隨機噪聲，圖4和圖5分別是含噪聲情況下常規FFT法和Tikhonov正則化迭代法所計算出的垂向二階以及垂向三階導數。其中，二階垂向導數的正則化參數為0.1，迭代8次；垂向三階導數的正則化參數為0.5，迭代8次。

圖3 模型重力異常及模型位置(a)、垂向二階導數(b) 和垂向三階導數(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垂向二階導數(a)和Tikhonov正則化迭代法垂向二階導數(b)Fig.4 Second vertical derivative using FFT method (a) and the second vertical derivative using Tikhonov regularization iteration method (b)

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

圖4和圖5分別是常規FFT法和Tikhonov正則化迭代法計算垂向二階和三階導數的對比圖。從圖中可以看出，由于噪聲的干擾，用常規FFT變換計算的垂向導數結果存在著明顯的波動性，且導數的階次越高，穩定性越差，無法識別有效異常，這是由于理論導數算子對高頻成分的放大作用。Tikhonov正則化迭代法計算垂向二階及三階導數結果穩定，對噪聲的壓制能力強，這是因為正則化迭代法的濾波特性，在低頻段較好地逼近理論因子，而在高頻段可以較好地壓制噪聲。對比圖3b和3c中不含噪聲計算得到的垂向二階導數和三階導數可以看出，隨著求導次數的增加，高頻部分的影響越來越大，即圖3c中曲線的震蕩幅度大于圖3b中。從圖4b和圖5b可以看出，Tikhonov正則化迭代法計算得到的導數，在數量級上和換算得到的導數(圖3b和3c)數量級是相同的，垂向三階導數的結果圖中并沒有出現明顯的震蕩，再次說明了Tikhonov正則化迭代法的準確性和穩定性。圖4b模型體1和2中，零值點對應的橫坐標分別為8.5和11.3、13.7和16.5，圖5b模型體1和2中，零值點對應的橫坐標分別為8.6和11.2、13.8和16.4，可以看出，隨著所計算導數的階次增加，曲線的零值點位置與地質體邊界的對應效果越來越好。

4 實際應用

為了試驗Tikhonov正則化迭代法求垂向導數對實際資料的處理效果，我們對老撾萬象盆地的布格重力異常Δg數據進行垂向二階導數的計算，其中，實際數據的點距為0.2 km。該研究區的地質情況復雜，如圖6。圖7a為該研究區的布格重力異常，比例尺為1∶50 000，局部細化為1∶25 000。

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

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

據文獻[16]修編。圖6 研究區地質圖 Fig.6 Geological map of the study area

圖7 原始重力異常(a)、常規FFT法垂向二階導數(b)、Tikhonov正則化迭代法得到的垂向二階導數(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 結論

1)本文將Tikhonov正則化與迭代法相結合，提出一種計算位場高階導數的新方法，通過分析該算子的濾波特性，表明該方法具有一定的穩定性及保幅性。

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

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

[1] 肖鋒，吳燕岡，孟令順. 重力異常圖中的邊界增強和提取技術[J]. 吉林大學學報(地球科學版)，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]. 吉林大學學報(地球科學版)，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]. 石油物探, 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] 汪炳柱. 用樣條函數法求重力異常二階垂向導數和向上延拓計算[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] 姚長利，黃衛寧，管志寧. 綜合利用位場及其垂直梯度的快速樣條曲化平方法[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] 曾小牛, 李夕海, 賈維敏,等. 位場各階垂向導數換算的新正則化方法[J]. 地球物理學報, 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] 王彥國, 張瑾. 位場高階導數的波數域迭代法[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] 卞光浪, 翟國君, 高金耀,等. 總強度磁場穩健向下延拓的改進泰勒級數法[J]. 測繪學報, 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.