胡鑫 ，王國權(quán)，劉俊洲 ，支麗霞，陳雙全

1 頁巖油氣富集機(jī)理與有效開發(fā)國家重點(diǎn)實(shí)驗室，北京 100083

2 中國石化彈性波理論與探測技術(shù)重點(diǎn)實(shí)驗室，北京 100083

3 中國石油大學(xué)(北京)CNPC物探重點(diǎn)實(shí)驗室，北京 102249

4 中國石油大學(xué)(北京)理學(xué)院，北京 102249

0 引言

目前，石油勘探開發(fā)正向深海、深地及非常規(guī)等油氣資源進(jìn)軍[1]。對于這類油氣資源的勘探開發(fā)，所面臨的地球物理問題變得越來越復(fù)雜，需要提供更為豐富和準(zhǔn)確的油藏儲層特征描述和相應(yīng)的油氣檢測資料。因此，面向疊前道集數(shù)據(jù)的地震反演方法已經(jīng)成為油藏儲層預(yù)測及流體檢測常用的技術(shù)手段[2]。此外，根據(jù)雙相介質(zhì)理論，地下巖石中儲層參數(shù)的微小變化會直接影響到地層的彈性模量，因此，結(jié)合地震巖石物理理論，疊前反演同樣可以實(shí)現(xiàn)對含油氣儲層的定量解釋[3]?？紤]疊前反演方法與各向異性理論的結(jié)合，反演得到的各向異性參數(shù)也有助于裂縫、縫洞型碳酸鹽巖和頁巖等復(fù)雜結(jié)構(gòu)儲層的裂縫預(yù)測[4]和地層壓力[5]、巖石脆性[6]以及地應(yīng)力的反演[7]。并且，隨著OVT(Offset Vector Tile)域五維數(shù)據(jù)采集和處理技術(shù)的發(fā)展，利用疊前地震反演方法開展儲層預(yù)測和流體解釋的精度將進(jìn)一步提升[8]。

疊前地震反演方法的理論基礎(chǔ)是Zoeppritz方程，該方程準(zhǔn)確地建立起了平面波入射到界面處產(chǎn)生的反射與透射波之間的關(guān)系式。但是，由于Zoeppritz方程的強(qiáng)非線性，將其直接應(yīng)用到實(shí)際地震數(shù)據(jù)反演中較為復(fù)雜，因此早期的地震反演一般采用Zoeppritz方程的近似公式。Aki和Richards在假設(shè)入射角角度較小且界面兩側(cè)彈性參數(shù)變化不大的情況下，利用射線參數(shù)與速度之間的關(guān)系，得到了一個新的用于描述反射系數(shù)與縱橫波速度及密度間關(guān)系的表達(dá)式，這個新的近似式被稱為Aki-Richards近似公式[9]。利用該近似式，相對簡化了疊前地震數(shù)據(jù)的反演過程，可以穩(wěn)定地獲得較為準(zhǔn)確的速度和密度等參數(shù)。因此，目前常規(guī)的疊前地震反演方法技術(shù)均是基于簡化的近似公式進(jìn)行的，并在許多商業(yè)化的疊前地震反演軟件中應(yīng)用廣泛[10]。此后，眾多學(xué)者對簡化方法進(jìn)行研究并形成了不同的表達(dá)式，如Shuey近似式[11]，Smith加權(quán)疊加方法[12]，F(xiàn)atti等反演縱橫波阻抗的AVO(Amplitude Variation with Offset)近似式[13]。Gary等建立的彈性參數(shù)反演公式，是目前另一種常用的近似公式[14]。

隨著復(fù)雜油氣儲層逐漸成為油氣勘探的主要目標(biāo)[15-21]，常規(guī)的AVO疊前反演方法已逐漸不能滿足其更高精度的要求，其中尤為突出的問題是密度反演的不穩(wěn)定性。Downton和Ursenbach很早即在研究中指出大偏移距數(shù)據(jù)的缺失會損失疊前反演中剪切模量和密度的精度[22]，而對大偏移距、大角度數(shù)據(jù)的處理能力正是Aki-Richards等一系列近似式所不具備的。此外考慮到單阻抗界面反射差異的問題[23]，降低近似方程因假設(shè)條件限制帶來的模型誤差，拓展現(xiàn)有的反演理論與方法同樣非常關(guān)鍵。因此，選擇精確的Zoeppritz方程直接構(gòu)建反演目標(biāo)函數(shù)，以提高反演的穩(wěn)定性和對實(shí)際數(shù)據(jù)的適用性，已經(jīng)成為當(dāng)前的研究重點(diǎn)和難點(diǎn)之一。早期的研究集中于將Zoeppritz方程改寫為一個由四個獨(dú)立的比值參數(shù)組成的形式[24-25]，新的形式顯著降低了原方程的強(qiáng)非線性，使得直接采用線性反演策略成為可能，但該方法并不能直接得到地下的彈性參數(shù)，因而引入了一定的不確定性。隨著正則化和貝葉斯理論等反演方法的進(jìn)步，基于Zoeppritz方程的疊前參數(shù)直接反演也逐漸得到發(fā)展。正則化方法通過在最小二乘方法的基礎(chǔ)上增加不同的約束項，有效增強(qiáng)反演算子的穩(wěn)定性。Zhi等通過采用IRLM(Iteratively Regularized Levenberg-Marquardt)方法，實(shí)現(xiàn)了疊前三參數(shù)的直接反演，算法整體表現(xiàn)出了很好的魯棒性[26]。Song等在此基礎(chǔ)上通過改寫精確形式的Zoeppritz方程，實(shí)現(xiàn)了楊氏模量和泊松比的高精度同時反演[27]。Ali等提出了約束非線性AVO反演方法，結(jié)合IRLM與Split-Bregman方法，能很好地反演出彈性參數(shù)[28]?；谪惾~斯理論的反演方法同樣應(yīng)用廣泛，考慮到柯西分布能夠促進(jìn)反演結(jié)果的稀疏性，張豐麒等提出了基于柯西分布的精確Zoeppritz方程線性反演方法，有效提升了反演的精度，取得了不錯的效果[29]。

考慮到橫波對巖石孔隙中流體的變化不敏感，更能反映巖石骨架的性質(zhì)，可以提高構(gòu)造解釋的精度，因此同時考慮縱橫波的聯(lián)合反演策略可以很大程度上減少反演結(jié)果的多解性。但考慮到橫波勘探成本的昂貴，且由于橫波能量衰減較快導(dǎo)致純橫波數(shù)據(jù)信噪比較低，因此目前依然廣泛采用縱波激發(fā)、三分量檢波器接收的轉(zhuǎn)換波勘探開展縱橫波聯(lián)合反演[30]。Stewart最先提出縱橫波聯(lián)合反演的方法，將加權(quán)疊加方法推廣到聯(lián)合反演進(jìn)而獲得縱波速度比和橫波速度比等參數(shù)[31]。陳天勝等圍繞P波和P-SV波的反射系數(shù)近似方程提出了一種基于方向加速度最優(yōu)化的速度比值同時掃描的縱橫波聯(lián)合反演方法[32]?？v橫波聯(lián)合反演也逐漸被應(yīng)用于流體識別、儲層刻畫等[33-35]。隨著基于Zoeppritz方程的縱波反演方法的興起，基于精確方程的聯(lián)合反演方法也逐漸成為學(xué)界研究的熱點(diǎn)，并取得了一定的效果[36-40]。

本文基于精確Zoeppritz方程，計算得到的縱波反射系數(shù)以及轉(zhuǎn)換波反射系數(shù)，解決了反射系數(shù)近似公式的各種假設(shè)條件的限制，降低了模型帶來的不確定性，同時利用多波資料進(jìn)行聯(lián)合反演，避免了反演過程中帶來的多解性問題。在反演過程中，采用IRLM方法，合成數(shù)據(jù)的測試結(jié)果表明該方法的抗噪性強(qiáng)且對于初始模型的依賴性較低。最后將該方法應(yīng)用于實(shí)際數(shù)據(jù)，井旁道的對比表明該方法能取得較為可靠的反演結(jié)果。

1 方法原理

1.1 正演模型

本文采用精確Zoeppritz方程來建立彈性參數(shù)與地震觀測數(shù)據(jù)之間的聯(lián)系，利用褶積模型進(jìn)行正演模擬。Zoeppritz方程描述了入射平面波在水平界面發(fā)生波的反射、透射與入射角之間的關(guān)系，反射系數(shù)及透射系數(shù)的大小可以表示成Zoeppritz方程。入射平面波為縱波時，Zoeppritz方程可以表示成：

其中，RPP，RPS，TPP，TPS分別表示縱波反射系數(shù)、橫波反射系數(shù)、縱波透射系數(shù)和橫波透射系數(shù)；VP1和VP2、VS1和VS2、ρ1和ρ2分別表示界面兩側(cè)地層介質(zhì)的縱波速度、橫波速度、密度，θ1和θ2分別為縱波反射角和透射角，φ1和φ2分別為橫波反射角和透射角，下標(biāo)數(shù)字“1”代表上覆地層，下標(biāo)數(shù)字“2”代表下覆地層。

根據(jù)Snell定律，射線參數(shù)與速度和角度滿足以下關(guān)系

式中，p表示射線參數(shù)。Zoeppritz方程精確描述了波的反射系數(shù)、透射系數(shù)與介質(zhì)參數(shù)VP1, VP2, VS1, VS2, ρ1, ρ2和入射角與透射角之間的關(guān)系式。為了能夠利用Zoeppritz方程進(jìn)行反演得到地層參數(shù)，需要將Zoeppritz方程(1)中的反射系數(shù)表示成解析表達(dá)式的形式。根據(jù)Aki和Richards給出的PP波與PS波的反射系數(shù)解析式，設(shè)模型參數(shù)的比值r=[r1,r2,r3,r4]T，其中

縱波反射系數(shù)RPP、轉(zhuǎn)換波反射系數(shù)RPS可以寫成：

其中，

公式(4)(5)給出了基于精確Zoeppritz方程得到的縱波與轉(zhuǎn)換波反射系數(shù)表達(dá)式。根據(jù)褶積模型，疊前不同角度的道集數(shù)據(jù)可以看成是地震子波與地層反射系數(shù)序列的褶積，即

1.2 縱波與轉(zhuǎn)換波聯(lián)合反演

地球物理反演是根據(jù)觀測數(shù)據(jù)與正演模型建立目標(biāo)函數(shù)進(jìn)行迭代求解的過程，本文采用Levenberg-Marquardt方法，并運(yùn)用正則化思想解決Zoeppritz方程求解的非線性和不適定問題。

1.2.1 反演目標(biāo)函數(shù)

首先根據(jù)褶積模型，將對應(yīng)的PP和PS波數(shù)據(jù)的正演過程寫成矩陣形式為：

因此，利用PP波和PS波疊前道集數(shù)據(jù)進(jìn)行聯(lián)合反演，其目標(biāo)函數(shù)

其 中，Σ=diag(λ1,λ2,… , λnp+ns)。聯(lián) 合 反 演 目 標(biāo) 函 數(shù)(10)中，λi(0 ＜λi＜1)代表加權(quán)因子，可以根據(jù)井旁道的反演測試結(jié)果定性地選擇相應(yīng)的權(quán)值。

1.2.2 IRLM反演方法

IRLM方法是求解非線性、不適定反問題的一個十分有效的算法。因此，采用IRLM方法來求解PP波和PS波聯(lián)合反演問題式(10)。對目標(biāo)函數(shù)式(10)進(jìn)行IRLM方法求解，

再根據(jù)迭代更新，得到當(dāng)前第k次迭代后的估計值

其中，α是步長。

上述向量函數(shù)f m()的Jacobi矩陣，可以對向量值函數(shù)f m()兩邊關(guān)于m求偏導(dǎo)，得

其中，WPP，WPS是PP波與PS波對應(yīng)子波的矩陣形式，APP(j)表示PP波反射系數(shù)相應(yīng)于入射角θj關(guān)于縱波速度的偏導(dǎo)數(shù)構(gòu)成的子矩陣；BPP(j)表示PP波反射系數(shù)相應(yīng)于入射角θj關(guān)于橫波速度的偏導(dǎo)數(shù)構(gòu)成的矩陣；CPP(j)表示PP波反射系數(shù)相應(yīng)于入射角θj關(guān)于密度的偏導(dǎo)數(shù)構(gòu)成的子矩陣；APS(j),BPS(j),CPS(j)分別表示PS波反射系數(shù)相應(yīng)于入射角θj關(guān)于縱波速度的偏導(dǎo)數(shù)、橫波速度的偏導(dǎo)數(shù)和密度的偏導(dǎo)數(shù)構(gòu)成的子矩陣，其具體的表達(dá)式見附錄A。具體算法流程如下：

1、給定模型參數(shù)的初值m0，ε，ρ，η，σ最大迭代次數(shù)kmax，令k = 0;

2、計算目標(biāo)函數(shù)f (m0)和梯度g0=?f (m0)，判斷是否滿足終止條件(15)

3、令k k=+1，如果k＜kmax，執(zhí)行以下步驟

1）解方程(16)

2）利用強(qiáng)Wolf條件(17)和(18)確定步長α

其中δ∈(0,1)，β δ∈(,1)，本文取δ=0.0001，β=0.9。

3）利用式(19)更新模型參數(shù)

4）計算目標(biāo)函數(shù)f (m0)和梯度gk+1=?f (mk+1)，判斷是否滿足終止條件式(18)和式(20)

5）根據(jù)式(21)和式(22)更新μk和λk，迭代停止前繼續(xù)步驟3

2 合成數(shù)據(jù)測試

為了驗證聯(lián)合PP和PS波數(shù)據(jù)進(jìn)行精確Zoeppritz方程的疊前反演方法的正確性及可行性，采用某實(shí)際測井?dāng)?shù)據(jù)合成疊前的PP和PS波道集數(shù)據(jù)。如圖1a所示為實(shí)際測井?dāng)?shù)據(jù)曲線，圖中包括縱波速度、橫波速度和密度曲線。合成道集數(shù)據(jù)是根據(jù)Zoepprtiz方程計算得到不同角度的PP波反射系數(shù)和PS波反射系數(shù)，再進(jìn)行與地震子波褶積得到，合成數(shù)據(jù)的不同角度范圍均采用相同的地震子波。如圖1b所示為最小相位子波[41]，子波的采樣間隔為2 ms，長度為200 ms。

圖1 模型數(shù)據(jù)Fig. 1 Model data (a) well logging curves; (b) seismic wavelet

由于縱波正演模擬時，使用近似方程作為正演算子時會產(chǎn)生理論誤差[42]?？紤]到轉(zhuǎn)換波正演模擬，也會產(chǎn)生理論誤差，如圖2所示，圖2a是精確Zoeppritz方程合成的PS地震記錄，圖2b是Aki-Richards近似方程合成的PS地震記錄，圖2c是兩者的理論誤差。據(jù)此分析，精確Zoeppritz方程與Zoeppritz近似方程在正演過程存在的理論誤差，可能會導(dǎo)致在反演過程中反演結(jié)果精度下降。因此直接采用精確Zoeppritz方程進(jìn)行反演，可以提高反演結(jié)果的精度，較Zoeppritz近似方程更有優(yōu)勢。

圖2 PS波合成角道集數(shù)據(jù)及誤差對比Fig. 2 Comparison of PS-wave synthetic gathers and data error (a) Zoeppritz equation; (b) Aki-Richards approximate equation; (c) gather data error

在反演過程中，小的擾動會造成反演求解產(chǎn)生很強(qiáng)的不適定問題，為了驗證本文方法具有較強(qiáng)的穩(wěn)定性，能很好地解決非線性問題，因此模型測試中加入不同信噪比的噪聲，合成相應(yīng)的地震記錄。合成地震記錄如圖3所示，圖3a、3b、3c分別為PP波角度域無噪聲記錄、信噪比為5與信噪比為2的含噪聲記錄；圖3d、3e、3f分別為PS波角度域無噪聲記錄、信噪比為5與信噪比為2的含噪聲記錄?？梢钥闯?，由于噪聲的加入，極大增加了合成地震記錄的擾動性，會導(dǎo)致反演求解產(chǎn)生很大的不適定問題。

圖3 不同信噪比PP波及PS波道集數(shù)據(jù)對比。PP波道集: (a)無噪聲; (b)信噪比為5；(c)信噪比為2；PS波道集: (d)無噪聲; (e)信噪比為5；(f)信噪比為2Fig. 3 Comparison of PP- and PS-waves synthetic gathers with different SNRs. PP-wave gathers: (a) no noise; (b) SNR=5; (c) SNR=2; PS-wave gathers: (d) no noise; (e) SNR=5; (f) SNR=2

圖4給出了利用本文方法對模型數(shù)據(jù)反演得到的結(jié)果，分別進(jìn)行了PP-PS波聯(lián)合反演與PP波單獨(dú)反演，反演中均采用了IRLM方法求解。其中，圖4a為無噪時的反演結(jié)果，圖4b為信噪比為5的反演結(jié)果，圖4c是信噪比為2的反演結(jié)果，圖中的黑色曲線代表初始模型，墨綠色曲線代表真實(shí)模型，紅色曲線代表反演結(jié)果。對比圖4的反演結(jié)果，在隨著數(shù)據(jù)噪聲加大時，不管是PP波與PS波聯(lián)合反演還是PP波單獨(dú)反演，整體反演結(jié)果仍然與真實(shí)結(jié)果吻合較好，信噪比為2時，亦能趨于穩(wěn)定，與信噪比為5時反演結(jié)果無明顯差別，特別是縱波速度與橫波速度。為了凸顯聯(lián)合反演的優(yōu)勢，計算了PP-PS波聯(lián)合反演以及PP波單獨(dú)反演的相關(guān)系數(shù)以及平均相對誤差，如表1所示，信噪比為5到信噪比為2時，相關(guān)系數(shù)與平均相對誤差逐漸趨于穩(wěn)定。但是密度反演結(jié)果，隨著噪聲的加強(qiáng)，與實(shí)際井?dāng)?shù)據(jù)有著一定的偏差，主要原因是由于密度曲線的反演需要更多的大角度范圍數(shù)據(jù)，而模型測試僅使用了三個角度的數(shù)據(jù)進(jìn)行反演，而且密度項受噪聲影響很大。由于整體的反演效果較好，可以很好地解決反演時精確Zoeppritz方程帶來的強(qiáng)非線性與不適定性，同時具有很好地抗噪性。因此，利用合成數(shù)據(jù)測試結(jié)果驗證了提出的反演方法的可行性。

表1 反演結(jié)果的相關(guān)系數(shù)與平均相對誤差對比Table 1 The comparisons of correlation coefficients and average relative error between the inversion results

3 實(shí)際數(shù)據(jù)應(yīng)用

最后，將該反演方法應(yīng)用到實(shí)際地震數(shù)據(jù)的PP波與PS波聯(lián)合反演中，驗證文中提出的反演方法在實(shí)際數(shù)據(jù)反演中的適用性。對于實(shí)際數(shù)據(jù)，疊前縱波與轉(zhuǎn)換波資料分別是縱波時間域與轉(zhuǎn)換波時間域的角度道集數(shù)據(jù)。由于轉(zhuǎn)換波的旅行時經(jīng)過縱波速度入射、橫波速度出射得到，而縱波的雙程旅行時均是縱波速度，因此同一反射層的反射PS波和PP波的雙程旅行時不同，因此，需要對其進(jìn)行時間域的匹配[43]。圖5為轉(zhuǎn)換波時域匹配前后的數(shù)據(jù)及縱波疊后剖面的展示，圖5a是轉(zhuǎn)換波匹配前的剖面，時間范圍是1.0～5.5 s，圖5b是時域匹配后的轉(zhuǎn)換波疊后剖面，時間范圍變到了0.7～3.8 s，圖5c是縱波疊后剖面?？梢钥吹轿雌ヅ淝暗霓D(zhuǎn)換波剖面與縱波剖面，在時間軸上差異很大，同相軸完全對應(yīng)不上。而在經(jīng)過時間域匹配之后，對比圖5b和圖5c，可以看到兩個剖面整體上同相軸比較一致。圖6為匹配后3個CDP道集對比結(jié)果，分別對應(yīng)CDP號為499～501，每個CDP道集有10個角度道(5°～32°)，圖中顯示的時間范圍為目的層段范圍(1.2～2.0 s)，盡管轉(zhuǎn)換波地震記錄信號偏弱，但是可以看到在1.3 s、1.6 s和1.9 s處，縱波與轉(zhuǎn)換波的主要同相軸可以大致對齊，波形也大致可以對齊。

圖5 PS波時間域匹配前后疊加剖面對比Fig. 5 Stacked sections comparsion of PS-wave time registration (a) PS-wave section before time registration; (b) PS-wave section after time registration; (c) the corresponding PP-wave section

圖6 時間域匹配后疊前道集數(shù)據(jù)對比Fig. 6 Comparison of gathers before and after time registration (a) PP-wave pre-stack gathers; (b) Time domain registration PS-wave pre-stack gathers

為了驗證時域匹配的效果，采用互相關(guān)時延譜進(jìn)行計算驗證，圖7為時域匹配后轉(zhuǎn)換波數(shù)據(jù)與縱波數(shù)據(jù)的互相關(guān)時延譜，在整個時間段上，最深的紅色集中在零延時處，其相關(guān)系數(shù)最大，說明轉(zhuǎn)換波數(shù)據(jù)很好地匹配到縱波時間域。

圖7 時域匹配后轉(zhuǎn)換波與縱波數(shù)據(jù)的互相關(guān)時延譜Fig. 7 Cross-correlation time-shift spectrum of the registration PS-wave and PP-wave data

圖8為匹配后的PP波與PS波分角度疊加的實(shí)際地震剖面，圖8a從上到下依次是分10°、20°、30°疊加的縱波疊后剖面，圖8b從上到下依次是分10°、20°、30°疊加的轉(zhuǎn)換波疊后剖面，由于轉(zhuǎn)換波本身數(shù)據(jù)質(zhì)量不如縱波數(shù)據(jù)，因此轉(zhuǎn)換波疊后剖面的成像質(zhì)量不如縱波疊后剖面，尤其是根據(jù)小角度疊加的剖面。但是，通過對比發(fā)現(xiàn)，盡管轉(zhuǎn)換波成像質(zhì)量不佳，但匹配后的效果來看，縱波剖面與轉(zhuǎn)換波剖面的主要同相軸是一致的，匹配后的數(shù)據(jù)能夠滿足聯(lián)合反演的要求。

圖8 三個角度部分疊加剖面對比，三個角度分別為10o，20o和30oFig. 8 Comparison of three different partial stacked sections of angle gathers, the angles are 10o, 20o and 30o, respectively PPwave sections and PS-wave sections

利用IRLM算法進(jìn)行了PP波單獨(dú)反演，分別反演得到相應(yīng)的縱波速度、橫波速度與密度，如圖9所示。通過將井?dāng)?shù)據(jù)插入到反演結(jié)果剖面中，反演結(jié)果與井?dāng)?shù)據(jù)的一致性較好。而且，通過計算井旁道反演結(jié)果與測井?dāng)?shù)據(jù)的縱波速度、橫波速度和密度的歸一化均方根誤差(Normalized Root Mean Square error, NRMSe)，其值分別為0.6%、0.62%、0.28%。說明本文提出的基于IRLM算法的演方法具有較強(qiáng)的抗噪性，有效保證了疊前反演的穩(wěn)定性。

圖9 PP波單獨(dú)反演結(jié)果對比，圖中分別插入了井曲線Fig. 9 The inverted results only using PP-wave gathers. The well-logging curves of velocity and density are inserted in the sections (a) P-wave velocity; (b)S-wave velocity; (c) density

圖10為利用本文方法采用聯(lián)合反演方法得到的縱波速度、橫波速度及密度剖面，反演結(jié)果穩(wěn)定，剖面特征與原始地震數(shù)據(jù)剖面一致，而且橫波速度剖面和密度剖面的結(jié)果較好。為了對比反演結(jié)果的準(zhǔn)確性，圖10中加入了實(shí)際井?dāng)?shù)據(jù)的縱波速度、橫波速度及密度曲線，將測井?dāng)?shù)據(jù)按相同的色標(biāo)進(jìn)行變密度顯示，井旁道反演結(jié)果與測井?dāng)?shù)據(jù)的縱波速度、橫波速度和密度的歸一化均方根誤差分別為0.69%、0.28%、0.12%。通過對比反演結(jié)果與井?dāng)?shù)據(jù)結(jié)果，可以看到反演結(jié)果與井的一致性較好，特別是密度曲線的反演結(jié)果。

圖10 采用IRLM方法反演結(jié)果，圖中分別插入了井曲線Fig. 10 The inverted results by IRLM method. The welllogging curves of velocity and density are inserted in the sections (a) P-wave ; (b) S-wave velocity; (c) density

為了進(jìn)一步驗證本文方法的有效性，對比了與商業(yè)軟件的常規(guī)聯(lián)合反演結(jié)果。如圖11所示，分別為用軟件反演得到對應(yīng)的縱波速度剖面、橫波速度剖面和密度剖面。從圖11可以看到，縱波速度反演結(jié)果較好，波阻抗特征橫向連續(xù)性清晰，與井?dāng)?shù)據(jù)也能很好地對應(yīng)上，井旁道反演結(jié)果與測井?dāng)?shù)據(jù)的縱波速度、橫波速度和密度的歸一化均方根誤差分別為0.61%、0.83%、0.31%。但是，反演得到的橫波剖面整體與縱波一致，具有一定的相關(guān)性，且井?dāng)?shù)據(jù)的對應(yīng)存在較大的誤差。而且，相對于縱橫波速度剖面，密度反演結(jié)果要更差一些，橫向連續(xù)性也不清晰。

圖11 常規(guī)反演方法反演結(jié)果Fig. 11 The inverted results by conventional inversion method (a) P-wave velocity; (b) S-wave velocity; (c) density

對比圖9、圖10以及圖11，反演結(jié)果中橫波速度剖面的黑色方框標(biāo)記處，可以看到本文方法反演得到的速度剖面與測井?dāng)?shù)據(jù)有較好的對應(yīng)關(guān)系，而常規(guī)方法反演的結(jié)果不理想，在典型的高速位與低速位與井?dāng)?shù)據(jù)對應(yīng)不上；而且對于密度剖面，本文方法反演得到的剖面特征橫向連續(xù)性清晰，與井?dāng)?shù)據(jù)的一致性較好，常規(guī)方法反演結(jié)果與井?dāng)?shù)據(jù)一致性低。而且，縱橫波聯(lián)合反演在精度和抗噪性上要高于PP波單獨(dú)反演，能更好地解決反演的不適定問題，尤其是密度項反演的抗噪能力較強(qiáng)。

4 結(jié)論

采用基于Zoeppritz方程的精確PP波與PS波反射系數(shù)表達(dá)式，建立起了基于PP波與PS波疊前道集數(shù)據(jù)的縱橫波聯(lián)合反演方法，推導(dǎo)了詳細(xì)的偏導(dǎo)數(shù)表達(dá)式。同時，在反演求解過程中，采用IRLM反演方法可以很好地解決縱橫波聯(lián)合反演的強(qiáng)非線性與不適定性問題。合成數(shù)據(jù)驗證了本文方法的有效性，并且通過含噪數(shù)據(jù)測試了算法的抗噪能力。在實(shí)際數(shù)據(jù)的反演應(yīng)用中，本文方法反演的縱波速度剖面、橫波速度剖面以及密度剖面，與原始地震剖面保持良好的一致性。相比于常規(guī)的方法，有著較好的優(yōu)勢，能彌補(bǔ)密度項本身受噪聲影響較大的缺陷，并且反演的橫波速度與密度剖面效果很好，與對應(yīng)井?dāng)?shù)據(jù)有著較好的一致性，驗證了本文提出的反演方法的適用性。

本文采用的PP波與PS波角度道集數(shù)據(jù)為不同時間域的數(shù)據(jù)，在進(jìn)行實(shí)際數(shù)據(jù)聯(lián)合反演前，需要利用準(zhǔn)確的速度比將PS時間域的數(shù)據(jù)匹配到PP時間域。因此，基于多波數(shù)據(jù)的時間域匹配是多波地震數(shù)據(jù)聯(lián)合解釋與反演的關(guān)鍵核心，而開展深度域的多波數(shù)據(jù)聯(lián)合直接反演將是后續(xù)研究的重要方向。

附錄A：偏導(dǎo)數(shù)計算表達(dá)式：

為了計算Jacobi矩陣（14），需要計算反射系數(shù)關(guān)于速度和密度的偏導(dǎo)數(shù)。假設(shè)模型參數(shù)為了方便，記（1）式的系數(shù)矩陣為A，反射系數(shù)和透射系數(shù)組成的向量為R，右端常數(shù)項構(gòu)成的向量為B，則有矩陣方程

令（A4）式的矩陣方程兩邊對的m每個分量求導(dǎo)可得

解方程組（A4）和（A5-A18）即可得反射系數(shù)對地層參數(shù)的偏導(dǎo)數(shù)，從而計算得到雅可比矩陣J (m)。