宋云紅 陳 浩 李 超 趙妙穎

(1 北華航天工業(yè)學(xué)院 廊坊 065000)

(2 北京市海洋深部鉆探測(cè)量工程技術(shù)研究中心 北京 100190)

(3 中國(guó)科學(xué)院聲學(xué)研究所 北京 100190)

(4 中國(guó)科學(xué)院大學(xué) 北京 100049)

0 引言

聲波測(cè)井作為主要的石油測(cè)井方法之一，是探測(cè)地層各向異性的重要手段[1-4]。在橫向各向同性(Vertical transversely isotropic,VTI)地層的聲波測(cè)井中，往往只能測(cè)量得到地層垂直方向的縱波速度VP和橫波速度VSV。Schoenberg 等[5]提出了一種簡(jiǎn)化的近似模型Annie 模型，隨后Thomsen[6]進(jìn)一步將VTI 介質(zhì)的5 個(gè)獨(dú)立彈性模量用3 個(gè)獨(dú)立模量C33、C44和C66代替。Xu等[7]利用經(jīng)典的5個(gè)模量表示的VTI模型測(cè)量的數(shù)據(jù)與采用3個(gè)模量的Annie 模型得到的數(shù)據(jù)對(duì)比，驗(yàn)證了Annie 模型的適用性。因此在測(cè)量獲得地層密度、垂直方向縱橫波速的基礎(chǔ)上，只需得知水平方向橫波速度就可近似表示VTI 地層。在對(duì)VTI 介質(zhì)的研究中，很早人們就發(fā)現(xiàn)了斯通利波對(duì)VTI 地層尤其是慢地層中的剪切彈性模量C66靈敏度較高，也就是說(shuō)斯通利波對(duì)VTI 地層的橫觀各向同性效應(yīng)比較敏感[8-9]。因此人們通常利用斯通利波來(lái)反演VTI地層的水平橫波速度VSH(對(duì)應(yīng)彈性模量C66)，結(jié)合測(cè)井中測(cè)量的垂直橫波速度VSV(對(duì)應(yīng)彈性模量C44)，就可以得到垂直方向與水平方向的橫波各向異性γ[6,10-11]。當(dāng)?shù)貙訖M波速度較快時(shí)，斯通利波對(duì)橫波各向異性的靈敏度下降，繼續(xù)采用斯通利波反演VTI 地層各向異性時(shí)結(jié)果會(huì)存在較大誤差[12]。為了更高精度地獲得快速VTI 地層的各向異性大小，人們轉(zhuǎn)而對(duì)VTI 介質(zhì)中的偶極聲場(chǎng)進(jìn)行了進(jìn)一步的研究。Sunaga 等[13]和Valero 等[14]發(fā)現(xiàn)VTI 介質(zhì)中的偶極彎曲波頻散在低頻段受C44影響較大，然而在高頻段對(duì)C66有較高的靈敏度。Fang 等[15]提出了利用寬頻帶的斯通利波和彎曲波來(lái)共同反演VTI地層的各向異性。Sinha等[16]提出利用一個(gè)等效的各向同性參考模型的彎曲波頻散與測(cè)量彎曲波陣列波形的頻散曲線在不同頻段對(duì)其高靈敏度彈性模量進(jìn)行反演，最終同時(shí)反演VTI 介質(zhì)的5個(gè)模量，這種方法中參考模型的參數(shù)選定會(huì)對(duì)反演結(jié)果產(chǎn)生比較大的影響。Xu等[7]提出將C44和C66兩個(gè)參數(shù)同時(shí)代入彎曲波頻散計(jì)算中，與測(cè)量波形頻散對(duì)比，同時(shí)反演出兩個(gè)橫波速度VSV和VSH以及各向異性大小γ。許松等[17]進(jìn)一步提出利用斯通利波和彎曲波分頻段比較對(duì)應(yīng)的測(cè)量波形的頻散，最終反演VTI 介質(zhì)的各向異性，這種方法在快地層和慢地層中同樣適用。

綜合來(lái)看，VTI 地層中的各向異性反演基本都是采用將參數(shù)代入模式波理論頻散計(jì)算，反復(fù)計(jì)算斯通利波或者彎曲波的理論頻散曲線，來(lái)最佳匹配測(cè)量波形頻散或者慢度。本文旨在對(duì)VTI 地層的模式波理論頻散進(jìn)行規(guī)律化研究，有效提高反演方法的效率，更加適用于大數(shù)據(jù)量的實(shí)際聲波測(cè)井反演。

1 VTI地層理論頻散分析

斯通利波是由井壁處固態(tài)-液態(tài)交界面的存在而產(chǎn)生的面波，主要在井中流體中傳播，所以其對(duì)井中流體彈性參數(shù)(流體速度Vf和密度ρf)的靈敏度在全頻段都保持比較高的水平。此外，前面提到斯通利波對(duì)TI 介質(zhì)的橫波相關(guān)模量(地層橫波速度VSV、VSH和密度ρs)靈敏度較高，對(duì)地層其他彈性模量靈敏度較低[1]。垂直橫波速度VSV在測(cè)井過(guò)程中可以測(cè)量得到，在計(jì)算斯通利波頻散曲線時(shí)可以作為已知量，在這里不多做研究。因此首先主要研究斯通利波頻散隨Vf、ρf、VSH和ρs的變化規(guī)律。

首先研究慢地層中當(dāng)VSV=1000 m/s 時(shí)，VSH在1000～1250 m/s 之間每隔10 m/s 變化，ρs在1900～2500 kg/m3之間每隔20 kg/m3變化，Vf在1400～1600 m/s 之間每隔10 m/s 變化，ρf在900～1300 kg/m3之間每隔10 kg/m3變化，跟隨4個(gè)參數(shù)的變化計(jì)算參數(shù)變化后的斯通利波理論頻散曲線。為了研究斯通利波理論頻散曲線跟隨每個(gè)參數(shù)變化的情況，下面采用控制變量法，對(duì)4個(gè)參數(shù)的變化影響逐個(gè)進(jìn)行分析。

首先分析當(dāng)Vf=1551 m/s、ρf=1000 kg/m3、ρs=2000 kg/m3時(shí)，VSH在1000～1250 m/s 之間每隔10 m/s 變化計(jì)算得到的斯通利波理論頻散曲線，如圖1 所示?？梢钥闯?，斯通利波理論頻散曲線隨VSH的增大而增大，其整體上看變化非常有規(guī)律性。

圖1 斯通利波理論頻散曲線隨VSH 變化Fig.1 Variation of Stoneley wave theory dispersion curve with VSH

斯通利波理論頻散曲線在ρs、Vf和ρf各自單獨(dú)改變時(shí)的變化在圖2 中給出。與斯通利波理論頻散曲線隨VSH的變化規(guī)律類似，參數(shù)變化時(shí)，斯通利波頻散曲線形狀整體而言變化比較規(guī)律。

圖2 斯通利波理論頻散曲線隨參數(shù)ρs、Vf、ρf 變化Fig.2 Variation of Stoneley wave theory dispersion curve with ρs, Vf, ρf

2 VTI地層頻散曲線的插值計(jì)算

通過(guò)以上分析已知，VTI 地層斯通利波理論頻散曲線隨VSH、ρs、Vf和ρf變化非常具有規(guī)律性。本文提出利用插值法計(jì)算其理論頻散曲線。首先采用兩個(gè)橫波速度VSH1和VSH2(VSH1＞VSH2)的理論頻散曲線VSH1_t和VSH2_t來(lái)插值計(jì)算位于VSH1和VSH2之間的VSH3處頻散曲線VSH3_i，并對(duì)VSH3的理論頻散曲線VSH3_t與插值頻散曲線VSH3_i進(jìn)行誤差分析。這里利用線性插值的方法，對(duì)每一個(gè)頻率點(diǎn)f處，有

當(dāng)VSH3位于二者中間位置時(shí)，VSH1_t和VSH2_t前面的系數(shù)相同，均為0.5，下文稱為對(duì)稱插值；當(dāng)VSH3位于其他位置時(shí)，VSH1_t和VSH2_t前面的系數(shù)不同，下文稱為非對(duì)稱插值。首先利用1140 m/s 和1160 m/s、1090 m/s 和1210 m/s 以及1050 m/s和1250 m/s三組VSH1_t和VSH2_t來(lái)對(duì)稱插值1150 m/s處的頻散曲線，插值結(jié)果如圖3(a)所示，可以看出3 組橫波速度的插值頻散曲線與理論頻散曲線都非常接近，利用誤差計(jì)算公式計(jì)算理論頻散與插值頻散之間的整體誤差和最大誤差分別為

圖3 斯通利波理論頻散與VSH 插值頻散曲線Fig.3 Theory dispersion and interpolation curve of different VSH of Stoneley wave

根據(jù)計(jì)算，3 組橫波速度的對(duì)稱插值的整體誤差分別為0.00%、0.08%和0.20%，最大誤差分別為0.04%、0.60%和1.04% (誤差精確到0.01%，因此小于0.005%的記為0.00%)，誤差都非常小，在可接受范圍。

當(dāng)需要計(jì)算頻散曲線的VSH3不位于已知理論頻散的兩個(gè)橫波速度的中間位置時(shí)，需要采用非對(duì)稱插值方式。下面利用1080 m/s 和1120 m/s、1090 m/s和1160 m/s以及1040 m/s和1250 m/s 三組VSH1_t和VSH2_t來(lái)非對(duì)稱插值1110 m/s 處的頻散曲線，插值結(jié)果如圖3(b)所示，可以看出3組橫波速度的插值頻散曲線與理論頻散曲線都非常接近，利用誤差計(jì)算公式計(jì)算理論頻散與非插值頻散之間的整體誤差分別為0.01%、0.02%和0.17%，最大誤差分別為0.09%、0.25%和0.94%，誤差非常小，在誤差允許范圍。通過(guò)上面的誤差分析，可以看出其他參數(shù)保持不變時(shí)，無(wú)論采用對(duì)稱插值還是非對(duì)稱插值方法，水平橫波速度每隔200 m/s 計(jì)算一次理論頻散曲線，就可以高精度地插值出位于其中不同速度的頻散曲線。

對(duì)ρs、Vf和ρf分別為唯一變量的頻散曲線采用對(duì)稱插值和非對(duì)稱插值方式，計(jì)算得到的插值頻散曲線與理論頻散曲線的對(duì)比情況在圖4 中給出。圖4 可以直觀地體現(xiàn)出插值頻散曲線與理論頻散曲線一致性很好，具體的誤差計(jì)算結(jié)果在表1 中給出?？梢钥闯隼眠@3 個(gè)參數(shù)對(duì)斯通利波理論頻散曲線插值計(jì)算的最大誤差的最大值為1.63%，在誤差允許范圍。

表1 斯通利波理論頻散與對(duì)應(yīng)ρs、Vf 和ρf 插值頻散曲線的誤差分析Table 1 Error analysis of theory and interpolation dispersion of Stoneley wave

圖4 斯通利波理論頻散與插值頻散曲線Fig.4 Theory dispersion and interpolation dispersion curve

上面計(jì)算的都是其他參數(shù)保持不變時(shí)，只對(duì)一個(gè)參數(shù)插值分析理論頻散曲線和插值頻散曲線的誤差，下面觀察多參數(shù)插值的結(jié)果。經(jīng)過(guò)上面的分析，已知慢地層中斯通利波理論頻散曲線對(duì)水平橫波速度VSH和流體速度Vf的靈敏度較高，所以對(duì)這兩個(gè)參數(shù)進(jìn)行插值。假設(shè)需要的是參數(shù)為ρf=1000 kg/m3、ρs=2000 kg/m3、Vf=1551 m/s、VSH=1150 m/s 的地層斯通利波理論頻散曲線，已知的是在ρf=1000 kg/m3、ρs=2000 kg/m3時(shí)，流體速度和水平橫波速度分別為Vf=1400 m/s、VSH=1050 m/s；Vf=1600 m/s、VSH=1050 m/s；Vf=1400 m/s、VSH=1250 m/s和Vf=1600 m/s、VSH=1250 m/s 四種地層的斯通利波理論頻散曲線。那么首先對(duì)Vf兩兩非對(duì)稱插值，得到Vf=1551 m/s、VSH=1050 m/s和Vf=1551 m/s、VSH=1250 m/s，再將插值結(jié)果對(duì)VSH對(duì)稱插值，最終得到需要的Vf=1551 m/s、VSH=1150 m/s 地層的斯通利波頻散曲線。圖5給出兩次插值得到的頻散曲線與直接計(jì)算的理論頻散曲線，圖5 中用VSH0表示1150 m/s，分別表示1250 m/s和1050 m/s，Vf0表示1551 m/s，分別表示1600 m/s 和1400 m/s。經(jīng)分析，VSH0Vf0處的理論頻散與兩參數(shù)插值頻散整體誤差為0.24%，最大誤差為0.96%，誤差較小，在誤差接受范圍，也就是說(shuō)對(duì)水平橫波速度VSH和流體速度Vf各自在200 m/s 的步長(zhǎng)范圍內(nèi)取值并進(jìn)行兩參數(shù)插值得到的插值頻散曲線與理論頻散曲線仍非常接近，在反演過(guò)程中可以用插值計(jì)算代替理論頻散曲線的計(jì)算，將提高VTI 各向異性反演的效率，這一方法將在下面快速VTI 各向異性反演方法研究中進(jìn)一步研究和探討。

圖5 斯通利波理論頻散與VSH-Vf 兩參數(shù)插值頻散曲線Fig.5 Theory dispersion and interpolation dispersion curve of two parameters VSH-Vf

3 VTI地層各向異性快速反演方法

在對(duì)垂直橫觀各向同性地層進(jìn)行反演方法研究前，首先對(duì)其進(jìn)行正演模擬，得到VTI 地層的陣列波形數(shù)據(jù)，用于后續(xù)的反演方法研究。由于VTI地層中地層對(duì)稱軸與井軸重合，具有解析性，可以利用實(shí)軸積分法計(jì)算單極子聲源得到的波形數(shù)據(jù)。井中激發(fā)聲源的頻率響應(yīng)譜為

其中，S(ω)是聲源頻譜，A(k,ω)是井對(duì)聲源的頻率-波數(shù)響應(yīng)函數(shù)。將頻率響應(yīng)譜沿k平面的實(shí)軸積分，疊加源輻射的直達(dá)波，就可以得到聲場(chǎng)函數(shù)：

其中，D(0)(ω)是單極子聲源輻射對(duì)波形的貢獻(xiàn)[1]。根據(jù)式(4)就可以計(jì)算得出單極子聲源在井中激發(fā)的波形數(shù)據(jù)。

3.1 傳統(tǒng)VTI各向異性反演方法

在VTI 地層中，尤其是慢速VTI 地層，斯通利波對(duì)其TI 效應(yīng)，也就是彈性參數(shù)C66靈敏度較高，所以人們希望利用斯通利波來(lái)反演VTI 地層的各向異性參數(shù)。利用斯通利波反演VTI地層各向異性利用了頻譜加權(quán)平均慢度定理[18]：

在快地層中，斯通利波對(duì)TI 效應(yīng)靈敏度下降，用其反演VTI各向異性的可靠性受到影響。人們通過(guò)研究快地層中橫波各向異性對(duì)彎曲波頻散特性的作用，提出了利用寬頻帶偶極測(cè)井?dāng)?shù)據(jù)頻散曲線與彎曲波理論頻散曲線對(duì)比來(lái)反演C66[7]：

式(7)中，Sm和Sd分別為計(jì)算得到的彎曲波理論頻散曲線和波形處理得到的頻散曲線，SSV為垂直橫波慢度。式(7)對(duì)垂直橫波慢度和各向異性大小γ同時(shí)進(jìn)行反演，為使目標(biāo)函數(shù)對(duì)兩個(gè)參數(shù)都保持較高的靈敏度，波形數(shù)據(jù)的頻帶應(yīng)包含高頻段、中頻段和低頻段。

由于慢地層中采用了斯通利波波形平均慢度，在進(jìn)行波形慢度反演時(shí)，通常需要確定波形處理時(shí)間窗口，如果確定的時(shí)間窗口不太合適，可能會(huì)導(dǎo)致計(jì)算的平均慢度較大程度地偏離斯通利波時(shí)差。而且需要用到理論頻散曲線與波形頻譜相乘后積分，計(jì)算較為繁瑣?；谏厦嫣岬降穆貙雍涂斓貙又袃煞N不同的反演VTI 地層各向異性的γ的方法，在慢地層中同樣利用理論頻散曲線與波形頻散直接比較來(lái)反演VTI 地層各向異性，給出與式(7)類似的反演目標(biāo)函數(shù)：

這里假設(shè)垂直橫波慢度已知，只對(duì)γ進(jìn)行反演?？梢酝ㄟ^(guò)一維最優(yōu)化算法尋找目標(biāo)函數(shù)最小值，如黃金分割法獲取最優(yōu)解，對(duì)應(yīng)反演得到地層各向異性參數(shù)。

3.2 快速VTI各向異性反演方法

前面提到了傳統(tǒng)VTI 各向異性反演過(guò)程中需要反復(fù)計(jì)算不同VSH時(shí)的斯通利波理論頻散曲線，與波形的頻散曲線對(duì)比，最終確定地層的水平橫波速度和各向異性大小。在實(shí)際測(cè)井?dāng)?shù)據(jù)處理中，數(shù)據(jù)量非常大，反復(fù)計(jì)算理論頻散曲線十分耗時(shí)，因此通過(guò)一種快速的反演方法來(lái)近似確定VTI 地層的各向異性大小。前面發(fā)現(xiàn)利用已知的斯通利波理論頻散曲線可以較高精度地插值計(jì)算出不同參數(shù)大小的頻散曲線?；谶@種高精度的插值頻散曲線，提出一種快速的VTI地層各向異性反演方法：

(1) 首先對(duì)流體速度Vf和地層水平橫波速度VSH兩參數(shù)大間隔取值，建立理論頻散表，本文針對(duì)慢地層，ρf=1000 kg/m3，ρs=2000 kg/m3，Vf取1400 m/s 和1600 m/s 兩個(gè)數(shù)值，VSH取1000 m/s、1200 m/s 和1400 m/s三個(gè)數(shù)值，共計(jì)算6 種組合參數(shù)下的理論頻散曲線，保存在理論頻散表ST中；

(2) 通過(guò)鉆井泥漿參數(shù)得到測(cè)量井孔的流體速度Vf0[1]，ST 中的6 組理論頻散數(shù)據(jù)針對(duì)Vf兩兩線型插值，得到3個(gè)VSH下Vf=Vf0的插值頻散；

(3) 計(jì)算測(cè)量波形的頻散，與3組插值頻散逐個(gè)比較，當(dāng)其處于兩個(gè)頻散中間區(qū)域時(shí)，得到VSH0的區(qū)間范圍；

(4) 對(duì)VSH0的區(qū)間的兩端和進(jìn)行高密度插值，如在200 m/s 的間隔內(nèi)每隔1 m/s 計(jì)算一次插值頻散，比較波形頻散與插值頻散，兩者最接近時(shí)得到VSH0≈，根據(jù)V′SH0 與地層垂直橫波速度VSV可以計(jì)算得到VTI地層各向異性大小。

經(jīng)過(guò)以上4 個(gè)步驟，就可以快速反演出VTI 地層的各向異性大小。

3.3 VTI地層各向異性反演結(jié)果分析

下面以ρf=1000 kg/m3、ρs=2000 kg/m3、Vs=1000 m/s、Vf=1551 m/s、VSH=1150 m/s的VTI 模型為例，根據(jù)Thomsen 定義的各向異性參數(shù)，可以得出該模型的各向異性γ=16.125%。這里假設(shè)其他參數(shù)已知，需要反演VSH的值。首先利用傳統(tǒng)的波形頻散與理論頻散反復(fù)對(duì)比的方法進(jìn)行反演，反演結(jié)果如圖6 所示，反演得到的γ=15.8321%，誤差δ=1.82%。根據(jù)反演參數(shù)計(jì)算的理論頻散與波形頻散一致性較好。在計(jì)算中進(jìn)行了16步尋找最優(yōu)解的過(guò)程，在個(gè)人計(jì)算機(jī)上耗時(shí)118.12 s。

圖6 傳統(tǒng)VTI 各向異性反演結(jié)果Fig.6 Traditional VTI anisotropy inversion result

通過(guò)本文提出的利用頻散插值快速反演VTI各向異性的方法進(jìn)行反演，反演結(jié)果如圖7 所示，反演得到的γ=16.7013%、δ=3.57%，相較傳統(tǒng)反演方法誤差略大，但仍在誤差允許范圍。與傳統(tǒng)反演方法計(jì)算時(shí)使用同一臺(tái)計(jì)算機(jī)設(shè)備，計(jì)算時(shí)間為0.78 s，相較于傳統(tǒng)反演方法所需的118.12 s 大大縮減，在保證反演精度的同時(shí)大幅度地提高了VTI各向異性的反演效率，這一點(diǎn)對(duì)于處理大數(shù)據(jù)量的實(shí)際測(cè)井?dāng)?shù)據(jù)的商業(yè)軟件來(lái)講非常重要，為后續(xù)的測(cè)井解釋商業(yè)軟件設(shè)計(jì)與改進(jìn)提供了理論基礎(chǔ)與指導(dǎo)。

圖7 兩參數(shù)插值VTI 各向異性反演結(jié)果Fig.7 Two-parameter interpolation VTI anisotropy inversion result

利用商業(yè)軟件處理實(shí)際測(cè)井資料時(shí)，除對(duì)于效率的要求外，由于現(xiàn)場(chǎng)測(cè)井的環(huán)境復(fù)雜，測(cè)量的數(shù)據(jù)包含很多因素導(dǎo)致的噪聲，信噪比較低。為了檢驗(yàn)該方法對(duì)于信噪比較低信號(hào)的可靠性，對(duì)數(shù)值模擬的數(shù)據(jù)人為添加噪聲，這里添加了信噪比為0 dB 的高斯白噪聲，添加噪聲后的波形信號(hào)的第一列如圖8(a)所示，可以看到噪聲受到的干擾較嚴(yán)重。對(duì)加噪后的陣列波形數(shù)據(jù)利用本節(jié)提出的快速反演方法進(jìn)行反演，結(jié)果如圖8(b)所示，反演誤差δ=5.00%，仍然比較可靠，證明該反演方法同時(shí)具備很好的抗噪性。

圖8 加噪波形與兩參數(shù)插值VTI 各向異性反演結(jié)果Fig.8 Waveform and two-parameter interpolation inversion result of noised data

4 結(jié)論

(1) 本文研究了斯通利波理論頻散曲線隨井孔和地層彈性模量相關(guān)參數(shù)的變化規(guī)律，發(fā)現(xiàn)理論頻散曲線隨參數(shù)變化整體來(lái)看曲線變化具有一定的規(guī)律性。

(2) 通過(guò)對(duì)斯通利波插值頻散曲線與理論頻散曲線的誤差分析，得出在較大間隔尺度內(nèi)計(jì)算的插值頻散與理論頻散之間的誤差均比較小，在誤差允許范圍。因此對(duì)相關(guān)參數(shù)大間隔取值建立理論頻散數(shù)值表，就可以高精度地、方便快捷地插值計(jì)算出不同參數(shù)值的近似頻散曲線。

(3) 基于以上分析，本文提出了通過(guò)對(duì)參數(shù)大間隔取值建立理論頻散數(shù)值表、高密度插值計(jì)算頻散曲線、比較波形頻散與插值頻散、快速反演VTI各向異性的方法。通過(guò)對(duì)比驗(yàn)證，快速反演方法在保證精度的同時(shí)有效提高了反演效率，且具有很好的抗噪性。本文的研究成果對(duì)于實(shí)際測(cè)井?dāng)?shù)據(jù)處理與解釋有一定的參考意義。