基于最優(yōu)輸運(yùn)原理的陸上單分量資料彈性波全波形反演

李青陽 吳國忱* 王玉梅 慎國強(qiáng) 段沛然

(①中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580；②海洋國家實(shí)驗(yàn)室海洋礦產(chǎn)資源評價與探測技術(shù)功能實(shí)驗(yàn)室，山東青島 266071；③中國石化勝利油田物探研究院，山東東營 257022)

0 引言

全波形反演(FWI)是一種將模擬和觀測數(shù)據(jù)之間的殘差最小化以獲得高分辨率介質(zhì)參數(shù)的方法[1]。Tarantola[2]從波動方程出發(fā)，通過波場殘差逆時傳播計(jì)算梯度方向，建立了基于最小二乘理論的時間域聲波波形反演(AFWI)。在聲波波形反演理論基礎(chǔ)上，Tarantola[3]進(jìn)一步提出了基于二階位移形式的彈性波動方程的波形反演(EFWI)理論，從此奠定了彈性波形反演理論的堅(jiān)實(shí)基礎(chǔ)。

在彈性波反演領(lǐng)域，Mora[4]應(yīng)用模型數(shù)據(jù)對 Tarantola的彈性波波形反演理論[3]進(jìn)行了測試并指出FWI對低頻信息依賴性強(qiáng)、計(jì)算量大且效率低。K?hn等[5]分析了不同參數(shù)化對彈性波反演結(jié)果串?dāng)_現(xiàn)象的影響。Jeong等[6]在頻率域利用超低頻信號實(shí)現(xiàn)了密度和速度的準(zhǔn)確反演。Ren等[7]利用彈性波波動方程解耦壓制參數(shù)之間的串?dāng)_。Wang等[8]提出P/S模式解耦梯度預(yù)條件方法，借助海森與分辨率矩陣“模式”分解建立了波模式解耦的EFWI方法，有效降低縱、橫波速度串?dāng)_效應(yīng)，改善了反演的收斂性。Raknes等[9]在Sleipner區(qū)域利用縱、橫波速度和密度間的經(jīng)驗(yàn)公式實(shí)施了彈性參數(shù)反演。相比縱、橫波速度參數(shù)，密度變化主要影響地震波場的振幅信息，而對諸如傳播時間和相位之類的波擴(kuò)散參數(shù)影響很小，因此密度很難準(zhǔn)確反演，由于其與速度的高度耦合性，通常反演中被設(shè)定為某個常數(shù)。李青陽等[10]利用速度和密度的高度耦合性進(jìn)行了密度建模和反演。

自20世紀(jì)90年代以來，海底電纜(OBC)系統(tǒng)的出現(xiàn)促進(jìn)了海洋地震勘探技術(shù)[11-13]取得了巨大的進(jìn)展。2000年后，陸上多波多分量勘探逐漸發(fā)展起來。對于多分量地震數(shù)據(jù)，EFWI比AFWI更適合于地下彈性參數(shù)的高精度重構(gòu)[14]。在復(fù)雜區(qū)域，AFWI無法描述彈性波的傳播過程[15]，例如P/S轉(zhuǎn)換波和振幅AVO效應(yīng)。彈性波多參數(shù)全波形反演還存在諸多問題，在工程方面主要面臨四方面困難。①計(jì)算量巨大。EFWI中彈性波數(shù)值模擬涉及彈性變量是常規(guī)標(biāo)量波數(shù)值模擬的五倍以上，巨大的內(nèi)存需求對設(shè)備提出極高要求，同時計(jì)算效率嚴(yán)重限制了其在實(shí)際中的應(yīng)用。②多分量采集數(shù)據(jù)不規(guī)則化。由于彈性波方程是矢量方程，高精度、準(zhǔn)確的彈性參數(shù)反演勢必需要高信噪比、高保真的多分量地震數(shù)據(jù)。③初始速度建模和低頻缺失問題。④地震子波的求取問題[16]。

目前大多數(shù)陸上勘探所采集的地震數(shù)據(jù)均是單分量資料，即從垂直地震檢波器接收到的粒子運(yùn)動響應(yīng)。當(dāng)實(shí)際工區(qū)地表存在低降速帶時，根據(jù)斯奈爾定理，這會導(dǎo)致傳播至檢波器的反射波波路徑幾乎與地表垂直，根據(jù)縱橫波傳播方向與偏振方向的關(guān)系，垂直檢波器接收的數(shù)據(jù)集可以稱為純縱波數(shù)據(jù)[17]。在此情況下，經(jīng)典的假設(shè)是將地下視為流體[9]，并使用聲波方程來匹配縱波數(shù)據(jù)，因此AFWI是目前陸上全波形反演應(yīng)用中最常用的方法。但AFWI無法反演橫波速度。盡管單分量采集的地震數(shù)據(jù)只包含縱波，不直接包含橫波，但包含了大量轉(zhuǎn)換縱波信息，例如P-S-P波，同時地下介質(zhì)是彈性的，縱波振幅也滿足彈性AVO特征。因此，用縱波資料反演彈性介質(zhì)中的彈性參數(shù)是可能的，Li等[17]利用標(biāo)量波信號實(shí)現(xiàn)了單分量資料的彈性參數(shù)反演。

FWI是強(qiáng)非線性反演過程[18]，局部極值問題一直困擾著該方法的發(fā)展。目標(biāo)函數(shù)的設(shè)計(jì)對FWI至關(guān)重要，采用不同目標(biāo)函數(shù)進(jìn)行反演能夠改變反演對初始速度和數(shù)據(jù)低頻的依賴程度，有可能解決周期性跳躍問題[19]。最小二乘L2范數(shù)由于其簡便性和高分辨率特性而被廣泛使用，但存在明顯的周期性跳躍現(xiàn)象。當(dāng)背景速度不準(zhǔn)確導(dǎo)致觀測記錄的波形與模擬記錄波形相位差超過半個周期時，F(xiàn)WI通常會收斂到局部最小值[20]。因此在FWI實(shí)現(xiàn)過程中，通常選取保證運(yùn)動學(xué)信息準(zhǔn)確的初始速度，或者是從相對較低的頻率開始反演。

全局優(yōu)化策略是解決FWI周期性跳躍問題的一種解決方案，例如引入蒙特卡洛算法、遺傳算法和模擬退火等方法。但全局優(yōu)化算法計(jì)算量巨大，與之相比，局部優(yōu)化策略的優(yōu)勢更加明顯。局部優(yōu)化策略中，多尺度反演[21]、包絡(luò)反演[22]、反射波全波形反演[23]、波動方程旅行時反演[20]、基于自適應(yīng)濾波器的波形反演[24]均一定程度上解決了FWI周期性跳躍問題。

近年來，最優(yōu)輸運(yùn)理論在地球物理界得到廣泛的關(guān)注，并成功應(yīng)用于FWI，解決局部極值問題有很好的效果。當(dāng)前有兩類基于最優(yōu)輸運(yùn)原理的目標(biāo)函數(shù)。第一類是二次Wasserstein距離(W2距離)[25]。另一類與Kantorovich-Rubinstein(KR)規(guī)范[26]有關(guān)。KR規(guī)范的優(yōu)點(diǎn)是它不需要數(shù)據(jù)滿足非負(fù)性或質(zhì)量守恒條件。

本文首先分析了地震勘探中單分量數(shù)據(jù)存在的AVO特征和轉(zhuǎn)換波信息，由此可知單分量數(shù)據(jù)具備指示地下介質(zhì)彈性參數(shù)(縱、橫波速度和密度)變化的潛力。再從彈性動力學(xué)理論出發(fā)，通過分解本構(gòu)方程為球張量和旋轉(zhuǎn)張量的形式，推導(dǎo)出彈性介質(zhì)中的偽壓力彈性波方程。然后介紹了基于L2范數(shù)構(gòu)建目標(biāo)函數(shù)的偽壓力彈性波波形反演方法(L2反演)，并將最優(yōu)輸運(yùn)理論下的W2范數(shù)引入反演，提出一種基于最優(yōu)輸運(yùn)原理的陸上單分量資料彈性波全波形反演方法(W2反演)，基于W2范數(shù)構(gòu)建的目標(biāo)函數(shù)具有更強(qiáng)的凸性。將基于W2范數(shù)的反演結(jié)果作為初始模型，再進(jìn)行L2反演(W2+L2反演)，克服了對初始模型的依賴性。最后應(yīng)用Marmousi模型模擬數(shù)據(jù)和實(shí)際資料驗(yàn)證了本文方法的正確性和實(shí)用性。

1 理論

1.1 彈性介質(zhì)中的偽壓力彈性波方程

1.1.1 單分量地震勘探AVO特征與轉(zhuǎn)換波

針對陸上單分量地震勘探，通常會用聲波方程為基礎(chǔ)的全波形反演預(yù)測地下的縱波速度參數(shù)，并結(jié)合測井得到的經(jīng)驗(yàn)關(guān)系計(jì)算橫波速度、密度參數(shù)[9]，進(jìn)而進(jìn)行儲層預(yù)測和流體識別[27]。該策略將地下介質(zhì)假設(shè)為流體，即縱波速度和密度參數(shù)化模式。采用這種處理方式主要有兩個原因，第一是傳統(tǒng)的地震信號處理技術(shù)聚焦于縱波的運(yùn)動學(xué)信息，通過聲波方程模擬結(jié)果可以完好地吻合其運(yùn)動學(xué)信息；第二是多分量地震資料進(jìn)行彈性波FWI存在較多的問題，例如帶寬不一致、子波不統(tǒng)一等問題亟待解決。目前陸上勘探采集仍以垂直檢波器為主，所采集地震信號為粒子在垂直地面方向的振動速度或者位移。

實(shí)際上地下介質(zhì)是彈性介質(zhì)，參數(shù)化除縱波速度和密度外，還應(yīng)包括橫波速度，因此采用彈性波方程模擬得到的數(shù)值解比聲波方程更加可靠[1]。真實(shí)的單分量地震數(shù)據(jù)包含大量轉(zhuǎn)換縱波，會對聲波FWI的結(jié)果產(chǎn)生較強(qiáng)的影響。幸運(yùn)的是近地表往往存在一個低降速帶，由斯奈爾定理可知，低降速帶會導(dǎo)致從地下波阻抗界面反射的地震波傳播到檢測器的角度幾乎垂直于地表(圖1)。縱波偏振方向與傳播方向相同，從垂直檢波器接收到的地震波可以視為純縱波，而偏振方向?yàn)樗椒较虻臋M波信號則不能被接收。因此可以假定縱波和橫波在陸上地震數(shù)據(jù)中是自然解耦的(基巖裸露和地表為強(qiáng)各向異性介質(zhì)除外)。針對該類介質(zhì)，本文提出應(yīng)用偽壓力彈性波方程模擬地震數(shù)據(jù)，通過偽壓力彈性波方程模擬結(jié)果與彈性波方程、聲波方程的模擬結(jié)果進(jìn)行對比，分析在單分量地震勘探中數(shù)據(jù)的AVO特征與轉(zhuǎn)換波信息。

圖1 近地表低降速帶彈性介質(zhì)模型地震波路徑示意圖

以一個四層彈性模型(圖2)為例，模型尺寸為2400m×1500m，常密度設(shè)置為1420kg/m3。震源是峰值頻率為20Hz的Ricker子波，位于(1200m，8m)，網(wǎng)格間距8m，時間采用間隔為0.5ms。

圖2 層狀彈性介質(zhì)模型示意圖

圖3a和圖3b分別為彈性波總位移波場和解耦后的彈性縱波位移波場的垂直分量波場快照。與聲波波場(圖3c)相比，彈性波總位移波場(圖3a)包含轉(zhuǎn)換縱波和橫波，它們在陸地地震數(shù)據(jù)中被天然解耦，即垂直檢波器接收到的波場(圖3b)為縱波位移波場的垂直分量。反射界面產(chǎn)生轉(zhuǎn)換波[28]，在第二界面處存在上行的轉(zhuǎn)換P-S-P波(圖3b黑框所示)，該轉(zhuǎn)換縱波偽壓力彈性波方程建模的偽壓力場(圖3d)中同樣存在，而在聲波方程建模的純壓力場(圖3c)中并不存在。

圖3 層狀模型不同方法正演的0.5s時刻波場快照

由于地表低降速帶對縱波和橫波的天然解耦作用，反射波垂直于地表射出，因此解耦彈性波方程的縱波位移分量與偽壓力彈性波方程的偽壓力都可以用于陸上彈性FWI中地震數(shù)據(jù)的匹配，但偽壓力數(shù)據(jù)是匹配觀測數(shù)據(jù)的更好選擇。因?yàn)槭褂每v橫波解耦彈性波方程的建模方法，計(jì)算成本比偽壓力彈性波方程高；而且常規(guī)的地震數(shù)據(jù)處理聚焦于縱波的運(yùn)動學(xué)信息，更適合采用標(biāo)量波波動方程建模。

圖4a和圖4c是彈性縱波垂直分量地震記錄和偽壓力地震記錄，對比可知在彈性介質(zhì)中二者反射波類型是一致的，其轉(zhuǎn)換縱波在聲波地震記錄(圖4b)中無法被模擬，這也表明了橫波速度反演的可能性。圖4c中紅圈中存在縱波的同相軸交叉，對比可知二者分別是淺中層轉(zhuǎn)換反射S-P波和深層反射P波。圖4d為中國東部M區(qū)采集的單分量地震數(shù)據(jù)，對比紅圈部分，該資料同樣存在同相軸交叉現(xiàn)象，驗(yàn)證了實(shí)際資料存在轉(zhuǎn)換縱波。

圖4 層狀模型不同方法正演記錄與實(shí)際地震記錄的對比

由上述分析可見，偽壓力彈性波方程建模不僅考慮了彈性AVO效應(yīng)，同時也包含轉(zhuǎn)換縱波，相比于聲波方程只能反演縱波速度和密度，偽壓力彈性波方程更加貼合實(shí)際介質(zhì)是彈性介質(zhì)的情況，同時還具備了反演橫波速度模型的潛能。關(guān)于壓力與速度數(shù)據(jù)存在的振幅及相位不匹配問題，在實(shí)際資料應(yīng)用中，通過處理(振幅補(bǔ)償、去噪、去面波處理)可用標(biāo)量波方程模擬去匹配實(shí)際數(shù)據(jù)；而子波相位的不匹配問題則可以通過提取實(shí)際數(shù)據(jù)的子波加以解決。

1.1.2 偽壓力彈性波方程的建立

經(jīng)典的彈性波方程是從彈性動力學(xué)中的幾何方程、本構(gòu)方程和牛頓運(yùn)動微分方程三個基本方程導(dǎo)出。表示彈性應(yīng)變和位移關(guān)系的幾何方程又稱柯西方程，為

(1)

式中：e表示應(yīng)變；u表示位移；下標(biāo)i、j的范圍均為1～3，并且滿足愛因斯坦求和約定。

各向同性介質(zhì)的本構(gòu)方程又稱廣義虎克定律，為

τij=Cijklekl=λθδij+2μeij

(2)

式中：Cijkl=λδijδkl+μ(δikδjl+δilδjk)表示剛度矩陣，其中λ和μ表示拉梅常數(shù)；θ為體應(yīng)變；δij為Kronecker符號；τ表示應(yīng)力。

牛頓運(yùn)動微分方程又稱納維爾方程，為

(3)

式(2)可以寫為等價脹縮與剪切張量形式的本構(gòu)方程[29]

(4)

式中：K為體積模量；P=τkk/3和θP=ekk/3分別為應(yīng)力和應(yīng)變張量的第一不變量；Pδij和Sij分別為應(yīng)力張量的球張量和偏斜張量，S可表示為

(5)

θPδij和Tij分別為應(yīng)變張量的球張量和偏斜張量，T可表示為

(6)

由此可得彈性介質(zhì)中不解耦的偽壓力彈性波方程(詳細(xì)推導(dǎo)見附錄A)

(7)

(8)

1.2 基于最優(yōu)輸運(yùn)的偽壓力彈性波波形反演

經(jīng)典全波形反演的目標(biāo)函數(shù)E通常由模擬記錄與觀測記錄之間的差異構(gòu)建

(9)

本文反演采用梯度類算法，其迭代格式為

(10)

為了避免直接計(jì)算海森矩陣，可以使用正向和反向傳播波場能量對梯度進(jìn)行加權(quán)。本文直接選取偽壓力波場互相關(guān)構(gòu)建海森矩陣

(11)

式中：xs為第s個炮點(diǎn)坐標(biāo)；tmax為記錄長度。

迭代步長計(jì)算一般采用線性搜索方法，在模型測試中常用的有非精確線搜索或者拋物線搜索方法。為降低計(jì)算空間復(fù)雜度，本文采用實(shí)際模型與梯度的最大值之比計(jì)算，即

(12)

式中：max(m)為模型參數(shù)m的最大值；max(Δm)為模型參數(shù)的最大變化量；β=0.02，對應(yīng)于2%的最大模型變化量。這種可變步長計(jì)算方法可以有效降低計(jì)算量。

1.2.1 基于L2范數(shù)偽壓力彈性波波形反演

Tarantola[2]應(yīng)用模擬記錄dcal與觀測記錄dobs之間的差異最小的L2范數(shù)為目標(biāo)函數(shù)

(13)

根據(jù)Plessix[30]的伴隨狀態(tài)法可得到目標(biāo)函數(shù)對拉梅參數(shù)和密度的梯度表達(dá)式(詳見附錄B)

(14)

式中上標(biāo)“*”表示反傳波場。

通過鏈?zhǔn)椒傻媚繕?biāo)函數(shù)對縱、橫波速度的梯度

(15)

常規(guī)L2范數(shù)的FWI是一個非線性反演問題，易發(fā)生周期性跳躍現(xiàn)象，進(jìn)而陷入局部極小值。因此對于初始模型精度有很高的要求，在實(shí)際應(yīng)用中往往難以滿足。為解決常規(guī)L2范數(shù)FWI局部極值問題，基于最優(yōu)輸運(yùn)理論，本文構(gòu)建模擬記錄與觀測記錄之間的Wasserstein距離的平方作為新的目標(biāo)函數(shù)，以提高目標(biāo)函數(shù)的凸性，從而提高反演結(jié)果的精度。

1.2.2 基于W2范數(shù)偽壓力彈性波波形反演

最優(yōu)輸運(yùn)理論是用于分析概率分布之間的關(guān)系，找尋一種最優(yōu)化的運(yùn)輸方案。Wasserstein距離表示兩種概率分布之間的差異，其衡量標(biāo)準(zhǔn)是將一種分布重新排列為另一種分布的最優(yōu)輸運(yùn)成本。將模擬和觀測地震數(shù)據(jù)集視為兩個概率密度函數(shù)的概率分布，可以將其想象為總體質(zhì)量相等(同時等于單位1)的兩堆沙子的分布。給定特定的成本函數(shù)，將一堆沙運(yùn)送到另一堆沙的不同方案會導(dǎo)致不同的運(yùn)輸成本，成本最低的方案就是最優(yōu)方案，而Wasserstein距離就是這一最低成本。本文采用二次Wasserstein范數(shù)(W2距離)，定義為

(16)

式中：b和q為概率密度函數(shù)，其對應(yīng)的概率分布函數(shù)分別被定義在X和Y空間；M為將b重新排列或搬運(yùn)為q的所有映射的集合；inf表示集合M的下確界；A表示運(yùn)輸方案。

推廣到地球物理問題，本文采用一維W2目標(biāo)函數(shù)

(17)

式中：Q-1與Q互為逆函數(shù)，為表示方便，后文用Q-1°B(x)代替Q-1[B(x)]；B和Q分別對應(yīng)于模擬數(shù)據(jù)和觀測數(shù)據(jù)的概率分布函數(shù)。地震數(shù)據(jù)不是嚴(yán)格的概率密度函數(shù)，因此本文采用Yang等[25]的地震數(shù)據(jù)預(yù)處理方案。

根據(jù)伴隨狀態(tài)法原理可知，梯度形式只和狀態(tài)方程(偽壓力彈性波方程)有關(guān)，因此梯度仍然是式(14)和式(15)，唯一不同的是反傳虛震源項(xiàng)，Yang等[25]給出了一維虛震源δη的表達(dá)形式

[t-Q-1°B(t)]2dt

(18)

式中dt表示地震數(shù)據(jù)的時間采樣間隔。

2 數(shù)值模擬結(jié)果及分析

2.1 目標(biāo)函數(shù)的凸性分析

選取完全滿足概率密度函數(shù)特征的Gauss子波(不含負(fù)極性)和Ricker子波(含負(fù)極性)測試L2范數(shù)和W2范數(shù)目標(biāo)函數(shù)的凸性。圖5a為Gauss子波模擬數(shù)據(jù)及觀測數(shù)據(jù)，圖5b為Gauss子波的模擬數(shù)據(jù)和觀測數(shù)據(jù)的基于L2范數(shù)與基于W2范數(shù)誤差量隨時移量變化的形態(tài)對比。圖5c為Ricker子波模擬數(shù)據(jù)及觀測數(shù)據(jù)，圖5d模擬數(shù)據(jù)和觀測數(shù)據(jù)的基于L2范數(shù)與基于W2范數(shù)誤差量隨時移量變化的形態(tài)對比。對于高斯子波數(shù)據(jù)，W2范數(shù)構(gòu)建了一個全局凸的目標(biāo)函數(shù)，而L2范數(shù)目標(biāo)函數(shù)雖然也是一個凸函數(shù)，但是在數(shù)據(jù)相位差超過一個波長的情況下，L2范數(shù)不能指示誤差下降。對于Ricker子波數(shù)據(jù)，兩種目標(biāo)函數(shù)曲線均存在局部極小值，但是W2范數(shù)的凸性比L2范數(shù)更好，在一定程度上降低了基于L2范數(shù)FWI的周期性跳躍問題的影響。可見，理論上W2范數(shù)比L2范數(shù)的凸性更好。

圖5 L2范數(shù)和W2范數(shù)目標(biāo)函數(shù)的凸性分析

2.2 Marmousi模型測試

利用重采樣后Marmousi縱、橫波速度模型(圖6)測試本文反演方法效果。模型網(wǎng)格數(shù)為300×150，網(wǎng)格間距均為20m。正演采用主頻為8Hz的Ricker子波作為震源，時間采樣間隔為2ms。共設(shè)置42炮，炮檢點(diǎn)均勻分布于地表，起始炮點(diǎn)位于0，炮間距為140m；每一炮的檢波點(diǎn)數(shù)均為300，起始檢波點(diǎn)位于0，檢波點(diǎn)距為20m。

圖6 Marmousi真實(shí)模型

本文采取兩種策略實(shí)現(xiàn)Marmousi模型彈性參數(shù)反演：一是基于L2范數(shù)的偽壓力彈性波方程縱、橫波速度參數(shù)反演；二是首先用基于W2范數(shù)進(jìn)行彈性參數(shù)反演得到修正的低波數(shù)縱、橫波速度模型，然后再用經(jīng)典L2范數(shù)進(jìn)行高波數(shù)擾動縱橫波速度模型反演。測試中，均不考慮密度參數(shù)的影響，密度參數(shù)設(shè)置為常值1500kg/m3。

為保證本文提出的方法能夠僅利用縱波資料準(zhǔn)確反演彈性參數(shù)模型，采用縱、橫波速度橫向不變、垂向線性增加的初始模型(圖7)。誤差大的初始速度會引起波形反演的周期性跳躍問題，進(jìn)而導(dǎo)致非線性反演落入局部極值。圖8為真實(shí)模型和初始模型采用偽壓力彈性波方程正演模擬得到的第21炮記錄(觀測記錄和模擬記錄)，可以看出模擬記錄不存在反射波信息。

圖7 Marmousi初始模型

圖8 Marmousi模型模擬地震記錄左為觀測記錄；中為初始模型模擬記錄；右為水平翻轉(zhuǎn)觀測記錄

圖9和圖10展示了在L2范數(shù)和W2范數(shù)目標(biāo)函數(shù)下的偽壓力彈性波波形反演的第21炮虛震源記錄。L2反演結(jié)果以中高波數(shù)為主，因此其虛震源記錄頻帶中的高低頻信息分布與觀測記錄頻帶特征基本吻合。由于線性初始速度無法提供可靠的運(yùn)動學(xué)信息，很容易導(dǎo)致FWI落入局部極值。W2的虛震源主要集中在低頻帶，記錄中分布了大量 “豎直線條狀”的低頻能量，這些能量主要更新梯度中的低波數(shù)更新量，它對于殘差模型(真實(shí)模型和初始模型的殘差)中的低波數(shù)分量恢復(fù)具有重要的貢獻(xiàn)，解決了波形反演由于初始速度精度低導(dǎo)致周期性跳躍和局部極值問題。

圖9 L2范數(shù)目標(biāo)函數(shù)的虛擬震源記錄

圖11為基于W2范數(shù)的縱、橫波速度第21次迭代反演結(jié)果。第22次迭代誤差上升，反演終止。反演模型的淺、中、深層的背景模型均得到修正，間斷面沒有被準(zhǔn)確刻畫，說明反演過程中主要以低波數(shù)模型的修正為主，驗(yàn)證了圖10的分析。圖12為用W2反演結(jié)果進(jìn)行偽壓力彈性波正演的模擬記錄和真實(shí)模型模擬的觀測記錄的對比，模擬記錄中強(qiáng)反射界面產(chǎn)生的反射波非常清晰，對比兩側(cè)的觀測記錄，銜接處的反射波位置基本一致，但是W2反演結(jié)果的模擬記錄缺失了高頻信息，箭頭處明顯與右邊觀測記錄不一致，缺失了一些細(xì)小的反射軸信息(高頻信息)。

圖10 W2范數(shù)目標(biāo)函數(shù)的虛擬震源記錄

圖11 Marmousi模型W2反演結(jié)果

基于L2范數(shù)的縱、橫波速度反演誤差在迭代57次后上升，反演終止，結(jié)果如圖13所示。由圖可見，基于L2范數(shù)的反演陷入局部極值且反演結(jié)果嚴(yán)重偏離了真實(shí)模型。為了更直觀的展示該結(jié)果，將真實(shí)縱波速度模型(圖6a)主要構(gòu)造界面映射到L2范數(shù)的縱波速度反演結(jié)果(圖13a黑線所示)，可見，基于L2范數(shù)的反演結(jié)果的地質(zhì)構(gòu)造整體上移，箭頭所指的四處彈性間斷面均與反演速度不一致，因此以線性初始速度模型為驅(qū)動的基于L2范數(shù)的反演失敗。

以W2反演結(jié)果作為初始模型，進(jìn)一步實(shí)施L2反演，W2+L2反演共迭代了78次，結(jié)果如圖14所示。對比圖6a、圖13a、圖14a可見，W2+L2反演結(jié)果的構(gòu)造界面與真實(shí)模型基本吻合，比L2反演結(jié)果更精確，說明W2反演的縱橫波速度模型準(zhǔn)確修正了模型中、低波數(shù)信息。圖15為用W2+L2最終反演結(jié)果進(jìn)行偽壓力彈性波正演的模擬記錄與真實(shí)模型模擬的觀測記錄的對比，可見，模擬記錄豐富了圖12數(shù)據(jù)的高頻信息，并且與觀測記錄吻合良好。

圖12 Marmousi模型W2反演結(jié)果模擬記錄與觀測記錄的對比

圖13 Marmousi模型L2反演結(jié)果

圖14 Marmousi模型W2+L2反演結(jié)果

圖15 Marmousi模型W2+L2反演結(jié)果模擬記錄與觀測記錄的對比

圖16為兩種反演過程中縱波速度模型梯度，可見，在W2反演過程中，速度更新量從淺層開始，逐步遞進(jìn)到中深層，并且更新量主要以低波數(shù)信息為主，尤其是圖16b的梯度剖面上中深層的低波數(shù)的更新量能量均衡，界面位置基本準(zhǔn)確。對比L2反演的第1次迭代與W2+L2反演的第21+1次迭代的梯度剖面(圖16c與16d)，L2反演梯度主要界面完全偏離了真實(shí)構(gòu)造位置(實(shí)線所示)；由于W2反演提供的準(zhǔn)確背景速度，W2+L2反演梯度主要界面與真實(shí)構(gòu)造位置(實(shí)線所示)一致。圖17和圖18分別為Marmousi模型反演結(jié)果與初始、真實(shí)模型的抽道對比(橫向4km處)，可見，W2+L2反演策略能夠準(zhǔn)確地恢復(fù)模型的低—中—高波數(shù)(在現(xiàn)有數(shù)據(jù)頻帶范圍內(nèi))成分。圖19為L2和W2+L2兩種反演策略的歸一化誤差曲線，W2+L2反演收斂速度更快，且誤差更小，結(jié)果更接近真值。

圖16 Marmousi模型反演過程中梯度剖面

圖17 Marmousi縱波速度反演結(jié)果與初值、真值的單道對比

圖18 Marmousi橫波速度反演結(jié)果與初值、真值的單道對比

圖19 Marmousi模型W2+L2和L2反演收斂曲線

2.3 實(shí)際資料應(yīng)用

采用中國東部M工區(qū)的二維地震數(shù)據(jù)驗(yàn)證本文的基于偽壓力彈性波方程的陸上單分量資料彈性波波形反演方法的實(shí)用性。該工區(qū)在近地表存在明顯的低降速帶，反射波路徑幾乎垂直于地面并被垂直檢波器接收，數(shù)據(jù)集可以視為縱波數(shù)據(jù)，符合本文方法的適用條件。

圖20為通過偏移速度分析得到的初始縱波速度模型，橫波速度模型是根據(jù)測井?dāng)?shù)據(jù)和巖石物理建模分析擬合出的經(jīng)驗(yàn)關(guān)系(vS=-728.168+0.748×vP)計(jì)算得到，密度參數(shù)設(shè)為常數(shù)。速度模型網(wǎng)格數(shù)為3151×500，水平和垂直方向的網(wǎng)格間距均為10m。地震數(shù)據(jù)采樣間隔為2ms，長度為4s。為滿足穩(wěn)定性條件，使用傅立葉變換方法進(jìn)行數(shù)據(jù)插值，重采樣間隔為0.8ms。震源和檢波器均在地面，總炮數(shù)為233炮，單炮最大炮檢距為4km，檢波器間距為10m。震源子波是由地震記錄提取，主頻約為20Hz。

圖20 縱波速度初始模型

由于橫波速度與縱波速度滿足經(jīng)驗(yàn)關(guān)系，因此在W2反演過程中，采用以下縱橫波速度梯度作為更新量

(19)

圖21為W2反演迭代5次的縱波速度模型，箭頭所指處更新較為明顯，橫波速度與其滿足經(jīng)驗(yàn)關(guān)系，因此沒有展示。然后將W2反演結(jié)果作為初始模型進(jìn)行L2反演，并采用式(16)的彈性波雙參數(shù)反演梯度進(jìn)行縱、橫波速度迭代更新。圖22為經(jīng)過11次L2反演迭代后縱、橫波速度模型，從圖中箭頭處可以看出反演速度的構(gòu)造細(xì)節(jié)，并且有較高的分辨率。

圖21 W2反演迭代5次的縱波速度模型

圖22 W2+L2反演迭代5+11次后的速度模型

圖23為觀測與用最終反演速度模型模擬的單炮記錄的對比，二者均主要包含縱波信息，由圖中箭頭處可見，模擬記錄與觀測記錄對應(yīng)關(guān)系良好。為驗(yàn)證反演結(jié)果的可靠性，應(yīng)用JASON軟件去除縱、橫波速度低波數(shù)背景信息，得到縱、橫波速度擾動模型，如圖24所示。然后計(jì)算縱橫波速度比(vP/vS)模型，如圖25所示。根據(jù)式(12)可知模型反演更新量與該模型的數(shù)值成比例，因此縱橫波速度的擾動量模型計(jì)算得到的vP/vS數(shù)值仍然穩(wěn)定，而去除背景干擾后，模型高頻細(xì)節(jié)更加突出。圖25白色方框部分的放大如圖26所示，其中兩條黑線為自然電位測井曲線。當(dāng)砂巖含量較高時，泊松比較低，對應(yīng)的vP/vS在測井曲線上表現(xiàn)為低值(圖中白色箭頭所指)，說明反演結(jié)果與測井解釋較吻合。

圖23 W2+L2反演迭代5+11次后速度模型模擬記錄與觀測記錄的對比

圖24 去除低波數(shù)背景的縱(a)、橫(b)波速度擾動模型

圖25 去除低波數(shù)背景的縱橫波速度比模型

圖26 去除低波數(shù)背景的縱橫波速度比模型局部放大顯示黑線為自然電位測井曲線

綜上所述，基于偽壓力彈性波方程的波形反演方法能夠適用于陸上單分量地震數(shù)據(jù)的縱、橫波速度彈性反演，有利于含油氣儲層的預(yù)測。

3 結(jié)論

本文針對陸上地震單分量勘探提出了一種高效、穩(wěn)定的彈性參數(shù)波形反演方法，并通過模型數(shù)據(jù)和實(shí)際數(shù)據(jù)驗(yàn)證了方法的正確性和實(shí)用性。

(1)由于地表低降速帶影響，地震反射波路徑幾乎與地表垂直，垂直分量檢波器采集接收的地震資料可視為縱波數(shù)據(jù)，這類數(shù)據(jù)包含大量轉(zhuǎn)換縱波信號和彈性AVO信息。偽壓力彈性波方程可以適應(yīng)陸上地震單分量勘探的需求，同時滿足轉(zhuǎn)換縱波和彈性AVO特征，具備在彈性介質(zhì)中用標(biāo)量波數(shù)據(jù)實(shí)現(xiàn)彈性參數(shù)反演的能力。

(2)W2范數(shù)相比于L2范數(shù)具有更優(yōu)的凸性。Marmousi模型數(shù)值實(shí)例說明L2范數(shù)目標(biāo)函數(shù)在波形反演應(yīng)用中存在周期性跳躍問題，W2范數(shù)避免了周期性跳躍現(xiàn)象，W2+L2反演策略減弱了波形反演對初始模型的依賴性。就收斂速度而言，W2反演更快。

(3)中國東部M工區(qū)二維實(shí)際資料反演應(yīng)用結(jié)果說明偽壓力彈性波方程反演方法能夠適用于實(shí)際工區(qū)的彈性參數(shù)求取，有利于含油氣儲層的預(yù)測。

本文方法介于縱波勘探和多分量勘探之間，在多分量勘探技術(shù)尚未成熟之前，為彈性參數(shù)反演提供了一個折中的思路。

感謝中國石化勝利油田物探研究院提供了數(shù)據(jù)支持，以及紐約大學(xué)Yunan YANG老師提供了有益的建議和幫助。

附錄A 偽壓力彈性波方程

將牛頓運(yùn)動微分方程(式(3))中應(yīng)力項(xiàng)替換為球張量和偏斜張量之和的形式，即將式(5)代入式(3)，可得

(A-1)

對該式兩邊同時求散度，有

(A-2)

由于P?Kuk,k，調(diào)整位移的時間導(dǎo)數(shù)項(xiàng)，有

(A-3)

可得張量表示的偽壓力彈性波方程

(A-4)

(A-5)

根據(jù)胡克定律，球應(yīng)力分量和偏斜應(yīng)力分量與應(yīng)變的關(guān)系為

(A-6)

將式(A-6)代入式(A-4)中，可得

(A-7)

(A-8)

附錄B 伴隨狀態(tài)法

Plessix[30]在非線性反演中提出伴隨狀態(tài)法，大大減化了非線性反演中梯度的計(jì)算難度。目標(biāo)函數(shù)E用泛函h定義為

E(m)=h[u(m),m]

(B-1)

其伴隨狀態(tài)方程為

G[u(m),m]=0

(B-2)

因此目標(biāo)函數(shù)/伴隨狀態(tài)方程的梯度由正傳波場(波動方程和狀態(tài)方程對模型參數(shù)的導(dǎo)數(shù))和反傳波場互相關(guān)計(jì)算得到，即

(B-3)

式中u*表示反傳波場，由伴隨狀態(tài)方程對模型參數(shù)的導(dǎo)數(shù)作為虛震源反傳計(jì)算

(B-4)

結(jié)合式(A-1)和式(A-6)可得偽壓力彈性波方程的狀態(tài)方程F(φ,m)，其中φ=(P,ux,uy,uz)。下面給出偽壓力彈性波方程

(B-5)

根據(jù)伴隨狀態(tài)法可得目標(biāo)函數(shù)對拉梅常數(shù)的梯度

(B-6)

化簡可得

(B-7)