基于Huber函數(shù)的頻率域全波形反演

呂曉春，顧漢明，成景旺

（1.中國(guó)地質(zhì)大學(xué)（武漢）地球物理與空間信息學(xué)院，湖北武漢430074；2.中國(guó)地質(zhì)大學(xué)（武漢）構(gòu)造與油氣資源教育部重點(diǎn)實(shí)驗(yàn)室，湖北武漢430074）

利用人工地震獲取的地震波信息反演地下介質(zhì)物性參數(shù)是地球物理學(xué)的重要研究?jī)?nèi)容。地震波場(chǎng)反演方法很多，如相位反演成像、振幅反演成像以及全波形反演成像等等。其中全波形反演不同于走時(shí)反演、振幅反演等傳統(tǒng)的成像方法，它利用地震記錄的整體信息來進(jìn)行反演，是一種高精度的成像方法。20世紀(jì)80年代，直接通過求取模型參數(shù)擾動(dòng)計(jì)算Frechet導(dǎo)數(shù)的做法工作量極為龐大，使全波形反演局限于一維模型或簡(jiǎn)單模型。Tarantola［1］借鑒逆時(shí)偏移的思想，利用正向傳播波場(chǎng)與殘差波場(chǎng)的互相關(guān)求取梯度方向，進(jìn)行了時(shí)間域的全波形反演研究。Pratt等［2－4］將全波形反演引入到頻率域中，根據(jù)炮點(diǎn)和接收點(diǎn)可互易性原理，導(dǎo)出了求解目標(biāo)函數(shù)梯度的高效算法，并利用高斯牛頓法對(duì)非均勻介質(zhì)的速度進(jìn)行了全波形反演。至此，地球物理學(xué)者開始對(duì)全波形反演開展更深入的研究。國(guó)外近十年中頻率域全波形反演已被廣 泛 應(yīng) 用 于 建 立 高 分 辨 率 的 速 度 模 型［5－7］。Virieux［8］總結(jié)了當(dāng)前全波形反演所取得的一些成果和認(rèn)識(shí)，指出了該方法所面對(duì)的難點(diǎn)以及下一步研究方向和趨勢(shì)。國(guó)內(nèi)很多學(xué)者也做了大量的研究工作，吳國(guó)忱等［9－10］以及殷文等［11］對(duì)頻率域的聲波和彈性波有限差分正演模擬做了很詳細(xì)的研究，為頻率域反演提供了前提手段；許焜等［12－13］基于頻域率有限元方法進(jìn)行了反演研究；龍桂華等［14］對(duì)粘彈聲波介質(zhì)進(jìn)行了頻率域反演；卞愛飛等［15］提出了自適應(yīng)剝層法全波形反演。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，以及地震處理要求的提高，全波形反演方法的高分辨率特性引起了人們的注意。然而，實(shí)際地震信號(hào)的低頻缺失以及噪聲的存在制約著全波形反演的應(yīng)用，其中噪聲的存在是一個(gè)非常關(guān)鍵的問題，因此有必要找到一個(gè)合適的目標(biāo)函數(shù)。L1范數(shù)相對(duì)于L2范數(shù)抗噪性強(qiáng)［16－17］，但數(shù)據(jù)誤差接近于零時(shí)，L1范數(shù)準(zhǔn)則求取的目標(biāo)函數(shù)梯度會(huì)出現(xiàn)不穩(wěn)定性，而L2范數(shù)則比較穩(wěn)定，不會(huì)出現(xiàn)這個(gè)問題。Guitton［18］在使用擬牛頓法進(jìn)行地震反演時(shí)，通過求取一個(gè)實(shí)數(shù)型的Huber目標(biāo)函數(shù)的梯度來實(shí)現(xiàn)。Ha等［19］基于Huber函數(shù)構(gòu)建反演目標(biāo)函數(shù)來提高算法的抗噪性和穩(wěn)定性。韓淼等［20］也進(jìn)行了反演目標(biāo)函數(shù)的幾種范數(shù)準(zhǔn)則的對(duì)比研究。

為了既滿足穩(wěn)定性，又滿足抗噪性，我們將L1范數(shù)和L2范數(shù)結(jié)合來建立目標(biāo)函數(shù)，提出一種基于Huber函數(shù)的頻率域全波形反演新方法。頻率域全波形反演是在復(fù)數(shù)域進(jìn)行的，因此引入了復(fù)數(shù)形式的Huber函數(shù)，并推導(dǎo)了基于Huber準(zhǔn)則的目標(biāo)函數(shù)的梯度表達(dá)式。通過在幾個(gè)簡(jiǎn)單的合成記錄中加入隨機(jī)脈沖噪聲、高斯噪聲和連續(xù)噪聲，驗(yàn)證了該方法的應(yīng)用效果。

1 復(fù)數(shù)形式的Huber函數(shù)

Huber［21］給出實(shí)數(shù)形式的Huber函數(shù)表達(dá)式：

其中，r為計(jì)算值與觀測(cè)值的誤差，是一個(gè)實(shí)數(shù)變量；ε是一個(gè)閾值，用來判定不同大小的誤差所對(duì)應(yīng)的目標(biāo)函數(shù)表達(dá)式（圖1）。類似實(shí)數(shù)域的Huber函數(shù)，對(duì)于復(fù)數(shù)r＝a＋bi，我們定義兩種不同的范數(shù)，其中L1范數(shù)計(jì)算表達(dá)式為，L2范數(shù)計(jì)算表達(dá)式為，就可以得到復(fù)數(shù)域的Huber函數(shù)表達(dá)式：

圖1 Huber函數(shù)示意

從公式（2）中可以看出，在誤差大于閾值ε時(shí)，目標(biāo)函數(shù)選取L1范數(shù)準(zhǔn)則。我們知道L1范數(shù)是所有范數(shù)中對(duì)數(shù)據(jù)加權(quán)最小的［22］，如果數(shù)據(jù)非常精確，那么某一預(yù)測(cè)值遠(yuǎn)離其觀測(cè)值的地方便是重要的影響因素，這時(shí)應(yīng)選擇L2范數(shù)，因?yàn)樗鼘?duì)較大的誤差加權(quán)也大。而對(duì)于包含噪聲的合成記錄來說，噪聲的存在使得數(shù)據(jù)不是很準(zhǔn)確，應(yīng)使用L1范數(shù)來減小這部分較大誤差的權(quán)值。由于L1范數(shù)總是大于L2范數(shù)，故（2）式還可以改寫成本文反演中采用公式（3）來建立目標(biāo)函數(shù)，其中的閾值ε的選取非常關(guān)鍵，我們采用Bube等［23］提出的方法，閾值ε的大小取為計(jì)算值與觀測(cè)值誤差數(shù)列的標(biāo)準(zhǔn)差的0.6倍。

2 頻域率全波形反演基本理論

2.1 頻率域波動(dòng)方程正演計(jì)算

頻率域聲波波動(dòng)方程可以表示為

式中：ρ（x）為密度；k（x）為體變模量；ω 為角頻率；u（x，ω）和s（x，ω）分別為波場(chǎng)壓力值與震源。將方程（4）加入PML邊界條件后采用有限差分算子將其離散，則（4）式可簡(jiǎn)單表示為

式中：A（ω）為與頻率和介質(zhì)速度有關(guān)的阻抗矩陣。假如模型空間網(wǎng)格數(shù)為（Nx，Nz），那么A（ω）是一個(gè)（Nx×Nz）×（Nx×Nz）的有限差分算子矩陣。該系數(shù)矩陣為大型稀疏矩陣，其非零元素對(duì)稱地分布在對(duì)角線上，但是它并非對(duì)稱矩陣，其值的大小是不一樣的。求解該大型稀疏矩陣方程組常用的方法有迭代法和直接分解法。對(duì)于多震源問題，我們采用直接分解（如LU分解）來求解，因?yàn)橄禂?shù)矩陣A（ω）只需分解一次，其LU分解結(jié)果就可以用于更多的新的震源求解，可以大量節(jié)省計(jì)算量，而且在全波形反演中的迭代反演上更為重要，因?yàn)槠溆写罅康恼鹪匆约疤撜鹪辞蠼?，每一次迭代過程中系數(shù)矩陣A（ω）不變，其LU分解可以多次重復(fù)使用。

2.2 目標(biāo)函數(shù)的建立以及梯度的求解

反問題的求解在數(shù)學(xué)上其實(shí)是一個(gè)最優(yōu)化求解問題，其中目標(biāo)函數(shù)的構(gòu)造是一個(gè)關(guān)鍵。根據(jù)上述復(fù)數(shù)形式的Huber表達(dá)式，我們構(gòu)造如下的目標(biāo)函數(shù)：

式中：p為模型變量；Ns為炮點(diǎn)數(shù)；Nr為接收點(diǎn)數(shù)；rs，r為計(jì)算波場(chǎng)與觀測(cè)波場(chǎng)的殘差：

結(jié)合（8）式和（9）式，我們將基于 Huber函數(shù)得到的目標(biāo)函數(shù)梯度表達(dá)式統(tǒng)一表示為

式中：r＝（r1，r2，…，rNr），

從（10）式中我們看到梯度求解的關(guān)鍵在于Jacobin矩陣的求解。對(duì)（5）式兩邊參數(shù)pk進(jìn)行求導(dǎo)，有

由于（13）式和（5）式的形式一樣，而且它的阻抗矩陣也一樣，故可以直接利用正演的LU分解結(jié)果求解偏導(dǎo)數(shù)波場(chǎng)。將（13）式代入（10）式可得

其中A［］－1Tr＊為殘差的反傳播場(chǎng)，我們只需將各個(gè)接收點(diǎn)處的殘差組成新的震源再正演一次便可得到。從（14）式可以看出目標(biāo)函數(shù)對(duì)速度變量的梯度可以通過計(jì)算虛震源與殘差反傳播場(chǎng)的零延遲卷積得到。

2.3 預(yù)梯度法反演

地球物理反演通過求取目標(biāo)函數(shù)的最小值，不斷對(duì)模型參數(shù)進(jìn)行修改，求取滿足條件的最佳模型。本文反演過程中采用近似Hessian矩陣的對(duì)角線來做梯度的預(yù)條件算子［14］，其模型更新量可以表示為：

式中：α為待定的步長(zhǎng)；diag（Ha）表示近似 Hessian矩陣Ha的對(duì)角線矩陣，用來對(duì)梯度進(jìn)行修正，以加快收斂速度；λ為可調(diào)節(jié)的正則化因子；I表示單位矩陣為目標(biāo)函數(shù)的梯度。其中Ha可以表示為

從（16）式不難發(fā)現(xiàn)Hessian矩陣可以表示成偏導(dǎo)數(shù)的互相關(guān)，根據(jù)（13）式可求出Hessian矩陣，令i＝j(luò)就可求出它的對(duì)角元素。

3 數(shù)值模型試算

為了驗(yàn)證基于Huber函數(shù)的頻率域全波形反演方法的實(shí)用性以及正確性，在原始合成記錄上分別加入了隨機(jī)脈沖噪聲、連續(xù)噪聲和高斯噪聲來進(jìn)行模型試算。

3.1 含隨機(jī)脈沖噪聲的數(shù)值試算

通過一個(gè)簡(jiǎn)單的地質(zhì)模型來進(jìn)行隨機(jī)脈沖噪聲的數(shù)值試算。模型計(jì)算范圍為5km×2km，一共有5層，中間有一個(gè)高速體（圖2）。原始合成記錄通過頻率域有限差分正演模擬得到，正演過程中采用主頻為10Hz的雷克子波，頻率間隔為0.24Hz，最大頻率為30Hz。共有51炮，炮間距100m。每炮有251道接收，道間距為20m。圖3a為炮點(diǎn)位置位于0處的原始合成記錄，圖3b為加了噪聲的記錄，其中有9個(gè)充零道和4個(gè)脈沖噪聲。

圖2 簡(jiǎn)單模型

圖3 震源位于0處的模擬炮集記錄

圖4 含隨機(jī)脈沖噪聲數(shù)據(jù)反演結(jié)果

分別采用基于L2范數(shù)的最小平方反演法和本文提出的Huber反演方法對(duì)不含噪聲和含噪聲的記錄進(jìn)行反演。反演過程中從2.4到20.0Hz共選取20個(gè)反演頻率，每個(gè)頻率反演40次。反演采用的初始模型為一個(gè)線性模型（圖4a），初始速度為1 500m／s，最深處的速度為3 000m／s，速度計(jì)算為v＝1 500＋kz，k為15m／s，z為網(wǎng)格點(diǎn)數(shù)。圖4b和4c分別為采用最小平方法對(duì)原始合成記錄（圖3a）和含隨機(jī)噪聲數(shù)據(jù)（圖3b）的反演結(jié)果。由圖4b可以看到，最小平方法對(duì)沒有噪聲的數(shù)據(jù)進(jìn)行反演可以得到較好的速度模型，根據(jù)Snell射線原理，高速體下方透射波較少，故該界面與理論有所差異，但模型的大部分層位都較好地反映出來。相對(duì)而言，對(duì)含有噪聲的數(shù)據(jù)最小平方法反演結(jié)果較差（圖4c），除了第一層較好地被反演出來，其他各層的分界面都反演較差。圖4d為采用本文方法對(duì)含隨機(jī)噪聲數(shù)據(jù)的反演結(jié)果，對(duì)比圖4c和圖4d可看出，基于Huber函數(shù)的全波形反演受噪聲影響小，反演結(jié)果較好。

3.2 含連續(xù)噪聲的數(shù)值試算

采用和3.1節(jié)同樣的模型（圖2）、觀測(cè)系統(tǒng)以及正、反演參數(shù)。如圖5所示，在原始合成單炮記錄（圖5a）上加入了2個(gè)不同斜率的連續(xù)噪聲（圖5b和圖5c）。為了和原始記錄的頻譜相一致，連續(xù)噪聲同樣采用主頻為10Hz的雷克子波，其中連續(xù)噪聲的振幅值隨著偏移距而改變［24］。

從2個(gè)含連續(xù)噪聲數(shù)據(jù)的反演結(jié)果（圖6）來看，連續(xù)噪聲對(duì)反演結(jié)果的影響較大。第一種斜率噪聲存在的情況下，利用最小平方反演法（圖6a）沒有將高速體反演出來，只反演出上界面，而且出現(xiàn)了假的層位；采用本文方法（圖6c）可以大致反演出高速體，而且下面也沒有出現(xiàn)假的層位。第二種噪聲存在的情況下，最小反演法反演（圖6b）的第一層的深度與理論不吻合，高速體同樣沒有反演出來；而采用本文方法得到的反演結(jié)果（圖6d）第一層的埋深有所改善，且高速體也被反演出來，但與不含噪聲反演結(jié)果（圖4b）相比，還有所差別。

圖5 震源位于0處的模擬炮集記錄

圖6 含連續(xù)噪聲數(shù)據(jù)反演結(jié)果

3.3 含高斯噪聲的數(shù)值試算

我們采用Marmousi模型（圖7）的合成記錄加入高斯噪聲來進(jìn)行反演測(cè)試。Marmousi模型計(jì)算范圍取為9.00km×3.18km，網(wǎng)格大小為451×159。原始合成記錄也是通過頻率域有限差分正演模擬得到的，正演過程中采用主頻為10Hz的雷克子波，頻率間隔為0.24Hz，正演最大頻率為30Hz。共有91炮，炮間距為100m。每一炮451道接收，道間距為20m。

圖8分別給出了沒有加噪聲（圖8a）和加了信噪比為5dB的高斯噪聲（圖8b）的合成記錄。高斯噪聲的產(chǎn)生采用Box－Muller方法［25］。

圖7 Marmousi模型

圖8 震源位于3 100m處的模擬炮集記錄

分別采用最小平方反演法和本文提出的Huber反演方法對(duì)圖8所示不含噪聲和含有噪聲的原始合成記錄進(jìn)行反演。反演過程中從2.4到20.0Hz共選取20個(gè)反演頻率，每個(gè)頻率反演40次。反演采用的初始模型為一個(gè)線性模型（圖9a），速度計(jì)算為v＝1 500＋kz，k為17m／s，z為網(wǎng)格點(diǎn)數(shù)。圖9b和圖9c分別為采用最小平方法對(duì)原始合成記錄和含高斯噪聲數(shù)據(jù)的反演結(jié)果，圖9d為本文方法對(duì)含高斯噪聲數(shù)據(jù)的反演結(jié)果。從3種反演結(jié)果來看，最小平方法對(duì)高斯噪聲敏感度較低，合成記錄中加入高斯噪聲也可以得到好的反演結(jié)果。由此可以看出最小平方法和Huber函數(shù)反演法受高斯噪聲的影響都較小，除了速度模型光滑度較差外，各個(gè)層的分界面均可以較好地反演出來。

圖9 含高斯噪聲數(shù)據(jù)反演結(jié)果

4 結(jié)束語

針對(duì)傳統(tǒng)的全波形反演采用L2范數(shù)準(zhǔn)則構(gòu)建目標(biāo)函數(shù)抗噪性差和L1范數(shù)則可能導(dǎo)致反演不穩(wěn)定的問題，我們引入復(fù)數(shù)形式的Huber函數(shù)準(zhǔn)則來建立目標(biāo)函數(shù)，同時(shí)考慮L2范數(shù)和L1范數(shù)的優(yōu)勢(shì)，提出了一種頻率域全波形反演新方法。通過以上理論研究以及加入不同噪聲的數(shù)值模型數(shù)據(jù)試算，可以得到以下認(rèn)識(shí)：

1）在頻率域全波形反演中，基于L2范數(shù)的最小平方法反演在無噪聲的情況下能夠得到精確的地下模型結(jié)構(gòu)。但是L2范數(shù)對(duì)噪聲的敏感度較高，在噪聲存在的情況下不能得到理想的反演結(jié)果。

2）基于Huber函數(shù)準(zhǔn)則建立的目標(biāo)函數(shù)，具有L1范數(shù)的抗噪性以及L2范數(shù)的穩(wěn)定性。在噪聲存在的情況下，依然能夠得到較好的反演結(jié)果。

3）L2范數(shù)對(duì)隨機(jī)脈沖噪聲和連續(xù)噪聲敏感度較高，而對(duì)高斯噪聲的敏感度相對(duì)來說較低。因此，對(duì)加入高斯噪聲的合成記錄，基于L2范數(shù)和基于Huber函數(shù)的方法反演效果均不錯(cuò)；而對(duì)加入隨機(jī)脈沖噪聲和連續(xù)噪聲的合成記錄，基于Huber函數(shù)反演要比基于L2范數(shù)反演效果好。

［1］Tarantola A.Inversion of seismic reflection data in the acoustic approximation［J］.Geophysics，1984，49（8）：1259－1266

［2］Pratt R G，Shin C S，Hicks G J.Guass－Newton and full Newton methods in frequency－space seismic waveform inversion［J］.Geophysics，1998，63（2）：341－362

［3］Pratt R G.Seismic waveform inversion in the frequency domain，Part 1：theory and verification in a physical scale model［J］.Geophysics，1999，64（3）：888－901

［4］Pratt R G.Seismic waveform inversion in the frequency domain，Part 2：fault delineation in sediments using crosshole data［J］.Geophysics，1999，64（3）：902－914

［5］Hicks G J，Pratt R G.Reflection waveform inversion using local decent methods：Estimating attenuation and velocity over a gas－sand deposit［J］.Geophysics，2001，66（2）：598－612

［6］Malinowski M，Operto S.Quantitative imaging of the Permo－Mesozoic complex and its basement by frequency domain waveform tomography of wide－aperture seismic data from the Polish Basin［J］.Geophysical Prospecting，2008，56（6）：805－825

［7］Ben－Hadj－Ali H，Operto S，Virieux J.Velocity model building by 3Dfrequency－domain，full－waveform inversion of wide－aperture seismic data［J］.Geophysics，2008，73（5）：VE101－VE117

［8］Virieux J，Operto S.An overview of full－waveform inversion in exploration geophysics［J］.Geophysics，2009，74（6）：WCC1－WCC26

［9］吳國(guó)忱，梁鍇.VTI介質(zhì)頻率－空間域準(zhǔn)P波正演模擬［J］.石油地球物理勘探，2005，40（5）：535－545 Wu G C，Liang K.Quasi P－wave forward modeling in frequency－space domain in VTI media［J］.Oil Geophysical Prospecting for Petroleum，2005，40（5）：535－545

［10］吳國(guó)忱，梁鍇.VTI介質(zhì)準(zhǔn)P波頻率空間域組合邊界條件研究［J］.石油物探，2005，44（4）：301－307 Wu G C，Liang K.Combined boundary conditions of quasi P－wave within frequency－space domain in VTI media，Oil Geophysical Prospecting for Petroleum，2005，44（4）：301－307

［11］殷文，印興耀，吳國(guó)忱，等.高精度頻率域彈性波方程有限差分方法及波場(chǎng)模擬［J］.地球物理學(xué)報(bào)，2006，49（2）：561－568 Yin W，Yin X Y，Wu G C，et al.The method of finite difference of high precision elastic wave equation in the frequency domain and wavefiled simulation［J］.Chinese Journal Geophysics（in Chinese），2006，49（2）：561－568

［12］許琨，王妙月.聲波方程頻率域有限元參數(shù)反演［J］.地球物理學(xué)報(bào)，2001，44（6）：852－864 Xu K，Wang M Y.Finite element inversion of the coefficients of acoustic equation in frequency domain［J］.Chinese Journal Geophysics（in Chinese），2001，44（6）：852－864

［13］許琨，王妙月.利用地質(zhì)規(guī)則塊體建模方法的頻率域有限元彈性波速度反演［J］.地球物理學(xué)報(bào)，2004，47（4）：708－717 Xu K，Wang M Y.Frequency－domain finite－element inversion of elastic－wave velocity using the geological regular－blocky－model method［J］.Chinese Journal Geophysics（in Chinese），2004，47（4）：708－717

［14］龍桂華，李小凡，張美根，等.頻率域粘彈性聲波透射波形速度反演［J］.地震學(xué)報(bào)，2009，31（1）：32－41 Long G H，Li X F，Zhang M G，et al.Visco－acoustic transmission waveform inversion for velocity structure in space－frequency domain［J］.Acta Seismologica Sinica，2009，31（1）：32－41

［15］Bian A F，Yu W H.Layer－Stripping full waveform inversion with damped seismic reflection data［J］.Journal of Earth Science，2011，22（2）：241－249

［16］Claerbout J，Muir J.Robust modeling with erratic data［J］.Geophysics，1973，38（5）：326－344

［17］Crase E，Pica A M，Noble M，et al.Robust elastic nonlinear waveform inversion：application to real data［J］.Geophysics，1990，55（5）：527－538

［18］Guitton A，Symes W.Robust inversion of seismic data using the Huber norm［J］.Geophysics，2003，68（4）：1310－1318

［19］Ha T，Chung W，Shin C.Waveform inversion using a back－propagation algorithm and a Huber function norm［J］.Geophysics，2009，74（3）：R15－R24

［20］韓淼，韓立國(guó)，崔杰，等.頻率域全波形反演目標(biāo)函數(shù)的幾種范數(shù)準(zhǔn)則的對(duì)比研究［J］.中國(guó)地球物理學(xué)會(huì)二十七屆年會(huì)論文集，2011，730 Han M，Han L G，Cui J，et al.The contrastive research of several norm－criterions of the objective function in FWI［J］.Chinese Geophysical Society 27thannual meeting proceedings，2011，730

［21］Huber J.Robust regression：Asymptotics，conjectures and Monte Carlo：annals of statistics［J］.Geophysics，1973，38（10）：799－821

［22］姚姚.地球物理反演基本原理與理論研究［M］.湖北武漢：中國(guó)地質(zhì)大學(xué)出版社，2002：17－18 Yao Y.Geophysical inversion of the basic theory and method of application［M］.Wuhan：China University of Geosciences Press，2002：17－18

［23］Bube K，Nemeth T.Fast line searches for the robust solution of linear systems in the hybrid L1／L2and Huber norms［J］.Geophysics，2007，72（2）：A13－A17

［24］Lu W，Zhang W，Liu D.Local linear coherent noise attenuation based on local polynomial approximation［J］.Geophysics，2006，71（6）：V163－V169

［25］Box G，Muller M.A note on the generation of random normal deviates［J］.Annals of Mathematical Statistics，1958，29（2）：610－611