1-Norm/2-Norm約束及小波多尺度2D體波走時成像

高 級，方洪健

(中國科學技術(shù)大學 地球和空間科學學院地震與地球內(nèi)部物理實驗室，合肥 230026)

1-Norm/2-Norm約束及小波多尺度2D體波走時成像

高 級，方洪健

(中國科學技術(shù)大學 地球和空間科學學院地震與地球內(nèi)部物理實驗室，合肥 230026)

在對淺地表復雜地下介質(zhì)層析成像時，由于觀測系統(tǒng)局限性及方法本身對特定巖石地球物理屬性的不敏感性等因素，造成地球物理反演存在多解性。在地震走時成像中，地震射線數(shù)量在低速區(qū)域分布較少形成陰影區(qū)，造成其對低速區(qū)域成像分辨率較低。同時因為觀測數(shù)據(jù)包含干擾信息以及不同反演方法的局限性等因素造成反演結(jié)果包含次生假異常。這里主要研究了空間域的L1-norm、L2-norm約束反演及小波域小波系數(shù)稀疏約束反演，對比不同反演方式及約束條件對地質(zhì)模型的分辨能力。通過測試孤立異常模型、層狀地層模型、傾斜地層模型，得出三種不同反演方式分別具有對不同特定模型的分辨特性。在實際資料成像中，可以通過使用多種反演方式，對比反演結(jié)果合理性以期達到更準確的地質(zhì)構(gòu)造解釋。

體波走時； 層析成像； 1-Norm； 2-Norm； 小波多尺度

0 引言

地震初至走時層析成像就是用地震數(shù)據(jù)來反演地下結(jié)構(gòu)的物質(zhì)屬性，并逐層剖析繪制其圖像的技術(shù)，其主要目的是確定地球內(nèi)部的精細結(jié)構(gòu)和局部不均勻性[1]。

地震層析成像主要包括：①模型的參數(shù)化；②正演計算地下介質(zhì)屬性的理論值；③反演及圖像重建；④反演結(jié)果的評價[2]。

模型參數(shù)化可分為兩類：①不分塊參數(shù)化，反演在泛函空間中進行，理論上可以計算任意位置的速度，結(jié)果不受離散化的影響，但其計算量較大，不利于實際利用；②離散化的模型參數(shù)化，其優(yōu)點是數(shù)學上容易處理，運算相對簡單，缺點是在一般方法中出現(xiàn)的某些簡化，在用離散時可能被掩蓋掉?，F(xiàn)在通用的大都是離散化的模型參數(shù)化，通常采用分塊網(wǎng)格或節(jié)點法。

分塊法將模型按某種準則劃分為塊體，假定每塊內(nèi)速度是常值。分塊法在進行射線追蹤和走時計算時簡便，但是在模型中人為的引入了塊體間的邊界，造成了速度的不連續(xù)性，并且結(jié)構(gòu)異常的輪廓只能表示成塊狀；節(jié)點法則是把節(jié)點上的速度值取為待求變量，其他任意點的速度值由于它相鄰的8個節(jié)點上的速度插值求得。與分塊法相比，節(jié)點法速度是連續(xù)變化的，沒有人為的邊界，可以獲得系數(shù)矩陣更密，節(jié)點法只能處理速度連續(xù)變化的情況，但無法處理存在速度間斷面的情況。

在實際工程勘探中，如地鐵施工中遇到的孤石、煤礦工作面中陷落柱、斷層等構(gòu)造異常情況，均勻地質(zhì)背景下會遇到速度突變的地質(zhì)問題。節(jié)點法在連續(xù)背景速度情況下對速度突變界面分辨率降低的情況，反演結(jié)果存在對初始模型依賴性強、易陷入局部極值等問題。在地震信號處理中，稀疏表示地震數(shù)據(jù)一直是近些年來發(fā)展的一個熱點問題[3-4]，采用盡可能少的非零元素表示信號的主要信息，從而簡化在地震信號處理中的求解過程[5-6]。

小波分析多尺度分解反演方法，將參數(shù)反演問題轉(zhuǎn)化到小波域中重要系數(shù)求解問題。利用多尺度之間的內(nèi)在聯(lián)系及小波域中小波系數(shù)的稀疏性，有效改進了局部極值、計算量等問題[7-8]。多尺度反演往往從較粗尺度開始，由于較粗尺度下目標函數(shù)呈現(xiàn)較強的凸性及較少的局部極值，有利于收斂至全局優(yōu)化解。進而有學者將小波變換為基礎(chǔ)的多尺度分析理論應(yīng)用于規(guī)則格點或方塊的模型參數(shù)化中[9-12]，并在反演過程中根據(jù)數(shù)據(jù)的分辨極限自動對模型進行多尺度再劃分。這種方法起到了不規(guī)則劃分網(wǎng)格的作用，先在粗尺度上迭代反演，得到一個比較好地參數(shù)估計，再將這個估計作為較小尺度的初始速度進行反演，直到反演出原問題的全局最優(yōu)解。其優(yōu)點是對初始模型選擇依賴性減小、更適合于大擾動異常體成像、圖像分辨率和質(zhì)量較高。

針對稀疏方程L0范數(shù)(向量中非零元素個數(shù))約束最小二乘解，但L0范數(shù)約束時方程求解為非凸函數(shù)求解，其解容易陷入局部極小值[13]。L1范數(shù)約束為L0與L2范數(shù)約束之間的一種折中且可行的方法，且L1范數(shù)約束時方程為凸函數(shù)求解，方程可求出全局最優(yōu)解。采用加權(quán)迭代最小二乘IRLS (Iteratively Reweighted Least Squares) 求解此類稀疏方程[14-16]，若求解模型尺度較大，則可使用迭代收縮算法ISA(Iterated Shrinkage Algorithms)[17-18]。

這里首先采用雙線性插值節(jié)點法離散速度模型；其次采用MSFM(Multi-Stencils Fast Marching Methods)正演計算走時場，并采用龍格庫塔方法求取路徑梯度求解最短射線路徑。在對速度結(jié)構(gòu)成像時，著重比較了L1全變差約束、L2范數(shù)約束、小波域多尺度稀疏約束，對比分析不同約束反演對不同模型分辨能力，指導實際地質(zhì)背景下的速度成像，為實際地質(zhì)問題提供更合理的解釋。

1 方法

1.1 地震走時方程建立

地震初值走時成像為通過射線的傳播時間來擬合地下介質(zhì)速度結(jié)構(gòu)的一種成像方法。在二維情況下，觀測地震初值到時可以表示為：

(1)

式中：Ti為第i條射線觀測到時;v(x,z)與s(x,z)為地下某一點處介質(zhì)速度慢度值;Γi為第i條射線路徑。

在做速度成像時，需要把速度模型離散成獨立的網(wǎng)格[19]，此時式(1)離散并寫成矩陣形式，見式(2)。

d=Gm

(2)

式中:d為觀測數(shù)據(jù)(N×1,N射線路徑數(shù));m為慢度(M×1，M為模型網(wǎng)格數(shù))，G為靈敏度矩陣(N×M)，即G的每一行為觀測數(shù)據(jù)對應(yīng)行對所有模型變量的偏導數(shù)。

若求解模型的最小長度,解為式(3)。

m=(GTG)-1GTd

(3)

1.2 模型與靈敏度矩陣的小波分解

小波是定義在有限間隔而且其平均值為零的一種函數(shù)，小波變換是信號f(t)與被縮放和平移的小波函數(shù)之積在信號存在的整個期間里求和，小波變換完成之后得到的系數(shù)是在不同的縮放因子下由信號的不同部分產(chǎn)生的。

對于連續(xù)小波變換，將任意L2(R)空間的函數(shù)f(t)在小波基下展開，成為函數(shù)的小波變換(小波分解)，其表達式為[20]：

WTf(a,τ)==

(4)

若小波變換滿足容許條件，則其逆變?yōu)槭?5)。

(5)

這里在對模型與靈敏度矩陣做小波分解時，選擇了Haar與Db2小波基。

1.3 小波域走時方程建立與求解

地震走時小波域多尺度反演成像步驟為，在建立靈敏度矩陣時采用基于Haar、Db2小波基的二維小波分解把靈敏度矩陣變換到小波域中，把在時域中求解速度模型的問題轉(zhuǎn)化到稀疏的多尺度小波系數(shù)的求取，并采用迭代過程中變權(quán)重法(IRLS)求解小波系數(shù)；最后采用小波逆變換把求解出的小波系數(shù)變量，變換到時域中的慢度更新量以達到對地下介質(zhì)速度結(jié)構(gòu)成像。

對于正交小波變換，有W-1=WT，WTW=I，I表示單位矩陣，式(2)可以寫為[7、14]：

Δd=GWTWΔm

(6)

(7)

1.4 二階吉洪諾夫與全變差正則化約束

井間地震波走時層析成像采用在孔中觀測的方式重建井間地層的速度結(jié)構(gòu)，因為我們需要根據(jù)孔中的觀測到時求解速度沿橫向與縱向的變化情況，導致式(2)嚴重欠定、條件數(shù)較大。必須采用正則化的方式減小解的假異常與不適定問題，一般情況下，為了權(quán)衡求解的穩(wěn)定性及假異常采用吉洪諾夫正則化與拉普拉斯變換約束。這兩種約束方法均為二范數(shù)約束。二范數(shù)約束目標函數(shù)為：

(8)

即折中模型平滑與因干擾異常產(chǎn)生的假異常，加上平滑約束參數(shù)相當于對模型求解方程增加了低通濾波器，使求解結(jié)果損失部分高頻信息，同時模糊了突變異常體的邊界，造成反演結(jié)果分辨率降低。隨著二范數(shù)(L2-norm)的廣泛應(yīng)用，最近開展了比較多的對Lp-norm約束的研究(0

(9)

在具體求解式(7)時，對其求導并令其一階導數(shù)為零可得[23]：

(2GTG+αLTWL)m=2GTd

(10)

式(8)寫成矩陣形式為：

(11)

采用IRLS迭代方式求解方程(11)，每次迭代分別重新計算W平滑權(quán)重系數(shù)[23]。

2 合成數(shù)據(jù)測試

我們設(shè)計了一個井-井觀測系統(tǒng)，模擬工程中經(jīng)常遇到的孤石或者溶洞情況。如圖1所示，孔距28 m、背景速度2 000 m/s，模型大小為32×64 m共計2 048個網(wǎng)格，兩個異常邊長為5 m×5 m，速度分別2 500 m/s、1 500 m/s。分別模擬充水充泥的低速溶洞與高速的孤石異常。從圖1中可以看出，速度模型變換到小波域中模型極為稀疏，2 048個小波系數(shù)中零值數(shù)目為1 937個，非零值111。即如果直接反演速度模型需要反演2 048個速度值，而在小波域中反演只需要反演111個小波系數(shù)值既可以完全恢復速度模型，使反演的目標模型較為稀疏。針對此種稀疏反演，這里采用IRLS求解稀疏模型方程，在迭代過程中動態(tài)選擇權(quán)重系數(shù)。

圖1 速度模型與其對應(yīng)的2階Haar小波系數(shù)Fig.1 Cross-borehole configuration and the 2rd Haar-Wavelet coefficient(a)速度模型；(b)速度模型的二階Haar小波變換后的小波系數(shù)

2.1 基于MSFM走時場計算

首先采用MSFM(Multi-Stencils Fast Marching Methods) 求解程函方程，正演計算初至到時。多模板快速推進算法是一種 FMM(Fast Marching Method) 的改進算法。該算法使用兩個模板分別用二階差分近似計算水平垂直方向以及對角線方向上網(wǎng)格節(jié)點的波前傳播時間，并選取滿足逆風條件的解，從而大大提高了 FMM 算法的計算精度和計算效率[24]。圖2為MSFM求解的波前走時場、反演迭代過程中射線路徑分布。

圖2 地震波波場與射線路徑Fig.2 Seismic time field and ray path(a)MSFM追蹤地震波走時場；(b)迭代過程中射線路徑

2.2 孤立異常模型

圖1中所示兩個孤立高低速異常，對應(yīng)于工程地質(zhì)中孤石或煤礦工作面溶洞構(gòu)造，其中高速對應(yīng)孤石，低速對應(yīng)溶洞。在具體求解過程中，采用四種反演策略，分別為在空間域中采用全變差L1-norm與二階吉洪諾夫正則化求解速度模型、小波域全變差求解基于Haar與Db2小波基的小波系數(shù)然后通過小波逆變換求解模型速度。初始速度模型均給定背景速度2 000 m/s。

圖3為不同約束條件下的速度模型反演結(jié)果。從反演結(jié)果可以看出，全變差約束與二階吉洪諾夫約束具有較為平滑的反演結(jié)果，為達到反演結(jié)果穩(wěn)定均產(chǎn)生了一定的假異常，但全變差約束具有更好的模型恢復及更少的假異常，對低速異常具有更好的分辨能力；基于Haar與Db2小波基的小波域反演比時空域反演具有更為干凈的背景，以及較少的干擾異常，但基于Haar小波基的反演產(chǎn)生了一定的拖尾異常，以及異常邊界過于尖銳，造成此情況的原因為Haar小波基為一個方波形態(tài)，過渡太迅速，高頻成分較多，針對此情況我們首先采取加大阻尼和平滑的權(quán)重，隨著權(quán)重的加大能一定程度壓制Haar小波的高頻成分，但也造成了對低速區(qū)域無法分辨的缺陷。從圖3(d)可以看出，一個方面減少了均勻背景上假異常的情況，同時具有較好的異常邊界分辨率。

圖3 地震波走時正反演Fig.3 Seismic travel time tomography inversion results(a) 全變差反演結(jié)果；(b) 吉洪諾夫二階正則化反演結(jié)果；(c) Haar小波域反演結(jié)果；(b) Db2小波域反演結(jié)果

圖4為不同約束條件下迭代殘差曲線圖，從圖4中可以看出，開始迭代時基于Haar小波基的反演具有最小的擬合誤差，二階吉洪諾夫具有最大的擬合誤差，隨著迭代次數(shù)的增加，4種不同的約束具有大致相同的擬合誤差，說明4種方式反演均具有隨著迭代次數(shù)增加擬合誤差減小及穩(wěn)定的收斂趨勢。同時也說明了在相同擬合誤差的情況下存在多種反演結(jié)果，即反演的不唯一性。

2.3 垂直與水平地層模型

圖5、圖6分別為豎直層狀與水平層狀速度模型反演結(jié)果，從圖5、圖6中可以看出，L1范數(shù)全變差與L2吉洪諾夫正則化約束下均具有較為平滑的結(jié)果，且全變差反演結(jié)果具有更少的假異常；小波域中反演的結(jié)果具有較好的邊界刻畫能力，其中基于Harr小波基的反演幾乎精確恢復速度模型，基于Db2小波基反演在異常邊界處也具有走向一致性。從反演結(jié)果數(shù)值恢復與邊界刻畫來看，空間域反演具有更好的數(shù)值恢復，小波域反演具有更好的邊界刻畫。

圖4 不同約束條件下迭代殘差曲線圖Fig.4 Iterative residual curves uneler different constraint conditions

圖5 豎直層狀模型反演Fig.5 Vertical layered velocity model(a)觀測系統(tǒng)；(b)速度模型；(c)全變差約束；(d)二階Tikhonov正則化約束；(e)稀疏約束-Haar小波基；(f)稀疏約束-Db2小波基

圖6 水平層狀模型反演Fig.6 Horizontal layered velocity model(a)觀測系統(tǒng)；(b)速度模型；(c)全變差約束；(d)二階Tikhonov正則化約束；(e)稀疏約束-Haar小波基；(f)稀疏約束-Db2小波基

2.4 傾斜斷層模型

圖7為傾斜斷層速度模型反演結(jié)果，從圖7中可以看出，L1范數(shù)全變差與L2吉洪諾夫正則化約束下均具有較為平滑的結(jié)果，邊界控制相對模糊，其中L1約束出現(xiàn)局部極值高速區(qū)域沒有L2約束下背景速度恢復的均勻平滑；針對傾斜速度結(jié)構(gòu)模型，小波域反演效果沒有孤立異常與層狀模型邊界刻畫的好，其反演結(jié)果均帶有一定的鋸齒狀形態(tài)存在，同時基于Db2小波基反演比基于Haar小波基具有更好的背景平滑度。

圖7 傾斜斷層模型反演Fig.7 Slope fault velocity modela)觀測系統(tǒng)；(b)速度模型；(c)全變差約束；(d)二階Tikhonov正則化約束；(e)稀疏約束-Haar小波基；(f)稀疏約束-Db2小波基

3 結(jié)論與討論

在井-井地震走時層析成像反演的基礎(chǔ)上，討論了小波域多尺度反演與空間域L1-norm、L2-norm約束反演下，對不同模型結(jié)構(gòu)的分辨能力，取得了以下幾點的認識：

1)對于塊狀速度結(jié)構(gòu)，在空間域欠定的迭代方程在小波域具有更好的完備性及穩(wěn)定性。

2)全變差L1-norm約束與二階L2-norm吉洪諾夫約束相比反演結(jié)果具有更少的假異常，在對孤立異常反演時，全變差約束具有更好的邊界識別及低速區(qū)域數(shù)值恢復。

3)針對層狀地層模型，基于Haar小波基的多尺度約束下成像具有準確的邊界分辨力。

4)L2-norm約束對不同速度結(jié)構(gòu)模型具有折中的成像效果。

5)針對不同的地質(zhì)模型，不同反演方法分別具有優(yōu)缺點，在對實際地質(zhì)資料反演成像時，因預先不知道地層結(jié)構(gòu)，一個推薦的方法為采用多種反演約束方法分別反演，對比之間的區(qū)別，綜合推斷一個更為合理的地層結(jié)構(gòu)解釋。

[1]楊文采,李幼銘.應(yīng)用地震層析成像[M].北京:地質(zhì)出版社,1993.YANG W C,LI Y M.Application of seismic tomography [M].Beijing:Geological publishing house,1993.(In Chinese)

[2]雷棟,胡祥云.地震層析成像方法綜述[J].地震研究,2006.29(4):418-425.LEI D,HU X Y.Seismic tomography methods review [J].seismic study,2006.29(4):418-425.(In Chinese)

[3]YUAN S,WANG S.Spectral sparse Bayesian learning reflectivity inversion[J].Geophysical Prospecting,2013,61(4):735-746.

[4]YUAN S,WANG S,LUO C,et al.Simultaneous multitrace impedance inversion with transform-domain sparsity promotion[J].Geophysics,2015,80(2):R71-R80.

[5]BRUCKSTEIN A M,DONOHO D L,ELAD M.From sparse solutions of systems of equations to sparse modeling of signals and images[J].SIAM review,2009,51(1):34-81.

[6]馮飛.結(jié)合稀疏變換的稀疏約束反演一次波估計研究[D].長春：吉林大學,2014.FENG F.Combined with a wave of sparse constraint inversion of the sparse transform estimation research [D].Changchun：Jilin University,2014.(In Chinese)

[7]馬堅偉,楊慧珠.地震波形多尺度反演的一點討論[J].地球物理學進展,2000,15(4):55-61.MA J W,YANG H Z.Multi-scale inversion of seismic waveform is discussed [J].Progress in Geophysics,2000,15(4):55-61.(In Chinese)

[8]馬堅偉,朱亞平,楊慧珠.二維地震波形小波多尺度反演[J].工程數(shù)學學報,2004,21(1):109-113.MA J W,ZHU Y P,YANG H Z.Two-dimensional wavelet multiscale seismic waveform inversion[J].Chinese Journal of Engineering Mathematics,2004,21(1):109-113.(In Chinese)

[9]裴正林,牟永光.復雜介質(zhì)小波多尺度井間地震層析成像方法研究[J].地球物理學報,2003,46(1):113-117.PEI Z L,MOU Y G.Complex medium wavelet multi-scale method of crosshole seismic tomography study[J].Chinese Journal of Geophysics,2003,46(1):113-117.(In Chinese)

[10]宋常瑜,裴正林,狄?guī)妥?等.小波多尺度井間地震衰減層析成像[J].石油地球物理勘探,2007,42(1):30-33.SONG C Y,PEI Z L,DI B R,et al.The wavelet multi-scale cross-hole seismic attenuation tomography[J].Geophysical prospecting for petroleum,2007,42(1):30-33.(In Chinese)

[11]NAZZAL,M.,& OZKARAMANLI,H.Wavelet domain dictionary learning-based single image superresolution.Signal,Image and Video Processing,2015,9(7):1491-1501.

[12]TIKHOTSKY,S.,& ACHAUER,U.Active Seismic Tomography Inversion with the Self Adaptive Wavelet Parameterization:Algorithm and its Application for the Vesuvius Volcano Structure[C].In Proc.7th Int.Conf.on Problems of Geocosmos,St.Petersburg,Russia ,2008:251-252.

[13] BRUCKSTEIN A M,DONOHO D L,ELAD M.From sparse solutions of systems of equations to sparse modeling of signals and images[J].SIAM review,2009,51(1):34-81.

[14]KARLOVITZ L A.Construction of nearest points in the L p,p even,and L∞ norms[J].Journal of Approximation Theory,1970,3(2):123-127.

[15]RAO B D,ENGAN K,COTTER S F,et al.Subset selection in noise based on diversity measure minimization[J].Signal Processing,IEEE Transactions on,2003,51(3):760-770.

[16]FANG H,ZHANG H.Wavelet-based double-difference seismic tomography with sparsity regularization[J].Geophysical Journal International,2014,199(2):944-955.

[17]KOLACZYK E D.A wavelet shrinkage approach to tomographic image reconstruction[J].Journal of the American Statistical Association,1996,91(435):1079-1090.

[18]FIGUEIREDO M A T,NOWAK R D.An EM algorithm for wavelet-based image restoration[J].Image Processing,IEEE Transactions on,2003,12(8):906-916.

[19]AKI K,LEE W H K.Determination of three-dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes:1.A homogeneous initial model[J].Journal of Geophysical research,1976,81(23):4381-4399.

[20]郝金來,姚振興.小波域反演確定區(qū)域地震震源機制[J].中國科學地球科學 (中文版),2012,42(2):191-201.HAO J L,YAO Z X.The wavelet domain inversion to determine the regional earthquake focal mechanism [J].Chinese science and earth sciences,2012,42(2):191-201.(In Chinese)

[21]BEHROOZ A,ZHOU H M,EFTEKHAR A A,et al.Total variation regularization for 3D reconstruction in fluorescence tomography:experimental phantom studies[J].Applied optics,2012,51(34):8216-8227.

[22]YUAN S,WANG S,LI G.Random noise reduction using Bayesian inversion[J].Journal of Geophysics and Engineering,2012,9(1):60-68.

[23]ASTER R C,BORCHERS B,THURBER C H.Parameter estimation and inverse problems[M].Academic Press,2011.

[24]王飛,曲昕馨,劉四新,等.一種新的基于多模板快速推進算法和最速下降法的射線追蹤方法[J].石油地球物理勘探,2014,49(6):1106-1114.WANG F,QU X X,LIU S X,et al..A new rapid advance algorithm based on template and the steepest descent method of ray tracing method [J].Geophysical prospecting for petroleum,2014,49(6):1106-1114.(In Chinese)

Compared 1-Norm / 2-Norm constraint and wavelet multi-scale 2D body wave velocity tomography

GAO Ji，F(xiàn)ANG Hong-jian

(Laboratory of Seismology and Physics of Earth's Interior,University of Science and Technology of China,Hefei 230026,China)

Because of the limited coverage of observation system and the complex physical relationship between physical parameters and observations,near surface geophysical method suffers issues of non-uniqueness and resolution limitation to some degree.For seismic travel time tomography,number distribution of ray paths in low speed region is less than the high speed region,and this phenomenon leading to the seismic travel time tomography have a low resolution for low speed region.At the same time,the measured data affected by many complicated environmental factors usually include strong noises and the limitations of different inversion method,which badly force the inversion result contain artifacts.We main studied tomography with different constraints of spatial domain of L1 norm,L2 norm and wavelet coefficients of wavelet domain.Compared the different inversion methods and constraint conditions on the resolution of geological model.Through test isolated anomaly model,layered stratum model,tilted strata model,three different inversion methods of different specific model distinguish characteristics respectively.In the actual data imaging,can through the use of multiple inversion methods,compare the inversion result rationality in order to achieve more accurate geological structure interpretation.

body wave travel time; tomography; 1-Norm; 2-Norm; wavelet multi-scale

2015-10-20 改回日期：2015-10-26

國家自然科學基金項目(41474039)

高級(1983-)，男，博士，主要從事地震資料處理，彈性波層析成像、電阻率層析成像、多屬性聯(lián)合反演研究,E-mail：gaoji617@mail.ustc.edu.cn。

1001-1749(2016)06-0765-08

P 631.4