楊 濤 張會星* 史才旺

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

0 引言

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

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

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

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

1 方法理論

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

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

(1)

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

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

(2)

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

(3)

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

(4)

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

(5)

(6)

(7)

(8)

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

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

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

(9)

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

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

(10)

選第k個檢波點作為特征道，在頻率域?qū)崪y波場數(shù)乘以模擬波場的特征道，模擬波場數(shù)乘以實測波場的特征道[18]，兩個乘積的差值為

(11)

(12)

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

(13)

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

(14)

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

(15)

將式(15)展開，可得

(16)

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

(17)

式中

(18)

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

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

2 模型試算

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

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

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

圖3 Overthrust縱、橫波速度模型

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

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

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

為進一步說明本文方法對子波的不依賴性，分別采用正確的子波和錯誤的子波進行不依賴子波的全波形反演，得到圖7和圖8所示最終反演結(jié)果，圖9是x=8.75km處的速度對比曲線?？梢钥闯?，用錯誤子波得到的結(jié)果與用正確子波得到的結(jié)果基本一致，證明了本文方法不依賴于子波。

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

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

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

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

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

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

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

3 結(jié)論與認識

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