利用物理模型三維縱波數據分析HTI介質的方位各向異性

尹志恒，李向陽，魏建新，狄幫讓，張四海，吳滿生

1 中國石油大學（北京）CNPC物探重點實驗室，北京 102249

2 中國石油化工集團國際石油工程有限公司，北京 100029

1 引 言

自然裂縫往往會導致地震波各向異性的出現，利用各向異性技術來研究裂縫型油氣藏已經成為國內外研究的熱點.對于三維縱波勘探來說，縱波走時、動校正速度、振幅以及振幅梯度都會因為裂縫的存在，呈現出方位各向異性的變化.Luo等[1］利用物理模型進行過裂縫識別.尹志恒等[2］利用3條2D測線數據研究了HTI（Traverse Isotropy with Horizontal axis，具有水平對稱軸的橫向各向同性）介質模型中品質因子隨方位角的變化.魏建新等[3-6］通過物理模型研究了不同參數裂縫對縱波的影響.Li[7］利用多分量地震數據進行了裂縫儲層描述.Ekanem[8］從走時和衰減特性兩個方面進行了裂縫探測研究.以往的研究多使用公式推導數值模擬或者實際野外資料來研究裂縫引起的各向異性變化，然而這些研究都存在本身的一些問題：數值方法需要實際數據的驗證，而直接對野外數據進行研究存在地層參數未知、研究難度大等問題.對此，物理模型的研究就突顯其優勢，一方面模型的各項參數已知，另一方面裂縫參數以及裂縫發育方向可控，對于尋求縱波地震響應與各向異性介質之間的關系更為方便有效.此外，過去的一些研究多采用2D數據，雖然可以看出存在各向異性的現象，但是利用3D數據進行的研究能使裂縫引起方位各向異性的認識更加深刻，對實際生產工作更有意義.

本文思路是使用物理模型的三維縱波數據，通過數值模擬結合物理模擬驗證的方法，來研究HTI介質引起的縱波動校正速度以及走時的方位各向異性變化.本文使用的公式出自以Tsvankin為代表的科羅拉多礦院CWP研究小組的研究[9-13］.首先，結合物理模型數據特點進行處理，抽出大面元數據，分方位角進行分析.其次，根據理論公式進行數值模擬.最后，對比數模與物模的結果，做出分析總結.

2 物理模型以及數據采集

如圖1所示，本文使用的物理模型由4個水平層組成.第一層是混合材料制作而成，第二層和第四層的材料是環氧樹脂，這三層都是各向同性介質.第一二層模擬儲層的上覆地層，第四層模擬整體地質體的底層.模型的第三層是本文的目的層，即人工石灰巖儲層，材料是由纖維織物和樹脂混合而成，模擬由垂直裂縫形成的HTI介質.在裂縫層中有兩個構造，一個是穹窿構造，另外一個是斷層構造，它們不是我們研究的重點.模型尺寸與實際地質體的比例是1∶10000，數據采集使用的換能器的頻率與實際情況比例為10000∶1，模型的速度與實際情況比例為1∶1，模型示意圖見圖1.實驗室通過測試不同入射角的速度信息計算得出了HTI介質的彈性參數，結果呈現弱的正交各向異性.雖然，我們是按照單組裂縫設計對目的層進行加工制造，但是由于制造工藝方面的問題，在進行材料壓實時，不可避免產生微裂縫.這一點也更符合野外的實際情況，造成正交各向異性的出現.通過計算XZ和YZ面內的Thomsen縱波各向異性參數ε，我們發現εYZ值為0.372，要遠大于εXZ的值0.076.這說明雖然裂縫層呈現弱正交各向異性，但是縱波各向異性的變化主要是由于YZ面發育的單組裂縫引起的，所以還是把它視為HTI介質進行分析（見表2）.

表1 模型的速度和密度參數Table 1 Velocity and density parameters of the model

圖1 穹窿模型示意圖（a）3D物理模型；（b）采集時沿X方向中線的2D切片示意圖.Fig.1 Sketch map of physical model（a）3Dphysical model；（b）2Dsection through middle line in Xdirection.

表2 HTI介質的彈性參數（單位109 N／m－2）以及Thomsen參數Table 2 Elastic constants（109 N／m－2）and Thomsen parameters of HTI media

進行3D數據采集時，將模型至于水槽之中，水面到模型頂面距離10mm.為保證后續分析的效果，我們設計了寬方位大偏移距的觀測系統，參數見表3.數據的最大偏移距／目的層深度的值為1.33，方位角和偏移距關系如圖2所示.數據crossline方向是Y方向，即平行于裂縫走向；inline方向是X方向，即垂直于裂縫走向.數據采集所使用的換能器是由本實驗室自己研制的，激發頻率是300kHz，頻帶范圍是120～470kHz，有負載時，主頻會降低.

圖2 滿覆蓋區域的方位角-偏移距關系圖Fig.2 Azimuth-offset in full-fold area

3 初步處理分析

數據的處理使用ProMAX軟件.物理模型數據信噪比高，處理流程相對簡單：建觀、加載、去噪、預測反褶積、速度分析、動校正、疊加、疊后Kirchhoff時間偏移.做最后兩步的目的是為了確定模型目的層同相軸的位置，為后續分析做準備以及觀察HTI介質對構造成像的影響.結果如圖3所示.可以看出，HTI介質對構造的成像沒有大的影響，成像都比較清晰準確.

由于模型都是水平層，所以未作疊前偏移，進行后續分析使用的數據是動校正后的道集數據.我們先提取了一個220m×220m的超道集，定性地分析一下HTI介質中的各向異性現象，這里以及以后的單位都已經經過換算，與實際情況一致.超道集的位置位于構造之外，避免了構造對分析結果的影響.在分析中，定義X方向為方位角0°方向，Y軸負方向為90°方向.

首先，將超道集里面偏移距0～1500m的方位角0°數據與方位角90°數據提出來，進行比較.可以看出在遠偏移距，走時會出現不同，如圖4方框所示，垂直于裂縫走向的方位角0°數據走時更長一些.由于觀測系統的原因，方位角0°附近的數據要更豐富.

然后，再把按照常規流程進行動校正之后的超道集數據，按照方位角90°～270°進行疊加顯示.采用這段數據，原因是圖2顯示這一范圍，信息豐富.如圖5顯示，對于近偏移距數據，各向同性的動校正方法誤差不大，軸基本都已經拉平.但是隨著偏移距的增加，就會出現校正不平的情況.隨著方位角的變化，目的層同相軸不再是水平的，而是出現波浪形狀.波谷出現在方位角180°的附近，說明動校正速度最小，縱波傳播的慢，這符合垂直于裂縫方向波的傳播特征.

由上述定性分析可以知道，裂縫介質的存在會對縱波走時以及動校正速度產生影響，要想進一步了解這種影響，則需要定量研究.

4 疊加速度分析

在各向異性的研究中，傳統彈性參數的分析方法存在參數個數多以及參數物理意義不明確的缺點 .鑒于此，Thomsen[14］（1986）提出了TI（TransverseIsotropic，橫向各向同性）介質中五個各向異性系數.Tsvankin[9］（1997）按照這個思想，提出了針對 HTI介質的五個分析參數，分別是：

表3 觀測系統參數Table 3 Geometry of model

這五個參數與Thomsen參數類似，無量綱，其中VPvert和VS⊥vert分別表示縱波和橫波的垂直速度，ε（V）描述了HTI介質中縱波各向異性差異，γ（V）描述了橫波的各向異性差異，δ（V）與波前的橢圓形狀有關.

此外，Tsvankin[9］（1997）還提出了 HTI介質中動校正速度的表達式：

圖5 超道集中按方位角疊加剖面（a）偏移距0～800m；（b）偏移距800～1600m.Fig.5 Stack profile at azimuth 90°～270°in super-gather（a）Offset 0～800m；（b）Offset 800～1600m.

其中α是測線與HTI介質對稱軸X方向的夾角，即所謂的方位角.對于縱波來說，A值為2δ（V）.將物理模型的參數代入公式，很容易得出單層的Vnmo（即（7）式中的Vnmoi）.要得到目的層的動校正速度還需要運用下列公式：

其中t0是反射層N的零偏移距的雙程走時，Vnmoi是每個單層i的動校正速度，Δti是層i的零偏移距的雙程走時.模型每層的厚度和速度已知，而且是各向同性層，所以代入公式很容易就得到目的層隨方位角變化的動校正速度.

對于物理模型來說，我們用ProMAX處理軟件來獲取Vnmo，軟件使用傳統速度掃描的方法.先將上邊提取的超道集按照方位角分為子道集.為了使得子道集中有足夠的數據進行速度分析，每個方位角道集包括方位角±15°的數據信息.例如，將75°～105°的數據都視為90°的方位角子道集的數據.對于HTI介質，由于對稱原因，不需要全部360°的數據進行分析，所以從0°到180°的每隔30°共分為了7個方位角子道集.然后對不同的方位角道集進行速度分析，記錄下目的層位置的動校正速度.為了保證結果的有效性，又拾取了另外兩個超道集，進行了相似的處理流程.然后記錄下ProMAX處理軟件在各個方位角道集按照各向同性分析得到的Vnmo平均值，進行分析的三個超道集拾取的動校正速度誤差小于±30m／s.同時，記錄了由公式（7）得到的理論數值，并與實驗數據進行對比，見表4.

表4 不同方位角的目的層動校正速度Table 4 Azimuthal different Vnmoof target zone

通過對比分析，可以看出兩者擬合的效果還是比較理想的.這是因為對于物理模型，各層介質的速度等參數都是已知的，所以與理論預測結果比較接近.從擬合的曲線可以看出，對于HTI介質，動校正速度隨方位角呈現橢圓形的變化（見圖6）.我們可以根據橢圓的長短軸來識別裂縫的走向，長軸指示裂縫的發育方向，短軸是其垂直方向，而從長短軸的比值可以看出介質速度各向異性的大小.

圖6 目的層動校正速度（m／s）數值模擬與實驗數據對比（a）直角坐標顯示；（b）極坐標顯示.Fig.6 Comparison of Vnmo（m／s）between numerical simulation and experimental data（a）Rectangular coordinates；（b）Polar coordinates.

5 走時分析

對于各向異性介質，時距曲線不再是標準的雙曲線形狀.Al-Dajani和 Tsvankin[10］（1998）提出了更為準確的各向異性介質中縱波走時的表達式：

其中h是層厚度，Vnmo用式（6）計算得出，即單層的動校正速度；t0是目的層的零偏移距的雙程走時；α是方位角；θ是入射角.通過上述公式的計算，就可以得到目的層單層的走時.將得到的單層時間加上上覆各向同性地層的走時，就可以得到隨方位角變化的目的位置（裂縫層底界面）全程反射時間曲線.從公式中可以看出走時同時受到入射角以及方位角的影響.

對于物理模型數據，也要按照不同入射角和方位角進行分析.同第4節一樣，分析是基于動校正之前的超道集數據.本來想要分為入射角10°，20°，30°三組數據進行觀察，對于更大的入射角數據，已經無法很好地區分有效波和干擾波.從式（13）得出，對應裂縫層底界面深度，入射角10°，20°，30°的位置分別對應著偏移距541m，1117m，1771m.每組包括偏移距±40m的信息.例如，將500～580m的數據都視為入射角10°的數據進行分析.

將不同入射角數據按照方位角進行疊加顯示.圖7是根據公式模擬出來的裂縫層底界面的走時方位變化曲線.從圖中可以看出，隨著入射角的增大，即偏移距的增大，各向異性的情況在變強.圖8實驗數據驗證了這一變化，我們使用其它2組數據也進行了分析，也呈現相同的變化趨勢，這里只展示了效果比較好的一組.說明對于近偏移距數據處理時，可以忽視各向異性的影響，但是偏移距越大影響越大，如果再忽略各向異性，就會出現前邊提到過的動校正不平的情況.可惜再大偏移距的走時變化，需要更有效的處理手段才能夠顯示出來，然而入射角30°的信息已經能夠說明問題.

公式（8）比傳統的雙曲線公式多了一個隨偏移距變化的四階項，它描述了遠偏移距的走時的變化.從實驗數據來看，公式（8）是準確的，可以為將來各向異性處理以及反演提供幫助.

圖7 隨方位角變化的不同入射角走時（數值模擬）Fig.7 Different azimuthal travel time with incident angle（numerical simulation）

圖8 隨方位角變化的不同入射角走時（實驗數據），從左至右依次是目的層深度入射角10°，20°，30°的走時顯示Fig.8 Different azimuthal travel time with incident angle（experimental data），from left to right，travel time in incident angle 10°，20°，30°

6 各向異性程度分析

為了方便進行各向異性差異的分析，首先我們模仿Thomsen（1986）分析縱波各向異性程度大小的方法，定義了一個用參數最大值和最小值來表征各向異性大小的量ΔZ：

其中Z就是用于描述的參數.

根據這個公式，我們計算得到動校正速度的實驗數據的ΔVnmo值為0.13，數值模擬的ΔVnmo值為0.16.由于時間剖面上的目的層走時難以準確拾取，只計算了數值模擬時30°入射角的走時ΔT值為0.10.而由已知的物理模型參數，即裂縫介質不同方向的縱波速度，計算得到的ΔVP值為0.18.由上邊的分析可以看出，裂縫引起的方位各向異性變化是一致的，即平行于裂縫走向與垂直于裂縫走向的方向上，各種分析參數的差異最大.和已知物理參數分析結果對比，實驗得出的動校正速度和走時的方位各向異性程度要略小，這可能是由于采集或者后續處理分析手段等影響因素引起的.

首先，我們認為最大偏移距／目標深度比值對各向異性分析影響較大.入射角信息，對于數據來說，是無法改變的，它和數據采集時的觀測系統定義有關.入射角反映在觀測系統上，就是最大偏移距與目標深度的比值，比值越大各向異性程度越大.走時各向異性差異大小為0.13時，入射角必須要達到37°，即為最大偏移距／目標深度比值等于1.5；當比值為2.5時，即入射角為51°，走時的ΔT各向異性值達到了0.18.這和圖7的變化趨勢是相同的，只是圖7中各向異性程度要小一些.

對于同一組數據，在用縱波傳播特征來進行裂縫識別等分析時，我們當然希望使用各向異性更加明顯的參數.因為各向異性程度越大變化越明顯，識別分析起來將更加容易，特別是對于實際數據，處理相對困難，有效波不是很明顯的變化，可能會被噪聲等其它干擾信息所掩蓋.所以大偏移距觀測系統的設計對于走時各向異性分析至關重要.其實，對于動校正速度分析也很重要.從上邊的分析可以看出，當只有近偏移距信息時，各向異性影響不明顯，不可能拾取得到準確的動校正速度.當遠偏移距信息足夠豐富，才能有效地進行速度拾取.

此外，在進行實驗研究的過程中，我們發現方位角信息分布不均勻對后續分析影響也比較大.從圖2可以看出，模型裂縫的發育方向的方位角為90°，但由于觀測系統的設計，造成這一區域的信息量最小，為了進行分析，必須加入此方位角周圍數據.如果方位角信息的分布更均勻一些的話，就可以進行更細致更準確的方位角研究.所以，在了解到裂縫大致發育方向的情況下，應該設計更合理的觀測系統，保證在此方向上數據量豐富.這一點對于實驗研究以及實際野外采集都很有指導意義.

當然，高質量的數據資料，尤其是好的遠道數據，也是進行方位各向異性分析所需要的，質量差的資料很可能造成結果的不準確，甚至研究失敗.

7 結 論

本文通過使用實驗室穹窿物理模型三維實驗數據研究了裂縫引起的HTI介質中動校正速度以及走時隨方位角的變化.由理論模擬和實驗分析，得到以下幾點認識：

（1）實驗數據顯示了HTI介質引起方位各向異性現象和理論公式預測結果耦合比較好，驗證了公式的準確性.文章中提到的理論公式可以用于HTI介質的處理流程或者速度反演分析等方面.

（2）理論值與實驗值的良好耦合，是由于物理模型各項參數已知，應用研究比較方便.這一方面說明在實際生產中，獲取準確地層參數對各向異性研究的重要意義，另外一方面也說明了物理模型在各向異性研究中的獨特優勢.

（3）通過分析研究，認識到最大偏移距與目標深度比值以及方位角信息的分布對方位各向異性的研究影響比較大，尤其是前者.這對三維觀測系統的設計有指導意義.

致 謝 感謝中國石油大學（北京）劉洋教授和陳雙全老師的指導.

（References）

[1]Luo M，Arihara N，Wang S G，et al.Abnormal transmission attenuation and its impact on seismic-fracture prediction—A physical modeling study.Geophysics，2006，71（1）：15-22.

[2]尹志恒，魏建新，狄幫讓等.利用Q值各向異性識別裂縫走向.石油地球物理勘探，2011，46（3）：429-433.Yin Z H，Wei J X，Di B R，et al.Fracture orientation detection using Q anisotropy.Oil Geophy.Prospec.（in Chinese），2011，46（3）：429-433.

[3]Wei J X，Di B R，Li X Y.Effects of fracture scale length and aperture on seismic waves：An experimental study.J.Seis.Explor.，2007，16（2）：265-280.

[4]魏建新，狄幫讓.裂隙密度對縱波傳播特性影響的實驗觀測.石油地球物理勘探，2007，42（5）：554-559.Wei J X，Di B R.Experimentally surveying influence of fractural density on P-wave propagating characters.Oil Geophy.Prospec.（in Chinese），2007，42（5）：554-559.

[5]魏建新，王椿鏞.各向異性介質中橫波測試技術實驗研究.石油物探，2004，43（5）：427-432.Wei J X，Wang C Y.Experiment on shear wave observation in anisotropy medium.Geophys.Prospec.Petrol.（in Chinese），2004，43（5）：427-432.

[6]魏建新，王椿鏞.橫波測試技術的實驗室研究.石油地球物理勘探，2003，38（6）：630-635.Wei J X，Wang C Y.Study of S-wave test and measurement technique in laboratory.Oil Geophys.Prospec.（in Chinese），2003，38（6）：630-635.

[7]Li X Y.Fracture reservoir delineating using multi-component seismic data.Geophy.Prospec.，1997，54（1）：39-64.

[8]Ekanem A M.Fracture detection using 2-D P-wave seismic data：A seismic physical modeling study.SEG Technical Program Expanded Abstracts，2009，28：2647-2651.

[9]Tsvankin I.Reflection moveout and parameter estimation for horizontal transverse isotropy.Geophysics，1997，62（2）：614-629.

[10]Al-Dajani A，Tsvankin I.Nonhyperbolic reflection moveout for horizontal transverse isotropy.Geophysics，1998，63（5）：1738-1753.

[11]Tsvanki I.Normal moveout from dipping reflectors in anisotropic media.Geophysics，1995，60（1）：268-284.

[12]Tsvanki I.P-wave signatures and notation for transversely isotropic media：An overview.Geophysics，1996，61（2）：467-483.

[13]Tsvankin I， Chesnokov E.Synthesis of body-wave seismograms from point sources in anisotropic media.J.Geophys.Res.，1990，95（B7）：11317-11331.

[14]Thomsen L.Weak elastic anisotropy.Geophysics，1986，51（10）：1954-1966.