敖瑞德， 董良國(guó)， 遲本鑫

同濟(jì)大學(xué)海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 上海 200092

?

不依賴子波、基于包絡(luò)的FWI初始模型建立方法研究

敖瑞德， 董良國(guó)*， 遲本鑫

同濟(jì)大學(xué)海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 上海 200092

地震全波形反演(FWI)從理論走向?qū)嶋H面臨著諸多難題，其中之一就是需要一個(gè)較高精度的初始模型，另一個(gè)難題就是需要一個(gè)較為精確的震源子波，初始模型和震源子波的準(zhǔn)確程度嚴(yán)重影響著全波形反演的最終結(jié)果.為此，本文提出了不依賴子波、基于包絡(luò)的FWI初始模型建立的方法，建立了相應(yīng)的目標(biāo)函數(shù)，推導(dǎo)出了反演的梯度，給出了伴隨震源的表達(dá)式，理論上分析了不依賴子波FWI的可行性.在數(shù)值試驗(yàn)中，討論了參考道的選取方式，通過(guò)分析歸一化目標(biāo)函數(shù)收斂速率，認(rèn)為近偏移距參考道優(yōu)于遠(yuǎn)偏移距參考道，在地震數(shù)據(jù)含干擾噪音時(shí)，平均道作為參考道要優(yōu)于最小偏移距參考道.通過(guò)包絡(luò)、包絡(luò)對(duì)數(shù)、包絡(luò)平方三種目標(biāo)函數(shù)反演結(jié)果的比較，發(fā)現(xiàn)包絡(luò)對(duì)數(shù)目標(biāo)函數(shù)對(duì)深層的反演效果最好.通過(guò)不同子波的試驗(yàn)進(jìn)一步驗(yàn)證了本方法的正確性.

全波形反演； 子波； 包絡(luò)； 初始模型； 參考道； 目標(biāo)函數(shù)

1 引言

全波形反演(Full waveform inversion，簡(jiǎn)稱FWI)作為一種高分辨率地震反演成像方法，是近些年來(lái)勘探地球物理界研究的熱點(diǎn).20世紀(jì)80年代，Tarantola(1984)從波動(dòng)方程出發(fā)，利用模擬波場(chǎng)與觀測(cè)數(shù)據(jù)殘差反傳波場(chǎng)的互相關(guān)構(gòu)造梯度方向，建立了一套時(shí)間域全波形反演的理論框架.20世紀(jì)90年代，Pratt等(1998)、 Pratt(1999)又將其發(fā)展到頻率域，利用炮點(diǎn)和檢波點(diǎn)的互易性原理推導(dǎo)出梯度的高效算法.這些研究是當(dāng)前地震全波形反演的理論基礎(chǔ).但是，F(xiàn)WI目前在實(shí)際中還沒(méi)得到很好的應(yīng)用，其中一個(gè)問(wèn)題就是如何得到一個(gè)較為精準(zhǔn)的初始速度模型，另一個(gè)問(wèn)題就是如何處理震源子波.

FWI要通過(guò)地震波正演模擬來(lái)實(shí)現(xiàn)，需要知道震源子波.但在實(shí)際地震數(shù)據(jù)中，子波往往是未知的，這成為全波形反演的又一大挑戰(zhàn).最直接的解決方法就是從地震數(shù)據(jù)中提取直達(dá)波作為震源子波，但直達(dá)波并不是真實(shí)的子波.Song等(1995)則將子波的計(jì)算統(tǒng)一在速度模型的反演中，通過(guò)迭代的方式不斷更新震源子波和速度模型，但如果速度模型與真實(shí)模型有較大的偏差，計(jì)算的子波也會(huì)不準(zhǔn)確，最終造成反演的失敗.另一種思路就是在反演過(guò)程中消除子波的影響.Lee和Kim(2003)、Zhou和Greenhalgh(2003)等提出了基于聲波方程的不依賴子波的FWI，這種方法通過(guò)在頻率域?qū)⒌卣饠?shù)據(jù)與參考道波場(chǎng)相除，從而獲得歸一化的波場(chǎng)來(lái)消除震源子波的影響.該方法的缺點(diǎn)在于需要直接計(jì)算Jacobian矩陣，計(jì)算量大.Choi等(2005)在此基礎(chǔ)上提出了利用地震數(shù)據(jù)的振幅譜消除子波的兩種方法，即反褶積法與褶積法：將觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)的振幅與各自參考道的振幅先相除后作差，以此構(gòu)建目標(biāo)函數(shù)的方法稱為反褶積法；將觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)的振幅與相互的參考道的振幅先相乘后作差，以此構(gòu)建目標(biāo)函數(shù)的方法稱為褶積法.兩種方法的梯度均可通過(guò)伴隨狀態(tài)法求取，減少了計(jì)算量，但這些方法都是在頻率域?qū)崿F(xiàn)的.之后，Choi和Alkhalifah(2011)又提出了時(shí)間域不依賴子波的FWI，他們?cè)跁r(shí)間域?qū)⒛M數(shù)據(jù)和觀測(cè)數(shù)據(jù)互相與參考道波場(chǎng)褶積消除子波的影響，同時(shí)驗(yàn)證了該方法在加噪數(shù)據(jù)上的可行性，但是該方法仍需要提供較準(zhǔn)確的初始模型.

本文在前人工作基礎(chǔ)上，提出不依賴子波、基于包絡(luò)的FWI初始速度模型的建立方法，該方法克服了FWI對(duì)子波的依賴性，降低了FWI的非線性程度，在觀測(cè)地震數(shù)據(jù)缺少低頻成分情況下也可以進(jìn)行有效的FWI反演，從而為下一步FWI提供較為準(zhǔn)確的初始模型.

2 方法原理

2.1 目標(biāo)函數(shù)

全波形反演的目標(biāo)函數(shù)一般定義為模擬數(shù)據(jù)與觀測(cè)數(shù)據(jù)殘差的L2泛函，形式如下：

(1)

其中s和r分別是炮號(hào)和道號(hào)，ui,j和di,j分別是第i炮的第j個(gè)檢波器處的模擬數(shù)據(jù)和觀測(cè)數(shù)據(jù).將ui,j和di,j寫(xiě)成格林函數(shù)和子波的褶積形式，目標(biāo)函數(shù)變?yōu)?/p>

(2)

其中，g是格林函數(shù)，w是子波，*代表褶積，上標(biāo)u和d分別代表模擬數(shù)據(jù)和觀測(cè)數(shù)據(jù).

為了避免在全波形反演中估算子波，可以通過(guò)交叉褶積波場(chǎng)的方法消除子波的影響(ChoiandAlkhalifah,2011).對(duì)于每一個(gè)炮集，模擬數(shù)據(jù)與同一炮觀測(cè)數(shù)據(jù)的第k個(gè)參考道進(jìn)行褶積，獲得模擬褶積波場(chǎng)；觀測(cè)數(shù)據(jù)與同一炮模擬數(shù)據(jù)的第k個(gè)參考道進(jìn)行褶積，獲得觀測(cè)褶積波場(chǎng).此時(shí)，新的目標(biāo)函數(shù)定義為

(3)

其中ui,k和di,k分別是第i炮的第k個(gè)參考道的模擬數(shù)據(jù)和觀測(cè)數(shù)據(jù).ui,j*di,k為模擬褶積波場(chǎng)，di,j*ui,k為觀測(cè)褶積波場(chǎng).

將(3)式中的模擬數(shù)據(jù)與觀測(cè)數(shù)據(jù)都表示為格林函數(shù)與震源子波的褶積形式：

(4)

從(4)式可以發(fā)現(xiàn)，其中的減號(hào)左右兩邊都與共同的新震源子波(wu*wd)進(jìn)行了褶積，這樣在目標(biāo)函數(shù)中就消除了實(shí)際子波與模擬子波差異的影響，而由(2)式可知，目標(biāo)函數(shù)(1)卻與子波差異有關(guān).因此，基于目標(biāo)函數(shù)(3)的全波形反演結(jié)果不再依賴于準(zhǔn)確的子波(圖2b)，而目標(biāo)函數(shù)(1)依賴于子波，采用不準(zhǔn)確的子波將嚴(yán)重影響反演結(jié)果的精度(圖2a).

進(jìn)一步地，利用褶積波場(chǎng)的包絡(luò)提取低頻信息：

(5)

(6)

(7)

除了(7)式中包絡(luò)誤差形式的目標(biāo)函數(shù)，常用的目標(biāo)函數(shù)還有下面的包絡(luò)平方目標(biāo)函數(shù)(Wuetal.,2014)和包絡(luò)對(duì)數(shù)目標(biāo)函數(shù)(Chietal.,2014)，它們都不受子波的影響：

(8)

(9)

2.2 梯度和伴隨震源

在目前基于局部?jī)?yōu)化的全波形反演過(guò)程中，目標(biāo)函數(shù)梯度的求取是關(guān)鍵.(10)式是包絡(luò)對(duì)數(shù)目標(biāo)函數(shù)(9)對(duì)模型參數(shù)p的梯度，本文中p代表縱波速度(詳細(xì)的推導(dǎo)過(guò)程請(qǐng)參考附錄A)，

(10)

(11)

其中?代表互相關(guān)運(yùn)算.將伴隨震源設(shè)為

當(dāng)取參考道di,k=ui,k=δ(t)時(shí)，褶積波場(chǎng)即為原始波場(chǎng)，參考道與空間位置無(wú)關(guān)，此時(shí)伴隨震源只剩r32.應(yīng)用δ函數(shù)的特殊性質(zhì)，伴隨震源進(jìn)一步退化為

當(dāng)采用其他兩種目標(biāo)函數(shù)時(shí)，梯度的推導(dǎo)過(guò)程與上面類似.

對(duì)包絡(luò)目標(biāo)函數(shù)(7)，在參考道k處的伴隨震源為

(15)

在任意檢波點(diǎn)j處的伴隨震源為

(16)

對(duì)包絡(luò)平方目標(biāo)函數(shù)(8)，在參考道k處的伴隨震源為

(17)

在任意檢波點(diǎn)j處的伴隨震源為：

(18)

2.3 初始模型建立過(guò)程

利用本文提出的方法進(jìn)行初始模型建立的步驟與時(shí)間域常規(guī)FWI反演的步驟基本一致，主要不同在于目標(biāo)函數(shù)和伴隨震源的計(jì)算，已在2.1和2.2節(jié)中分別給出.本文中的具體過(guò)程如下：(1)給定一個(gè)模型，正演計(jì)算每一炮的模擬波場(chǎng)；(2)選定合適的參考道，求取模擬褶積波場(chǎng)和觀測(cè)褶積波場(chǎng)；(3)計(jì)算目標(biāo)函數(shù)和伴隨震源；(4)利用伴隨震源計(jì)算伴隨波場(chǎng)，并求取單炮梯度；(5)累加各單炮的梯度，從而獲得全局梯度；(6)通過(guò)L-BFGS算法對(duì)梯度進(jìn)行校正并求取步長(zhǎng)，進(jìn)行模型更新，直到滿足終止條件，最終輸出反演結(jié)果.

3 數(shù)值試驗(yàn)

采用的Marmousi速度模型(圖1a)參數(shù)如下：水平網(wǎng)格數(shù)nx=501，垂直網(wǎng)格數(shù)nz=140，水平網(wǎng)格大小dx=15 m，垂直網(wǎng)格大小dz=15 m，水層深30 m，水平總長(zhǎng)度7500 m，垂直總長(zhǎng)度2085 m.震源和檢波器都位于地表處，共有80炮數(shù)據(jù)，第一炮位于195 m處，炮間距為90 m.第一個(gè)檢波器位于0 m處，檢波器間隔15 m.地震記錄的采樣間隔為0.0015 s，采樣點(diǎn)數(shù)為2000，總記錄時(shí)間為3 s.用高階有限差分法(董良國(guó)等，2000)進(jìn)行二維聲波的正演.理論子波是主頻為8 Hz的Ricker子波，延時(shí)0.15 s，子波最大振幅為1.0.

即使采用相對(duì)精確的Gauss平滑初始模型(圖1b)，利用錯(cuò)誤的子波(主頻3 Hz的Ricker子波，延時(shí)0.5 s，最大振幅2.0)，常規(guī)FWI也會(huì)得到很差的反演結(jié)果(圖2a).而不依賴子波的FWI的反演結(jié)果對(duì)模型的高頻成份進(jìn)行了很好的恢復(fù)(圖2b).

如果以梯度模型作為FWI的初始模型，常規(guī)FWI將不僅面對(duì)錯(cuò)誤子波的問(wèn)題，較差的初始模型還會(huì)使反演極易陷入局部極值，常規(guī)FWI的結(jié)果將變得更加糟糕.本文提出的不依賴子波、基于包絡(luò)的FWI方法，就是要在不知道震源子波以及沒(méi)有任何模型先驗(yàn)信息的條件下，為下一步的常規(guī)FWI反演提供一個(gè)含有可靠的低波數(shù)成分的初始模型.其中參考道的選擇、方法的抗噪性、不同目標(biāo)函數(shù)的比較以及不同子波的對(duì)比，是本文提出的方法中涉及的幾個(gè)關(guān)鍵問(wèn)題，下面通過(guò)數(shù)值試驗(yàn)和理論分析對(duì)這幾個(gè)問(wèn)題進(jìn)行探討.為了表達(dá)方便，將依賴于正確子波、基于包絡(luò)的FWI稱為傳統(tǒng)的基于包絡(luò)的FWI，在下面的試驗(yàn)中，它與本文提出的不依賴子波、基于包絡(luò)的FWI均是在梯度模型(圖3)的基礎(chǔ)上進(jìn)行的.

3.1 參考道的選擇

本文提出的不依賴子波、基于包絡(luò)的FWI方法中的一個(gè)關(guān)鍵問(wèn)題，就是參考道應(yīng)該如何選取.理論上，參考道的選擇并不重要，然而實(shí)際中，參考道的選擇對(duì)于反演結(jié)果有較大影響.關(guān)于參考道的選取，前人已經(jīng)做了比較詳細(xì)的分析.劉寶龍等(2014)比較了每個(gè)炮集中不同參考道的全波形反演結(jié)果，認(rèn)為最小偏移距參考道為最優(yōu)選擇方案.但這種選取方案的FWI并沒(méi)有涉及包絡(luò)目標(biāo)函數(shù)的問(wèn)題.

為了比較不同參考道對(duì)本方法反演結(jié)果的影響，我們進(jìn)行了兩組試驗(yàn)，分別利用不含噪數(shù)據(jù)和含有噪聲的數(shù)據(jù)進(jìn)行測(cè)試和對(duì)比.選取的子波都為主頻3 Hz的Ricker子波，延時(shí)0.5 s，最大振幅2.0，顯然與觀測(cè)數(shù)據(jù)的子波(8 Hz的Ricker子波，延時(shí)0.15 s，子波最大振幅為1.0)差別很大.

首先比較數(shù)據(jù)不含噪聲時(shí)不同偏移距的參考道.理論上，最小偏移距參考道的能量最大；考慮實(shí)際波場(chǎng)的傳播，子波具有時(shí)變和空變特征，偏移距越小，畸變?cè)叫?，因此?dāng)采用不依賴子波、基于包絡(luò)的FWI建立初始模型時(shí)，最小偏移距參考道應(yīng)優(yōu)于遠(yuǎn)偏移距參考道.但直接比較不同偏移距的反演結(jié)果難以看出差別(圖4)，而從圖5的歸一化的目標(biāo)函數(shù)隨迭代次數(shù)的變化過(guò)程可以看出，隨著參考道的偏移距增加，收斂速度減慢.這個(gè)試驗(yàn)說(shuō)明：近偏移距參考道優(yōu)于遠(yuǎn)偏移距參考道.因此在下面的試驗(yàn)中選取零偏移距參考道.

圖1 Marmousi準(zhǔn)確模型(a)以及高斯平滑結(jié)果(b)Fig.1 The Marmousi velocity model (a) and Gauss-window smoothed model (b)

圖2 子波錯(cuò)誤時(shí)不同目標(biāo)函數(shù)的FWI反演結(jié)果(a) 目標(biāo)函數(shù)(1)的常規(guī)FWI反演結(jié)果；(b) 目標(biāo)函數(shù)(3)的不依賴子波的FWI反演結(jié)果.Fig.2 FWI results of different misfit functions with an incorrect wavelet(a) The traditional FWI result corresponding to misfit function (1); (b) Source-independent FWI result corresponding to misfit function (3).

圖4 無(wú)噪時(shí)選取不同偏移距參考道迭代16次的反演結(jié)果(a) 偏移距750 m； (b) 偏移距600 m； (b) 偏移距0 m.Fig.4 The inversion results of different offset reference traces after 16 iterations for noise-free data(a) 750 m-offset； (b) 600 m-offset； (c) 0 m-offset.

圖3 梯度初始模型Fig.3 The gradient initial model

圖5 無(wú)噪時(shí)不同偏移距參考道反演過(guò)程中目標(biāo)函數(shù)收斂速度對(duì)比紅色：偏移距750 m參考道；綠色：偏移距600 m參考道；藍(lán)色：0 m偏移距參考道.Fig.5 The history of normalized misfit functions of 750 m-offset reference trace(red), 600 m-offset reference trace (green) and 0 m-offset reference trace(blue) for noise-free data.

圖6 加入噪聲的一炮地震記錄(信噪比10 dB)Fig.6 One shot record contaminated by noise with 10 dB S/N ratio

對(duì)于含有噪聲的數(shù)據(jù)，如果只選取某一道作為參考道，噪聲的影響會(huì)很大.而取平均道作參考道后，不同地震道的噪聲可以相互抵消，從而可以更好地保留信號(hào)中的有效信息.所以在第二組試驗(yàn)中，對(duì)原始地震記錄加入信噪比10 dB的噪聲(圖6)，其中的噪聲是產(chǎn)生的隨機(jī)噪音與正演模擬的子波褶積的結(jié)果，以保證噪音與有效信號(hào)頻帶一致.選取每個(gè)炮集炮點(diǎn)附近100道的平均道作為參考道.對(duì)比在無(wú)噪數(shù)據(jù)和加噪數(shù)據(jù)情況下，零偏移距參考道和平均道作參考道時(shí)不依賴子波、基于包絡(luò)的FWI初始模型建立結(jié)果(圖7)，可以發(fā)現(xiàn)，本方法在初始模型建立過(guò)程中對(duì)噪聲并不敏感，均很好地反演出了模型的背景成份.從歸一化的目標(biāo)函數(shù)(圖8)中也可以看出，當(dāng)數(shù)據(jù)中含有噪聲后，平均道作參考道的收斂速度要優(yōu)于零偏移距參考道的收斂速度.

由上可知，對(duì)于無(wú)噪數(shù)據(jù)，近偏移距參考道的目標(biāo)函數(shù)收斂情況優(yōu)于遠(yuǎn)偏移距參考道；對(duì)于加噪數(shù)據(jù)，平均道作參考道的目標(biāo)函數(shù)收斂情況優(yōu)于最小偏移距參考道.利用不依賴子波、基于包絡(luò)的FWI方法，建立的初始模型效果對(duì)于不同參考道和數(shù)據(jù)噪聲并不敏感，這是因?yàn)槔冒j(luò)目標(biāo)函數(shù)的反演方法僅對(duì)模型的低頻成份敏感.最終的反演結(jié)果得到的是宏觀的背景速度，可以為接下來(lái)的常規(guī)FWI提供一個(gè)較好的初始速度模型.

3.2 不同目標(biāo)函數(shù)的比較

由2.1節(jié)可知，不依賴子波、基于包絡(luò)的FWI方法可以選擇三種不同的目標(biāo)函數(shù).雖然三種目標(biāo)函數(shù)均可以反演背景速度，但結(jié)果差別較大.圖9是選用同樣的錯(cuò)誤子波、在梯度模型的基礎(chǔ)上、采用不同目標(biāo)函數(shù)的反演結(jié)果，它們均是不依賴子波的，都對(duì)模型的背景成份進(jìn)行了一定的恢復(fù).但是，明顯可以看出，包絡(luò)對(duì)數(shù)目標(biāo)函數(shù)對(duì)深層的恢復(fù)更好(圖9a)，其次為包絡(luò)目標(biāo)函數(shù)(圖9b)，而包絡(luò)平方目標(biāo)函數(shù)的效果最差(圖9c).在同樣的迭代次數(shù)下，不依賴子波、基于包絡(luò)平方目標(biāo)函數(shù)的FWI方法對(duì)模型的深部速度基本不更新，收斂速度很慢.

最后，選取最小偏移距參考道和包絡(luò)對(duì)數(shù)目標(biāo)函數(shù)進(jìn)行不依賴子波、基于包絡(luò)的FWI建立初始速度模型，在此基礎(chǔ)上進(jìn)行常規(guī)的FWI，并與傳統(tǒng)的基于包絡(luò)的FWI反演結(jié)果對(duì)比(圖12).可以發(fā)現(xiàn)，錯(cuò)誤的子波對(duì)于傳統(tǒng)的基于包絡(luò)的FWI反演影響是巨大的(圖12a)，而以圖12a為初始模型的常規(guī)FWI結(jié)果當(dāng)然非常糟糕(圖12b).而當(dāng)采用了不依賴子波、基于包絡(luò)的FWI方法后，反演的初始模型不再受錯(cuò)誤子波的影響，較好地反映了模型速度變化的低波數(shù)成分(圖12c)，而以圖12c為初始模型的常規(guī)FWI反演結(jié)果(圖12d)與準(zhǔn)確速度模型(圖1a)非常接近，這一點(diǎn)也可以從沿深度抽線結(jié)果中清楚地看出(圖13).

圖7 不同參考道迭代16次反演結(jié)果(a) 無(wú)噪時(shí)零偏移距參考道反演結(jié)果； (b) 加噪時(shí)零偏移距參考道反演結(jié)果；(c) 無(wú)噪時(shí)平均道作參考道反演結(jié)果； (d) 加噪時(shí)平均道作參考道反演結(jié)果.Fig.7 The inversion results with different reference traces after 16 iterations(a) 0 m-offset reference trace for noise-free data; (b) 0 m-offset reference trace for noisy data; (c) Average trace as reference trace for noise-free data; (d) Average trace as reference trace for noisy data.

圖8 不同參考道反演時(shí)，歸一化目標(biāo)函數(shù)收斂曲線對(duì)比紅色：平均道作參考道；藍(lán)色：零偏移距參考道. (a) 無(wú)噪時(shí)歸一化目標(biāo)函數(shù)收斂情況；(b)加噪時(shí)歸一化目標(biāo)函數(shù)收斂情況.Fig.8 The history of normalized misfit functions for different reference traces(a) Noise-free data; (b) Noisy data. Red: average trace as reference trace; Blue: 0 m-offset reference trace.

圖9 不依賴子波、基于不同包絡(luò)目標(biāo)函數(shù)的FWI反演結(jié)果對(duì)比(a) 包絡(luò)對(duì)數(shù)目標(biāo)函數(shù)； (b) 包絡(luò)目標(biāo)函數(shù)； (c) 包絡(luò)平方目標(biāo)函數(shù).Fig.9 The inversion results of source-independent envelope-based FWI with different misfit functions(a) Logarithmic envelope misfit function; (b) Envelope misfit function; (c) Squared envelope misfit function.

圖10 第一次迭代后，不同目標(biāo)函數(shù)反演的波場(chǎng)殘差(a) 包絡(luò)對(duì)數(shù)目標(biāo)函數(shù)； (b) 包絡(luò)目標(biāo)函數(shù)； (c) 包絡(luò)平方目標(biāo)函數(shù).Fig.10 The residual wavefields for different misfit functions after 1st iteration(a) Logarithmic envelope misfit function； (b) Envelope misfit function; (c) Squared envelope misfit function.

圖11 第一次迭代后，不同目標(biāo)函數(shù)反演的伴隨震源(a) 包絡(luò)對(duì)數(shù)目標(biāo)函數(shù)；(b) 包絡(luò)目標(biāo)函數(shù)；(c) 包絡(luò)平方目標(biāo)函數(shù).Fig.11 The adjoint sources for different misfit functions after 1st iteration(a) Logarithmic envelope misfit function; (b) Envelope misfit function; (c) Squared envelope misfit function.

3.3 子波試驗(yàn)

為了進(jìn)一步說(shuō)明本文提出的方法在為常規(guī)FWI建立初始模型時(shí)并不依賴于子波，我們選取了不同參數(shù)的子波(均不同于原始觀測(cè)數(shù)據(jù)中的子波)來(lái)進(jìn)行測(cè)試，選用的子波參數(shù)如下：(1)主頻20 Hz、濾去5 Hz以下低頻成分的Ricker子波，延時(shí)0.5 s，最大振幅2.0；(2)主頻為12 Hz的高斯子波，延時(shí)0.15 s，最大振幅1.0.這里需要再次說(shuō)明的是，合成地震數(shù)據(jù)時(shí)所用的是主頻為8 Hz、延時(shí)0.15 s、最大振幅為1.0的Ricker子波.

圖14a和圖14c是本文方法采用上述兩種錯(cuò)誤子波時(shí)的反演結(jié)果，可以看出，即使采用錯(cuò)誤的子波，甚至在地震數(shù)據(jù)缺少低頻信息的情況下，反演結(jié)果也都較好地恢復(fù)了模型中的背景成分.用圖14a和圖14c作為初始速度模型，常規(guī)FWI的反演結(jié)果都比較好地刻畫(huà)了真實(shí)模型的高低波數(shù)成分(圖14b和圖14d).

為了驗(yàn)證本文方法反演的背景速度的可靠性，我們分別用四種速度模型對(duì)地震數(shù)據(jù)進(jìn)行了逆時(shí)偏移，結(jié)果如圖15所示.偏移所用的四種速度模型分別是：準(zhǔn)確模型的Gauss平滑結(jié)果(圖1b)、上述兩種錯(cuò)誤子波情況下本文方法的反演結(jié)果(圖14a和圖14c)以及梯度模型(圖3).相比于梯度模型的偏移結(jié)果(圖15d)，本文提出的方法反演的模型的成像結(jié)果(圖15b、c)大大改善，與Gauss平滑模型的偏移結(jié)果(圖15a)基本一致，這說(shuō)明了本方法在子波未知情況下，仍對(duì)模型中的低波數(shù)成分進(jìn)行了很好的恢復(fù)，這樣的反演結(jié)果當(dāng)然就為接下來(lái)的常規(guī)FWI提供了很好的初始模型.

圖12 子波錯(cuò)誤時(shí)FWI反演結(jié)果對(duì)比(a) 傳統(tǒng)基于包絡(luò)的FWI建立的初始模型；(b) 以(a)為初始模型的常規(guī)FWI反演結(jié)果；(c)不依賴子波、基于包絡(luò)的FWI建立的初始模型；(d) 以(c)為初始模型的常規(guī)FWI反演結(jié)果.Fig.12 The comparison of different FWI results with an incorrect wavelet(a) The initial model built by traditional envelope-based FWI; (b) The inversion result of traditional FWI using (a) as the initial model； (c) The initial model built by source-independent envelope-based FWI; (d) The inversion result of traditional FWI using (c) as the initial model.

圖13 本文方法的反演抽線結(jié)果(a) 1000 m處速度抽線； (b) 5800 m處速度抽線. 綠色：準(zhǔn)確模型；藍(lán)色：本文方法建立的初始模型；紅色：最終FWI反演結(jié)果.Fig.13 Depth velocity slices of inversion results(a) X=1000 m; (b) X=5800 m. Green: true model; Blue: the inverted initial model by our method; Red: the final inversion result.

圖14 本文方法采用兩個(gè)錯(cuò)誤子波時(shí)建立的初始模型以及相應(yīng)的常規(guī)FWI反演結(jié)果(a) 采用主頻為20 Hz濾去5 Hz以下的Ricker子波建立的初始模型；(b)以(a)為初始模型的常規(guī)FWI反演結(jié)果；(c) 采用主頻為12 Hz的高斯子波建立的初始模型；(d) 以(c)為初始模型的常規(guī)FWI反演結(jié)果.Fig.14 The inverted initial models using two incorrect wavelets by our method and the corresponding traditional FWI results(a) The inverted initial model using 20 Hz Ricker wavelet without low frequency below 5 Hz; (b) The corresponding traditional FWI result using (a) as the initial model; (c) The inverted initial model using 12 Hz Gauss wavelet; (d) The corresponding traditional FWI result using (c) as the initial model.

圖15 不同速度模型的逆時(shí)偏移結(jié)果對(duì)比(a) 高斯平滑初始模型；(b) 主頻為20 Hz濾去5 Hz以下低頻Ricker子波建立的初始模型；(c) 12 Hz高斯子波初始模型；(d) 梯度初始模型.Fig.15 The comparison of RTM results with different velocities(a) Gauss-window smoothed initial model; (b) The inverted initial model using 20 Hz Ricker wavelet without low frequency below 5 Hz; (c) The inverted initial model using 12 Hz Gauss wavelet; (d) Gradient initial model.

4 結(jié)論

(1) 為進(jìn)一步推動(dòng)FWI的實(shí)用性，本文在前人研究的基礎(chǔ)上，提出了不依賴子波、基于包絡(luò)的FWI反演方法.在子波未知、沒(méi)有任何模型先驗(yàn)信息、以及地震數(shù)據(jù)缺少低頻信息的情況下，這種方法可以為常規(guī)FWI提供一個(gè)可靠的初始速度模型. ，也可以作為疊前深度偏移的速度模型.

(2) 比較了不同參考道的反演結(jié)果和歸一化目標(biāo)函數(shù)的收斂速度，認(rèn)為近偏移距參考道優(yōu)于遠(yuǎn)偏移距參考道.當(dāng)數(shù)據(jù)含有噪音時(shí)，平均道作為參考道優(yōu)于最小偏移距參考道.

(3) 在不依賴子波、基于包絡(luò)的FWI的不同目標(biāo)函數(shù)中，包絡(luò)對(duì)數(shù)目標(biāo)函數(shù)為最優(yōu)選擇，對(duì)深層的反演效果最好.

(4) 不依賴子波、基于包絡(luò)的FWI反演方法能夠克服傳統(tǒng)FWI方法對(duì)震源子波的依賴.

附錄A 不依賴子波、基于包絡(luò)對(duì)數(shù)目標(biāo)函數(shù)的梯度公式推導(dǎo)

不依賴子波、基于包絡(luò)對(duì)數(shù)目標(biāo)函數(shù)為公式(9)，它對(duì)模型參數(shù)p求導(dǎo)，得

(A1)

其中

(A2)

(A3)

h為希爾伯特變換的濾波因子.

將(A2)和(A3)帶入方程(A1)中并展開(kāi)，得到如下四項(xiàng)：

(A4)

(A5)

(A6)

(A7)

其中第一項(xiàng)(A4)和第三項(xiàng)(A6)，第二項(xiàng)(A5)和第四項(xiàng)(A7)推導(dǎo)過(guò)程類似.

下面推導(dǎo)(A4). (A4)可變?yōu)?/p>

(A8)

令

(A9)

(A8)進(jìn)一步變換為

(A10)

(A11)

令

(A12)

(A11)進(jìn)一步變換為

(A13)

作類似的變換，第二項(xiàng)(A5)化為

(A14)

對(duì)(A14)括號(hào)內(nèi)的部分應(yīng)用相關(guān)褶積定理，同時(shí)注意濾波因子h為奇函數(shù)，得

對(duì)于包絡(luò)目標(biāo)函數(shù)、包絡(luò)平方目標(biāo)函數(shù)，它們對(duì)參數(shù)的梯度的推導(dǎo)過(guò)程與此類似，這里不再贅述.

Brossier R, Operto S, Virieux J. 2009. Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion.Geophysics, 74(6): WCC105-WCC118. Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion.Geophysics, 60(5): 1457-1473.

Chi B X, Dong L G, Liu Y Z. 2013. Full waveform inversion based on envelope objective function. ∥ 75th Annual International Conference and Exhibition, EAGE, Extended Abstracts.

Chi B X, Dong L G, Liu Y Z. 2014.Full waveform inversion method using envelope objective function without low frequency data.JournalofAppliedGeophysics, 109: 36-46.

Choi Y, Shin C, Min D J, et al. 2005. Efficient calculation of the steepest descent direction for source-independent seismic waveform inversion: An amplitude approach.JournalofComputationalPhysics, 208(2): 455-468. Choi Y, Alkhalifah T. 2011. Source-independent time-domain waveform inversion using convolved wavefields: Application to the encoded multisource waveform inversion.Geophysics, 76(5): R125-R134.

Dong L G, Ma Z T, Cao J Z, et al. 2000. A staggered-grid high-order difference method of one-order elasticwave equation.ChineseJournalofGeophysics(in Chinese), 43(3): 411-419, doi: 10.3321/j.issn:0001-5733.2000.03.015.

Lee K H, Kim H J. 2003. Source-independent full-waveform inversion of seismic data.Geophysics, 68(6): 2010-2015.

Liu B L, Dong L G, Chi B X. 2014. Method of Source-independent Full Waveform Inversion[Ph. D. thesis] (in Chinese). Shanghai: Tongji University. Pratt R G, Shin C, Hick G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.GeophysicalJournalInternational, 133(2): 341-362.Pratt R G. 1999. Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model.Geophysics, 64(3): 888-901.

Sirgue L, Pratt R G. 2004. Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies.Geophysics, 69(1): 231-248.Song Z M, Williamson P R, Pratt R G. 1995. Frequency-domain acoustic-wave modeling and inversion of crosshole data: Part II-Inversion method, synthetic experiments and real-data results.Geophysics, 60(3): 796-809.

Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation.Geophysics, 49(8): 1259-1266.

Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics.Geophysics, 74(6): WCC1-WCC26.

Wu R S, Luo J R, Wu B Y. 2014. Seismic envelope inversion and modulation signal model.Geophysics, 79(3): WA13-WA24.

Zhou B, Greenhalgh S A. 2003. Crosshole seismic inversion with normalized full-waveform amplitude data.Geophysics, 68(4): 1320-1330.

附中文參考文獻(xiàn)

董良國(guó), 馬在田, 曹景忠等. 2000. 一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法. 地球物理學(xué)報(bào), 43(3): 411-419, doi: 10.3321/j.issn:0001-5733.2000.03.015.

劉寶龍, 董良國(guó), 遲本鑫. 2014. 不依賴于子波的全波形反演方法研究[碩士論文]. 上海: 同濟(jì)大學(xué).

(本文編輯 胡素芳)

Source-independent envelope-based FWI to build an initial model

AO Rui-De, DONG Liang-Guo*, CHI Ben-Xin

StateKeyLaboratoryofMarineGeology,TongjiUniversity,Shanghai200092,China

Full waveform inversion (FWI) is facing many difficulties in practice. One of the issues is how to construct a relatively precise initial model. Another one is the treatment of source wavelet. Both seriously affect the final inversion result. To tackle this problem, we propose a source-independent envelope-based FWI approach to build an initial model for the traditional FWI.On one hand, the cross-convolved wavefields, including the convolution of observed wavefields with a reference trace from the modeled wavefields and the convolution of modeled wavefields with a reference trace from the observed wavefields, can eliminate source wavelet effect in the misfit function of waveform inversion. On the other hand, the envelope of wavefields contains abundant low-frequency components, which are crucial for initial model inversion. Therefore, we use the envelope of cross-convolved wavefields to establish misfit functions, the feasibility of which is proved by theoretical analysis. The expressions of the corresponding gradient and adjoint source are derived. Furthermore, we use a Marmousi synthetic data example to test our method.The proposed method is tested by using different reference traces, misfit functions and source wavelets. (1) The normalized misfit of long offset reference trace converges more slowly than that of short-offset reference trace for noise-free data and the average trace as reference trace is better than the minimum-offset reference trace for noisy data. (2) The logarithmic envelope is the best and the squared envelope is the worst misfit function to invert the deep part of the model. The initial model built by our method choosing the logarithmic misfit function and minimum-offset reference trace is much better than that of the traditional envelope-based FWI when the source wavelet is incorrect. (3) Using other kinds of source wavelets, such as 20 Hz Ricker wavelet without low frequency below 5 Hz and 12 Hz Gauss wavelet, both of which are quite different from the true 8 Hz Ricker wavelet, the inversion result of our method can provide a good initial model for the traditional FWI. In addition, the initial model can be used in seismic depth migration, and the imaging result is close to the RTM result with true velocity model and much better than that with gradient initial model.In conclusion, first, we propose a source-independent envelope-based FWI approach to build an initial model for the traditional FWI without any prior information when the source wavelet is unknown and even the data lack low frequency. The initial model built by our method can also be used as the velocity model for pre-stack depth migration. Second, the short-offset reference trace is better than long-offset reference trace for noise-free data and average trace as reference trace is better than the minimum-offset reference trace for noisy data. Third, the logarithmic envelope misfit function is the optimal misfit function. Finally, by changing the source wavelet, it is further proved that our method is source-independent.

Full waveform inversion; Source wavelet; Envelope; Initial model; Reference trace; Misfit function

10.6038/cjg20150615.

國(guó)家油氣重大專項(xiàng)(2011ZX05005-005-007HZ)資助.

敖瑞德，1991年生，男，2014年畢業(yè)于同濟(jì)大學(xué)，理學(xué)學(xué)士，主要從事地震全波形反演的相關(guān)研究，現(xiàn)為清華大學(xué)碩士研究生. E-mail:ruideao@163.com

*通訊作者 董良國(guó)，同濟(jì)大學(xué)教授，主要從事地震波傳播與數(shù)值模擬、地震波波形反演與層析成像等方面的研究.E-mail:dlg@#edu.cn

10.6038/cjg20150615

P631

2014-12-05，2015-05-25收修定稿

敖瑞德， 董良國(guó)， 遲本鑫.2015.不依賴子波、基于包絡(luò)的FWI初始模型建立方法研究.地球物理學(xué)報(bào)，58(6):1998-2010,

Ao R D, Dong L G, Chi M X. 2015. Source-independent envelope-based FWI to build an initial model.ChineseJ.Geophys. (in Chinese),58(6):1998-2010,doi:10.6038/cjg20150615.