基于梯度投影法的全變差正則化全波形反演

姚含，徐海

(貴州省地質(zhì)環(huán)境監(jiān)測(cè)院，貴州 貴陽 550001)

0 引言

我國(guó)油氣資源勘探從構(gòu)造勘探階段逐步走向巖性勘探階段，勘探難度逐漸加大。地震勘探具有勘探深度大、精度高、數(shù)據(jù)量大等特點(diǎn)，是目前油氣勘探最有效的方法，應(yīng)用廣泛。為進(jìn)一步提高地震勘探精度，地震全波形反演方法迅速發(fā)展起來，成為國(guó)內(nèi)外學(xué)者研究的熱點(diǎn)。

全波形反演(full waveform inversion，F(xiàn)WI)是一種高精度、高分辨率的速度反演方法[1-2]。FWI利用地震資料中的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)信息，通過最小化模擬數(shù)據(jù)與實(shí)際觀測(cè)數(shù)據(jù)之間的誤差來重構(gòu)地下介質(zhì)模型[3-4]，能得到地下介質(zhì)中更多細(xì)節(jié)的參數(shù)變化特征[5]。由于實(shí)際觀測(cè)中，子波未知，數(shù)據(jù)頻帶有限，且通常存在一定程度的噪聲干擾，波場(chǎng)與介質(zhì)參數(shù)之間呈現(xiàn)較強(qiáng)的非線性關(guān)系，F(xiàn)WI體現(xiàn)為高度欠定問題，對(duì)初始模型的精度要求較高，導(dǎo)致FWI在實(shí)際應(yīng)用中達(dá)不到預(yù)期的效果。

為了減少FWI的不適定性，加入先驗(yàn)?zāi)Ｐ图s束的正則化項(xiàng)，能使反演具有更好的穩(wěn)定性[6-7]。全變差(total variation，TV)約束方法是一種非光滑約束，由Rudin等[8]在1992年提出，主要用于圖像的去噪處理，能有效處理不連續(xù)性解的求解問題，重構(gòu)圖像中的一些間斷點(diǎn)，保留邊緣信息。近幾年來，全變差約束已經(jīng)應(yīng)用在電法成像[9]、圖像去噪[10-11]、偏移速度分析[12-13]、混合數(shù)據(jù)偏移[14]等領(lǐng)域。地下地層通常存在較多復(fù)雜的構(gòu)造變化帶及巖性突變帶，導(dǎo)致地層物理參數(shù)在這些地方產(chǎn)生不連續(xù)變化，應(yīng)用基于TV約束的FWI可以有效保留這些不連續(xù)特征[15]。

本文采用梯度投影法[16]，構(gòu)建基于全變差約束的全波形反演目標(biāo)函數(shù)，實(shí)現(xiàn)多約束集下的聲波方程全波形反演，逐步搜索到全局極值點(diǎn)，避免陷入局部極值，提高全波形速度反演的精度與速度。

1 方法原理

FWI是一個(gè)最優(yōu)化問題，通過對(duì)觀測(cè)波場(chǎng)與理論波場(chǎng)的殘差求取極小值來構(gòu)造目標(biāo)函數(shù)進(jìn)行迭代更新，最終重構(gòu)地下介質(zhì)的參數(shù)模型。當(dāng)前常用的目標(biāo)函數(shù)主要是通過數(shù)據(jù)殘差的l2范數(shù)來構(gòu)建，即：

(1)

我們考慮模型參數(shù)m，它對(duì)應(yīng)于速度平方的倒數(shù)，是一個(gè)實(shí)向量m∈N，其中N是空間離散化中的點(diǎn)數(shù)。對(duì)于每個(gè)源和頻率，波場(chǎng)和觀測(cè)數(shù)據(jù)分別用usv∈N和dobs∈Nr表示，其中，s=1,…,Ns指源，v=1,…,Nv指頻率，Nr為接收器的數(shù)量。dcal=Pusv，P是將波場(chǎng)投射到接收器位置的算子。

目標(biāo)函數(shù)對(duì)m求一階導(dǎo)數(shù)可以得到梯度矩陣，公式為：

(2)

其中，J為Frechet導(dǎo)數(shù)矩陣，上標(biāo)T和H分別表示矩陣的轉(zhuǎn)置和共軛。對(duì)目標(biāo)函數(shù)求二階導(dǎo)數(shù)可以得到近似Hessian矩陣，公式為：

Ha=Re[JTJH] 。

(3)

根據(jù)高斯牛頓法可以得到模型的更新量為：

(4)

利用梯度下降法[17]最小化目標(biāo)函數(shù)G，模型更新量可以寫成：

2 全變差正則化

如果把m表示成N1×N2的速度模型，網(wǎng)格間距為h，用(i,j)代表各個(gè)網(wǎng)格點(diǎn)坐標(biāo)，其中i=1,2,…,N1，j=1,2,…,N2，則在縱向和橫向兩個(gè)方向的梯度分別為：

(6)

全變差定義式如下：

式(7)代表離散模型中每一點(diǎn)的離散梯度的l2范數(shù)之和。假設(shè)Neumann邊界條件使這些差值在邊界處為零。我們可以通過定義一個(gè)差分算子D來更緊湊地表示‖m‖TV，這樣Dm就是離散梯度的一個(gè)連接，(Dm)n表示在n索引位置對(duì)應(yīng)的離散梯度向量，n=1,…,N,N=N1×N2，可以定義：

(8)

將全變差約束項(xiàng)應(yīng)用于常規(guī)全波形反演算法中，可以得到基于全變差約束的目標(biāo)函數(shù)，公式為：

對(duì)于式(9)，如果添加約束m∈[B1,B2]和‖m‖TV≤τ，則式(9)具有等價(jià)的優(yōu)化形式：

mn+Δm∈[B1,B2],‖m‖TV≤τ，

(10)

其中，τ是與λ有關(guān)的一個(gè)正常數(shù)，τ值越小，全變差約束項(xiàng)在目標(biāo)函數(shù)中的權(quán)重作用越大。此時(shí)反問題可以表示為：

mn+1=mn+Δm

(11a)

Δm∈[B1,B2],‖m‖TV≤τ。

(11b)

根據(jù)改進(jìn)的原始對(duì)偶混合梯度(PDHG)方法，求解式(11)的全變差約束問題相當(dāng)于求解以下問題的鞍點(diǎn)(Δm,p)：

(12)

其中，p是拉格朗日乘子，mn+Δm∈[B1,B2]，‖·‖∞,2為‖·‖1,2的對(duì)偶模，相當(dāng)于求最大值，而不是l2模的和。

修改后的PDHG方法需要迭代：

mn+Δm∈[B1,B2]

迭代公式可以更明確地寫成：

pk+1=pk+δD(mn+Δmk)-

∏‖·‖1,2≤τδ[pk+δD(mn+Δmk)]

(14)

其中，δ和α為迭代步長(zhǎng)，滿足0<αδ<1/‖DTD‖，∏‖·‖1,2≤τδ(z)表示z在半徑為τδ的球上‖·‖1,2范數(shù)的正交投影。終止條件為：

(15)

3 數(shù)值模擬

考慮一個(gè)2D合成實(shí)驗(yàn)，模擬區(qū)域網(wǎng)格點(diǎn)大小為200×200，網(wǎng)格寬度h為10 m。圖1所示的合成速度模型有一個(gè)由較慢的平滑背景包圍的恒定的高速區(qū)域。將該平滑背景作為反演的初始模型m0。類似于van Leeuwen和Herrmann[18]中的例1，將Ns=34個(gè)源置于矩形高速異常體左側(cè)，Nr=81個(gè)接收器置于矩形高速異常體右側(cè)。源qsv對(duì)應(yīng)于主頻為30 Hz的Ricker小波。

圖1 真實(shí)速度模型

正則化參數(shù)τ有3種不同的選擇：①τopt，即慢度平方的總變化量；②τlarge=1 000τopt， 此時(shí)τlarge足夠大，全變差約束沒有影響；③τsmall=0.875τopt，比理想的選擇稍微小一點(diǎn)。注意，通過使用從式(5)開始的高斯牛頓步驟作為初始模型，式(11)中的凸子問題在τlarge情況下立即收斂。對(duì)于所有實(shí)驗(yàn)，懲罰參數(shù)λ固定為1。從5 Hz數(shù)據(jù)和6 Hz數(shù)據(jù)開始，使用計(jì)算出的m作為反演6 Hz和7 Hz數(shù)據(jù)的初始模型，按照由低頻到高頻的順序，以此類推。每一個(gè)頻帶都進(jìn)行50次外部迭代，以確保凸子問題收斂，達(dá)到式(15)所示終止條件時(shí)，迭代停止。

圖2為開展的6次數(shù)值實(shí)驗(yàn)結(jié)果。在無噪聲和有噪聲的情況下，TV正則化均減少了反演模型中的振蕩現(xiàn)象，并使得高速異常體得到更好的恢復(fù)。

圖2 不同正則化參數(shù)τ在有噪聲和無噪聲情況下反演的速度模型

從噪聲數(shù)據(jù)中反演得到的不連續(xù)點(diǎn)的位置要比從無噪聲數(shù)據(jù)中反演的更準(zhǔn)確。這是因?yàn)樵肼曉诮邮掌魃显斐闪烁蟮牟贿B續(xù)，從而加強(qiáng)了其他地方TV正則化的效果。τsmall實(shí)驗(yàn)證明了這一點(diǎn)，增加TV正則化的權(quán)重會(huì)導(dǎo)致較大不連續(xù)點(diǎn)的分辨率更高。

選取國(guó)際標(biāo)準(zhǔn)復(fù)雜地質(zhì)模型開展全波形反演實(shí)驗(yàn)。真實(shí)速度模型如圖3a所示，取自墨西哥灣典型復(fù)雜巖體的合成數(shù)據(jù)——SEAM模型。模擬區(qū)域網(wǎng)格點(diǎn)大小為267×384，網(wǎng)格寬度h為20 m。在靠近地面處放置總共72個(gè)源和191個(gè)接收器。初始模型采用圖3b所示平滑模型。圖3c及圖3d為無噪聲數(shù)據(jù)下的反演結(jié)果，圖3e及圖3f為噪聲數(shù)據(jù)下的反演結(jié)果。圖3d及圖3f為選取τsmall正則化參數(shù)下的反演結(jié)果，圖3c及圖3e為選取τlarge正則化參數(shù)下的反演結(jié)果。由圖可知，全變差正則化有助于消除偽影，主要是減少鹽下的噪聲。

圖3 Seam模型的真實(shí)速度模型、初始速度模型，以及不同正則化參數(shù)τ在有噪聲和無噪聲情況下反演的速度模型

4 結(jié)論及討論

本文提出了一種在計(jì)算上可行的梯度投影算法，該算法可將附加凸約束條件下的全波形反演二次懲罰公式最小化。本文展示了在模型上添加全變差約束和邊界約束時(shí)，全波形反演中凸子問題的求解過程。數(shù)值實(shí)驗(yàn)表明，在數(shù)據(jù)有限或含有噪聲的情況下，TV正則化可以通過消除偽影和精確識(shí)別不連續(xù)點(diǎn)的位置，提高模型反演效果。

本文圍繞的梯度投影法，其本質(zhì)依賴于目標(biāo)函數(shù)的一階偏微分信息，這使得收斂速度受到一定制約。因此，理論上，我們可以采用二階方法進(jìn)行非線性優(yōu)化，以探索收斂速度是否有提升。