葉麗霞，張玉潔

(1.中國地質(zhì)大學(xué) 數(shù)學(xué)與物理學(xué)院，湖北 武漢 430074；2.中國地質(zhì)大學(xué) 地質(zhì)探測(cè)與評(píng)估教育部重點(diǎn)實(shí)驗(yàn)室, 湖北 武漢 430074)

1 引 言

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

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

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

本文提出基于基追蹤的重磁場(chǎng)物性稀疏反演方法，用壓縮感知技術(shù)代替?zhèn)鹘y(tǒng)的反演方法，直接從反演模型出發(fā)，并且引入物性上下界約束，構(gòu)建目標(biāo)函數(shù)，將其轉(zhuǎn)化為有約束的優(yōu)化問題，然后利用基追蹤對(duì)數(shù)障礙規(guī)劃算法進(jìn)行求解。

2 基本原理

2.1 正演

重磁異常的正演計(jì)算表達(dá)式可表示為：

Gm+e=d

(1)

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

2.2 反演

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

y=Ax

(2)

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

由于M

(3)

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

(4)

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

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

(5)

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

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

(6)

其中，σ為觀測(cè)數(shù)據(jù)的標(biāo)準(zhǔn)差；I為單位矩陣。

為方便采用上述算法進(jìn)行求解，令

(7)

則式(5)可變?yōu)椋?/p>

(8)

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

(9)

其中,λ為拉格朗日參數(shù)。

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

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

(10)

式(10)可以寫成下式：

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

(11)

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

(12)

式(12)的求解可以采用對(duì)數(shù)障礙規(guī)劃算法。

由于反演的物性參數(shù)會(huì)出現(xiàn)負(fù)值和超出范圍的現(xiàn)象，對(duì)求解后的結(jié)果進(jìn)行如下約束：

0≤m≤mmax

(13)

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

3 實(shí) 驗(yàn)

本節(jié)將在模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)上驗(yàn)證該算法的性能，并與傳統(tǒng)的基于L2范數(shù)的共軛梯度法[21]進(jìn)行對(duì)比。

3.1 模擬數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

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

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

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

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

3.2 實(shí)際數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

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

4 結(jié) 論

本文突破傳統(tǒng)的反演理論，將基于對(duì)數(shù)障礙規(guī)劃的算法用于位場(chǎng)數(shù)據(jù)反演。模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)實(shí)驗(yàn)結(jié)果表明，與傳統(tǒng)反演方法相比，該方法不僅可以快速反演位場(chǎng)的稀疏分布，而且可以獲得分辨率高、具有尖銳邊界的反演結(jié)果。下一步的工作將進(jìn)一步從實(shí)際問題出發(fā)，增加聯(lián)合反演的工作，從而更好地解釋反演結(jié)果，并改進(jìn)算法對(duì)深層區(qū)域加以改善。