基于曲波變換的大地電磁二維稀疏正則化反演

蘇揚， 殷長春*， 劉云鶴， 任秀艷， 張博， 邱長凱， 熊彬

1 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026 2 中國地質(zhì)調(diào)查局發(fā)展研究中心， 北京 100037 3 桂林理工大學(xué)地球科學(xué)學(xué)院， 桂林 541006

0 引言

大地電磁(MT)方法廣泛用于石油、天然氣、地?zé)豳Y源勘探以及地球深部電性結(jié)構(gòu)研究(Sheard et al.，2005；Bedrosian，2007； 趙國澤等，2007).過去幾十年中，大地電磁正反演算法發(fā)展迅速，經(jīng)歷了從20世紀80年代發(fā)展的一維反演(Constable et al.，1987)，到十幾年前的二維反演(de Groot-Hedlin and Constable，2004；Mehanee and Zhdanov，2002)，以及現(xiàn)在主流的三維反演(鄧琰等，2019；阮帥等，2020).然而目前主流的三維大地電磁反演代碼主要是基于L2范數(shù)的光滑約束實現(xiàn)的(Rodi and Mackie，2001；Newman and Alumbaugh，2002；Siripunvaraporn et al.，2005；Egbert et al.，2014)，其優(yōu)勢是可保證反演穩(wěn)定收斂，但其主要缺點是過度光滑會嚴重損失邊界的分辨率.為了克服這一缺點，一系列提高邊界刻畫能力的正則化反演方法被提出，de Groot-Hedlin和Constable(2004)提出了一種基于Occam反演方法的尖銳邊界反演(sharp boundary inversion)，提高了邊界的分辨率，但是反演效果嚴重依賴于初始模型.Mehanee和Zhdanov(2002)提出了二維聚焦反演方法，并將其與非線性共軛梯度(NLCG)和快速松弛反演(RRI)方法進行比較，驗證聚焦反演可以得到更為清晰的邊界.Farquharson(2008)提出了基于L1范數(shù)的最小結(jié)構(gòu)二維反演，該反演不僅考慮了水平和垂直方向的模型粗糙度，還考慮傾斜方向的模型粗糙度，由此可以恢復(fù)異常體的傾斜邊界.張羅磊等(2010)將最小梯度支撐泛函用于二維大地電磁正則化反演，得到了具有尖銳邊界的反演結(jié)果.Jahandari和Farquharson(2017)采用非結(jié)構(gòu)有限元法實現(xiàn)最小結(jié)構(gòu)反演，并將其成功應(yīng)用于三維大地電磁數(shù)據(jù)反演中.雖然上述基于L1范數(shù)的正則化反演在刻畫異常體邊界上具有一定的優(yōu)勢，但是其具體參數(shù)選取以及收斂穩(wěn)定性仍需進一步研究.

近年來，地球物理學(xué)家采用一種新的系數(shù)正則化項替代傳統(tǒng)的模型粗糙度，以獲得更可靠的反演結(jié)果.Abubakar等(2012)使用離散傅里葉變換和離散小波變換壓縮模型空間，并將該方法應(yīng)用于挪威TrollWest Gas Province地區(qū)的可控源電磁(CSEM)數(shù)據(jù)反演.Nittinger和Becken(2016)提出基于復(fù)小波變換的二維大地電磁反演，獲得了與L2范數(shù)相似的反演結(jié)果.蘇揚等(2018)將哈爾(Haar)小波變換應(yīng)用于一維航空電磁數(shù)據(jù)反演，并提出廣義模型約束反演算法，提高對地層界面位置的刻畫.Liu等(2018)引入小波變換建立三維航空電磁稀疏正則化反演方法，并對挪威Bynest地區(qū)的實測數(shù)據(jù)進行測試，結(jié)果表明基于小波變換的稀疏正則化反演可以很好地提高模型分辨率.然而小波變換的固有缺陷是缺少方向性，不能最佳地表示高維數(shù)據(jù).與小波變換相比，Candès等(2006)提出的第2代曲波變換不僅具有多尺度性，還具有多方向特征，可以很好地提取目標體的邊緣特征.曲波變換最初主要應(yīng)用于圖像壓縮和去噪(Starck et al.，2002).Herrmann和Verschuur(2004)首先采用曲波變換進行地震數(shù)據(jù)處理.Shan等(2009)提出了小波和曲波的組合方案，通過求解L1范數(shù)優(yōu)化問題對地震數(shù)據(jù)進行去噪，取得了良好的效果.此外，Herrmann等(2008)還將曲波變換應(yīng)用于一次波和多次波分離，Wang等(2008)將貝葉斯概率分析與曲波變換相結(jié)合，建立一種稀疏正則化反演方法進行一次波和多次波分離，Herrmann和Hennenfent(2008)將曲波變換用于地震數(shù)據(jù)恢復(fù).

截止目前，曲波變換還沒有像小波變換那樣被廣泛應(yīng)用于電磁數(shù)據(jù)反演中.在本文中，我們從傳統(tǒng)的空間域正則化反演出發(fā)，利用曲波變換將空間域中的電阻率模型轉(zhuǎn)換為曲波域中的稀疏系數(shù)來構(gòu)建新的正則化項，構(gòu)建迭代加權(quán)最小二乘(IRLS)算法，并通過共軛梯度(CG)方法求解高斯-牛頓(GN)方程.首先介紹大地電磁的反演策略，然后分析曲波變換的多尺度和多方向特征并介紹實施步驟，進而通過設(shè)計理論合成算例，分析不同尺度系數(shù)和不同正則化項對反演結(jié)果的影響，驗證基于曲波變換的正則化反演比傳統(tǒng)L1和L2范數(shù)反演方法的優(yōu)越性.最后，通過對南非SAMTEX項目中的二維大地電磁實測數(shù)據(jù)反演進一步驗證本文反演算法的有效性.

1 方法原理

1.1 反演策略

我們先給出空間域中基于Lp范數(shù)正則化反演的目標函數(shù)：

Φ=φd(p)+λφm(q),

(1)

其中，φd(p)是數(shù)據(jù)擬合項，φm(q)是模型約束項，λ是正則化因子.p和q的具體形式為

p=Wd(dobs-dprd),

(2)

q=Wm(m-mref),

(3)

式中，Wd是一個對角矩陣，其元素為觀測電磁響應(yīng)中噪聲的倒數(shù)，dobs為數(shù)據(jù)向量，dprd為反演模型m經(jīng)正演計算得到的響應(yīng)數(shù)據(jù)，Wm是平滑算子矩陣，mref是包含地質(zhì)信息的參考模型.公式(2)中dprd可以表示為

dprd=F(m),

(4)

其中，F(xiàn)是正演算子.

公式(1)中φd和φm的一般形式為

(5)

其中，x表示公式(2)和(3)中的p或者q，i表示向量x中元素的索引.ρ(xi)可有多種選擇，具體選擇在下文中給出.我們對數(shù)據(jù)擬合項采用L2范數(shù)約束，而對正則化項采用L1范數(shù)約束，給出迭代過程中公式(2)和(3)的形式

p=Wd(dobs-dn-1-Jδm),

(6)

q=Wm(mn-1+δm-mref),

(7)

其中，dn-1是上次迭代的模型mn-1經(jīng)正演計算得到的響應(yīng)數(shù)據(jù)，δm=mn-mn-1，J是空間域模型參數(shù)的靈敏度矩陣，可以通過伴隨正演技術(shù)計算得到(Liu and Yin，2016).將公式(5)對模型更新量δm求導(dǎo)可得(Farquharson，2008)

(8)

(8)式可表示為

(9)

式中，?φ/?δm=(?φ/?δm1,…,?φ/?δmN)T，Bij=?xi/?δmj，u=(ρ′(x1),…,ρ′(xN))T.引入廣義約束矩陣R，將公式(9)重寫為

(10)

其中，R為一個對角矩陣，R=diag{ρ′(x1)/x1,…,ρ′(xN)/xN}.對于數(shù)據(jù)擬合項，B=-WdJ，對于模型約束項，B=Wm.

對于公式(5)中的ρ(xi)，我們使用Ekblom(1987)提出Lp范數(shù)的形式：

ρ(x)=(x2+ε2)p/2，

(11)

式中，ε是一個很小的正實數(shù)，用于確保導(dǎo)數(shù)在x=0處存在.由此，對角矩陣R的非零元素表示為

(12)

由于我們對數(shù)據(jù)擬合項使用L2范數(shù)約束，因此根據(jù)公式(12)，矩陣Rd的對角元素為2；對模型約束項使用L1范數(shù)約束，廣義約束矩陣用Rm表示.

對目標函數(shù)求導(dǎo)并令其為零，利用Gauss-Newton法求解最優(yōu)化問題，可以得到反演迭代中求解的線性方程組為

(13)

(14)

(15)

(16)

同時，根據(jù)曲波基函數(shù)的正交性，將(16)式改寫為

(17)

這里不再需要給出正則化項的平滑算子矩陣Wm，因為我們將模型轉(zhuǎn)換為曲波系數(shù)本身就是對模型的一種整體約束.同樣的我們對系數(shù)正則化項采用L1范數(shù)約束，Donoho(2006)，Loris等(2007)和Simons等(2011)證明了反演中對稀疏域系數(shù)采用L1范數(shù)約束更傾向于得到稀疏解.(17)式中對角矩陣R以及靈敏度矩陣J都取決于模型，它們在每次迭代時更新，這種方法被稱為迭代加權(quán)最小二乘(IRLS)算法.對于每次迭代，我們用共軛梯度方法求解上述方程.在利用方程(17)求解出曲波系數(shù)改正量后，空間域中的模型更新可表示為

(18)

其中s是通過線性搜索確定的比例因子(Haber，2014).

1.2 曲波變換

作為多尺度分析方法，曲波變換目前已廣泛應(yīng)用于圖像處理、醫(yī)學(xué)成像和地震勘探等諸多領(lǐng)域.曲波變換克服了傳統(tǒng)的多尺度分析方法(比如小波變換)缺乏方向性的弱點.為了更好地理解基于曲波變換的稀疏正則化反演的優(yōu)越性，我們在對頻率域和空間域中曲波進行數(shù)學(xué)描述的基礎(chǔ)上，重點介紹如何實現(xiàn)曲波變換并分析曲波變換的稀疏性特征.

1.2.1 頻率域中的曲波

根據(jù)Herrmann和Hennenfent(2008)，二維曲波窗函數(shù)U是由徑向窗函數(shù)W和角度窗函數(shù)V共同定義的多尺度和多方向的楔形基(見圖1a).這里我們僅給出笛卡爾坐標系中二維離散曲波的定義形式.為此，首先引入徑向窗函數(shù)(Candès et al.，2006)

(19)

其中，ω= (ω1,ω2)是頻率變量，Φ是兩個一維低通窗函數(shù)φ的乘積，即

Φj(ω1,ω2)=φ(2-jω1)φ(2-jω2).

(20)

同樣，我們引入角度窗函數(shù)

Vj(ω)=V(2?j/2_|ω2/ω1)，

(21)

其中，?j/2_|是j/2的整數(shù)部分.由此，笛卡爾坐標系下的曲波窗函數(shù)可以定義為

Uj(ω)∶=Wj(ω)Vj(ω).

(22)

進而，我們引入一組等間距的斜率tanθl∶=l·2-?j/2_|,l=-2-?j/2_|,…,2?j/2_|-1，則由(22)式，可以得到

Uj,l(ω)=Wj(ω)Vj(Sθlω),

(23)

其中，Sθ是剪切矩陣，具體形式為

(24)

圖1a給出頻率域(又稱曲波域)中的楔形曲波基與空間域中曲波之間的對應(yīng)關(guān)系.其中，實線表示徑向窗函數(shù)W劃分的同心正方格，而虛線表示角度窗函數(shù)V又將同心正方格劃分為多個曲波基.k1和k2為頻率域參數(shù)，曲波基是長度為2j，寬度為2j/2的各向異性楔形基.其中，最內(nèi)側(cè)的第1尺度曲波基是無方向的，第1尺度的系數(shù)是粗尺度系數(shù)，而其他尺度的系數(shù)為細節(jié)系數(shù).隨著尺度增加，楔形基的方向數(shù)也隨之增加(Herrmann and Hennenfent, 2008).

圖1 (a) 頻率域中的楔形曲波基對應(yīng)于空間域中的曲波； (b) 通過曲波逆變換可從圖(a)中的三個楔形基獲得空間域 中的三個不同曲波.兩個域中曲波之間的方向正交性是由傅里葉變換旋轉(zhuǎn)90°引起的(參考Herrmann et al.，2007)Fig.1 (a) Curvelets in frequency domain with each wedge corresponding to the frequency support of a curvelet in the space domain; (b) Three different curvelets in the space domain are obtained from 3 polar wedges in (a) by inverse curvelet transform. The directional orthogonality between the curvelets in two domains is resulted from 90° rotation of Fourier transform (after Herrmann et al., 2007)

1.2.2 空間域中的曲波

(25)

其中，x=(x1,x2)是空間域參數(shù)，0≤θl<2π,k=(k1,k2)∈Z2是旋轉(zhuǎn)參數(shù), 而Rθl是弧度θl的旋轉(zhuǎn)矩陣，即

(26)

曲波系數(shù)可以通過信號f∈L2(R2)與曲波φj,l,k(x)進行內(nèi)積得到，即

(27)

其中，上劃線表示復(fù)共軛.空間域中的曲波如圖1b所示.藍色和紅色箭頭指示第3尺度上不同方向和位置的曲波，黃色箭頭指示第4尺度的曲波.空間域中曲波的形狀是狹長的，其包絡(luò)的有效長度為2-j/2，有效寬度為2-j，滿足各向異性尺度關(guān)系.此外，曲波沿脊方向平滑衰減，而在垂直脊方向上振蕩衰減(Ying et al.，2005).隨著尺度的增加，曲波變得越來越細小，對圖像中彎曲邊緣的敏感性更高.從公式(27)可以看出曲波系數(shù)是圖像和曲波在空間域中的內(nèi)積，如果曲波的方向與目標體的邊界方向一致，則曲波系數(shù)很大；否則，曲波系數(shù)接近于0.由此，目標體邊緣在每個精細尺度下都對應(yīng)于大的系數(shù).換句話說，通過對各尺度下較大的系數(shù)進行反演可以實現(xiàn)對目標體邊界的精細刻畫.

1.2.3 曲波變換的實現(xiàn)

Candès等(2006)提出了兩種用于實現(xiàn)快速離散曲波變換(FDCT)的算法，分別是Unequispaced FFTs (USFFT)和Wrapping 算法.本文中我們采用Curvelab工具箱(www.curvelet.org)提供的Wrapping算法開發(fā)我們的反演方法.該方法比USFFT算法更易于實現(xiàn)且運行速度快.利用Wrapping算法的快速離散曲波變換包括以下步驟 (Ying et al.，2005；Candès et al.，2006)：

(1)將二維快速傅里葉變換(FFT)用于模型參數(shù)，獲得如下傅里葉樣本集：

其中，[n1,n2]是二維離散函數(shù)的數(shù)組索引，n是數(shù)組的維數(shù).

(2)將傅里葉樣本與每個尺度j和角度l的Uj,l相乘，形成如下乘積：

(3)對乘積進行周期性重復(fù)，并以直角坐標系原點為中心進行數(shù)據(jù)提取，得到

(28)

(4)對每個fj,l使用二維逆FFT，獲得離散曲波系數(shù).

該算法的關(guān)鍵是將楔形基環(huán)繞在直角坐標系原點周圍，以形成可用于二維逆FFT的標準矩形.楔形基數(shù)據(jù)包含在大小為2j×2j/2的平行四邊形中(圖2中的黑色框).我們對平行四邊形內(nèi)的數(shù)據(jù)進行分塊，然后將其Wrapping到圍繞坐標系原點的標準矩形區(qū)域內(nèi)，這樣就可將原始平行四邊形包含的數(shù)據(jù)集中到大小為2j×2j/2的標準矩形中并完成數(shù)據(jù)提取.

圖2 用于楔形基數(shù)據(jù)提取的Wrapping技術(shù) 黑色平行四邊形包含楔形基數(shù)據(jù)，而灰色平行四邊形是周期性平滑移動產(chǎn)生的復(fù)制結(jié)果.楔形基數(shù)據(jù)通過平滑移動復(fù)制到環(huán)繞原點的紅色平行四邊形中(參照Candès et al.，2006).Fig.2 Wrapping wedge data The black parallelogram contains the polar wedge data, whereas the gray parallelograms are the replicas resulting from periodization. The parallelogram data are wrapped around the origin and contained in the red rectangle (after Candès et al., 2006).

在本文中，為方便讀者理解，我們定義曲波正變換算子為Wc(張華等，2017)，

Wc=TF,

(29)

其中，F(xiàn)表示二維快速傅里葉變換，T表示曲波拼接算子，即將頻率域轉(zhuǎn)換到曲波系數(shù)的過程.與小波變換類似，曲波基函數(shù)也滿足正交性，因此可以定義曲波逆變換算子為

(30)

其中，F(xiàn)-1表示二維快速傅里葉逆變換，T-1表示曲波逆拼接算子，即將曲波系數(shù)轉(zhuǎn)換到頻率域的過程.

1.2.4 目標體邊緣的最佳稀疏表示

參見圖3，由于曲波的多方向性和遵循拋物線比例關(guān)系(各向異性尺度關(guān)系)，因此與小波相比，我們利用數(shù)量較少的曲波就可以擬合目標體邊緣.這意味著曲波變換的一個突出優(yōu)點是可以對目標體的邊緣進行最佳稀疏表示.

圖3 (a) 小波擬合彎曲邊緣； (b) 曲波擬合彎曲邊緣 (參看Alzubi et al.，2011)Fig.3 Comparison of wavelets (a) and curvelets (b) for curved edges (after Alzubi et al.，2011)

2 反演算例

本文中我們提出的新方法是使用曲波變換將電阻率模型轉(zhuǎn)換為曲波系數(shù)進行求解的，這與基于相鄰單元之間電阻率差異的模型約束反演不同，基于曲波變換的反演是對頻域中曲波系數(shù)的整體約束.第1尺度的系數(shù)恢復(fù)了電阻率模型的整體結(jié)構(gòu)，而其余尺度中的較大系數(shù)恢復(fù)了目標體的邊緣細節(jié).因此，以曲波系數(shù)為正則化項的反演算法是一種多尺度、多分辨率反演方法.

2.1 二維大地電磁理論算例

在下面的反演算例中，我們均采用如下規(guī)則的矩形網(wǎng)格對模型進行剖分：垂直方向上共有64層.第1層的厚度為500 m，隨后各層的厚度為前一層的1.07倍，模型總體深度范圍約240 km.水平方向上共有128個網(wǎng)格單元，計算區(qū)域的寬度為5000 m，擴邊區(qū)域的網(wǎng)格寬度為前一個網(wǎng)格寬度的1.2倍，總體延伸約400 km.我們在計算區(qū)域內(nèi)設(shè)計31個均勻分布的測點，間距為15316 m.進而，利用ModEM軟件(Kelbert et al., 2014)計算了0.001～110 Hz范圍內(nèi)10個頻率的TE和TM阻抗，并對所有數(shù)據(jù)添加阻抗模長的3%高斯噪聲后進行反演.對于基于曲波變換的稀疏正則化反演，我們使用具有4個尺度的曲波變換，其中第2尺度有16個方向，第3和第4尺度有32個方向，這樣選擇可以在有限的系數(shù)下就很好地表達模型特征.反演中首先給定初始的正則化因子，在每步迭代中正則化因子以一個固定的衰減率下降，通過對比不同初始正則化因子的反演結(jié)果，發(fā)現(xiàn)采用這種方法反演結(jié)果相差不大，我們選擇迭代次數(shù)最少的初始正則化因子.對于實測數(shù)據(jù)，為提高曲波系數(shù)稀疏正則化項的權(quán)重，反演初始階段我們選擇一個較大的正則化因子，隨著反演進程逐漸趨近于收斂，正則化因子逐漸減小.

2.1.1 算例1：不同尺度系數(shù)反演結(jié)果對比

本文的第1個算例是比較不同尺度系數(shù)的反演結(jié)果.參見圖4a，背景電阻率為100m.中間紅色區(qū)域的電阻率為10m，兩側(cè)藍色區(qū)域的電阻率為1000m.設(shè)定初始模型為100m的均勻半空間，設(shè)定正則化因子為10，并且每次迭代減小十分之九.

圖4b—e給出不同尺度曲波系數(shù)的反演結(jié)果.從圖4b可以看出，第1尺度系數(shù)的反演結(jié)果很粗糙，存在較多的假異常.隨著我們在反演中逐漸增加更多尺度的系數(shù)，反演結(jié)果越來越接近真實模型(參見圖4c—e).從圖4e中可以看出，當所有尺度的系數(shù)參與反演時，反演結(jié)果越可靠.由于本文算法中曲波變換具有多尺度特征，不同尺度的系數(shù)包含了模型的不同信息，粗尺度系數(shù)表征模型輪廓的低頻信息；而細尺度系數(shù)則表征模型精細特征(比如模型邊界)的高頻信息.因此，如果僅將粗尺度系數(shù)用于反演，結(jié)果將是粗糙的并且存在假異常，只有將所有的曲波系數(shù)都用于反演，才能獲得最可靠的反演結(jié)果.

圖5給出數(shù)據(jù)擬合差、目標函數(shù)和正則化項數(shù)值.從圖5a中可以看出，所有的反演都具有相似的收斂速度，但僅使用第1尺度系數(shù)的反演沒有很好地收斂.相比之下，其他3套組合系數(shù)反演的數(shù)據(jù)擬合差經(jīng)過20次迭代后均降至1，說明在恢復(fù)了模型的主要特征后，對模型的細節(jié)特征進行刻畫不會改變反演的收斂性.

2.1.2 算例2：正則化項對反演結(jié)果的影響

本文的第2個算例是比較L1和L2范數(shù)，以及基于小波變換和曲波變換的稀疏正則化反演方法對反演結(jié)果的影響.參見圖6a，背景半空間的電阻率為100m，藍色起伏層的電阻率為1000m，而紅色侵入層的電阻率為5m.所有4種反演算法的正則化因子均設(shè)定為10，并且每次迭代將其降低十分之九.初始模型設(shè)定為100m的均勻半空間.

圖4 不同尺度系數(shù)的反演結(jié)果 (a) 理論模型； (b) 第1尺度系數(shù)反演結(jié)果； (c) 第1和第2尺度系數(shù)反演結(jié)果； (d) 第1、第2和第3尺度系數(shù)反演結(jié)果； (e) 所有4個尺度系數(shù)的反演結(jié)果.Fig.4 Inversion results of different scale coefficients (a) Original theoretical model; (b) Inversion result of 1st scale coefficients; (c) Inversion result of 1st and 2nd scale coefficients; (d) Inversion result of 1st, 2nd and 3rd scale coefficients; (e) Inversion result of all 4 scales coefficients.

圖5 (a)數(shù)據(jù)擬合差； (b) 目標函數(shù)Φ； (c) 正則化項Fig.5 (a) Data misfit; (b) Objective functional Φ; (c) Regularization term

圖6 (a) 真實模型； (b) L1范數(shù)反演模型； (c) L2范數(shù)反演模型； (d) 基于小波變換的稀疏正則化反演模型； (e) 基于曲波變換的稀疏正則化反演模型Fig.6 (a) True model; (b) Model for L1-norm inversion; (c) Model for L2-norm inversion; (d) Model for sparse inversion based on wavelet transform; (e) Model for sparse inversion based on curvelet transform

圖6b—e分別給出L1范數(shù)，L2范數(shù)，基于小波變換和基于曲波變換的稀疏正則化反演結(jié)果.從圖中可以看出，L1范數(shù)反演的結(jié)果是塊狀的，無法恢復(fù)異常體的真實結(jié)構(gòu)；L2范數(shù)反演結(jié)果過于光滑，目標體的邊界沒有被清楚地恢復(fù).對于基于小波變換的稀疏正則化反演，我們選擇db20小波，這是因為通過使用較大消失矩的小波可以更好地表征模型(Liu et al.，2018).由圖3可以看出小波基函數(shù)是塊狀，而曲波基函數(shù)是楔形，因此與小波變換相比，用較少的曲波基函數(shù)就可以擬合曲線邊界，所以曲波變換的系數(shù)比小波變換的系數(shù)更加稀疏.對于稀疏正則化反演，正則化項采用L1范數(shù)約束，系數(shù)越稀疏，越容易在反演過程中恢復(fù)出關(guān)鍵系數(shù).同時，由于小波變換的振蕩特性，反演結(jié)果會出現(xiàn)一些冗余構(gòu)造.從圖6d給出的結(jié)果可以看出，基于小波變換的稀疏正則化反演結(jié)果存在許多假異常，目標體的邊界也沒有準確地恢復(fù).然而，從圖6e給出的曲波系數(shù)正則化反演可以看出，反演結(jié)果較好地恢復(fù)了異常體的邊界.實際上，當我們將模型參數(shù)轉(zhuǎn)換為曲波系數(shù)時，精細尺度中較大系數(shù)對應(yīng)于異常體的邊界，對這些大系數(shù)進行反演可恢復(fù)模型的邊界特征.因此，較之于其他算法，基于曲波變換的稀疏正則化反演可以更好地刻畫地下結(jié)構(gòu)的邊界特征.

圖7展示了所有4種算法的數(shù)據(jù)擬合差，目標函數(shù)和正則化項數(shù)值.由圖可以看出，所有算法都具有良好的收斂性(圖7a)，反演迭代20次后數(shù)據(jù)擬合差都降至1.在后期迭代中，正則化項接近一個恒定值(圖7c).另外，所有反演中基于L2范數(shù)的反演收斂最快，基于曲波變換的稀疏正則化反演收斂最慢，而L1范數(shù)和基于小波變換的稀疏正則化反演居于其中.這種現(xiàn)象很容易進行解釋.事實上，L2范數(shù)反演的模型粗糙度項是凸函數(shù)，易于收斂；相比之下，基于曲波變換的稀疏正則化反演方法在反演出模型的總體特征(輪廓)后，還需要搜索精細的邊界信息，導(dǎo)致反演速度較慢.

從計算成本的角度來看，我們得出相同的結(jié)論.對于圖6中的模型，L1范數(shù)反演需要70 min，L2范數(shù)反演需要60 min，基于小波變換的稀疏正則化反演需要90 min，而基于曲波變換的稀疏正則化反演需要76 min.這說明基于曲波的反演雖然對地下電阻率結(jié)構(gòu)具有更好的分辨率，但要獲得模型的細節(jié)特征，需花費的時間要比L1范數(shù)和L2范數(shù)反演長.需要說明的是，由于小波變換不具有方向性，為刻畫目標體的邊界等細節(jié)特征，基于小波變換的反演需要花費更長的時間用于搜索模型空間.

2.1.3 算例3：深部目標體的分辨能力

為了研究本文曲波變換稀疏正則化反演對深層目標體的探測能力，我們設(shè)計了一個如圖8a所示的模型：背景電阻率為100m，淺層異常體的電阻率為10m，深藍色異常體的電阻率為1000m，紅色異常體的電阻率為2m.我們分別采用L2范數(shù)，L1范數(shù)和基于曲波變換的稀疏正則化反演對理論數(shù)據(jù)進行反演.反演中初始正則化因子均設(shè)為100，每次迭代將其降低十分之九，初始模型均設(shè)定為100m的均勻半空間.

圖7 (a) 數(shù)據(jù)擬合差； (b) 目標函數(shù)Φ； (c) 正則化項Fig.7 (a) Data misfit; (b) Objective functional Φ; (c) Regularization term

圖8 (a) 真實模型； (b) L1范數(shù)反演模型； (c) L2范數(shù)反演模型； (d) 基于曲波變換的稀疏正則化反演模型Fig.8 (a) True model; (b) Model for L1-norm inversion; (c) Model for L2-norm inversion; (d) Model for sparse inversion based on curvelet transform

圖8給出L1范數(shù)，L2范數(shù)和基于曲波變換的稀疏正則化反演的對比結(jié)果.從圖8b和8c可以看出，L1和L2范數(shù)反演都恢復(fù)了淺層地下模型，異常體的淺層邊界較為清晰.然而，對于50 km以下的三個異常體，L1范數(shù)和L2范數(shù)反演都無法準確地恢復(fù)其電阻率和邊界位置.從圖8d中可以看出，除了對淺部異常體具有較高的分辨率外，基于曲波變換的稀疏正則化反演對深部異常體的分辨也較L1和L2范數(shù)反演效果好.

圖9展示了所有三種算法的數(shù)據(jù)擬合差，目標函數(shù)和正則化項數(shù)值.由圖可以看出，所有算法都具有良好的收斂性(圖9a)，反演迭代50次后數(shù)據(jù)擬合差都降至1.

2.2 實測數(shù)據(jù)反演

為了進一步測試基于曲波變換的稀疏正則化反演的有效性，我們對來自非洲南部大地電磁實驗項目SAMTEX的ZIM數(shù)據(jù)進行反演(Miensopust et al.，2011).參見圖10，二維剖面跨越津巴布韋(ZIM)克拉通，Magondi帶和Ghanzi-Chobe帶，全長約500 km，共計31個測點(圖中紅點標識).大地電磁觀測周期為0.008986～856.091 s.本文著重于分析不同反演方法的差異以及噪聲估計對反演結(jié)果的影響，相關(guān)的詳細地質(zhì)解釋參見Miensopust等(2011).

我們建立了一個128×64的電阻率網(wǎng)格模型.計算區(qū)域中沿水平方向的網(wǎng)格單元大小為5000 m，而垂直方向上頂部網(wǎng)格單元厚度為500 m，下部網(wǎng)格單元厚度以1.07的比率增加.背景半空間初始電阻率設(shè)置為1000m.對于L2范數(shù)、L1范數(shù)及曲波變換稀疏正則化三種反演算法，初始正則化因子均設(shè)定為1000，每次迭代將其降低至先前值的十分之九，直到降為1時保持不變.我們首先假設(shè)所有測點的噪聲水平為5%.圖11b—d展示了L1范數(shù)，L2范數(shù)和基于曲波變換的稀疏正則化的反演結(jié)果對比.由圖可以看出，三種方法給出了相似的反演結(jié)果，地下主要電性結(jié)構(gòu)特征得到刻畫.在斷面南部107～115測點下有一個近地表的高電阻率層，即Okavango Dike Swarm(ODS)，同時在橫向范圍內(nèi)有若干個中低阻異常體(de Beer et al.，1975)，在剖面北部位于125～131測點處的高電阻率結(jié)構(gòu)對應(yīng)于Ghanzi-Chobe帶(GCB).然而，從圖11a可以看到，104～108測點和119～124測點的數(shù)據(jù)擬合差明顯大于其他測點.為此，我們將這些數(shù)據(jù)擬合差較高測點的噪聲水平從5%更改為20%，并重新進行反演.圖11e給出基于新噪聲估計的反演結(jié)果.由圖可以看出，這些測點的數(shù)據(jù)擬合差顯著降低.對比圖11b—d和10f—h中的結(jié)果，我們發(fā)現(xiàn)當這些測點噪聲水平被改變后，L2范數(shù)和L1范數(shù)反演結(jié)果發(fā)生很大變化.相比之下，基于曲波變換的稀疏正則化反演結(jié)果幾乎沒有發(fā)生變化.這意味著與L2范數(shù)和L1范數(shù)反演相比，基于曲波變換的稀疏正則化反演對噪聲估計不太敏感.

圖9 (a) 數(shù)據(jù)擬合差； (b) 目標函數(shù)Φ； (c) 正則化項Fig.9 (a) Data misfit; (b) Objective functional Φ; (c) Regularization term

圖10 非洲南部SAMTEX項目的ZIM數(shù)據(jù)集對應(yīng)的 測點位置(根據(jù)Miensopust et al.，2011修改)Fig.10 Locations of ZIM subset data from SAMTEX survey in Southern Africa (modified from Miensopust et al., 2011)

圖12展示了所有三種算法在兩種噪聲水平類型下的數(shù)據(jù)擬合差，目標函數(shù)和正則化項，可以看出所有算法都正常收斂.

3 討論

L1和L2范數(shù)反演的模型約束是空間域中相鄰單元之間的電阻率差異.對于L1范數(shù)反演，允許相鄰單元之間存在較大差異，并且模型差分沿水平和垂直方向.從圖6b和圖8b可以看出，異常體的反演結(jié)果在水平和垂直方向均是塊狀的；對于L2范數(shù)反演，模型要求相鄰單元之間的差異不大.因此，從圖6c和圖8c給出的反演結(jié)果中未觀察到異常體之間的明顯邊界，電阻率在水平和垂直方向上均呈現(xiàn)光滑過度的特征.從圖6e和圖8d可以看出，基于曲波變換的稀疏正則化反演的結(jié)果更接近真實模型(既沒有出現(xiàn)塊狀又沒有出現(xiàn)過度光滑的反演結(jié)果).基于曲波變換的稀疏正則化反演的這種優(yōu)點可以得到很好的解釋.事實上，與傳統(tǒng)的空間域反演不同，我們將空間域模型轉(zhuǎn)換為曲波域中的稀疏系數(shù)，并對曲波系數(shù)實施L1范數(shù)正則化來保證系數(shù)的稀疏性.曲波系數(shù)包含大量的模型信息.根據(jù)曲波變換的多尺度和多方向性，第1尺度的系數(shù)相對較大，包含模型的整體概貌，細尺度中較大系數(shù)對應(yīng)于模型的細節(jié)信息(比如異常體邊界).通過對稀疏系數(shù)應(yīng)用L1范數(shù)正則化，我們可以求解出第1尺度的系數(shù)和細尺度的關(guān)鍵系數(shù).這些反演獲得的稀疏系數(shù)既包含了模型的整體概貌，又包含與邊界相關(guān)的細節(jié)信息.因此，基于曲波變換的稀疏正則化反演更有利于恢復(fù)地下復(fù)雜的電性結(jié)構(gòu).

圖11 非洲南部SAMTEX項目的ZIM數(shù)據(jù)集反演結(jié)果 (a) L1 / L2范數(shù)反演和基于曲波變換的稀疏正則化反演的數(shù)據(jù)擬合差； (b) L1范數(shù)反演結(jié)果； (c) L2范數(shù)反演結(jié)果； (d) 基于曲波變換的稀疏正則化反演結(jié)果； (e) L1 / L2范數(shù)反演和基于曲波變換的稀疏正則化反演的數(shù)據(jù)擬合差； (f) L1范數(shù)反演結(jié)果； (g) L2范數(shù)反演結(jié)果； (h) 基于曲波變換的稀疏正則化反演結(jié)果.對于(a)—(d)，所有測點的噪聲設(shè)定為5%，而對于(e)— (h)，在測點104～108和119～124的噪聲設(shè)定為20%，其他測點的噪聲仍為5%.圖中僅展示奇數(shù)測點.Fig.11 Inversion results of ZIM subset data from SAMTEX survey in Southern Africa (a) RMS for L1-/L2-norm and the curvelet-based inversions; (b) Model for L1-norm inversion; (c) Model for L2-norm inversion; (d) Model for curvelet-based inversion; (e) RMS for L1-/L2-norm and the curvelet-based inversions; (f) Model for L1-norm inversion; (g) Model for L2-norm inversion; (h) Model for curvelet-based inversion. For (a)—(d), the noise is assumed to be 5% at all survey sites, while for (e)—(h), the noise is assumed to be 20% at survey sites 104～108 and 119～124 and 5% at other sites. For simplicity, we only show the odd stations.

圖12 (a，d)數(shù)據(jù)擬合差；(b，e)目標函數(shù)Φ；(c，f)正則化項 對于(a)—(c)，所有測點的噪聲水平設(shè)定為5%，而對于(d)—(f)，在測點104～108和119～124的噪聲水平設(shè)定為20%，其他測點的噪聲水平仍為5%.Fig.12 (a,d) Data misfit; (b,e) Objective functional Φ; (c,f) Regularization term For (a)—(c), the noise is assumed to be 5% at all survey sites, while for (d)—(f), the noise is assumed to be 20% at survey sites 104～108 and 119～124 and 5% at other sites.

4 結(jié)論

本文將曲波變換引入到電磁數(shù)據(jù)反演中，成功實現(xiàn)了二維大地電磁稀疏正則化反演.與傳統(tǒng)的將模型粗糙度作為正則化約束項不同，我們的反演結(jié)合了曲波系數(shù)的L1范數(shù)約束和數(shù)據(jù)擬合項的L2范數(shù)約束以獲得最佳解，保證反演結(jié)果具有良好的分辨率.理論模型反演結(jié)果表明，當所有尺度的曲波系數(shù)參與反演時，基于曲波變換的稀疏正則化反演能獲得最佳分辨率.另外，對比L1范數(shù)和L2范數(shù)反演及基于小波變換的稀疏正則化反演結(jié)果，我們進一步證實基于曲波變換的稀疏正則化反演可以更好地恢復(fù)地下異常電阻率邊界特征的結(jié)論.對非洲南部SAMTEX項目實測大地電磁數(shù)據(jù)的反演表明，基于曲波變換的稀疏正則化反演不僅可以獲得較理想的反演結(jié)果，反演結(jié)果受數(shù)據(jù)噪聲估計的影響也較小.

致謝我們衷心感謝Curvelab的作者提供曲波變換開源代碼，感謝Mod2DMT代碼的開發(fā)人員Gary Egbert教授和Anna Kelbert教授，同時也感謝SAMTEX聯(lián)盟提供的大地電磁實測數(shù)據(jù).