基于凸集投影的重力數據擴充下延一體化方法

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

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

0 引言

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

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

1 原理方法

1.1 基于POCS的重力插值擴邊方法

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

設含缺失數據的重力網格數據為g(x,y)，則其POCS插值擴邊計算公式為

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

k=1,2,…,K

(1)

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

(2)

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

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

(3)

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

(4)

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

1.2 重力數據向下延拓原理

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

(5)

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

(6)

針對該不適定性問題，Tikhonov正則化和各類迭代法是常用的處理方法。這些方法將對應的算子變換到波數域，其實質是進行低通濾波。以Tikhonov正則化為例，其對應的波數域正則化算子為[6,9]

(7)

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

(8)

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

式(8)可稱為譜截斷正則化，其作用等同于截斷奇異值分解正則化(Truncated Singular Value Decomposition，TSVD)方法[24]。TDC(u,v)的閾值cDC對應于Tikhonov的正則參數α。用于下延的理想低通濾波器TDC(u,v)只需要濾除高頻噪聲即可確保向下延拓穩定，因此，正則參數cDC一般比較大且主要與數據的噪聲水平相關。上一節的插值擴邊迭代中，最大截止波數cK也與原始數據的噪聲水平相關，由此，將正則參數cDC設置為最大截止波數cK具有合理性，這是因為插值和擴邊同樣不需要噪聲譜。此外，現有確定正則參數的方法比較成熟，故將正則參數cDC設置為cK，即利用正則參數cDC確定式(3)中插值和擴邊迭代過程的最大截止波數cK。因此經插值和擴邊后的下延場為

(9)

這樣，通過理想低通濾波器以及截止波數與正則參數的關系實現了重力數據的插值、擴邊和下延一體化。

1.3 正則化參數的選取

正則化的關鍵在于最優正則參數的選取，其關系到正則化的精度和計算效率。常用的正則參數選取方法有L—曲線法、GCV法、偏差準則等[15]，本文采用L—曲線法。

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

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

(10)

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

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

2 理論模型實驗

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

表1 理論模型球體參數

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

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

(11)

式中：Dc(i)和Dt(i)分別表示計算值和理論值；Q表示參與誤差統計的數據個數。

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

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

對比本文方法與經典組合方法的插值、擴邊、下延結果及其殘差和均方根誤差可知，本文方法的插值、擴邊和下延結果與理論值更接近，殘差和均方根誤差更小。

圖2 理論模型及重力數據計算結果

圖3 理論模型本文方法插值、擴邊、下延的結果及誤差

圖4 理論模型經典組合方法插值、擴邊、下延的結果及誤差

3 實測航空重力數據實驗

為檢驗本文方法在實際數據處理中的實用性，對美國地質勘探局(United States Geological Survey,USGS)2006年在阿富汗實測的航空重力數據進行計算。

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

圖5 實測航空重力數據實驗

圖6 實測數據本文方法插值、擴邊、下延結果和誤差

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

圖7 實測數據經典方法插值、擴邊、下延結果和誤差

對比圖6與圖7可知，本文方法的插值、擴邊和下延結果遠優于經典組合方法，其計算殘差和均方誤差也遠低于后者。

4 結束語

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

目前該方法的局限在于只能處理規則網格數據，基于該方法的不均勻采樣重力數據的處理是后續進一步研究的內容。