劉 帆, 李正文, 種 洋

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

0 引言

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

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

1 正則化參數(shù)選取方法

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

(1)

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

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

(2)

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

(3)

(4)

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

1.1 頻率域的L曲線法

(5)

(6)

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

(7)

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

(8)

(9)

1.2 借鑒GCV法的構(gòu)造泛函法

(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)的計算流程少，計算更加簡潔、易于實現(xiàn)。式(17)既滿足了頻域計算要求,又兼顧了方便快速計算的優(yōu)點,筆者使用式(17)計算得到的V(α)最小值對應(yīng)的α為最優(yōu)正則化參數(shù)，為了驗證此方法的正確性，分別使用仿真磁異常數(shù)據(jù)和實測磁異常數(shù)據(jù)進行了向下延拓實驗(本文提出的參數(shù)選取方法為構(gòu)造泛函法)。構(gòu)造泛函法選取正則化參數(shù)求解位場向下延拓問題的流程圖如圖1所示。

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

2 仿真數(shù)值實驗

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

3x1z1sin 2IcosD+

3x1y1cos2Isin 2D-

3y1z1sin 2IsinD]

(18)

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

表1 球體磁異常模型參數(shù)

表2 不同高度面上理論磁異常數(shù)據(jù)統(tǒng)計

圖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含噪聲理論磁異常

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

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

表3 仿真數(shù)值向下延拓實驗精度統(tǒng)計

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

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

迭代方法實際最優(yōu)參數(shù)RMSE/nTMRE/%Tikhonov正則化迭代法0.39632.63126.76Landweber正則化迭代法1.94792.66886.85

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

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

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

3 實測數(shù)據(jù)實驗

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

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

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

圖7 原始磁異常數(shù)據(jù)格網(wǎng)Fig.7 The original magnetic anomaly data grid

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

表6 實測航磁數(shù)據(jù)向下延拓實驗精度統(tǒng)計

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

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

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

4 結(jié)論

通過對仿真數(shù)值實驗和實測數(shù)據(jù)實驗的統(tǒng)計結(jié)果進行分析，得到以下結(jié)論：

1)在頻域內(nèi)使用正則化迭代法進行向下延拓時，使用構(gòu)造泛函法選取的參數(shù)優(yōu)于使用L曲線法選取的參數(shù)，使用Tikhonov正則化迭代法時，構(gòu)造泛函法的優(yōu)勢更加明顯。

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

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

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

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

[3] 劉東甲, 洪天求, 賈志海,等．位場向下延拓的波數(shù)域迭代法及其收斂性 [J].地球物理學(xué)報，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].地球物理學(xué)報,1989,32(5):600-608.LIANGJW.Downwardcontinuationofregularizationmethodsforpotentialfields[J].ChineseJ.Geophys., 1989, 32(5):600-608 . (InChinese)

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

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

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

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

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

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

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

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