西部山前帶黃土礫石區靜校正方法

崔慶輝 尚新民 騰厚華 金昌昆 趙勝天 宋桂橋

(①中國石油化工股份有限公司勝利油田分公司物探研究院，山東東營257022；②中國石油化工股份有限公司油田勘探開發事業部，北京100728)

0 引言

中國西部復雜山前帶地區地表高程和近地表速度結構變化劇烈，給近地表建模帶來巨大挑戰，且嚴重影響激發井深設計、靜校正等[1-2]，也吸引許多學者針對該類地區特點開展各種近地表建模及靜校正方法研究[3-5]。馬晶晶等[6]利用折射靜校正的中、長波長靜校正優勢，在沙丘為主的準噶爾盆地山前帶取得較好效果。藍益軍[7]提出基于連續介質理論的空變時深曲線靜校正方法，并成功應用于近地表以沙丘和巨厚礫石為主的山前帶地區。方勇等[8]將菲涅耳層析靜校正技術應用于庫車地區，郭振波等[9]將回轉波層析靜校正應用于西部復雜地區，但層析靜校正不能解決激發井深設計問題，且需耗費大量時間拾取初至，不適用于地震數據快速處理。

本次研究區位于烏魯木齊市東北部、博格達山西北緣，區內總體地勢呈西北低、東南高，地形高差較大，最大可達約3000m。該區近地表主要為黃土、礫石二元結構，低速層為表層黃土，厚度較小，降速層為泥質礫石或黏土質礫石，厚度一般較大。這類近地表分布廣泛，且為該區各種近地表中最復雜的類型，構建精確的近地表黃土和礫石厚度模型及查清近地表地球物理特征，是解決該區近地表問題的關鍵。尚新民等[10]曾在該區聯合利用露頭、巖性錄井、微測井與淺層二維地震等信息建立了全區近地表模型，基于該模型逐點設計激發井深。

已有的研究成果為解決山前帶黃土礫石區靜校正問題奠定了基礎。本文基于山前帶黃土礫石區近地表特征研究，應用地質統計學方法建立了黃土層底界面和礫石層底界面，分別根據黃土層和礫石層的時深曲線計算靜校正量，從而實現了一種高效的山前帶黃土礫石區靜校正方法，為該類地區激發井深設計及靜校正提供了新思路。

1 方法原理

1.1 近地表調查方法

為了查清工區近地表巖性分布及速度結構，共布設了173口微測井和3條二維淺層地震測線。微測井采用井中激發、井口接收，以1km×1km 間隔均勻布設，鉆井過程中根據上返巖屑做巖性錄井，井深大于50m 的微測井比例超過25%。全區還部署了3條針對近地表淺層調查的地震測線，采用中間激發雙邊接收觀測方式，兩邊各有501 道，道距為2m，最大炮檢距為1000m，炮點距為6m，激發井深為2.5m。藉此二維測線層析反演以全面了解該工區近地表速度結構。

1.2 基于微測井的近地表特征研究

近地表特征主要包括近地表巖性分布、垂向速度變化規律、近地表介質與速度的關系、近地表介質時深對應關系及近地表介質空間分布特征等。

1.2.1 近地表巖性分布

沿南北向抽取一條線上的微測井(數據)，繪制出巖性剖面(圖1)。可見在工區北部和中部地表高程較低，黃土層和礫石層均較厚，為巨厚連續沉積；向南隨著地表高程增加，黃土、礫石厚度逐漸變薄，礫石分布逐漸消失。

圖1 近地表巖性探測南北向剖面

1.2.2 近地表速度特征

該區近地表速度分層性較差，常規微測井資料解釋方法存在較大人為誤差。利用微測井初至信息做層析反演可獲得可靠性更高的井筒速度模型[11]。

圖2為區內一口微測井初至、巖性剖面及層析反演的井周速度模型，該微測井初至層析反演的速度與巖性錄井數據較吻合。統計全區微測井層析反演結果后得知：表層黃土的速度在1000m/s以下，礫石層速度分布于1000～2000m/s 范圍，高于2000m/s的介質大多為基巖。據此可將近地表的速度模型轉化為巖性模型。

圖2 典型微測井初至、巖性剖面(a)及初至層析反演的速度曲線(b)

1.2.3 近地表時深關系

已有研究成果表明，與大沙漠區近地表相似，巨厚礫石區近地表呈連續介質特征，利用時深曲線靜校正方法可有效解決靜校正問題，尤其對于中、長波長靜校正具有一定的優勢[12，13]。

為了驗證該類方法是否同樣適用于山前帶黃土礫石區，進一步分析了該區近地表微測井時深關系。圖3a是全區所有微測井的時深散點圖，可見時深關系規律性差。分別按照0～1000m/s 和1000～2000m/s速度區間，將每一口微測井分成黃土層段(圖3b，深度0點為地表)和礫石層段(圖3c，深度0點為黃土層與礫石層的分界面)，在相同介質層內時深散點規律性較強，都呈現明顯的連續介質特征。可見黃土層和礫石層內部速度隨深度變化的規律性很強，這一結果為時深曲線靜校正方法的應用提供了依據，也提示近地表建模的關鍵是建立分層結構模型，亦即黃土層與礫石層的厚度模型。

為了從其他資料進一步佐證以上觀點，選取一條二維淺層地震測線的大炮初至進行層析反演得到近地表模型(圖4)。在該測線上均勻選取多個樣點，提取這些樣點的近地表速度曲線，將速度曲線轉為時深曲線，分別按照0～1000m/s 和1000～2000m/s速度區間提取數據，繪制時間—深度散點圖(圖5)。層析反演結果轉換的時深關系具有明顯的連續介質特征，與微測井結果基本一致。

1.2.4 黃土、礫石分布的地質統計特征

近地表建模是以微測井測量結果作為控制點，基于空間插值得到近地表黃土和礫石底界面模型。由于已知炮點和檢波點高程，巖性界面模型和厚度模型本質上是相同的，知道其中一個模型即可得到另一個模型。在地表起伏較大、近地表變化劇烈地區，巖性界面高程和厚度的插值結果可能會有很大差異，需從中選取最適合的。

圖3 不同層段范圍微測井時深散點圖(a)全井段； (b)黃土層段； (c)礫石層段

圖4 工區二維測線大炮初至層析反演結果

圖5 二維層析反演結果轉換的時深散點(a)0～1000m/s速度區間；(b)1000～2000m/s速度區間

所有空間插值算法都基于待插值數據的空間自相關，它表示變量在同一分布區內觀測數據之間潛在的相互依賴性，若變量空間不相關，則插值結果可靠性無法保證。變量的空間自相關一般用莫蘭(Moran)指數、P值和z得分三個參數表征[14]。莫蘭指數是一個歸一化的有理數，取值范圍是[-1.0，1.0]。該指數表征空間相關性大小： 大于0 時，表示數據呈現空間正相關，其值越大空間相關性越明顯；小于0時，表示數據呈空間負相關，其值越小空間差異越大；等于0 時，空間呈隨機性。P值表示觀測值呈隨機分布的概率，P值越小說明觀測值隨機分布的概率越小，如P=0.05 表示觀測值只有5%的可能性是隨機生成的結果。z得分表示觀測值標準差的倍數，z得分越高表示數據聚集性越強。

分別按0～1000m/s和1000～2000m/s速度區間提取每口微測井處黃土層、礫石層的厚度和界面高程，計算兩種介質層高程和厚度分布的莫蘭指數、P值和z得分(表1)。結合圖3 可知：山前帶近地表黃土、礫石層厚度和界面高程的莫蘭指數均大于0.2、P值均小于0.05、z得分均大于3.5，具有明顯的空間正相關和高聚集特征，故可通過空間插值建立全區模型。黃土、礫石層厚度的莫蘭指數和z得分均明顯大于界面高程的，P值明顯小于界面高程的(表1)，說明對厚度的空間插值結果更可靠。

表1 空間自相關分析結果

1.3 山前帶近地表黃土、礫石厚度建模

1.3.1 空間插值方法基本原理

空間插值方法主要有以反距離加權插值為代表的確定性插值和基于數據測量統計特征的地質統計學插值法。通過設定待測點周圍控制點的權重，進行加權求和，得到待測點的值

(1)

式中：Z(x0)為待測點x0的值；Z(xi)為控制點i的測量值；m為控制點個數；λ犻為控制點i的權系數

(2)

式中：d(x0,xi)為控制點犻到待測點的距離； α 為正實數，該值越大插值結果受最近點影響越大，一般取2。顯然，此算法簡捷、實用，當采樣點較密集且均勻分布時精度更高。

克里金插值的原理是利用區域化變量的已知測量數據及其變異函數，對待測點進行線性無偏最優估計，也可用式(1)表示，但其權系數的計算不同，除與距離有關外，還與區域化變量的方向變化有關。以無偏為前提，克里金方差為最小可得到求解待定權系數λi的方程組[12]

(3)

式中犆C(xi,xj)是Z(xi)與Z(xj)的協方差函數。

1.3.2 不同插值方法效果對比

分別利用反距離加權和普通克里金插值建立工區近地表黃土和礫石層厚度模型(圖6、圖7)。結果表明： 反距離加權插值法簡捷高效，但易產生“牛眼”現象；普通克里金插值法考慮了測量數據相關性在不同方向的差異，結果更平滑、自然。

圖6 黃土層厚度模型(a)反距離加權插值； (b)普通克里金插值

進一步對插值結果做交叉驗證，圖8中以普通克里金插值法所得黃土和礫石厚度誤差總體上均明顯小于反距離加權插值結果，克里金插值法所得黃土厚度誤差約為-5～5m、礫石厚度誤差大致為-10～10m，克里金插值法效果優于反距離加權插值法。

1.3.3 基于地表相似性的模型校正方法

研究區地形起伏較大，控制點相對稀少且呈均勻分布，起伏較大的局部區域無控制點或控制點密度不夠，降低了黃土和礫石層厚度插值結果的可靠性。根據已有的研究成果[15-17]，利用地表與近地表巖性界面之間的相似性可進一步提高建模的精度。根據給定半徑內的控制點信息計算出每個插值點巖性界面與地表的相似系數

(4)

式中：K(S,N)為地表高程S與界面(黃土層底或礫石層底)高程N的相似系數；Cov(S，N)為S與N的協方差；Var(·)為方差。

從利用式(4)計算的巖性界面與地表高程的相似系數(圖9a、圖9c)看，整體上工區南部黃土界面和礫石界面與地表高程相似系數較小，北部的較大，這是由于東南部屬于地表起伏劇烈的陡坡區，近地表物性不穩定，而北部屬于沖積扇，地表起伏較小，近地表沉積物規律性較好。

對圖6和圖7模型進行校正，得到校正后的模型(圖9b、圖9d)，再結合圖10 可知，本區黃土層厚度和礫石層厚度均呈現出由東南向西北部增加的趨勢，與工區地表地形和高程分布基本一致。進一步通過交叉驗證對校正效果進行評價

H′ =H+ (S-N)(1-K)

(5)

式中：H′為校正后厚度；H為插值后所得厚度；N為插值后所得界面高程。可見經過相似系數校正后，黃土和礫石層厚度誤差基本分布在-5～5m(圖11)，即誤差進一步減小，模型精度有所提高。

圖7 礫石層厚度模型(a)反距離加權插值； (b)普通克里金插值

圖8 交叉驗證對比(a)黃土厚度模型； (b)礫石層厚度模型

圖9 地表相似系數及校正后厚度模型

1.4 分層時深曲線靜校正

時深曲線靜校正主要用于大沙漠區，基本原理是對低降速帶內的微測井時深散點進行擬合，得到一條自變量為深度、因變量為時間的時深關系曲線，再根據建立的低降速帶厚度模型，利用時深關系公式求取地震波在低降速帶的傳播時間，進一步計算靜校正量[7，12，18]。與其他靜校正方法相比，此方法不需反演詳細的近地表速度結構，也不需拾取大炮初至信息，在近地表滿足連續介質假設條件下，能高效地解決靜校正問題，特別是對中、長波長靜校正具有一定優勢，在大沙漠和巨厚礫石區應用較多。

圖10 校正后模型的二維剖面(a)東西向5000m 處南北向剖面； (b)南北向10000m 處東西向剖面

圖11 校正后黃土層(a)和礫石層(b)厚度模型的交叉驗證

據該區微測井資料，時深散點在黃土層內和礫石層內呈現明顯的規律性(圖3)，在此提出分層時深曲線靜校正方法。該方法采取“兩步法”計算：地表到高速層頂界的時間根據黃土厚度和礫石厚度分別用黃土層和礫石層時深曲線計算后相加得到，此為“剝去量”；再通過替換速度計算高速層頂界到基準面的時間，此為“充填量”。任意炮點或檢波點靜校正量依據不同高程/厚度從下式選擇計算

(6)

式中：f1、f2分別為待計算點黃土層和礫石層時深曲線；h1、h2分別為該點黃土層和礫石層厚度；EstationEsurface、Ebase、Ehvlt分別為待計算點高程、地表高程、基準面高程、高速層頂高程；sdepth為待計算點深度；vH、vR分別為待計算點高速層速度和替換速度。檢波點一般位于地表，采用式(6)中的第一種情況進行計算。

2 實際應用效果

針對圖3中黃土層和礫石層段時深散點，用最小二乘法擬合出黃土層和礫石層時深關系分別為

f1(x)= -0.0459x2+2.6958x

(7)

f2(x)= -0.0012x2+0.6607x

(8)

利用式(6)計算校正后厚度模型(圖9b、圖9d)的靜校正量(圖12c)，與高程靜校正、某商業軟件的折射靜校正結果進行對比(圖12)，可見不同方法計算的靜校正量的總體趨勢較一致。對比沿圖12 中藍色直線抽取的一條檢波線上不同方法的靜校正量(圖13)，看到本文方法與折射靜校正法計算的靜校正量總體分布一致，但在局部有明顯差別。

進一步對比疊加剖面(圖14)效果，且(三種方法)采用同一疊加速度以消除其影響。可見本文方法明顯優于高程靜校正，總體上也優于折射波靜校正，同時在計算效率上具有明顯優勢，證明本文的近地表建模及靜校正方法的正確性。

3 結論與討論

通過對中國西部山前帶近地表黃土、礫石結構展開系統研究，得出以下認識和結論。

圖12 不同靜校正方法計算的靜校正量平面圖(a)高程靜校正； (b)折射靜校正； (c)本文方法

圖13 不同靜校正方法的檢波點靜校正量對比

圖14 不同靜校正方法局部疊加剖面對比(a)未做靜校正 (b)高程靜校正 (c)折射靜校正 (d)本文方法

(1)研究區地表黃土層速度在1000m/s以下，礫石層速度為1000～2000m/s，山前帶近地表黃土層和礫石層內時深關系具有較好規律性，近似為連續介質，可以用時深關系曲線描述。

(2)黃土層和礫石層厚度分布均具有較強空間相關性，可通過對微測井測量的黃土層和礫石層厚度做空間插值建立整個工區厚度模型；普通克里金插值效果好于反距離加權插值法，利用介質界面與地表面的相似性調整厚度模型可進一步提高精度。

(3)本文的近地表建模與分層時深關系曲線靜校正方法具有較高計算效率，避免了基于初至的靜校正方法需耗費大量時間做初至拾取，更適合當前高密度地震勘探現場處理的需求，同時建立的巖性界面還可用于激發井深設計。

(4)受限于施工成本及難度，本區僅有幾口微測井井深超過100m，在礫石層較厚區域未能探測到礫石底界位置。建議在該區盡量采用更多數量的超深微測井，以更準確地落實礫石底界分布。條件允許情況下，也可采用大排列的折射法加以補充，這樣將使本文方法取得更好效果。