曾小牛 李夕海 劉繼昊 侯維君

(火箭軍工程大學(xué),陜西西安 710025)

0 引言

重力場(chǎng)向下延拓能夠突出局部或淺部異常，對(duì)重力數(shù)據(jù)的定量解釋起重要作用。重力場(chǎng)下延是經(jīng)典的不適定問(wèn)題。當(dāng)前，該問(wèn)題的求解方法除了邊界單元法[1]、樣條函數(shù)法[2-3]、等效源法[4]等少數(shù)空間域方法以外，絕大部分下延方法需要在波數(shù)域加速[5-11]，即便是基于泰勒級(jí)數(shù)[12]、中值定理[13]、Milne法[14]等原理的下延方法，也需要利用波數(shù)域方法快速獲得垂向?qū)?shù)。換言之，目前絕大部分下延方法均基于Fourier變換。Fourier變換要求輸入的重力數(shù)據(jù)無(wú)空白區(qū)，且數(shù)據(jù)長(zhǎng)度滿(mǎn)足快速Fourier變換(FFT)對(duì)數(shù)據(jù)長(zhǎng)度(2的冪次方)的要求。顯然，實(shí)測(cè)重力數(shù)據(jù)很難完全滿(mǎn)足上述兩個(gè)條件。因此，對(duì)實(shí)測(cè)重力數(shù)據(jù)進(jìn)行下延處理前，必須對(duì)數(shù)據(jù)空白部分進(jìn)行插值，且波數(shù)域下延方法還需對(duì)數(shù)據(jù)進(jìn)行擴(kuò)邊，這將直接關(guān)系到重力數(shù)據(jù)下延的精度。重力數(shù)據(jù)的插值和擴(kuò)邊對(duì)應(yīng)于信號(hào)的重構(gòu)和外推。信號(hào)重構(gòu)和外推是信號(hào)處理的逆問(wèn)題，也是不適定的[15]。因此，可以將重力數(shù)據(jù)的插值、擴(kuò)邊和向下延拓統(tǒng)一考慮。

凸集投影(Projection onto Convex Sets，POCS)是帶限信號(hào)重構(gòu)的重要理論方法之一，在圖像重建[16-17]和地震數(shù)據(jù)插值[18-21]等領(lǐng)域應(yīng)用廣泛，具有原理簡(jiǎn)單、易于執(zhí)行和精度較高的特點(diǎn)。文獻(xiàn)[22]基于凸集投影方法進(jìn)行了重磁網(wǎng)格數(shù)據(jù)插值重建，結(jié)果表明該方法的插值精度優(yōu)于克里金方法和反距離加權(quán)法，與最小曲率法相當(dāng)。實(shí)際工程往往是計(jì)算實(shí)測(cè)重力數(shù)據(jù)與模型正演數(shù)據(jù)的差獲得目標(biāo)體的重力異常。因此，實(shí)測(cè)重力異常數(shù)據(jù)可以視為帶限信號(hào)[23]，滿(mǎn)足采用凸集投影方法處理的條件。本文基于凸集投影原理，將重力數(shù)據(jù)的插值、擴(kuò)邊與下延這三個(gè)不適定問(wèn)題統(tǒng)一考慮，提出重力數(shù)據(jù)同時(shí)填充、擴(kuò)邊和下延的一體化方法。

1 原理方法

1.1 基于POCS的重力插值擴(kuò)邊方法

POCS對(duì)網(wǎng)格化后的重力數(shù)據(jù)缺失部分的插值擴(kuò)邊，其原理與圖像缺失重建方法類(lèi)似，主要基于Gerchberg-Saxton迭代計(jì)算，每次迭代計(jì)算的核心是二維Fourier變換。通過(guò)二維Fourier正變換將數(shù)據(jù)從空間域轉(zhuǎn)換到波數(shù)域，給定一個(gè)閾值，保留大于和等于該閾值的譜分量，其余譜分量置零；然后，對(duì)保留的譜分量做二維Fourier反變換，并將原始已知數(shù)據(jù)重置回反變換后的數(shù)據(jù)中；照此過(guò)程逐次迭代，并逐漸減小每次迭代的閾值，以便恢復(fù)缺失數(shù)據(jù)和邊界的更多細(xì)節(jié)信息，直至達(dá)到設(shè)置的最大迭代次數(shù)即完成插值擴(kuò)邊。從以上的方法實(shí)現(xiàn)過(guò)程描述可以看出，POCS方法是利用整個(gè)數(shù)據(jù)的譜(能量)進(jìn)行插值和擴(kuò)邊的。

設(shè)含缺失數(shù)據(jù)的重力網(wǎng)格數(shù)據(jù)為g(x,y)，則其POCS插值擴(kuò)邊計(jì)算公式為

gk(x,y)=g(x,y)+MF-1TkFgk-1(x,y)

k=1,2,…,K

(1)

式中：gk(x,y)表示第k次迭代插值擴(kuò)邊后的數(shù)據(jù)；M為采樣矩陣，矩陣元素為0或1，0代表該點(diǎn)有數(shù)據(jù)，無(wú)需插值，為1則代表該點(diǎn)無(wú)數(shù)據(jù)，需要進(jìn)行插值； F和F-1分別表示Fourier正、反變換；K表示最大迭代次數(shù)；Tk表示閾值矩陣，其元素滿(mǎn)足

(2)

式中：Sk-1(u,v)表示第k-1次迭代得到的插值擴(kuò)邊數(shù)據(jù)的頻譜， 滿(mǎn)足Sk-1(u,v)=Fgk-1(x,y)，其中u和v分別是x和y方向的波數(shù)；pk∈{p1,p2,…,pK}表示第k次迭代的閾值，滿(mǎn)足p1>p2>…>pK和min{|Sk-1(u,v)|}

Tk(u,v)實(shí)質(zhì)等同于理想低通濾波器。為減少參數(shù)個(gè)數(shù)，將Tk(u,v)設(shè)置為

(3)

式中：ck表示第k次迭代的截止波數(shù)，滿(mǎn)足c1≤c2≤…≤cK和1

(4)

顯然，式(3)與式(2)等價(jià)。式(3)只是將式(2)的閾值pk轉(zhuǎn)化成截止波數(shù)ck。同理，最大截止波數(shù)cK的大小與原始數(shù)據(jù)的噪聲水平相關(guān)。

1.2 重力數(shù)據(jù)向下延拓原理

若觀測(cè)平面(z=0,z軸向下為正)上的重力場(chǎng)g(x,y)已知，則由g(x,y)求z=d(d>0,場(chǎng)源深度大于d)平面上的重力數(shù)據(jù)f(x,y)稱(chēng)為重力場(chǎng)的向下延拓

(5)

G(u,v)=H(u,v)F(u,v)

(6)

針對(duì)該不適定性問(wèn)題，Tikhonov正則化和各類(lèi)迭代法是常用的處理方法。這些方法將對(duì)應(yīng)的算子變換到波數(shù)域，其實(shí)質(zhì)是進(jìn)行低通濾波。以Tikhonov正則化為例，其對(duì)應(yīng)的波數(shù)域正則化算子為[6,9]

(7)

式中α為正則參數(shù)。該正則化方法的實(shí)質(zhì)是利用正則化低通濾波算子Tα壓制高頻噪聲的影響。同理，也可以利用類(lèi)似式(3)的理想低通濾波器實(shí)現(xiàn)穩(wěn)定的向下延拓

(8)

式中的TDC(u,v)與式(3)中的理想濾波器的區(qū)別在于其閾值為cDC，而式(3)濾波器的閾值為cK。

式(8)可稱(chēng)為譜截?cái)嗾齽t化，其作用等同于截?cái)嗥娈愔捣纸庹齽t化(Truncated Singular Value Decomposition，TSVD)方法[24]。TDC(u,v)的閾值cDC對(duì)應(yīng)于Tikhonov的正則參數(shù)α。用于下延的理想低通濾波器TDC(u,v)只需要濾除高頻噪聲即可確保向下延拓穩(wěn)定，因此，正則參數(shù)cDC一般比較大且主要與數(shù)據(jù)的噪聲水平相關(guān)。上一節(jié)的插值擴(kuò)邊迭代中，最大截止波數(shù)cK也與原始數(shù)據(jù)的噪聲水平相關(guān)，由此，將正則參數(shù)cDC設(shè)置為最大截止波數(shù)cK具有合理性，這是因?yàn)椴逯岛蛿U(kuò)邊同樣不需要噪聲譜。此外，現(xiàn)有確定正則參數(shù)的方法比較成熟，故將正則參數(shù)cDC設(shè)置為cK，即利用正則參數(shù)cDC確定式(3)中插值和擴(kuò)邊迭代過(guò)程的最大截止波數(shù)cK。因此經(jīng)插值和擴(kuò)邊后的下延場(chǎng)為

(9)

這樣，通過(guò)理想低通濾波器以及截止波數(shù)與正則參數(shù)的關(guān)系實(shí)現(xiàn)了重力數(shù)據(jù)的插值、擴(kuò)邊和下延一體化。

1.3 正則化參數(shù)的選取

正則化的關(guān)鍵在于最優(yōu)正則參數(shù)的選取，其關(guān)系到正則化的精度和計(jì)算效率。常用的正則參數(shù)選取方法有L—曲線(xiàn)法、GCV法、偏差準(zhǔn)則等[15]，本文采用L—曲線(xiàn)法。

常規(guī)L—曲線(xiàn)法使用對(duì)數(shù)尺度描述正則解的范數(shù)‖fcK(x,y)‖和殘差范數(shù)‖g(x,y)-h(x,y)fcK(x,y)‖。Hansen定義L—曲線(xiàn)的隅角為最大曲率，但最大曲率的求解相對(duì)復(fù)雜，本文采用極小化泛函[25]

h(x,y)fck(x,y)‖×‖fck(x,y)‖]

(10)

確定最優(yōu)正則參數(shù)cK。然后結(jié)合迭代次數(shù)K，閾值ck按線(xiàn)性方式從c1遞增到最大值cK。本文所述的一體化方法流程可歸納為圖1。

圖1 頻率域POCS迭代插值、擴(kuò)邊、下延一體化流程

2 理論模型實(shí)驗(yàn)

構(gòu)建一個(gè)由位于不同深度、大小不同的5個(gè)球體組合而成的理論球體模型(圖2a)，各球體的具體參數(shù)見(jiàn)表1。

表1 理論模型球體參數(shù)

設(shè)觀測(cè)點(diǎn)、線(xiàn)距均為50m，正演計(jì)算地面(0高度，201×201個(gè)觀測(cè)點(diǎn)，10km×10km觀測(cè)范圍)的重力場(chǎng)，如圖2b所示。以地面之上1km高度處(256×256個(gè)觀測(cè)點(diǎn)，12.8km×12.8km觀測(cè)范圍)的重力數(shù)據(jù)為觀測(cè)數(shù)據(jù)，并加入均值為零、方差為0.01mGal2的高斯白噪聲(信噪比為40.51dB)模擬實(shí)際情況，其結(jié)果如圖2c所示。為檢驗(yàn)方法的插值和擴(kuò)邊能力，將圖2c所示數(shù)據(jù)的內(nèi)部部分?jǐn)?shù)據(jù)(40×30個(gè)數(shù)據(jù)點(diǎn))和邊界(四個(gè)邊界各28個(gè)數(shù)據(jù)點(diǎn)，共25536個(gè)數(shù)據(jù)點(diǎn))“挖去”構(gòu)成空白區(qū)(圖2d)。采用本文一體化方法，對(duì)圖2d所示的重力數(shù)據(jù)進(jìn)行填充、擴(kuò)邊后，下延1km(20倍網(wǎng)格距)到地面，并用1km高度處的理論重力數(shù)據(jù)檢驗(yàn)本文方法的填充和擴(kuò)邊精度，以圖2b所示的理論地面重力數(shù)據(jù)檢驗(yàn)方法的下延精度。

為定量分析算法的精度，填充、擴(kuò)邊和下延精度都以均方根誤差(Root Mean Square Error,RMSE)進(jìn)行定量分析

(11)

式中：Dc(i)和Dt(i)分別表示計(jì)算值和理論值；Q表示參與誤差統(tǒng)計(jì)的數(shù)據(jù)個(gè)數(shù)。

根據(jù)圖1的流程，首先利用式(10)求解截止波數(shù)cK(圖2e)，最小值cK=7。選定迭代次數(shù)K=100，然后按照算法步驟進(jìn)行迭代計(jì)算。圖2f為填充、擴(kuò)邊和向下延拓均方根誤差隨迭代次數(shù)k的變化曲線(xiàn)。由圖可知，本文算法具有收斂性，最終的填充、擴(kuò)邊和下延的均方根誤差分別為0.04、0.36、1.43mGal。最終的插值和擴(kuò)邊結(jié)果見(jiàn)圖3a，與真實(shí)值在數(shù)據(jù)缺失處的殘差見(jiàn)圖3b。圖3a的下延結(jié)果見(jiàn)圖3c，其與真實(shí)值(圖2b)的殘差見(jiàn)圖3d。

為了對(duì)比分析，采用經(jīng)典的最小曲率插值法、余弦函數(shù)擴(kuò)邊法和Tikhonov正則化下延法組合方法對(duì)圖2d所示的數(shù)據(jù)分步進(jìn)行插值、擴(kuò)邊和下延，并與本文方法的處理結(jié)果進(jìn)行對(duì)比。最小曲率法的插值在軟件Golden Surfer上實(shí)現(xiàn)，其最大殘差和最大迭代次數(shù)分別選定為1.0×10-5和1.0×10-5； 余弦函數(shù)擴(kuò)邊法采用文獻(xiàn)[26]提供的“taper2d”程序；Tikhonov正則化參數(shù)選定為使下延均方根誤差最小的最優(yōu)正則化參數(shù)，經(jīng)對(duì)比計(jì)算選定為0.005。插值擴(kuò)邊結(jié)果見(jiàn)圖4a，與真實(shí)值在數(shù)據(jù)缺失處的殘差見(jiàn)圖4b。圖4a的下延(下延1km)結(jié)果見(jiàn)圖4c，與真實(shí)值(圖2b)的殘差見(jiàn)圖4d。經(jīng)統(tǒng)計(jì)，采用經(jīng)典組合方法進(jìn)行插值、擴(kuò)邊和下延結(jié)果的均方根誤差分別為0.15、1.88、2.23mGal。

對(duì)比本文方法與經(jīng)典組合方法的插值、擴(kuò)邊、下延結(jié)果及其殘差和均方根誤差可知，本文方法的插值、擴(kuò)邊和下延結(jié)果與理論值更接近，殘差和均方根誤差更小。

圖2 理論模型及重力數(shù)據(jù)計(jì)算結(jié)果

圖3 理論模型本文方法插值、擴(kuò)邊、下延的結(jié)果及誤差

圖4 理論模型經(jīng)典組合方法插值、擴(kuò)邊、下延的結(jié)果及誤差

3 實(shí)測(cè)航空重力數(shù)據(jù)實(shí)驗(yàn)

為檢驗(yàn)本文方法在實(shí)際數(shù)據(jù)處理中的實(shí)用性，對(duì)美國(guó)地質(zhì)勘探局(United States Geological Survey,USGS)2006年在阿富汗實(shí)測(cè)的航空重力數(shù)據(jù)進(jìn)行計(jì)算。

該航空重力測(cè)量數(shù)據(jù)經(jīng)各項(xiàng)處理并最終歸算至離地7000m高度，形成網(wǎng)格距為1000m的布格重力異常(圖5a)。由圖可知，實(shí)際航空重力測(cè)區(qū)范圍很不規(guī)則。將圖5a的布格重力異常采用本文方法進(jìn)行填充、擴(kuò)邊(擴(kuò)邊至1024×1024個(gè)點(diǎn))并下延至地面。首先，采用式(10)求得截止波數(shù)cK=70，如圖5b所示。選定迭代次數(shù)K=100，按照?qǐng)D1所示的步驟，圖5a數(shù)據(jù)經(jīng)填充和擴(kuò)邊的結(jié)果見(jiàn)圖6a，圖6a下延至地面的結(jié)果見(jiàn)圖6b。將圖6b的布格重力異常場(chǎng)上延7000m恢復(fù)至原數(shù)據(jù)高度，結(jié)果如圖6c所示(僅展示有原始數(shù)據(jù)的區(qū)域)。圖6c與圖5a所示真實(shí)值的殘差見(jiàn)圖6d，其均方根誤差僅0.05mGal，驗(yàn)證了本文方法的有效性和高精度。

圖5 實(shí)測(cè)航空重力數(shù)據(jù)實(shí)驗(yàn)

圖6 實(shí)測(cè)數(shù)據(jù)本文方法插值、擴(kuò)邊、下延結(jié)果和誤差

同理論模型實(shí)驗(yàn)一樣，同樣采用經(jīng)典的最小曲率插值法、余弦函數(shù)擴(kuò)邊法和Tikhonov正則化下延法對(duì)圖5a所示數(shù)據(jù)分步進(jìn)行插值、擴(kuò)邊和下延。最小曲率插值法采用與理論模型計(jì)算部分一致的參數(shù)；Tikhonov正則化下延法的正則化參數(shù)采用使下延結(jié)果經(jīng)上延后均方根誤差最小的正則參數(shù)，經(jīng)計(jì)算確定為0.001。采用經(jīng)典組合方法對(duì)圖5a數(shù)據(jù)進(jìn)行插值和擴(kuò)邊，其結(jié)果見(jiàn)圖7a～圖7d，其中圖7c僅展示有原始數(shù)據(jù)的區(qū)域。圖7c與圖5a的均方根誤差達(dá)5.94mGal(圖7d)。

圖7 實(shí)測(cè)數(shù)據(jù)經(jīng)典方法插值、擴(kuò)邊、下延結(jié)果和誤差

對(duì)比圖6與圖7可知，本文方法的插值、擴(kuò)邊和下延結(jié)果遠(yuǎn)優(yōu)于經(jīng)典組合方法，其計(jì)算殘差和均方誤差也遠(yuǎn)低于后者。

4 結(jié)束語(yǔ)

本文針對(duì)實(shí)測(cè)重力數(shù)據(jù)不規(guī)則且存在空白區(qū)的實(shí)際情況，結(jié)合重力數(shù)據(jù)填充、擴(kuò)邊和向下延拓同屬不適定反問(wèn)題的理論基礎(chǔ)，基于帶限信號(hào)重建的經(jīng)典算法——凸集投影算法，提出了重力插值、擴(kuò)邊、下延的一體化方法，并利用理論模型和實(shí)測(cè)數(shù)據(jù)驗(yàn)證了方法的收斂性和有效性。從算法流程來(lái)看，方法采用FFT運(yùn)算，在波數(shù)域僅涉及理想低通濾波器，且只需要預(yù)先設(shè)定迭代次數(shù)，具有原理簡(jiǎn)單、實(shí)際操作方便和運(yùn)算高效的特點(diǎn)。理論模型數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)處理結(jié)果顯示，本文方法的擴(kuò)邊效果光滑、無(wú)畸變且插值和向下延拓結(jié)果精度較高，優(yōu)于經(jīng)典的組合方法。因此，本文方法適宜于實(shí)測(cè)規(guī)則網(wǎng)格重力數(shù)據(jù)的插值、擴(kuò)邊和向下延拓。

目前該方法的局限在于只能處理規(guī)則網(wǎng)格數(shù)據(jù)，基于該方法的不均勻采樣重力數(shù)據(jù)的處理是后續(xù)進(jìn)一步研究的內(nèi)容。