一種基于GCV法改進的頻域內位場正則化向下延拓的參數選取方法

劉 帆, 李正文, 種 洋

(1.61365部隊，天津 300140；2.61206部隊，北京 100042;3.信息工程大學 地理空間信息學院，鄭州 450001;4.地理信息工程國家重點實驗室，西安 710054 5.海軍工程大學 導航工程系，武漢 430033)

0 引言

位場向下延拓是航空重力和磁力測量數據處理的重要技術，向下延拓能夠突出局部場源，以此提高異常解釋的可靠性；向下延拓也可用于構建高精度的空間地磁數據庫以滿足高精度地磁導航的需求。通常在空間域上的位場延拓算法面臨著大型矩陣存儲和運算困難的問題，因為借助FFT(快速傅里葉變換)的頻域延拓算法具有計算快速和處理格網數據方便的優勢[1]，所以頻率域上的向下延拓方法一直是研究熱點。2006年徐世浙[2]院士提出了積分迭代法，此法實現了大深度的位場向下延拓；劉東甲[3]將徐世浙的空間域迭代方法推廣至波數域，得到了波數域迭代法。由于向下延拓問題是一個不適定問題，下延算子會放大觀測數據中的高頻噪聲； Tikhonov[4]于上世紀60年代使用正則化方法給出了穩定的向下延拓頻率響應；梁錦文[5]從四種頻率響應的分析中得出了合適的頻率響應公式。如今，反問題的正則化解法被廣泛應用于位場的向下延拓，其延拓精度要優于絕大多數的延拓方法，但是正則化方法的實施難點在于正則化參數的確定，并且正則化參數對延拓精度有很大的影響。目前有關正則化參數確定方法的研究較少，馬濤等[6]研究了基于L曲線法的位場下延正則化參數的選取問題；劉曉剛等[7]研究了重磁數據下延中的最優正則化參數確定方法；曾小牛[8]提出了基于徑向譜的位場向下延拓正則參數選取方法。

正則化參數的選取方法大致分為先驗選取和后驗選取兩大類，先驗策略難以在實際應用中檢驗其成立的條件，因而確定正則化參數的后驗策略研究具有很大的實用意義，目前在位場延拓問題的正則化求解中廣泛使用的后驗策略有L曲線法、GCV法(廣義交叉驗證法)、擬最優準則等。L曲線法是被普遍認為比較可靠的方法，但是基于頻域的L曲線算法十分復雜；GCV法的計算比較方便快速，但由于將空間域中矩陣的求跡運算轉換至頻域內求解十分困難，所以GCV法在頻域中難以實施，其他的方法幾乎難以轉換至頻域內實現。筆者借鑒GCV函數值的計算形式，提出一種在頻域內可實施的參數選取方法，并分別通過仿真數據和實測數據的向下延拓實驗來進行驗證，實驗結果表明這種方法實用可行，且具有良好的延拓精度。

1 正則化參數選取方法

在向下延拓中u0(x,y)表示觀測面上的位場數據，uh(ε,η)表示延拓面上的位場數據，h為延拓深度，由位場向上延拓的表達式[9]，可以得到：

(1)

式(1)可表示為卷積方程式：

u0(x,y)=Kuh(x,y)

(2)

其中:K為第一類Fredholm積分算子，使用Tikhonov正則化迭代法時，式(2)的解為[10]：

(3)

(4)

在式(3)、式(4)中α取值較大會使擬合誤差過大，但α取值較小會使正則解過于接近原問題的不適定解而發生震蕩，筆者使用L曲線法和一種新的方法來確定正則化參數，并比較使用這兩種參數選取方法的延拓精度。

1.1 頻率域的L曲線法

(5)

(6)

通常取L曲線的曲率最大點對應的α作為最優正則化參數，曲率計算公式如式(7)所示。

(7)

用一階差分和二階差分來代替函數ρ(α)、θ(α)的一階導數與二階導數。為了快速計算出c(α)序列，把ρ(α)、θ(α)轉換到頻域進行數值計算，根據Parseval等式，這兩個函數的頻域計算式如下：

(8)

(9)

1.2 借鑒GCV法的構造泛函法

(10)

式中：

Hα=K(KTK+αI)-1KT

(11)

(12)

則有：

(13)

=K(x,y)*K(x,y)

(14)

對式(14)兩端同時進行傅里葉變換，得到：

Φ(u,v)Φ(u,v)Gα(u,v)+αGα(u,v)=

Φ(u,v)Φ(u,v)

(15)

(16)

(17)

對比L曲線法，從式(16)、式(17)中可以看到，由于Φ(u,v)很容易求得,所以Gα是容易得到的,并且式(17)遠比式(7)的計算流程少，計算更加簡潔、易于實現。式(17)既滿足了頻域計算要求,又兼顧了方便快速計算的優點,筆者使用式(17)計算得到的V(α)最小值對應的α為最優正則化參數，為了驗證此方法的正確性，分別使用仿真磁異常數據和實測磁異常數據進行了向下延拓實驗(本文提出的參數選取方法為構造泛函法)。構造泛函法選取正則化參數求解位場向下延拓問題的流程圖如圖1所示。

圖1 構造泛函法用于正則化迭代法求解下延位場的流程Fig.1 The process for downward continuation using regularized iteration method with the constructing functional method

2 仿真數值實驗

筆者使用球體磁異常正演模型生成仿真磁力測量數據[11]：

3x1z1sin 2IcosD+

3x1y1cos2Isin 2D-

3y1z1sin 2IsinD]

(18)

式中：取Z軸向上為正；μ0=4π×10-7H/m為真空導磁率；M為磁化強度；v為球體的體積；I為磁化傾角；D為磁化偏角；x1=x-x0、y1=y-y0、z1=z-z0；(x,y,z)是計算點坐標；(x0,y0,z0)為模型球體的球心坐標，實驗所用球體磁異常正演模型參數見表1；分別生成高度為0m和100m的理論磁異常平面格網數據，每個平面上的格網間距為20m×20m,測點數88×88個，格網數據的統計特性見表2，實驗中在高度為100m的格網數據中加入了標準差為3nT的高斯白噪聲，然后將含噪聲的數據向下延拓5倍點距，最后將延拓結果與高度為0m的磁異常理論值相比較，用平均相對誤差(MRE)和均方根誤差(RMSE)來評價不同參數選取方法的延拓精度，各觀測面數據等值線見圖2。為了減小邊界效應，在進行延拓實驗前，先對高度為100m含有噪聲的格網數據使用了區域場[13]擴邊法進行擴邊處理，擴邊點數為各方向20個，計算結束后再縮邊20點得到向下延拓結果。實驗中分別使用了Tikhonov正則化迭代法和Landweber正則化迭代法進行向下延拓，兩種迭代法的迭代次數均設置為20次，Tikhonov正則化迭代法中參數α是從首項α1為10-8、末項αN為1的等比數列中選取的，N=200；Landweber正則化迭代法中參數α是從首項α1為0.7、末項αN為2的等比數列中選取的，N=200；兩種迭代法都分別使用L曲線法和本文提出的方法選取最優參數以做對照。

表1 球體磁異常模型參數

表2 不同高度面上理論磁異常數據統計

圖2 各觀測面上的磁異常理論值等值線圖Fig.2 The contour map of magnetic anomaly on each observation plane(a)z=0 m理論磁異常；(b)z=100 m理論磁異常；(c)z=100 m含噪聲理論磁異常

向下延拓精度統計見表3，Tikhonov正則化迭代法延拓結果等值線圖的對比見圖3，Landweber正則化迭代法延拓結果等值線圖的對比見圖4。Tikhonov正則化迭代法延拓結果過異常中心的剖面對比見圖5(a)，Landweber正則化迭代法延拓結果過異常中心的剖面對比見圖5(b)。

為了進一步比較和說明各個方法所選參數的優劣，統計了試驗中兩種迭代法迭代20次所對應的最優參數及其延拓精度(表4)。兩種迭代法中根據式(17)計算出的V(α)值與參數α對應的關系曲線見圖6。

表3 仿真數值向下延拓實驗精度統計

圖3 頻域Tikhonov正則化迭代法向下延拓結果等值線圖Fig.3 The contour map of Tikhonov regularization iterative method’s downward continuation results in frequency domain(a)最優參數向下延拓結果；(b)L曲線法所選參數向下延拓結果;(c)構造泛函數所選參數向下延拓結果

圖4 頻域Landweber正則化迭代法向下延拓結果等值線圖Fig.4 The contour map of Landweber regularization iterative method’s downward continuation results in frequency domain(a)最優參數向下延拓結果；(b)L曲線法所選參數向下延拓結果;(c)構造泛函數所選參數向下延拓結果

迭代方法實際最優參數RMSE/nTMRE/%Tikhonov正則化迭代法0.39632.63126.76Landweber正則化迭代法1.94792.66886.85

圖5 頻域內兩種正則化迭代法向下延拓的剖面對比圖Fig.5 Sectional comparison of two regularized iteration methods in the frequency domain(a)Tikhonov正則化迭代法延拓結果過異常中心的剖面對比;(b)Landweber正則化迭代法延拓結果過異常中心的剖面對比

圖6 頻域內兩種正則化迭代法向下延拓的α-V(α)曲線Fig.6 The α-V(α) curves of two regularization iterative methods’ downward continuation in frequency domain(a)頻域Tikhonov正則化迭代法的v倍曲線；(b)頻域Landweber正則化迭代法的v倍曲線

對比表3和表4，再結合圖1～圖4，可以發現在兩種迭代法中，使用構造泛函法均比使用L曲線法得到的參數更優，均不同程度地提高了延拓精度；并且在這兩種正則化迭代法中，使用構造泛函法選取的參數更接近實際最優參數，延拓精度也更接近最優延拓精度。對于Tikhonov正則化迭代法而言，使用構造泛函法選取參數的延拓精度比使用L曲線法的延拓精度有大幅度提高；對于Landweber正則化迭代法而言，使用構造泛函法選取參數的延拓精度比使用L曲線法的延拓精度有略微提高。從圖6中可以看出，在Tikhonov正則化迭代法中，V(α)曲線在α∈(0.15,0.6)時較為平緩，曲線在極小值點附近變化很小；而在Landweber正則化迭代法中V(α)曲線是明顯的先下降再上升的形態，在極小值點附近變化較為明顯。

3 實測數據實驗

筆者使用美國地質勘探局(USGS)海軍科學研究實驗所(NRL)于2008年在阿富汗實測的航空磁力資料進行數據處理檢驗。從官方網站上所下載到的數據，已經過數據預處理并被歸算到離地面5 000m的高度。我們截取了東經64°～64.814°、北緯31°～31.814°之間的經過交叉測線改正后的磁異常觀測數據進行延拓實驗。試驗中先將該區域內的數據處理成點距為1km的88×88點的格網(延拓前的原始磁異常格網)，然后將此格網向上延拓5倍點距后再向下延拓5倍點距至原高度，用向下延拓后的結果與原始實測格網異常數據進行比較并統計精度。試驗中使用區域場法對原始格網數據進行擴邊，擴邊點數為20個，對計算得到的結果再縮邊20個點得到向下延拓結果，向下延拓的實驗方法與仿真數值試驗方法相同。圖7為用于延拓實驗的原始磁異常格網，圖8為分別使用L曲線法和構造泛函法選取參數的Tikhonov正則化迭代法的延拓結果。圖9為分別使用L曲線法和構造泛函法選取參數的Landweber正則化迭代法的延拓結果。Tikhonov正則化迭代法延拓結果過異常中心的剖面對比見圖10(a)，Landweber正則化迭代法延拓結果過異常中心的剖面對比見圖10(b)。表5為延拓前的原始磁異常格網數據的統計特性，表6為實測航磁數據延拓實驗的精度統計。

從表6中可以看出，在兩種迭代法中使用構造泛函法選取參數的延拓精度都要優于L曲線法。在實測數據實驗中，基于L曲線法選擇參數的Tikhonov正則化迭代法的延拓精度明顯好于仿真實驗，因為上延是低通濾波過程，平滑掉了原始數據中的高頻干擾，所以這一現象也說明了Tikhonov迭代法抗干擾能力比Landweber迭代法差。

對比表3和表6，可以看到無論是仿真實驗還是實測數據實驗，在Tikhonov正則化迭代法中，使用構造泛函法和L曲線法選出的參數相差較大，但是在Landweber正則化迭代法中，使用兩種方法選出的參數相差不大，這說明在Tikhonov正則化迭代法中，構造函數法選取參數的優勢更加明顯，在實測數據實驗中，所有方法的延拓精度都十分理想，平均相對誤差均小于7%，相比圖7，可以看到圖8和圖9中的延拓結果幾乎沒有明顯的畸變和邊界效應。

圖7 原始磁異常數據格網Fig.7 The original magnetic anomaly data grid

觀測面高度/m最小值/nT最大值/nT平均值/nT標準差/nT5000(原始異常)-282.690859.8296-151.906166.6123

表6 實測航磁數據向下延拓實驗精度統計

圖8 頻域Tikhonov正則化迭代法向下延拓結果等值線圖Fig.8 The contour map of Tikhonov regularization iterative method’s downward continuation results in frequency domain(a)L曲線法所選參數向下延拓結果；(b)構造泛函數所選參數向下延拓結果

圖9 頻域Landweber正則化迭代法向下延拓結果等值線圖Fig.9 The contour map of Landweber regularization iterative method’s downward continuation results in frequency domain(a)L曲線法所選參數向下延拓結果；(b)構造泛函數所選參數向下延拓結果

圖10 頻域內兩種正則化迭代法向下延拓的剖面對比圖Fig.10 Sectional comparison of two regularized iteration methods in the frequency domain(a)Tikhonov正則化迭代法延拓結果過異常中心的剖面對比;(b)Landweber正則化迭代法延拓結果過異常中心的剖面對比

4 結論

通過對仿真數值實驗和實測數據實驗的統計結果進行分析，得到以下結論：

1)在頻域內使用正則化迭代法進行向下延拓時，使用構造泛函法選取的參數優于使用L曲線法選取的參數，使用Tikhonov正則化迭代法時，構造泛函法的優勢更加明顯。

2)相比L曲線法而言,構造泛函法的計算過程不涉及到求導運算，計算過程更加方便快速，并且具有更好的延拓精度。

3)構造泛函法選出的參數與實際最優參數十分接近，且在兩種迭代法中選取參數的效果比較穩定，不同數據的延拓精度處于同一個水平。

[1] 周波陽,羅志才,許闖,等．航空重力數據向下延拓的FFT快速算法比較[J].大地測量與地球動力學,2013,33(1):64-73.ZHOUBY,LUOZC,XUC,etal．Comparisonamongfastfouriertransformalgorithmsfordownwardcontinuationofairbornegravitydata[J].Journalofgeodesyandgeodynamics,2013,33(1):64-73.(InChinese)

[2] 徐世浙. 位場延拓的積分-迭代法[J].地球物理學報,2006,49(4):1176-1182.XUSZ.Theintegral-iterationmethodforcontinuationofpotentialfields[J].ChineseJ.Geophys,2006,49(4):1176-1182.(InChinese)

[3] 劉東甲, 洪天求, 賈志海,等．位場向下延拓的波數域迭代法及其收斂性 [J].地球物理學報，2009,52(6):1599-1605.LIUDJ,HONGTQ,JIAZH,etal．WavenumberdomainiterationmethodfordownwardcontinuationofPotentialfieldsanditsconvergenee. [J].ChineseJ.Geophys,2009,52(6):1599-1605. (InChinese)

[4]TIKHONOVAN,ARSENINVY.Solutionofcertainintegralequationsofthefirstkind[J].JournalofAssociateComputerMachine,1962(9):84-97.

[5] 梁錦文.位場向下延拓的正則化方法[J].地球物理學報,1989,32(5):600-608.LIANGJW.Downwardcontinuationofregularizationmethodsforpotentialfields[J].ChineseJ.Geophys., 1989, 32(5):600-608 . (InChinese)

[6] 馬濤,陳偉龍,吳美平,等．基于L曲線法的位場向下延拓正則化參數選擇[J] .地球物理學進展,2013,28(5):2485-2494.MAT，CHENLW，WUMP，eta1．TheselectionofregularizationparameterindownwardcontinuationofpotentialfieldbasedonL-curvemethod[J].ProgressinGeophys,2013,28(5):2485-2494. (InChinese)

[7] 劉曉剛,李迎春,肖云,等．重力與磁力測量數據向下延拓中最優正則化參數確定方法[J] .測繪學報,2014,43(9):881-887.LIUXG,LIYC,XIAOY,etal．Optimalregularizationparameterdeterminationmethodindownwardcontinuationofgravimetricandgeomagneticdata[J].ActaGeodaeticaetCartographicaSinica,2014,43(9):881-887. (InChinese)

[8] 曾小牛,李夕海，陳鼎新，等．基于徑向譜的位場向下延拓正則參數選取方法[J].石油地球物理勘探，2015,50(4):749-755.ZENGXN,LIXH，CHENDX,etal．Newregularizationparameterselectionforpotentialfielddownwardcontinuationbasedontheradialspectrum[J].OilGeophysicalProspecting，2015,50(4):749-755. (InChinese)

[9] 曾小牛,李夕海,韓紹卿,等．位場向下延拓三種迭代方法之比較[J].地球物理學進展,2011 ,26(3):908-915.ZENGXN,LIXH,HANSQ,etal.Acomparisonofthreeiterationmethodsfordownwardcontinuationofpotentialfields[J].ProgressinGeophys,2011 ,26(3):908-915. (InChinese)

[10]徐新禹, 李建成, 王正濤,等．Tikhonov正則化方法在GOCE重力場求解中的模擬研究[J].測繪學報,2010,39(5):465-470.XUXY,LIJC,WANGZT,etal．ThesimulationresearchontheTikhonovregularizationappliedingravityfielddeterminationofGOCEsatellitemission[J].ActaGeodaeticaetCartographicaSinica,2010,39(5):465-470. (InChinese)

[11]駱遙,吳美平．位場向下延拓的最小曲率方法[J].地球物理學報，2016,59(1):240-251.LUOY,WUMP．Minimumcurvaturemethodfordownwardcontinuationofpotentialfielddata[J].ChineseJ.Geophys,2016,59(1):240-251. (InChinese)

[12]胡小平．水下地磁導航技術[M].北京:國防工業出版社,2013.HUXP．Technologiesonunderwatergeomagneticfieldnavigation[M].Beijing:NationalDefenceIndustryPress,2013. (InChinese)

[13]段本春, 徐世浙．磁(重力)異常局部場與區域場分離處理中的擴邊方法研究[J].物探化探計算技術 ,1997,19(4):298-304.DUANBC,XUSZ．Astudyoftheschemeofextendingedgeintheprocessingofseparatinglocalfieldfromregionalfieldformagnetic/gravityanomaly[J].ComputingTechniquesForGeophysicalAndGeochemicalExploration,1997,19(4):298-304. (InChinese)