楊 濤 張會(huì)星* 史才旺

(①中國(guó)海洋大學(xué)，山東青島 266100； ②青島海洋科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室，山東青島 266071； ③海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東青島 266100)

0 引言

速度是地震波在地下傳播過(guò)程中的重要物性參數(shù)之一，速度建模是高精度地震成像和資料解釋的重要基礎(chǔ)。全波形反演綜合利用了地震波場(chǎng)中的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)信息，可進(jìn)行高精度、多參數(shù)建模，滿足復(fù)雜構(gòu)造的精度要求，是勘探地球物理領(lǐng)域的研究熱點(diǎn)[1-4]。

全波形反演方法具有完備的理論基礎(chǔ)[5-7]，該方法基于最小二乘最佳擬合原理，利用正傳波場(chǎng)與伴隨震源的反傳波場(chǎng)的互相關(guān)構(gòu)造梯度，以此作為速度模型的更新方向。由于反演是一個(gè)非線性過(guò)程，常規(guī)全波形反演容易陷入局部極小值[8]，為解決這個(gè)問(wèn)題，Bunks等[9]提出了時(shí)間域多尺度分頻處理的理念，通過(guò)對(duì)資料做分頻處理，從低頻信息開(kāi)始反演，把反演更新的速度模型作為初始模型逐步向高頻反演。Pica等[10]和Pratt[11-12]提出了頻率域全波形反演的思想，相比時(shí)間域反演，其計(jì)算速度更快。Brossier等[13]在頻率域?qū)崿F(xiàn)了彈性波全波形反演，隨著模型的增大，頻率域正演需要海量?jī)?nèi)存[14]，在很大程度上限制了頻率域反演在實(shí)際三維地震勘探中的應(yīng)用。結(jié)合時(shí)間域正演和頻率域反演的優(yōu)勢(shì)，Sirgue等[15]提出了時(shí)間域正演及頻率域反演的混合反演策略，避免了時(shí)間域波場(chǎng)反傳的計(jì)算，提高了計(jì)算效率。

對(duì)于實(shí)際資料，全波形反演的建模精度受震源子波的影響。在野外采集環(huán)節(jié)很難獲得準(zhǔn)確的震源子波信息，且反演過(guò)程中震源子波的微小誤差即可能造成反演結(jié)果和實(shí)際模型不匹配。目前子波求取方法大致包括三類。第一類方法主要基于地震反褶積模型，分為確定性和統(tǒng)計(jì)性兩種方法[16-17]。提取確定性子波需要利用測(cè)井資料求解反射系數(shù)，因此在沒(méi)有測(cè)井資料時(shí)該方法很難實(shí)現(xiàn)；統(tǒng)計(jì)性子波方法僅使用實(shí)測(cè)地震記錄就可以提取地震子波，但如果沒(méi)有明確的地層反射系數(shù)信息，精度很難保證[18]。第二類方法為反演子波法，Song等[19]提出在反演過(guò)程中統(tǒng)一計(jì)算子波，子波隨模型的更新而更新。此種方法對(duì)初始模型的依賴性大，在沒(méi)有精確初始模型的情況下很難取得理想效果。胡勇等[20]利用直達(dá)波信息構(gòu)造反演目標(biāo)震源函數(shù)，用全波形反演方法反演子波，取得了較好效果。第三類方法為不依賴子波法，即在反演過(guò)程中去除子波的影響。Choi等[21]在頻率域提出利用褶積法和反褶積法消除子波對(duì)反演結(jié)果的影響，將觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)的振幅與各自參考道的振幅先相除后計(jì)算差值，構(gòu)建目標(biāo)函數(shù)，此方法稱為反褶積法； 之后Choi等[22]提出了時(shí)間域去除子波影響的全波形反演，其實(shí)質(zhì)是把波場(chǎng)和特征道在時(shí)間域做卷積處理，取得了較好的效果。敖瑞德等[23]把該方法應(yīng)用到基于包絡(luò)的全波形反演。王毓偉等[24]把地震道包絡(luò)的全波形反演推廣到彈性波，并結(jié)合時(shí)間域正演、頻率域反演實(shí)現(xiàn)了對(duì)彈性參數(shù)的高精度重建。

以上去除子波影響的方法都是基于聲波方程實(shí)現(xiàn)的，為了更接近地震波在地下介質(zhì)傳播的真實(shí)情況，滿足全波形反演在內(nèi)存和計(jì)算效率上的需求，本文提出不依賴子波的彈性波混合域全波形反演。文中給出了不依賴子波的彈性波混合域全波形的目標(biāo)函數(shù)，并推導(dǎo)了反傳震源的公式。相比于常規(guī)的全波形反演，在震源子波不準(zhǔn)確的情況下，該方法依然能取得較好的反演結(jié)果。

1 方法理論

1.1 常規(guī)彈性波混合域全波形反演

正演模擬采用二階彈性波應(yīng)力—位移方程

(1)

式中：ux、uz分別表示x、z分量的位移；τxx、τzz表示正應(yīng)力；τxz表示剪應(yīng)力；ρ為密度；λ、μ為拉梅常數(shù)。

全波形反演的目標(biāo)函數(shù)有多種求法，本文使用最常用的L2范數(shù)下的目標(biāo)函數(shù)

(2)

式中：v表示模型參數(shù)；ui表示第i炮頻率域模擬記錄；di表示第i炮頻率域觀測(cè)記錄； 上標(biāo)“*”表示復(fù)共軛。頻率域梯度可表示為

(3)

式中： Re(·)表示取實(shí)部；ns為總炮數(shù)。頻率域波動(dòng)方程可以表示為S(ω)u(ω)=f(ω)，這里S是阻抗矩陣，u為頻率域波場(chǎng)列向量，f為震源矢量。等式兩邊對(duì)v求導(dǎo)并代入式(3)，可得[6,25-27]

(4)

在彈性波混合域反演中，一般先計(jì)算目標(biāo)函數(shù)關(guān)于λ、μ的導(dǎo)數(shù)，然后通過(guò)鏈?zhǔn)椒▌t得到關(guān)于縱波速度vP、橫波速度vS的導(dǎo)數(shù)。由式(5)～式(8)可推導(dǎo)彈性波頻率域梯度

(5)

(6)

(7)

(8)

1.2 不依賴子波的全波形反演

地震波場(chǎng)記錄d(t) 可以用震源子波s(t)與脈沖響應(yīng)I(t)的卷積表示

d(t)=s(t)?I(t)

(9)

將式(9)變換到頻率域,可得

d(ω)=s(ω)I(ω)

(10)

選第k個(gè)檢波點(diǎn)作為特征道，在頻率域?qū)?shí)測(cè)波場(chǎng)數(shù)乘以模擬波場(chǎng)的特征道，模擬波場(chǎng)數(shù)乘以實(shí)測(cè)波場(chǎng)的特征道[18]，兩個(gè)乘積的差值為

(11)

(12)

式中j表示檢波點(diǎn)號(hào)。把式(11)、式(12)代入式(10)可得x、z分量的殘差統(tǒng)一表達(dá)式

(13)

依據(jù)以上理論，構(gòu)造彈性波的L2范數(shù)的目標(biāo)函數(shù)

(14)

式中：nr表示檢波點(diǎn)總數(shù)。用目標(biāo)函數(shù)對(duì)模型參數(shù)v的偏導(dǎo)數(shù)作為模型更新的梯度方向，可得

(15)

將式(15)展開(kāi)，可得

(16)

方程S(ω)u(ω)=f(ω) 兩邊同時(shí)對(duì)模型參數(shù)v求導(dǎo)并代入式(16)，可以得到

(17)

式中

(18)

式(17)即為不依賴子波的彈性波頻率域梯度表達(dá)式，其形式與式(4)基本相同，唯一不同之處在于伴隨震源的求取方式，其中β的非零項(xiàng)所在的行為第k行。還可以看出，在頻率域以(α-β)*為震源正演模擬得到的波場(chǎng)，或在時(shí)間域?qū)?α-β)反傳波場(chǎng)與正演波場(chǎng)做互相關(guān)就可以求出單炮梯度，最終梯度為所有炮梯度的累加。從圖1可以看出，不依賴子波的彈性波混合域全波形反演的計(jì)算過(guò)程相比于時(shí)間域每一炮的計(jì)算都減少了一次波場(chǎng)反傳的過(guò)程，提高了計(jì)算效率，額外增加的離散傅里葉變換計(jì)算過(guò)程對(duì)整個(gè)計(jì)算效率的影響可以忽略不計(jì)。

圖1 不依賴子波的彈性波混合域全波形反演流程

2 模型試算

對(duì)Overthrust模型進(jìn)行基于彈性波混合域全波形反演四種方案的試算：錯(cuò)誤子波的常規(guī)全波形反演、正確子波的常規(guī)全波形反演和分別采用正確子波、錯(cuò)誤子波的不依賴子波全波形反演。試算范圍為12.5km(x)×4.675km(y)，剖分網(wǎng)格距為25m。觀測(cè)數(shù)據(jù)所用子波為雷克子波，主頻為8.0Hz(圖2a)，模擬數(shù)據(jù)所用子波為原始子波經(jīng)過(guò)Hilbert變換后的子波(圖2b)，即錯(cuò)誤的子波。從子波波形圖可以看出模擬子波與實(shí)際子波在振幅和相位上均有很大區(qū)別。圖3分別為真實(shí)縱(圖3a)、橫(圖3c)波速度模型和經(jīng)過(guò)高斯平滑后的反演初始縱(圖3b)、橫(圖3d)波模型。本文采用全排列接收，采樣間隔為2ms，共500個(gè)檢波器，間隔為25m，炮間距為250m，共49炮。使用二階應(yīng)力—位移方程的交錯(cuò)網(wǎng)格高階有限差分方法進(jìn)行彈性波正演模擬，反演頻率從低頻3.0Hz到高頻16.5Hz，頻率間隔為1.5Hz，反演特征道選擇最小炮檢距道[23]，每個(gè)頻率點(diǎn)迭代20次，共200次。

為驗(yàn)證本文方法的有效性，分別用經(jīng)過(guò)Hilbert變換后的錯(cuò)誤子波進(jìn)行常規(guī)的全波形反演和本文不依賴子波的全波形反演。由于使用的子波與真實(shí)子波差別太大，常規(guī)反演方法合成的地震記錄不僅受模型影響，還受子波影響，結(jié)果與實(shí)際記錄相差太大，無(wú)法獲得正確更新方向，致使計(jì)算誤差越來(lái)越大，最終在第68次迭代時(shí)程序終止。因此，選取第67次迭代更新得到的速度模型反演結(jié)果(圖4),可以看出，用錯(cuò)誤的子波進(jìn)行常規(guī)反演更新得到的模型與實(shí)際速度模型相差很大，說(shuō)明用錯(cuò)誤子波求取的速度更新方向本身是錯(cuò)誤的。圖5為用本文方法得到的第67次迭代結(jié)果，與實(shí)際模型對(duì)比來(lái)看，不依賴子波的方法已經(jīng)對(duì)模型的正確更新初見(jiàn)成效，反演出了模型的大致構(gòu)造形態(tài)。

圖2 模型試算采用的兩種子波

圖3 Overthrust縱、橫波速度模型

圖4 錯(cuò)誤子波常規(guī)方法第67次迭代縱波(a)、橫波(b)速度更新結(jié)果

圖5 錯(cuò)誤子波本文方法第67次迭代縱波(a)、橫波(b)速度更新結(jié)果

為對(duì)比本文方法與常規(guī)方法反傳震源的區(qū)別，選取了頻率為3Hz第15次迭代第21炮x分量的記錄殘差(圖6)，從記錄殘差上可以看到，用不依賴子波的全波形反演得到的記錄殘差與用常規(guī)方法得到的記錄殘差的能量不在同一個(gè)數(shù)量級(jí)上，這主要是由于不依賴子波的反傳震源求取中使用了記錄與特征道相乘的方法。從圖6b可以看出，本文方法的第21炮記錄殘差的最小炮檢距道出現(xiàn)異常，其原因是式(18)中β的特征道是非零項(xiàng)。由此可見(jiàn)，用本文方法得到的反傳震源與式(18)是吻合的。

為進(jìn)一步說(shuō)明本文方法對(duì)子波的不依賴性，分別采用正確的子波和錯(cuò)誤的子波進(jìn)行不依賴子波的全波形反演，得到圖7和圖8所示最終反演結(jié)果，圖9是x=8.75km處的速度對(duì)比曲線。可以看出，用錯(cuò)誤子波得到的結(jié)果與用正確子波得到的結(jié)果基本一致，證明了本文方法不依賴于子波。

圖10是模擬子波與實(shí)際子波一致的情況下常規(guī)正演波形反演得到的彈性波反演結(jié)果。在實(shí)際應(yīng)用中，很難得到準(zhǔn)確的震源子波。從圖8可見(jiàn)，在模型參數(shù)全部一致的情況下，不依賴子波的反演結(jié)果與常規(guī)正演波形反演結(jié)果基本一致。這說(shuō)明即使使用錯(cuò)誤的子波，利用不依賴子波的正演波形反演依然能夠得到正確速度模型。為了對(duì)比兩種方法的模型速度變化，抽取模型x=8.75km處的速度值進(jìn)行對(duì)比(圖11)，可以看出兩種方法的速度誤差極小，速度曲線基本重合。

圖6 第21炮x分量3Hz記錄殘差

圖7 正確子波本文方法的縱波(a)、橫波(b)速度反演結(jié)果

圖8 錯(cuò)誤子波本文方法的縱波(a)、橫波(b)速度反演結(jié)果

圖9 不依賴子波反演x=8.75km處縱波(a)、橫波(b)速度曲線對(duì)比

圖10 正確子波常規(guī)方法的縱波(a)、橫波(b)速度反演結(jié)果

圖11 x=8.75km處縱波(a)、橫波(b)速度不同方法反演結(jié)果對(duì)比

3 結(jié)論與認(rèn)識(shí)

混合域全波形反演結(jié)合了時(shí)間域和頻率域的優(yōu)勢(shì)，不僅提高了計(jì)算效率，還解決了頻率域正演需要海量?jī)?nèi)存的困難。相比聲波全波形反演，彈性波全波形反演具有更豐富的物性參數(shù)，更趨近地下真實(shí)情況，也更有利于復(fù)雜構(gòu)造的多參數(shù)精確建模。用不準(zhǔn)確的子波進(jìn)行常規(guī)全波形反演得到的結(jié)果與真實(shí)模型相差很大，而用不依賴子波的全波形反演在子波不準(zhǔn)確的情況下依然能夠得到較好的反演效果，且不增加額外的正演運(yùn)算過(guò)程。將此方法沿用到彈性波混合域，對(duì)三維實(shí)際資料的全波形反演具有重要意義。