施俊杰，張建中

（1.廈門大學(xué)信息科學(xué)與技術(shù)學(xué)院，福建 廈門361005；2.中國海洋大學(xué)海洋地球科學(xué)學(xué)院，山東 青島266100）

井間地震層析成像是估算井間介質(zhì)速度最有力的技術(shù)之一。自McMechan［1］提出井間地震層析成像以來，人們提出了不少井間地震走時(shí)層析成像方法。但這類方法幾乎都是使用直達(dá)波走時(shí)數(shù)據(jù)，往往只能給出地層速度，不能直接準(zhǔn)確給出地層界面展布，且反演速度的分辨率不能滿足勘探和開發(fā)的需要。

反射波走時(shí)比直達(dá)波走時(shí)對地層分界面的變化敏感，而直達(dá)波走時(shí)對地層速度的變化更敏感，因此聯(lián)合兩者同時(shí)反演分界面和速度，將會發(fā)揮各自的優(yōu)勢，并增加觀測數(shù)據(jù)個(gè)數(shù)，減少層析反演問題的多解性。Calnan和Schuster［2］利用垂直井、水平地層二維模型，在已知水平反射界面和直射線路徑的假設(shè)下對地層速度進(jìn)行聯(lián)合層析反演，結(jié)果表明聯(lián)合走時(shí)層析成像比直達(dá)波走時(shí)層析成像具有更高的速度分辨率。Wu等［3］用二維矩形單元離散地層，用斜線表示反射界面，提出了1種二維井間曲射線反射和透射層析方法，提高了層析成像的分辨率，可以分辨出薄層。Van Schaack［4］用規(guī)則矩形網(wǎng)格離散地層，用一元深度函數(shù)表示反射界面，建立了地層與反射界面不耦合的井間二維離散模型，以及與該離散模型對應(yīng)的直達(dá)波走時(shí)與層速度、反射波走時(shí)與界面深度的層析反演方程，提出了直達(dá)波和反射波走時(shí)聯(lián)合層析成像方法。Spetzler［5］利用物理實(shí)驗(yàn)數(shù)據(jù)，研究了有限頻率直達(dá)波和反射波走時(shí)聯(lián)合層析反演問題。Bube and Langan［6］利用兩層均勻介質(zhì)的簡單井間模型，從理論上研究了井間直達(dá)波和反射波走時(shí)聯(lián)合層析反演的分辨率問題，認(rèn)為聯(lián)合層析反演可以很好地確定反射層深度，提高慢度分辨率。這些研究使用簡單的或者地層與反射界面不耦合的二維井間模型，分別用直達(dá)波走時(shí)約束層速度，用反射波走時(shí)推斷反射界面，還沒有達(dá)到聯(lián)合使用2種走時(shí)，同時(shí)反演層速度和反射界面的更符合實(shí)際情形的程度。

實(shí)際井孔軌跡一般不是過井口的鉛垂線，而是1條曲線，地震波傳播路徑也是在1個(gè)三維空間內(nèi)展布，因此在三維空間進(jìn)行射線追蹤和層析反演，更符合實(shí)際情況。另一方面，由于激發(fā)點(diǎn)和接收點(diǎn)只能放在井孔內(nèi)，觀測孔徑有限，地震波射線難以覆蓋對應(yīng)的三維空間，也難以為三維反演問題提供足夠的旅行時(shí)數(shù)據(jù)，從而增大了反演問題的不適定性。綜合考慮這些情況，本文建立了基于不規(guī)則單元的地層與反射界面完全耦合的井間三維離散模型，并假設(shè)在垂直于連井截面方向上模型參數(shù)不變，即在正演時(shí)采用三維模型，以模擬三維空間的射線分布和走時(shí)數(shù)據(jù)；反演時(shí)相當(dāng)于采用二維模型，以減小反演問題的不適定性。在此不規(guī)則單元離散模型基礎(chǔ)上，建立了聯(lián)合使用直達(dá)波和反射波走時(shí)同時(shí)反演地層速度和反射界面的層析反演方程，實(shí)現(xiàn)了井間直達(dá)波和反射波走時(shí)的聯(lián)合層析反演。本方法已在文獻(xiàn)［7］中得到了初步應(yīng)用，本文將詳細(xì)討論本方法并在模型應(yīng)用上展開分析。

1 方法原理

采用張建中等［8］的方法建立井間離散模型，圖1所示是三維離散模型的xoz截面圖。

模型中地層速度和分界面在連接左右兩井軌跡中線的xoz平面內(nèi)變化，而在垂直于該平面的y軸方向上保持不變。對于正演計(jì)算，模型在y軸方向上被等間距地離散成幾個(gè)單元，對于反演問題，在y軸方向上相當(dāng)于只有1個(gè)單元。由于在反射界面上速度發(fā)生突變，每個(gè)節(jié)點(diǎn)需設(shè)置地層上下2個(gè)速度值。這種離散模型的地層與界面完全耦合。

圖1 井間離散模型截面示意圖Fig.1 Sketched section of crosshole discretized model

考慮到走時(shí)數(shù)據(jù)存在測量和拾取誤差，反問題的不適定性等因素，構(gòu)造了下列聯(lián)合層析反演問題的目標(biāo)函數(shù)：

對上述非線性目標(biāo)函數(shù)在初始模型鄰域用泰勒級數(shù)展開至線性項(xiàng)，利用函數(shù)極小點(diǎn)滿足的關(guān)系式，得到反演的迭代方程：

其中：k代表迭代次數(shù)；A表示Frechet矩陣；Δb表示觀測走時(shí)與當(dāng)前迭代模型計(jì)算走時(shí)之差；x表示模型參數(shù)；Δx為x的修正量，即線性方程中的未知量。各量的具體表達(dá)式可以寫為：

式中：si是地層第i個(gè)節(jié)點(diǎn)上的慢度；zj是界面上第j個(gè)節(jié)點(diǎn)上的深度。

直達(dá)波走時(shí)和反射波走時(shí)對地層單元節(jié)點(diǎn)慢度si的導(dǎo)數(shù)計(jì)算式可以統(tǒng)一為：

其中：K為射線穿過的以節(jié)點(diǎn)i為頂點(diǎn)的單元個(gè)數(shù)；lk為第k個(gè)單元內(nèi)的射線長度，由于網(wǎng)格慢度由4個(gè)頂點(diǎn)慢度平均所得，所以節(jié)點(diǎn)慢度的系數(shù)需除以4。

直達(dá)波和反射波走時(shí)對界面單元兩端節(jié)點(diǎn)深度的導(dǎo)數(shù)式也可以統(tǒng)一為：

其中：xc和zc是直達(dá)波或反射波射線與界面單元交點(diǎn)的坐標(biāo)；xl，zl，xr，zr表示該 界 面 單 元左、右端 點(diǎn) 的 坐標(biāo)；sin和sout分別表示入射射線和出射射線所在地層單元的慢度；α和β分別是入射射線和出射射線與z軸的夾角

實(shí)際計(jì)算中，（4）式和（5）式中的方差均折算到權(quán)重系數(shù)。權(quán)重系數(shù)的設(shè)置應(yīng)基于模型預(yù)判和經(jīng)驗(yàn)值：ρtd可作為基準(zhǔn)值設(shè)置為1，ρtr一般不超過1；慢度平滑系數(shù)λs，1分為x軸系數(shù)和z軸系數(shù)，其x軸系數(shù)在各種模型中變化不大，但z軸系數(shù)可因模型層次疏密程度進(jìn)行寬幅變化；深度平滑系數(shù)λz，1一般設(shè)置為1左右即可得到良好效果；考慮到先驗(yàn)信息的可靠性問題，先驗(yàn)系數(shù)一般較小，在缺少先驗(yàn)信息的模型中甚至可以不設(shè)置先驗(yàn)方程，下文的模擬模型使用了測井先驗(yàn)信息，實(shí)例模型則未使用任何先驗(yàn)信息。在迭代反演過程中，應(yīng)使后四項(xiàng)系數(shù)的值隨著迭代次數(shù)逐次降低，以期在迭代初期獲得穩(wěn)定且有意義的收斂解，而在迭代后期獲得與走時(shí)數(shù)據(jù)更匹配的高分辨率解。在求解線性方程（3）時(shí)，本文選擇阻尼LSQR方法。

對于本文不規(guī)則單元的離散模型，使用基于波前擴(kuò)展和走時(shí)插值的三維射線追蹤方法［9］確定射線路徑，該方法適應(yīng)性強(qiáng)，計(jì)算精度和效率高。

2 模型實(shí)驗(yàn)

圖2是含有起伏反射界面的層狀介質(zhì)井間模型試驗(yàn)。圖2a是模型的三維尺寸，藍(lán)色線為炮點(diǎn)陣列，綠色線為檢波點(diǎn)陣列；圖2b是理論模型的xoz截面。模型尺寸為200m×10m×960m；含有5個(gè)勻速層，4個(gè)反射界面；炮點(diǎn)和檢波點(diǎn)各300個(gè)，分別均勻布置在33～930m深度范圍內(nèi)。層析反演時(shí)，選取了射線與水平面的夾角在15°～45°范圍內(nèi)共計(jì)15 544個(gè)直達(dá)波和10 972個(gè)反射波走時(shí)數(shù)據(jù)。初始界面取為左右井理論界面深度的連線，即平面或斜面，各層的初始速度均取為2.0km／s。地層離散單元起始尺寸為8m×5m×5m，且高度隨著迭代過程中反射界面的變化而變化。圖2c是直達(dá)波走時(shí)層析反演第六次迭代結(jié)果，圖2d是聯(lián)合層析反演的第六次迭代結(jié)果?？梢钥闯?，單直達(dá)波走時(shí)反演給出了較好的速度分布，但只能以速度過渡帶的形式反映界面，特別是2個(gè)斜截面和下凹界面的中段反演精度不高；聯(lián)合反演既給出了清晰的反射界面，也給出了精度更高的地層速度，整體結(jié)果與理論模型很接近。

3 實(shí)例

本文對勝利油田某區(qū)J41－108井的實(shí)際井間地震數(shù)據(jù)進(jìn)行了層析成像。炮點(diǎn)329個(gè)，深度從869～1 853m，檢波點(diǎn)300個(gè)，深度從866.1～1 763.1m，炮（檢）深度間距均為3m；經(jīng)坐標(biāo)軸旋轉(zhuǎn)后，炮點(diǎn)陣列在左井，且為直井，檢波點(diǎn)陣列在右井，為斜井；地震有效道數(shù)43 568個(gè)，拾取了36 693道直達(dá)波走時(shí)數(shù)據(jù)，并參考其他地質(zhì)資料，選擇了6個(gè)反射界面及其6 102個(gè)反射波走時(shí)數(shù)據(jù)，6個(gè)反射面分別為1 412、1 480、1 570、1 614、1 697和1 726m。通 過 射 線 出 射 角 限 定 后，有11 210個(gè)直達(dá)波和1 203個(gè)反射波走時(shí)數(shù)據(jù)參與了反演計(jì)算。圖3是J41－108第163炮的初始拾取結(jié)果，其中反射波皆為上行波。

圖2 含非水平界面模型合成走時(shí)數(shù)據(jù)的層析反演結(jié)果Fig.2 Tomographic inversion models using synthetic traveltimes of a crosshole model with nonlevel interfaces

圖3 勝利油田J41－108第163炮部分拾取結(jié)果Fig.3 Part of the picking result of 145th shot of 41－108wells in Shengli Oilfield

根據(jù)井斜資料，建立194m×15m×1 000m的模型，深度范圍從860～1 860m；設(shè)置6個(gè)水平初始反射面，深度如上文所述，各地層初始速度均為2.0km／s；初始離散單元為5m×5m×3m。圖4a為只用直達(dá)波走時(shí)反演的第6次迭代結(jié)果，圖4b為聯(lián)合層析反演第6次迭代結(jié)果，反演中沒有施加任何慢度和深度的先驗(yàn)信息約束。為了直觀表示反演結(jié)果的可靠程度，速度圖兩側(cè)附加了聲波測井曲線和邊界反演速度的對比，其中，青色為聲波測井速度，黑色為對應(yīng)的邊界反演速度。由于聲波測井曲線的記錄深度有限，為了與測井速度進(jìn)行對比，圖4只顯示出深度范圍為1 000～1 860m的反演速度結(jié)果。

由圖3可見，該井間反射記錄復(fù)雜，能夠清晰拾取的反射走時(shí)數(shù)據(jù)很少，參與反演計(jì)算的反射波走時(shí)數(shù)據(jù)比直達(dá)波走時(shí)數(shù)據(jù)少得多，使得圖4a和圖4b的結(jié)果整體上比較接近，但圖4b中聯(lián)合層析反演模型的速度分布更精細(xì)，垂向上速度變化層次更多，層間變化更劇烈，此現(xiàn)象在反射數(shù)據(jù)走時(shí)最集中的1 480、1 570和1 697m的3個(gè)界面附近最為明顯。兩種方案沿井徑的層析速度與測井聲波速度的變化趨勢也基本吻合，而聯(lián)合層析的速度曲線與測井聲波速度曲線的吻合度更好。實(shí)例結(jié)果說明，井間地震直達(dá)波和反射波聯(lián)合層析反演方法是可靠的，在層析成像中加入反射波走時(shí)，不僅能夠確定反射界面的起伏，而且能夠提高模型速度的分辨率。

圖4 41－108井實(shí)際資料層析反演結(jié)果Fig.4 Tomographic model section from field seismic data between wells 41and 108

4 結(jié)語

本文提出了1種基于不規(guī)則六面體（和五面體）單元的井間直達(dá)波和反射波走時(shí)聯(lián)合層析反演方法。建立的地層與界面完全耦合的三維井間離散模型更合理，可以很好的模擬連井截面內(nèi)反射界面起伏和每個(gè)地層內(nèi)的速度變化；使用的三維射線追蹤算法，能同時(shí)獲得井間三維空間的直達(dá)波和反射波射線路徑和走時(shí)數(shù)據(jù)；使用二維反演方程，大大減小了反演問題的不適定性。與井間直達(dá)波層析反演相比，井間直達(dá)波和反射波走時(shí)聯(lián)合層析反演，既可以確定反射界面，也提高了井間層析成像的分辨率和可靠性。

［1］ McMechan G A.Seismic tomography in boreholes［J］.Geophys J Roy Astr Soc，1983，74：601－612.

［2］ Calnan C，Schuster G T.Reflection ＋ transmission crosswell tomography［C］.［s.l.］：59th Annual International Meeting，SEG，Expanded Abstracts，1989，8：908－911.

［3］ Wu L，Song W，Zhang M.Seismic crosshole curved ray reflection＆transmission tomography（CCRTT）［J］.CT Theory and applications，1994，3（1）：42－47.

［4］ Van Schaack M A.Velocity estimation for crosswell reflection imaging using combined direct and reflected arrival traveltime tomography［D］.Ph.D.dissertation.USA：Stanford University：1997.

［5］ Spetzler J.Finite－frequency wavefield theory for high－resolution velocity estimation using transmission and reflection data［C］.［s.l.］：73th Annual International Meeting，SEG，Expanded Abstracts，2003，22：2315－2318.

［6］ Bube K P，Langan R T.Resolution of slowness and reflectors in crosswell tomography with transmission and reflection traveltimes［J］.Geophysics，2008，73（5）：321－335.

［7］ 左建軍，林松輝，孔慶豐，等.井間地震直達(dá)波和反射波聯(lián)合層析成像及應(yīng)用［J］.石油地球物理勘探，2011，46（2）：226－231.

［8］ 張建中，丁興號.一種2.5維井間透射波層析成像方法［J］.聲學(xué)學(xué)報(bào)，2007，32（1）：91－96.

［9］ Huang Y Q，Zhang J Z，Liu Q H.Three－dimensional GPR ray tracing based on wavefront expansion with irregular cells［J］.IEEE Transactions on Geoscience and Remote Sensing，2011，49（2）：676－686.

［10］ 張建中，陳世軍，徐初偉.動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤［J］.地球物理學(xué)報(bào)，2004，47（5）：899－904.