金昌昆 王延光 尚新民 崔慶輝 王東凱

(①中國(guó)石化勝利油田物探研究院，山東東營(yíng) 257022； ②中國(guó)石化勝利油田博士后科研工作站，山東東營(yíng) 257002； ③中國(guó)石化勝利油田分公司，山東東營(yíng) 257000)

0 引言

近地表問題一直是陸上地震勘探的瓶頸。復(fù)雜近地表會(huì)嚴(yán)重影響地震波的激發(fā)和接收，并容易產(chǎn)生散射、次生干擾等噪聲，淹沒有效信號(hào)。近地表速度變化容易引起反射波同相軸畸變、振幅畸變，嚴(yán)重影響地震成像[1-3]。因此，研究近地表結(jié)構(gòu)，尤其是近地表速度，一直是勘探領(lǐng)域的熱點(diǎn)和難點(diǎn)。

近地表速度信息在地震數(shù)據(jù)采集前可通過近地表調(diào)查觀測(cè)，在數(shù)據(jù)采集后可應(yīng)用折射波解釋[4-6]或初至層析反演[7-9]求取?，F(xiàn)今常用近地表調(diào)查方法[10]有微測(cè)井調(diào)查、小折射法、淺層反射法、面波頻散法[11-13]、地質(zhì)雷達(dá)[14-15]、重磁電勘探[16]等。其中前兩者最為常用，這兩種方法各有特點(diǎn)，應(yīng)用條件也有所不同。

小折射法主要適用于地形平坦、變化平緩的層狀介質(zhì)地區(qū)。而在地表地形起伏劇烈，低降速帶速度和厚度橫向變化大，高速層頂面不穩(wěn)定，甚至速度發(fā)生倒轉(zhuǎn)的復(fù)雜地區(qū)，小折射法往往難以適用。與小折射法不同，微測(cè)井的觀測(cè)方式不受地形變化限制，施工中直接鉆穿低降速帶進(jìn)行激發(fā)，接收的地震波信噪比較高，便于初至拾取。顯然，微測(cè)井是一種相對(duì)成熟、可靠的近地表調(diào)查方法，即使是在巨厚礫石區(qū)，也能取得較滿意成果。該方法的不足之處是成本較高，施工效率較低，通常僅在有限、離散的點(diǎn)展開[17]。若要得到完整的近地表模型，則要對(duì)離散數(shù)據(jù)進(jìn)行內(nèi)插，插值結(jié)果的精度取決于離散點(diǎn)的數(shù)量、分布和插值方法。插值的表層速度模型對(duì)地震數(shù)據(jù)現(xiàn)場(chǎng)采集、后期數(shù)據(jù)處理中的靜校正、振幅補(bǔ)償、深度成像[18]等具有十分重要的意義。

為構(gòu)建精細(xì)近地表速度模型，本文提出一種基于微測(cè)井層析的方位加權(quán)插值建模方法，并在實(shí)際應(yīng)用中與常規(guī)方法做了對(duì)比，結(jié)果令人滿意。

1 微測(cè)井調(diào)查及插值建模簡(jiǎn)介

對(duì)于山前帶等復(fù)雜地區(qū)的地震勘探，確定表層速度、厚度和最佳激發(fā)參數(shù)是十分重要的課題，而微測(cè)井調(diào)查發(fā)揮了舉足輕重的作用。常規(guī)微測(cè)井調(diào)查首先鉆穿低降速帶，并從井底到井口依次布設(shè)激發(fā)震源，再由地面檢波器接收透射波信號(hào)； 然后通過拾取初至走時(shí)并做垂直校正，得到垂向走時(shí)信息，繪制時(shí)深曲線。在地表速度呈層狀分布的情況下，垂向走時(shí)分布具有一定規(guī)律性，同一地層內(nèi)的垂向走時(shí)呈線性分布，而直線的斜率反映了地層慢度。根據(jù)這一規(guī)律，將時(shí)深曲線劃分成不同地層并做擬合，計(jì)算出各層速度和厚度[19]。該解釋過程可有效刻畫局部速度結(jié)構(gòu)，即使是在速度反轉(zhuǎn)的情況下。但層位的劃分通常依賴于處理人員的視覺和經(jīng)驗(yàn)，參數(shù)選取的正確性和可靠性受人為因素影響，難以達(dá)到絕對(duì)準(zhǔn)確、客觀，尤其是在表層速度連續(xù)變化情況下，無法刻畫層內(nèi)速度變化，更易引入人為誤差。

微測(cè)井調(diào)查獲得的只是各點(diǎn)的垂向速度分布，若要描述完整的近地表速度分布，還必須進(jìn)行數(shù)據(jù)內(nèi)插。構(gòu)建表層速度模型的常用方法有反距離加權(quán)法、徑向基函數(shù)法、克里金插值法等[20]。反距離加權(quán)法效率高，能便捷地?cái)U(kuò)展到各向異性和高維情況； 但該方法對(duì)數(shù)據(jù)的分布較敏感，易出現(xiàn)“牛眼”效應(yīng)。徑向基函數(shù)法根據(jù)確定的基函數(shù)，構(gòu)建線性方程組并求解加權(quán)系數(shù)。該方法插值效果好，但隨著數(shù)據(jù)量增大，其計(jì)算量和儲(chǔ)存量也顯著增加，且求解過程可能出現(xiàn)不穩(wěn)定問題?？死锝鸩逯凳悄壳白畛Ｓ玫牟逯捣椒ㄖ?，該方法統(tǒng)計(jì)分析變差函數(shù)，并利用最小二乘算法求解插值系數(shù)。在插值過程中，采用的變差函數(shù)對(duì)插值結(jié)果的影響很大[18]。實(shí)際應(yīng)用中，微測(cè)井觀測(cè)點(diǎn)之間有較大間隔，實(shí)測(cè)數(shù)據(jù)欠充分均勻[21]，直接擬合得到的變差函數(shù)可靠性較差。

為克服上述問題，通?；跀?shù)據(jù)統(tǒng)計(jì)特征和影響因素，選取適當(dāng)?shù)淖儾詈瘮?shù)模型，再擬合確定其變程等參數(shù)。常用的變差函數(shù)模型有：高斯模型、指數(shù)模型、球形模型、冪函數(shù)模型等。若要描述更復(fù)雜的空間變異性，則需考慮尺度效應(yīng)和各向異性情況，適宜的變差函數(shù)理論模型可描述不同尺度、不同方向的空間變異性。這需要大量的統(tǒng)計(jì)分析工作作為支撐，實(shí)際應(yīng)用難度較大。此外，數(shù)據(jù)量較大時(shí)，求解克里金方程組也將變得困難。

2 微測(cè)井層析成像與方位加權(quán)插值

2.1 微測(cè)井層析成像原理

與微測(cè)井調(diào)查方法相同，初至波層析反演也是獲取近地表速度分布的常用手段。它不依賴于層狀模型假設(shè)，適用于復(fù)雜近地表情形； 對(duì)初始模型依賴性小，穩(wěn)定性高，所得結(jié)果為地下各點(diǎn)的速度分布。微測(cè)井層析成像是一種將微測(cè)井的垂向走時(shí)、解釋結(jié)果(地層厚度和速度)與初至波走時(shí)層析原理相結(jié)合，反演井壁速度分布的方法。

在圖1中，假設(shè)微測(cè)井附近的地層速度橫向變化可忽略，地震射線近似垂向傳播，則可將實(shí)際三維問題簡(jiǎn)化為一維問題。以垂直校正后的初至走時(shí)作為觀測(cè)走時(shí)T； 根據(jù)激發(fā)點(diǎn)、接收點(diǎn)的相對(duì)位置，計(jì)算垂向射線路徑，并以各網(wǎng)格中的路徑長(zhǎng)度為元素，構(gòu)建敏感度矩陣A； 設(shè)定井壁各點(diǎn)慢度值未知量為s，將微測(cè)井層位解釋速度信息網(wǎng)格化，并以其倒數(shù)作為先驗(yàn)慢度值sH，則微測(cè)井層析的目標(biāo)函數(shù)可表示為

(1)

式中：L為由一維二階拉普拉斯算子構(gòu)成的平滑矩陣；ε1為平滑約束權(quán)重；ε2為先驗(yàn)約束權(quán)重。引入單位矩陣I，與目標(biāo)函數(shù)(式(1))對(duì)應(yīng)的方程組為

(2)

圖1 微測(cè)井層析成像示意圖

應(yīng)用高效、穩(wěn)定的聯(lián)合迭代重建算法(SIRT)求解式(2)，可得反演問題的最小二乘解?；诼戎蹬c速度值的倒數(shù)關(guān)系，慢度解可轉(zhuǎn)化為井壁速度值。與常規(guī)初至走時(shí)層析反演不同，本方法射線路徑及敏感度矩陣不隨速度變化而改變，反演變量為井壁慢度值，而非慢度更新量。

以上處理將原問題轉(zhuǎn)化為對(duì)線性反演問題求解，無須多次迭代，計(jì)算效率和反演穩(wěn)定性高，可在微測(cè)井觀測(cè)中進(jìn)行現(xiàn)場(chǎng)處理。

2.2 方位加權(quán)插值

前已述及，微測(cè)井層析只能獲得觀測(cè)點(diǎn)處垂向速度分布，若要取得整個(gè)工區(qū)精細(xì)近地表速度場(chǎng)，就必須采用高精度橫向插值方法。因此，本文提出一種基于微測(cè)井層析和層位信息并考慮方位因素的局部插值方法，以獲得符合地表層趨勢(shì)的精細(xì)模型。該方法主要涉及三個(gè)方面：基于層位的深度變換、構(gòu)建基函數(shù)及求解權(quán)重方程組。

2.2.1 基于層位的深度變換

常規(guī)插值方法只考慮距離因素，且結(jié)果帶有明顯的平滑效應(yīng)，而近地表往往是速度變化最劇烈的區(qū)域，橫縱向變化十分明顯。在地表層位起伏情況下，做橫向插值所用數(shù)據(jù)可能來自不同速度的地層，若不做區(qū)分直接進(jìn)行插值，其平滑效應(yīng)會(huì)使插值結(jié)果產(chǎn)生橫向拉伸現(xiàn)象。

為避免此情形，本文在插值前對(duì)微測(cè)井?dāng)?shù)據(jù)進(jìn)行基于層位信息的深度變換，將沿起伏層位的數(shù)據(jù)轉(zhuǎn)化為沿水平分布的數(shù)據(jù)，橫向插值采用來自同一地層、同一深度的數(shù)據(jù)，插值完成后將沿水平分布數(shù)據(jù)反變換為沿起伏層位數(shù)據(jù)，可減少拉伸效應(yīng)，并使插值結(jié)果符合層位趨勢(shì)。有如下具體步驟。

(1)基于微測(cè)井解釋的地層信息，通過反距離加權(quán)計(jì)算得到地下界面。

(2)以各界面的平均深度為基準(zhǔn)，對(duì)微測(cè)井速度所對(duì)應(yīng)的深度進(jìn)行標(biāo)準(zhǔn)化，變換公式為

(3)

2.2.2 構(gòu)建基函數(shù)

數(shù)據(jù)的空間相關(guān)性是空間插值研究的基礎(chǔ)。Tobler定理認(rèn)為空間上距離較近的數(shù)據(jù)比距離較遠(yuǎn)的數(shù)據(jù)具有更大的相關(guān)性[22]，反距離加權(quán)、克里金插值等方法均服從這一定理。另外，數(shù)據(jù)在不同方向上可呈現(xiàn)出不同的變化，所以數(shù)據(jù)的空間相關(guān)性也具有各向異性特征[23]。因此，不同方向、相同距離上的數(shù)據(jù)對(duì)插值點(diǎn)的權(quán)重也應(yīng)有所不同。為表征不同距離、不同方向的權(quán)重，本文構(gòu)建的權(quán)重基函數(shù)包括徑向基函數(shù)和方位基函數(shù)。

徑向基函數(shù)：與全局徑向基函數(shù)法不同，本文方法為一種局部插值方法，將凸起的余弦函數(shù)(漢寧窗)作為徑向基函數(shù)，以降低計(jì)算量和存儲(chǔ)量，并保證插值精度?；瘮?shù)的具體形式[24]為

(4)

方位基函數(shù)：本文以八個(gè)主方位系數(shù)(正北、東北、正東、東南、正南、西南、正西、西北)加權(quán)結(jié)果表示各個(gè)方位權(quán)重。每個(gè)主方位的加權(quán)函數(shù)為二階三角B樣條函數(shù)，其分段表示形式[25]如下

(5)

(6)

圖2b顯示八個(gè)主方位的B樣條函數(shù)曲線在極坐標(biāo)系下的分布，各主方位系數(shù)[0.8，0.8，0.8，0.8，0.4，0.4，0.8，0.8]的插值曲線(圖2c)展示出相應(yīng)的各向異性特征。

式(6)中B(θ)無法給出(θ，r=0)處的權(quán)重值。為此，將插值函數(shù)B(θ)進(jìn)一步修改，作為方位基函數(shù)

(7)

圖2 方位基函數(shù)顯示圖

(a)主方位B樣條曲線； (b)極坐標(biāo)系中八個(gè)主方位B樣條曲線分布； (c)方位系數(shù)[0.8，0.8，0.8，0.8，0.4，0.4，0.8，0.8]在全方位上的插值結(jié)果(方位順序依次為正北、東北、正東、東南、正南、西南、正西、西北)； (d)與圖c對(duì)應(yīng)的方位基函數(shù)在定義域中分布顯示(c=0.1)

基于式(4)與式(7)，本文方法采用兩者乘積的形式作為權(quán)重基函數(shù)

(8)

通過以上方法構(gòu)建的權(quán)重基函數(shù)具有緊支撐性和各向異性特征。

2.2.3 求解權(quán)重方程組

(9)

式中：vi為第i口井速度；W(θij，rij)為權(quán)重基函數(shù)；B(θij)為基于方位角θij的插值函數(shù)；θij、rij、Rij分別為第i口井到第j口井的方位角、距離參數(shù)和徑向權(quán)重；wik為第i口井的第k個(gè)主方位系數(shù)。

將式(6)代入式(9)得

(10)

基于式(10)，建立如下插值目標(biāo)函數(shù)

(11)

式中：敏感度矩陣G的元素為由式(10)中與主方位系數(shù)相乘的系數(shù)值； 微測(cè)井的主方位系數(shù)w為未知量；V為微測(cè)井層析速度； 平滑矩陣L由微測(cè)井相鄰主方位系數(shù)之間的平滑約束項(xiàng)構(gòu)成；ε為平滑權(quán)重系數(shù)。與式(11)對(duì)應(yīng)的反演方程組為

(12)

應(yīng)用最小平方QR分解算法(LSQR)求解上述方程組，得到各微測(cè)井的主方位系數(shù)。與各向異性克里金插值統(tǒng)計(jì)區(qū)域各向異性不同，本文方法通過反演計(jì)算每口微測(cè)井相關(guān)性的各向異性特征。基于求解的主方位系數(shù)，插值點(diǎn)p的計(jì)算結(jié)果可表示為

(13)

3 實(shí)際應(yīng)用

3.1 準(zhǔn)南緣地區(qū)激發(fā)井深設(shè)計(jì)

中國(guó)西部探區(qū)大多具有地表及地下地質(zhì)構(gòu)造“雙重”復(fù)雜介質(zhì)特征[26]。準(zhǔn)南緣米泉區(qū)塊位于博格達(dá)山西北緣，地勢(shì)西北低、東南高，地表高差達(dá)3km，為典型的山前帶地貌[27-28]。工區(qū)北部低緩丘陵區(qū)分布有黃土、巨厚礫石層，南部山區(qū)出露二疊系灰質(zhì)細(xì)砂巖或石炭系灰?guī)r，中部高大丘陵區(qū)為侏羅系細(xì)砂巖和煤層。此次近地表調(diào)查采用了微測(cè)井、淺層地震等多種技術(shù)手段，其中微測(cè)井調(diào)查共173口，有5口為100m以上深井微測(cè)井，用以查明巨厚礫石層的底界。

圖3展示了該地區(qū)微測(cè)井的解釋結(jié)果和層析結(jié)果，包括微測(cè)井的時(shí)深曲線及層位解釋(圖3a)、微測(cè)井層析成像反演結(jié)果(圖3b)。對(duì)比兩者可見，微測(cè)井解釋的斜率與時(shí)深曲線并不完全吻合，而層析結(jié)果卻能清晰地展示地下速度變化細(xì)節(jié)，包括地表速度反轉(zhuǎn)。

三維地震施工前，為確定不同地表?xiàng)l件下激發(fā)井深等設(shè)計(jì)參數(shù)，在工區(qū)內(nèi)設(shè)計(jì)了自西北到東南的4個(gè)激發(fā)試驗(yàn)點(diǎn)。其中，S1試驗(yàn)點(diǎn)位于黃土礫石區(qū)，鉆井深度達(dá)100m，未鉆穿礫石層。圖4a為S1點(diǎn)的微測(cè)井層析速度曲線，可見井壁速度整體在2.5km/s以下。S2試驗(yàn)點(diǎn)位于工區(qū)中北部，鉆井井深為26m，從該點(diǎn)的微測(cè)井層析速度曲線(圖4b)可見：在井深20m之下為3～4km/s的高速層，實(shí)際鉆遇地層為侏羅系細(xì)砂巖。S3試驗(yàn)點(diǎn)位于工區(qū)中部，鉆井井深為50m，從該點(diǎn)的微測(cè)井層析速度曲線(圖4c)可見，深部高速層為侏羅系細(xì)砂巖，淺部低速為表層土、黃土及砂礫層。S4試驗(yàn)點(diǎn)位于工區(qū)中南部，鉆井井深為22m，從該點(diǎn)的微測(cè)井層析速度曲線(圖4d)可知：淺層速度低，主要成分為表層土和被風(fēng)化的松散巖層； 深層速度在3～4km/s之間，為二疊系灰質(zhì)細(xì)砂巖。米泉勘探通過各試驗(yàn)點(diǎn)的微測(cè)井調(diào)查以及大量的激發(fā)接收試驗(yàn)，明確了激發(fā)井深與激發(fā)巖性、激發(fā)速度的關(guān)系，確定了不同地表巖性條件下的最佳激發(fā)參數(shù)和激發(fā)原則，進(jìn)而完成了逐點(diǎn)激發(fā)井深設(shè)計(jì)任務(wù)[29]。

圖3 米泉微測(cè)井解釋結(jié)果(a)與層析結(jié)果(b)對(duì)比

圖4 試驗(yàn)點(diǎn)的微測(cè)井層析反演結(jié)果

3.2 陳官莊地區(qū)精細(xì)表層建模

陳官莊地區(qū)位于濟(jì)陽(yáng)坳陷，地表為第四系沉積物，潛水面淺，低降速帶穩(wěn)定，變化平緩，但仍存在異常區(qū)域。該區(qū)共做微測(cè)井129口，平面分布如圖5a所示，微測(cè)井間距約為2km，鉆井深度不超過30m，井壁速度不超過2km/s，鉆遇地層以膠泥層為主，微測(cè)井解釋結(jié)果(圖5b)和層析速度曲線(圖5c)均具有清晰的層位特征，層內(nèi)速度較穩(wěn)定。

根據(jù)速度曲線特征，工區(qū)近地表可劃分為三層：第一層速度在0.6km/s 以下；第二層速度約為1km/s；第三層速度在1.4km/s以上?；谖y(cè)井的解釋層位，插值得到全區(qū)層位分布(圖6a)所示。對(duì)微測(cè)井?dāng)?shù)據(jù)做深度轉(zhuǎn)換，并以4km作為最大井間距做方位加權(quán)插值，獲得最終反變換后結(jié)果(圖6b)，可見該插值速度與層位趨勢(shì)(圖6a)相符合。

作為對(duì)比，采用普通克里金法進(jìn)行插值。由于微測(cè)井分布稀疏，計(jì)算變差函數(shù)時(shí)采用各向同性模型表示。通過分析、擬合，以趨勢(shì)相符的線性模型作為理論變差函數(shù)模型并得到變程、基臺(tái)值、塊金常數(shù)等參數(shù)。之后基于無偏性和最小二乘原理，利用拉格朗日乘法構(gòu)造目標(biāo)函數(shù)，并求導(dǎo)得到克里金方程組，求解得到插值權(quán)重，進(jìn)而得到插值模型。

圖7為兩種插值結(jié)果的速度切片(深度為6m)?？傮w而言，兩種方法的插值效果基本相當(dāng)； 但對(duì)比黑色方框部分可知，本文方法結(jié)果(圖7b)比普通克里金法(圖7a)具有更高的橫向分辨率。

為測(cè)試插值結(jié)果的準(zhǔn)確性，以橫向距離為14.09km、縱向距離為6.113km處的微測(cè)井?dāng)?shù)據(jù)作為驗(yàn)證數(shù)據(jù)，得到如圖8a所示的解釋結(jié)果。基于其他微測(cè)井的層位及速度信息，應(yīng)用本文方法和普通克里金法對(duì)待驗(yàn)證位置做插值處理。從插值結(jié)果與驗(yàn)證數(shù)據(jù)的對(duì)比(圖8b)可見：當(dāng)深度小于10m時(shí)，兩種結(jié)果與驗(yàn)證數(shù)據(jù)具有較高的吻合度，而在2.16m和8.34m處速度界面附近，本文方法結(jié)果的縱向分辨率高于克里金插值結(jié)果，更符合驗(yàn)證數(shù)據(jù)； 當(dāng)深度大于10m時(shí)，兩種插值結(jié)果的變化平緩，雖未體現(xiàn)出驗(yàn)證數(shù)據(jù)的局部異常，但在趨勢(shì)上與驗(yàn)證數(shù)據(jù)一致，證明了插值結(jié)果的準(zhǔn)確性。

圖5 陳官莊地區(qū)微測(cè)井調(diào)查情況

圖6 陳官莊地區(qū)近地表層位分布(a)及插值結(jié)果(b)

圖7 普通克里金法(a)與本文方法(b)的速度切片對(duì)比

圖8 驗(yàn)證數(shù)據(jù)與插值結(jié)果顯示圖

4 結(jié)束語

本文提出一種基于微測(cè)井層析的方位加權(quán)插值建模方法，實(shí)際數(shù)據(jù)測(cè)試結(jié)果表明：

(1)與常規(guī)微測(cè)井解釋相比，文中的微測(cè)井層析成像可獲得更精細(xì)、連續(xù)的井壁速度分布；

(2)該方法反演微測(cè)井的主方位權(quán)重，無須人工進(jìn)行各方位的數(shù)據(jù)統(tǒng)計(jì)與分析，方便快捷； 所得速度數(shù)據(jù)符合層位趨勢(shì)，且縱橫向分辨率局部?jī)?yōu)于普通克里金插值的結(jié)果。