基于基追蹤的位場數據稀疏反演

葉麗霞，張玉潔

(1.中國地質大學 數學與物理學院，湖北 武漢 430074；2.中國地質大學 地質探測與評估教育部重點實驗室, 湖北 武漢 430074)

1 引 言

重磁場物性反演作為了解地下密度和磁化率分布的重要手段，在礦產、石油、工程勘探等很多領域已經有了廣泛的應用。相比于參數反演，物性反演具有其特定的優勢，它能夠較好地提供異常體的位置、形狀等信息，從而對于重磁數據的解釋具有重要意義，也是該領域的一個熱點。重磁勘探的主要目的是通過測量重力場源和地磁場源的數據異常，并對測量的數據進行反演和解釋，從而得到地下異常源的密度分布和磁化率分布[1]。磁場物性反演指的是將地下空間的異常體看作是一系列較強或較弱的磁化強度所組成的異常區域，根據地面觀測的磁異常反演地下磁化強度分布時，將地下空間介質劃分為一系列塊體單元，當塊體單元足夠小時，假設每個塊體單元內部磁化強度均勻分布，允許不同塊體單元具有不同的磁化強度。在反演的過程中，磁化強度是唯一需要反演的參數。這種反演方法與層析成像方法相似，因而也稱磁化強度成像。

由于觀測數據數目遠小于模型參數數目，會導致反演問題的多解性，從而造成分辨率低的現象?；谶@種現象，人們提出了許多反演理論，其中包括聯合反演和約束反演[2]。Pilkington[3]將預優共軛梯度算法用來求解磁數據離散反演問題，減少了計算量，提高了反演效率。Li等[4]在目標函數中加入了深度加權函數，以此來克服反演結果的“趨膚效應”。這套反演理論雖然穩定性好，但是得到的反演結果比較發散。此外，他們還將地表全磁數據與鉆孔三分量磁數據進行反演，這種聯合反演方法有效地提高了反演的分辨率[5]。Sun等[6]在反演算法中加入了已知的物性信息，進一步改進了反演模型。此外，反演的方法還包括基于L2范數的正則化光滑反演方法[7]、基于L0范數的正則化聚焦反演方法[8]、二元和多元聚焦反演方法[9]等?？偠灾?，這些方法都是為了解決重磁場反演的多解性問題以及分辨率低的問題。

傳統重磁反演可以采用觀測數據和模型參數的L2范數極小化，L2范數易于處理，容易求解。但L2范數獲得的是非稀疏解，在稀疏重構方面所得的解稀疏性差，會增大非零值的范圍，降低了反演的分辨率。壓縮感知[10]理論能很好地刻畫稀疏性，該理論指出，只要信號是稀疏的或者能夠被稀疏表示，那么就可以無失真地重構出原始信號。在壓縮感知的推動下，Lp范數(0≤p≤1)優化在理論和算法上有了巨大的發展[11,12]，這為求解Lp范數稀疏反演問題提供了理論基礎。Li等[13]將物性上下界不等式約束添加到目標函數中，用內點法處理優化問題，對三維位場數據進行Lp范數稀疏反演。此外，零范數稀疏恢復在地球物理領域也得到了廣泛應用[14,15]。Meng等[16]提出了一種基于稀疏恢復的三維重力反演方法，選取L0范數作為目標函數，然后用近似零范數函數迭代求解。在過去的十幾年里，壓縮感知理論得到了較為成熟地發展，在很多領域有了廣泛應用。壓縮感知用于地震數據反演已經取得了較好的效果[17,18]，越來越多的研究者嘗試將壓縮感知用于重磁數據反演。

本文提出基于基追蹤的重磁場物性稀疏反演方法，用壓縮感知技術代替傳統的反演方法，直接從反演模型出發，并且引入物性上下界約束，構建目標函數，將其轉化為有約束的優化問題，然后利用基追蹤對數障礙規劃算法進行求解。

2 基本原理

2.1 正演

重磁異常的正演計算表達式可表示為：

Gm+e=d

(1)

其中，G是核矩陣，維度為M×N(M

2.2 反演

重磁場物性反演的過程就是利用觀測數據向量和核矩陣來估計模型參數向量。傳統的反演方法得到的是分辨率比較低的非稀疏解，而大多數情況下，物性分布本身就是稀疏的，因而本文采用壓縮感知理論對重磁場物性進行反演。壓縮感知模型如下：

y=Ax

(2)

其中，y是M維的觀測向量；A表示測量矩陣，維度為M×N(M

由于M

(3)

Li等[19]指出零范數的求解是NP難問題，但是當滿足一定條件的時候，L1范數最小化和L0范數最小化是等價的，因而式(3)可以轉化為下式：

(4)

式(3)的求解比較常用的是凸優化算法和貪婪算法。貪婪算法雖然速度快且精度相當，但是貪婪算法大多數需要知道信號的稀疏度。而凸優化算法是將式(3)轉化為凸優化問題，并求解目標函數得到全局最優解。比較典型的L1范數最小化求解方式是基追蹤(Basis Pursuit, BP)算法，BP算法的精度高，本文選用的是基于對數障礙規劃算法的基追蹤反演方法。

考慮到L1范數對稀疏度的刻畫比較容易求解，本文將常規的基追蹤方法對噪聲的壓制的L2范數改為L1范數的約束，那么基于壓縮感知理論的基追蹤反演問題表示為：

(5)

其中，m=(m1,m2,…,mN)表示異常體的磁異常;ε表示噪聲強度。

式(5)所示的目標函數沒有考慮深度加權，容易使反演結果趨于地表。因此，考慮加入深度加權Wm，其形式與文獻[4]中基本一致，只是參數設置不同。同時我們也引入數據加權矩陣Wd，無噪聲情況下，其形式為：

(6)

其中，σ為觀測數據的標準差；I為單位矩陣。

為方便采用上述算法進行求解，令

(7)

則式(5)可變為：

(8)

由拉格朗日乘子法，式(8)可轉化為求解下式優化問題：

(9)

其中,λ為拉格朗日參數。

式(9)的求解可以根據文獻[20]，基于改進的基追蹤反演方法的目標函數，令m=x-y，dnew-Gnewm=z-w，其中xi=(mi)+，yi=(-mi)+，zi=(dnew-Gnewm)+，wi=(-dnew+Gnewm)+；可以推導出下式：

z-w+Gnew(x-y)=dnew

(10)

式(10)可以寫成下式：

(I,-I,Gnew,-Gnew)(z,w,a,b)=dnew

(11)

其中,I為單位矩陣，進一步令(I,-I,Gnew,-Gnew)=D，(z,w,a,b)=c，dnew=b,k=(1,…,1,λ,…,λ)T，則式(9)的求解與下式的優化問題等價：

(12)

式(12)的求解可以采用對數障礙規劃算法。

由于反演的物性參數會出現負值和超出范圍的現象，對求解后的結果進行如下約束：

0≤m≤mmax

(13)

其中,mmax指的是物性參數的下界，一般是已知的。

3 實 驗

本節將在模擬數據和實際數據上驗證該算法的性能，并與傳統的基于L2范數的共軛梯度法[21]進行對比。

3.1 模擬數據實驗

模擬數據實驗選取的是6種二維磁性棱柱體模型分別是直立板狀體模型、平行豎直板狀體模型、傾斜板狀體模型、向斜模型、斷層切割模型和垂向尖滅模型。模擬實驗假設異常體都是被均勻磁化的，磁化強度大小為M=100 A/m，磁化傾角為I=45°，測線磁方位角為A=0°，即有效磁化傾角為45°。觀測數據的點間距為20 m，觀測點個數為51個。地下異常體研究空間被剖分為800個網格單元，即20行40列，網格單元是邊長為25 m的正方形。用共軛梯度算法和基于對數障礙規劃算法進行二維磁化強度反演，選取的上下界范圍在0～100 A/m內，反演結果見圖1～圖6，圖中色標柱為磁化強度M，單位為A/m。

圖6 垂向尖滅模型反演結果Fig.6 Vertical pinch model inversion result image

以圖1進行說明， 圖1(a)和圖1(c)分別表示對數障礙規劃算法和共軛梯度算法觀測數據和預測數據的擬合曲線，可以看出兩個算法擬合都非常好；圖1(b)和圖1(d)分別表示直立板狀體模型對數障礙規劃算法和共軛梯度算法的反演結果，白色邊框表示理論模型的輪廓，從中可以看出共軛梯度法的成像結果大致體現了地下異常體的位置，但是邊界的反演卻很模糊，得到的磁化強度大小為50～70 A/m。而用壓縮感知技術的對數障礙規劃算法得到的磁化強度大小更接近真實值100 A/m。其他模型的反演結果類似，因而對數障礙規劃算法能夠獲得物性的稀疏分布，相比傳統的物性反演方法，稀疏重構的結果分辨率得到提高，可以獲得具有尖銳邊界的反演結果，對目標地質體的定位更加準確。

圖1 直立板狀體模型反演結果Fig.1 Upright plate model inversion result image

圖2 平行豎直板狀體模型反演結果Fig.2 Parallel vertical plate model inversion result image

圖3 傾斜板狀體模型反演結果Fig.3 Inclined plate model inversion result image

圖4 向斜模型反演結果Fig.4 Syncline model inversion result image

圖5 斷層切割模型反演結果Fig.5 Fault cutting model inversion result image

3.2 實際數據實驗

實際數據選擇的是青海省尕林格鐵礦保護區，將本文方法應用于該地區[22]，對鐵礦區兩條沿測線的剖面數據進行反演解釋。圖7顯示的是青海省尕林格礦區磁異常平面等值線圖，該區各磁異常分布明顯，磁異常的幅值最高為1 600 nT。該礦區的主磁鐵礦被大面積的砂礫層圍巖所覆蓋，反演結果圖中Mt為圍巖中出露的磁鐵礦體[23]。該礦區磁鐵礦的平均磁化強度是40 A/m，因此設置磁化強度上下界范圍0～40 A/m。對于該地區的磁力數據反演，前人已經做過很多工作[24,25]。因此，本方法只是用來驗證該算法在實際數據反演中的有效性。212和196兩條平行測線穿過主要的磁異常礦體暴露區域而且兩條平行測線均經過很多鉆孔點，對于后續反演測線結果的驗證具有重要意義。每條平行測線的總長為1 200 m，點間距設置為20 m，即觀測點數據為61個，將地下空間均勻剖分為20行48列，網格單元是邊長為25 m的正方體，反演結果見圖8和9。

圖7 青海省尕林格礦區磁異常平面等值線[22]Fig.7 Plane contour map of magnetic anomaly in Galinge mining area, Qinghai province

以圖8進行說明，圖8(a)表示數據擬合曲線圖，圖8(b)表示212測線反演結果圖。從圖8中可以看出，212測線反演出來的異常體區域的位置大致與鉆孔剖面反映出來的信息相符合，但是在一定程度上未能夠反映礦體的傾斜方向，而且更深層的位置沒有反演出來。圖9中196測線右邊反演出來的異常體區域大體上也符合鉆孔信息，但是也存在212測線深層區域沒有反演出來的問題。另外，196測線左邊的反演結果是未經過鉆孔驗證的，只能推測測線左邊的礦體位于300～400 m左右的位置，且礦體規模小于測線右邊。

圖8 尕林格鐵礦區212測線成像效果示意圖Fig.8 Schematic diagram of imaging effect of 212 survey line in Galinge iron ore area

圖9 尕林格鐵礦區196測線成像效果示意圖Fig.9 Schematic diagram of imaging effect of 196 survey line in Galinge iron ore area

4 結 論

本文突破傳統的反演理論，將基于對數障礙規劃的算法用于位場數據反演。模擬數據和實際數據實驗結果表明，與傳統反演方法相比，該方法不僅可以快速反演位場的稀疏分布，而且可以獲得分辨率高、具有尖銳邊界的反演結果。下一步的工作將進一步從實際問題出發，增加聯合反演的工作，從而更好地解釋反演結果，并改進算法對深層區域加以改善。