基于凸集投影的重力同時(shí)填充擴(kuò)邊和去噪方法

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

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

0 引言

實(shí)測重力數(shù)據(jù)不僅含有大量噪聲，且可能由于各種因素的限制，常常會存在部分缺失。波數(shù)域重力數(shù)據(jù)處理過程中，要求輸入的重力數(shù)據(jù)無空白，且數(shù)據(jù)長度也要盡量滿足快速Fourier變換對數(shù)據(jù)長度(2的冪次方)的要求。顯然，實(shí)測數(shù)據(jù)很難滿足上述條件。因此，在進(jìn)行波數(shù)域數(shù)據(jù)處理前，必須對重力數(shù)據(jù)空白部分進(jìn)行填充、對邊界進(jìn)行擴(kuò)邊，這將直接關(guān)系到數(shù)據(jù)的處理精度[1-2]。

常用的重力數(shù)據(jù)去噪方法有Fourier域?yàn)V波法、優(yōu)化濾波法、經(jīng)驗(yàn)?zāi)B(tài)分解法、等效源法、小波分解法、滑動窗口平均法等。Fourier域?yàn)V波法的基本原理是在波數(shù)域提取并分離干擾信號，其主要問題是信號和噪聲成份不易區(qū)分。由此，基于重力位場物理特性的等效層法[3]、等效源法[4]和基于奇異值分解、經(jīng)驗(yàn)?zāi)B(tài)分解[5]等的信號處理方法被引入重力數(shù)據(jù)去噪中，其目的就是有效區(qū)分有用信號和噪聲。小波分解法[6]在各類信號的去噪中都得到了廣泛應(yīng)用，但小波去噪的效果與小波基函數(shù)以及閾值的選擇有較大關(guān)系。

對缺失數(shù)據(jù)進(jìn)行填充重建，理論上所有的數(shù)據(jù)插值方法都可以使用，經(jīng)典的插值方法有線性插值、樣條插值、 Kriging法、最小曲率法等。這些方法均有各自的理論適應(yīng)性和計(jì)算特點(diǎn)，它們的操作共同點(diǎn)是： 首先將測量數(shù)據(jù)與坐標(biāo)結(jié)合，以網(wǎng)格中已知數(shù)據(jù)點(diǎn)為型值點(diǎn)，空缺位置坐標(biāo)為待插值點(diǎn)，進(jìn)行插值[7]。

擴(kuò)邊實(shí)質(zhì)上是信號的外推過程，常用的擴(kuò)邊方法有補(bǔ)零擴(kuò)邊法、對折擴(kuò)邊法、區(qū)域場擴(kuò)邊法、三方向擴(kuò)邊法、余弦擴(kuò)邊法和最小曲率法等[1-2]。余弦擴(kuò)邊方法不能很好地反映實(shí)際數(shù)據(jù)的區(qū)域變化特征，并在實(shí)際數(shù)據(jù)和擴(kuò)邊數(shù)據(jù)的銜接處導(dǎo)數(shù)不連續(xù)，從而產(chǎn)生較嚴(yán)重的Gibbs效應(yīng)[1]。最小曲率法是目前重力數(shù)據(jù)插值和擴(kuò)邊中最受歡迎的方法，但該方法主要考慮曲面的光滑性(在沒有數(shù)據(jù)點(diǎn)的情況下具有最小曲率表面的是球形)，因此，不能得到精確的插值結(jié)果。當(dāng)網(wǎng)格化間距小于實(shí)際的網(wǎng)格間距時(shí)，稀疏數(shù)據(jù)控制區(qū)域的插值精度會快速下降[8]。

凸集投影(Projection onto Convex Sets，POCS)方法是帶限信號重構(gòu)的重要理論方法之一，在圖像重建[9]和地震數(shù)據(jù)插值[10-15]等領(lǐng)域應(yīng)用廣泛，且具有原理簡單、易于執(zhí)行和精度較高的特點(diǎn)。在實(shí)際工程應(yīng)用中，往往是通過實(shí)際測量獲得的重力數(shù)據(jù)減去由模型計(jì)算的正常場數(shù)據(jù)而獲得重力異常數(shù)據(jù)。因此，實(shí)測重力異常數(shù)據(jù)可以視為帶限信號[16]，滿足采用凸集投影方法處理的條件。在重力數(shù)據(jù)處理方面，閆浩飛等[7]利用POCS方法實(shí)現(xiàn)了重力數(shù)據(jù)插值重建；曾小牛等[17]提出一種基于凸集投影原理的重力數(shù)據(jù)填充、擴(kuò)邊、下延一體化方法。Ji等[18]結(jié)合預(yù)條件共軛梯度迭代法和POCS方法實(shí)現(xiàn)了磁場參數(shù)的正則化反演。本文基于凸集投影原理，將重力數(shù)據(jù)的填充、去噪和擴(kuò)邊這三個(gè)問題統(tǒng)一考慮，提出基于凸集投影的規(guī)則網(wǎng)格重力數(shù)據(jù)同時(shí)填充、擴(kuò)邊和去噪方法。

1 方法原理

1.1 基于POCS的重力數(shù)據(jù)同時(shí)擴(kuò)充和去噪

設(shè)含噪且存在數(shù)據(jù)缺失的觀測重力網(wǎng)格數(shù)據(jù)為gobs，則基于凸集投影的填充和擴(kuò)邊計(jì)算可表示為

gk=gobs+(I-S)F-1TkFgk-1

k=1,2,3,…,K

(1)

式中：gk表示第k次迭代擴(kuò)邊填充后的數(shù)據(jù)；I為單位矩陣；S為采樣矩陣，由0和1組成，1代表該點(diǎn)有數(shù)據(jù)，無需插值，反之為0，則代表該點(diǎn)需要進(jìn)行插值； F和F-1分別表示Fourier正、反變換；K表示最大迭代次數(shù)；Tk表示閾值矩陣，其元素滿足

(2)

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

顯然，式(1)的每次迭代都將含噪聲的原始數(shù)據(jù)gobs帶入了插值結(jié)果，因此，傳統(tǒng)的POCS方法對含噪數(shù)據(jù)的插值效果不理想。為此，曹靜杰等[12]提出了一種用于地震數(shù)據(jù)重建的改進(jìn)凸集投影方法，本文將其應(yīng)用到重力數(shù)據(jù)波數(shù)域的同時(shí)去噪和擴(kuò)充，可表示為

gk=F-1TkF[gobs+(I-S)gk-1]

k=1,2,3,…,K

(3)

由式(1)和式(3)可知，該類迭代方法的核心是確定波數(shù)域閾值pk。在基于POCS方法的地震和圖像信號去噪重建中，一般先確定閾值pk的最大值和最小值，再由線性[10]、指數(shù)或數(shù)據(jù)驅(qū)動[11]模式完成從最大值到最小值的下降過程。實(shí)際應(yīng)用中，pk的最大值和最小值的確定與原始數(shù)據(jù)的噪聲水平有關(guān)，對結(jié)果影響較大，需要謹(jǐn)慎選擇[10-12]。相對地震和圖像信號來說，重力數(shù)據(jù)變化較平緩，連續(xù)性較好，頻譜Sk-1(u,v)相對集中。Tk(u,v)的實(shí)質(zhì)是理想低通濾波器，針對重力數(shù)據(jù)頻譜特點(diǎn)及頻譜的物理特征，本文將Tk(u,v)修改為

(4)

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

(5)

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

1.2 基于徑向平均功率譜的最大截止波數(shù)cK確定方法

在基于POCS方法的地震和圖像信號去噪、重建過程中，閾值最大值和最小值的確定比較隨意，缺乏一定的準(zhǔn)則[10-12]。對應(yīng)于式(4)的理想低通濾波器，截止波數(shù)cK顯然最好為有效信號與噪聲的分界點(diǎn)，這樣迭代過程既能保留有效信號應(yīng)用于擴(kuò)充和填充，又能避免噪聲的干擾。

Spector等[19]提出的徑向平均功率譜能夠反映重磁位場的頻譜特征，并被廣泛應(yīng)用于重磁位場數(shù)據(jù)的處理[3,20-22]。一個(gè)假想的徑向平均功率譜如圖1所示，則垂直疊加的深源和淺源產(chǎn)生重力異常的徑向平均功率譜可以采用下式進(jìn)行擬合

(6)

式中：A為能譜幅值；h為源深度。如果觀測噪聲為白噪，則與重力信號不相關(guān)的白噪聲功率譜應(yīng)該為常數(shù)，即噪聲大致對應(yīng)于徑向平均功率譜的水平部分。這樣，對于重力異常的徑向平均功率譜，存在一個(gè)截止波數(shù)將重力異常的信號譜與噪聲譜大致分開。由此，可將迭代法的最大截止波數(shù)cK設(shè)置為由徑向平均功率譜擬合所確定的劃分信號譜與噪聲譜的截止波數(shù)。顯然，這樣的設(shè)定具有明確的物理意義，即只采用有效信號的頻譜進(jìn)行擴(kuò)充，同時(shí)又能消除噪聲譜。

通過徑向平均功率譜可獲得截止波數(shù)cK，結(jié)合迭代次數(shù)K，迭代過程的截止波數(shù)ck依線性遞增的模式從1增大到最終的截止波數(shù)cK。本文所提的迭代方法可歸納為圖2。

圖1 基于徑向平均功率譜確定最大截止波數(shù)cK示意圖

圖2 本文填充、擴(kuò)邊及去噪方法流程

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

理論模型由位于不同深度、不同大小的4個(gè)長方體組合而成(圖3a)，各長方體的參數(shù)見表1。

表1 理論模型中的異常體參數(shù)

點(diǎn)、線距均設(shè)為50m，平面剖分網(wǎng)格大小為256×256，正演計(jì)算地面(z=0)的重力響應(yīng)如圖3b所示。為檢驗(yàn)方法的去噪能力，給圖3b所示的重力數(shù)據(jù)增加零均值、標(biāo)準(zhǔn)差為1.00mGal的高斯白噪聲(圖3c)。為檢驗(yàn)方法的填充和擴(kuò)邊能力，首先，隨機(jī)去掉圖3c所示數(shù)據(jù)總量的1/4，并在內(nèi)部和邊界“挖去”部分?jǐn)?shù)據(jù)，得到缺失數(shù)據(jù)(圖3d)。圖3c內(nèi)部“挖去”的數(shù)據(jù)選取在異常梯度帶及邊部異常平緩處，這樣的選擇具有一定的典型性。采用本文提出的擴(kuò)邊、填充及去噪法，對圖3d所示的加噪且有缺失的重力數(shù)據(jù)進(jìn)行處理，并與圖3b所示的真實(shí)重力數(shù)據(jù)進(jìn)行對比。采用均方根誤差(Root Mean Square Error，RMSE)進(jìn)行定量分析

(7)

式中：Dc(i)和Dt(i)分別表示計(jì)算值和理論值；Q表示參與誤差計(jì)算的總數(shù)據(jù)量。

根據(jù)圖2的算法步驟，選定最大迭代次數(shù)K，先計(jì)算圖3d重力數(shù)據(jù)的徑向平均功率譜(圖3e)。根據(jù)形狀特征將功率譜分為三段： 1～5(頻段1)、5～18(頻段2)和18～129(頻段3)。頻段1主要對應(yīng)低頻異常信息，頻段2主要對應(yīng)中高頻信息，而頻段3的功率譜主要對應(yīng)高斯白噪聲。因此，確定截止波數(shù)cK為頻段2和頻段3的分界點(diǎn)17。最終的填充、去噪和擴(kuò)邊結(jié)果見圖4a中的上圖。對比圖4a中的上圖與原始重力數(shù)據(jù)(圖3b)可知，去噪、填充的結(jié)果同理論值非常接近，擴(kuò)邊的結(jié)果稍差，但形態(tài)特征非常一致，且已知數(shù)據(jù)與缺失數(shù)據(jù)的連接處光滑無畸變。

為了對比分析，分別采用經(jīng)典的最小曲率法和Kriging插值法分別與小波分解去噪、余弦函數(shù)擴(kuò)邊法組成兩組組合方法對圖3d所示數(shù)據(jù)分步進(jìn)行插值、去噪和擴(kuò)邊，并與本文方法的處理結(jié)果進(jìn)行對比。最小曲率和Kriging插值法均可在Golden Surfer 12軟件中實(shí)現(xiàn)：最小曲率插值法的最大殘差和最大迭代次數(shù)參數(shù)分別選定為1.0×10-5和1.0×105； Kri-ging插值的變異函數(shù)模型選定為Gaussian函數(shù)。小波分解采用Matlab提供的小波分析工具，分別試驗(yàn)了常用的Symlets小波系中的sym1～sym8小波和Daubechies小波系中db1～db8小波的去噪效果，選定去噪結(jié)果與理論值均方根誤差最小的sym4小波進(jìn)行最終的去噪實(shí)驗(yàn)；余弦函數(shù)擴(kuò)邊法采用Cooper等[23]提供的“taper2d”程序。

圖3 理論模型重力數(shù)據(jù)

圖4 理論模型重力數(shù)據(jù)的填充、擴(kuò)邊和去噪結(jié)果(上)及填充、去噪(中)和擴(kuò)邊(下)殘差

采用這兩組組合方法對圖3d所示數(shù)據(jù)按填充、去噪、擴(kuò)邊的順序進(jìn)行試驗(yàn)。圖4展示了本文方法填充、去噪和擴(kuò)邊殘差。最終結(jié)果與理論值(圖3b)相比較的填充、去噪和擴(kuò)邊的均方根誤差分別為： 組合方法1為0.86、0.03、4.80mGal； 組合方法2為0.22、0.03、4.80mGal。圖5是填充、去噪和擴(kuò)邊均方根誤差隨迭代次數(shù)k的變化曲線。由圖可知，算法具有收斂性，最終的填充、去噪和擴(kuò)邊均方根誤差分別為0.08、0.01、0.82mGal。由圖4和圖5可知，本文方法的填充、去噪和擴(kuò)邊結(jié)果與理論值更一致，殘差和均方根誤差也更小。

圖5 理論模型重力數(shù)據(jù)填充、擴(kuò)邊和去噪均方根誤差隨迭代次數(shù)的變化

3 實(shí)測重力數(shù)據(jù)對比實(shí)驗(yàn)

為檢驗(yàn)本文方法在實(shí)際數(shù)據(jù)處理中的實(shí)用性，采用美國地質(zhì)調(diào)查局(United States Geological Survey，USGS)公布的阿富汗均衡重力異常數(shù)據(jù)進(jìn)行對比實(shí)驗(yàn)。該重力數(shù)據(jù)融合了2006年的航空實(shí)測數(shù)據(jù)以及地面已有的重力數(shù)據(jù)，經(jīng)各項(xiàng)處理并最終歸算至離地7km的高度、數(shù)據(jù)網(wǎng)格距為1km的重力異常(圖6a)。由圖可見，融合后的重力數(shù)據(jù)存在較大面積的空白區(qū)。

3.1 局部完整數(shù)據(jù)對比實(shí)驗(yàn)

首先從原始數(shù)據(jù)中截取256km×256km大小的完整數(shù)據(jù)(圖6a黑色方框區(qū)域)進(jìn)行方法驗(yàn)證，圖6b所示即是所選區(qū)域的數(shù)據(jù)。同時(shí)，因原始數(shù)據(jù)的高頻噪聲經(jīng)過低通濾波，為模擬實(shí)際情況，對圖6a數(shù)據(jù)添加零均值、標(biāo)準(zhǔn)差為1.00mGal的高斯白噪聲，結(jié)果如圖6c。與理論模型實(shí)驗(yàn)一樣，為了檢驗(yàn)本文所提方法的填充和擴(kuò)邊能力，將圖6c所示方框內(nèi)數(shù)據(jù)隨機(jī)去掉數(shù)據(jù)總量的1/4，并在邊界和內(nèi)部“挖去”部分?jǐn)?shù)據(jù)構(gòu)成空白區(qū)(圖6d)。

采用本文提出的擴(kuò)充去噪迭代法對圖6d數(shù)據(jù)進(jìn)行處理，并對比圖6b所示的真實(shí)重力數(shù)據(jù)以檢驗(yàn)方法的填充、去噪和擴(kuò)邊精度。根據(jù)圖2所示算法步驟，首先計(jì)算圖6d所示數(shù)據(jù)的徑向平均功率譜，如圖7所示。根據(jù)功率譜形狀，將其分為1～6(頻段1)、6～18(頻段2)、18～49(頻段3)和49～129(頻段4)四個(gè)頻段。頻段1主要對應(yīng)低頻異常信息，頻段2主要對應(yīng)中高頻信息，而頻段3和頻段4主要對應(yīng)高斯白噪聲。因此，確定截止波數(shù)cK為頻段2和頻段3的分界點(diǎn)18。選定迭代次數(shù)K=200，然后按照圖2所示的算法步驟進(jìn)行迭代計(jì)算。

與理論模型實(shí)驗(yàn)過程一樣，采用經(jīng)典的最小曲率法和Kriging插值法分別與小波分解去噪、余弦函數(shù)擴(kuò)邊法組成兩組組合方法，對圖6d所示的數(shù)據(jù)分步進(jìn)行填充、去噪和擴(kuò)邊，并與本文方法的處理結(jié)果(圖8a上圖)進(jìn)行對比。經(jīng)對比Symlets小波系和Daubechies小波系試驗(yàn)，小波分解采用使去噪均方根誤差最小的sym8小波進(jìn)行去噪；兩種插值法和余弦函數(shù)擴(kuò)邊的實(shí)現(xiàn)和參數(shù)均與理論模型一致。組合方法1和2的填充、擴(kuò)邊和去噪結(jié)果分別見圖8b上圖和圖8c上圖。組合方法1處理結(jié)果(圖8b上圖)相比于圖6b的填充、去噪和擴(kuò)邊均方根誤差分別為3.36、0.20、11.75mGal； 組合方法2處理結(jié)果(圖8c上圖)相比于圖6b的填充、去噪和擴(kuò)邊均方根誤差分別為1.97、0.19、11.56mGal。圖8中圖展示了圖8上圖與圖6b相比的填充和去噪殘差，圖8下圖展示了圖8上圖與圖6b相比的擴(kuò)邊殘差。填充、去噪和擴(kuò)邊的均方根誤差隨迭代次數(shù)k的變化見圖9。由圖可知，本文方法具有收斂性，最終的填充、去噪和擴(kuò)邊的均方根誤差分別為1.22、0.15、5.37mGal。對比圖8所顯示的本文方法與經(jīng)典組合方法的填充、去噪和擴(kuò)邊結(jié)果，以及殘差和均方根誤差，可以看出，同理論數(shù)據(jù)計(jì)算結(jié)果一樣，本文方法的填充、去噪和擴(kuò)邊結(jié)果優(yōu)于組合方法，其殘差和均方根誤差也均小于組合方法。

圖6 實(shí)測重力數(shù)據(jù)及計(jì)算結(jié)果

圖7 徑向平均功率譜及擬合結(jié)果示意圖

3.2 完整數(shù)據(jù)對比實(shí)驗(yàn)

為更進(jìn)一步對比方法在數(shù)據(jù)缺失較多情況下的效果，采用提出的迭代法和兩組組合方法對圖6c所示的重力數(shù)據(jù)進(jìn)行對比實(shí)驗(yàn)。首先計(jì)算圖6c所示數(shù)據(jù)的徑向平均功率譜(圖10)?；趯β首V形狀的分析，將其分為1～12(頻段1)、12～82(頻段2)、82～270(頻段3)和270～512(頻段4)四段。頻段1主要對應(yīng)低頻異常信息，頻段2和頻段3主要對應(yīng)中高頻信息，而頻段4主要對應(yīng)高斯白噪聲。因此，確定截止波數(shù)cK為頻段3與頻段4的分界點(diǎn)270。選定迭代次數(shù)K=200。本文方法的填充、去噪和擴(kuò)邊結(jié)果見圖11a。

與前兩次實(shí)驗(yàn)過程類似，采用經(jīng)典的最小曲率法和Kriging插值法與小波分解去噪對圖6c所示的數(shù)據(jù)分步進(jìn)行填充、擴(kuò)邊和去噪。為保證收斂，采用最小曲率法插值時(shí)，插值和擴(kuò)邊的參數(shù)分別確定為1.0×10-2和1.0×105； Kriging插值法的變異函數(shù)依然為Gaussian函數(shù)。去噪方面，經(jīng)對比Symlets小波系和Daubechies小波系試驗(yàn)結(jié)果，小波分解均采用使去噪均方根誤差最小的sym8小波進(jìn)行去噪。組合方法1填充、擴(kuò)邊和去噪的結(jié)果見圖11b，小波去噪的均方根誤差為0.30mGal。組合方法2填充、擴(kuò)邊和去噪的結(jié)果見圖11c，小波去噪的均方根誤差為0.28mGal。

圖8 局部重力數(shù)據(jù)的填充、擴(kuò)邊和去噪結(jié)果(上)及填充、去噪(中)和擴(kuò)邊(下)殘差

圖9 局部重力數(shù)據(jù)去噪、填充和擴(kuò)邊均方根

對比圖11a與圖11b、圖11c可知，雖然采用最優(yōu)小波去噪的均方根誤差略小于本文方法，但本文方法的填充和擴(kuò)邊結(jié)果基本保持了真實(shí)數(shù)據(jù)的形態(tài)特征，已知數(shù)據(jù)和空白數(shù)據(jù)的連接處光滑無畸變，而最小曲率法的填充和擴(kuò)邊結(jié)果僅在靠近已知數(shù)據(jù)的邊緣形態(tài)較好，在離已知數(shù)據(jù)較遠(yuǎn)的區(qū)域則效果很差。Kriging插值法的填充結(jié)果較好，但在離已知數(shù)據(jù)較遠(yuǎn)的區(qū)域擴(kuò)邊效果同樣不甚理想。

圖10 徑向平均功率譜及擬合結(jié)果示意圖

圖11 全區(qū)重力數(shù)據(jù)實(shí)驗(yàn)結(jié)果對比

圖12為去噪均方根誤差隨迭代次數(shù)k的變化，最終的去噪均方根誤差為0.32mGal，可見算法具有收斂性。

圖12 去噪均方根誤差隨迭代次數(shù)的變化曲線

4 結(jié)論

本文針對實(shí)際重力測量獲得的重力異常數(shù)據(jù)往往存在噪聲且有空缺的實(shí)際情況，基于帶限信號重建的經(jīng)典算法——凸集投影算法，提出了一種重力數(shù)據(jù)同時(shí)填充、擴(kuò)邊和去噪的迭代法，并利用理論模型和實(shí)測數(shù)據(jù)驗(yàn)證了方法的收斂性和有效性。從算法流程來看，迭代法采用FFT運(yùn)算、在波數(shù)域僅涉及理想低通濾波器，因此只需要預(yù)先設(shè)定迭代次數(shù)以及依據(jù)擬合徑向平均功率確定截止波數(shù)即可完成迭代，具有原理簡單、實(shí)際操作方便和運(yùn)算高效的特點(diǎn)。從理論和實(shí)測數(shù)據(jù)處理的實(shí)驗(yàn)效果來看，在擴(kuò)邊和補(bǔ)空的銜接處，數(shù)據(jù)光滑無畸變，方法取得了較好的填充、擴(kuò)邊和去噪效果。本文方法的處理結(jié)果優(yōu)于由最小曲率法和Kriging插值與小波分解去噪、余弦擴(kuò)邊方法組成的組合方法。