基于空間克里金插值的高精度多波匹配方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于地質(zhì)解釋領(lǐng)域中多波匹配問題，特別涉及一種基于空間克里金插值的 高精度多波匹配方法。
【背景技術(shù)】
[0002] 多波地震勘探是進行巖性油氣藏和隱蔽油氣藏勘探的一種非常有潛力的手段，但 是，由于諸多原因，多波多分量理論研究和油氣田實際勘探地質(zhì)需求的結(jié)合問題、復(fù)雜條件 下的轉(zhuǎn)換波地震資料處理問題和多波綜合解釋、全波屬性的地質(zhì)應(yīng)用等問題一直沒有取得 顯著進展，并且已經(jīng)成為制約多波地震勘探技術(shù)進一步發(fā)展的"瓶頸"。而解決這些問題 的基礎(chǔ)是做好多波多分量資料處理，提供高質(zhì)量的各向同性和各向異性處理成果。其中， 多波傳播機理的基礎(chǔ)研究、多波資料中的縱橫波匹配方法研究是目前多波地震資料后續(xù)處 理的重點和難點，是多波精確成像和疊前縱橫波聯(lián)合反演及巖性識別、儲層預(yù)測和含氣性 識別的重要基礎(chǔ)，是體現(xiàn)多波多分量地震勘探技術(shù)實際勘探開發(fā)應(yīng)用價值的關(guān)鍵。因此， 基于多波傳播機理，研究縱橫波高精度的匹配新方法，有利于充分利用多波多分量地震資 料、準確認識多波地質(zhì)響應(yīng)特征，突出多波多分量地震資料解決地質(zhì)問題的能力，具有重 大的意義。
[0003] 目前多波匹配有基于反射特征的匹配方法和基于多波層位的匹配方法，前者通過 PP波和PS波地震數(shù)據(jù)的波形和波組特征進行對比生成γ。值，然后基于該γ。值實現(xiàn)兩者 的時間域匹配。后者首先分別基于ΡΡ波和PS波地震數(shù)據(jù)追蹤解釋出對應(yīng)的層位，然后通 過標定對應(yīng)的層位產(chǎn)生時移體，最后將時移體應(yīng)用于PS地震數(shù)據(jù)，實現(xiàn)PS與ΡΡ地震數(shù)據(jù) 的匹配。
[0004] 隨著勘探目標要求的提高，多波匹配技術(shù)的研究越來越受到人們的重視，縱橫波 匹配技術(shù)已經(jīng)成為地球物理學的研究熱點。JamesE.G(1996)詳細介紹了縱橫波速度比 (γ。）的求取方法，并用最大相關(guān)法求取γ<：、平均γ<：、層間Yc等，使用VSP資料驗證從PP 波和PS波剖面中求取γ。，且γ?？梢杂枚滩ㄩL振幅反演。Wai-KinChan等（1997)利用常 數(shù)值多次試算法，在時間對數(shù)域?qū)崿F(xiàn)縱、橫波匹配。由于淺層和深層的γc值不同，因此 該匹配方法只能適應(yīng)于特定的目的層。RichardV.D和JamesG等（2001)利用最大相似 性原理，掃描PP波和PS波的γ。譜，拾取平均γ。值，將二者在時間域進行匹配。VanDok 和Kristiansen(2003)采用人工手動拾取強同相軸和其他明顯的地質(zhì)特征，實現(xiàn)了ΡΡ波和 PS波在時間域的粗略匹配。MichaelV.D等（2003)從ΡΡ波和PS波中求得縱橫波速度比 和泊松比，基于建立的橫波模型，在深度域?qū)崿F(xiàn)縱橫波匹配。就縱橫波匹配技術(shù)軟件現(xiàn)狀而 言，Transform軟件的匹配功能相對完善，提供的基于反射特征的匹配方法和基于多波層位 的匹配方法都能實現(xiàn)PP數(shù)據(jù)和PS數(shù)據(jù)的匹配。縱觀國內(nèi)外縱橫波匹配技術(shù)的發(fā)展，幾乎 所有的方法均是直接以PP波和PS波同一地層反射同相軸相似性最大的原則求取縱橫波的 速度比，直接對PP波和PS波在時間域進行匹配，其匹配程度是比較粗略的。

【發(fā)明內(nèi)容】

[0005] 本發(fā)明的目的在于克服現(xiàn)有技術(shù)，提供一種通過空間克里金插值算法對初始的 ga_a數(shù)據(jù)體進行插值，采用重采樣算法對PS數(shù)據(jù)進行初始匹配；然后通過控制層位點以 及劃定時窗的方式完成PS數(shù)據(jù)和PP數(shù)據(jù)的精細匹配，完成最終多波匹配，提高了多波匹配 的精度以及匹配的效率的基于空間克里金插值的高精度多波匹配方法。
[0006] 本發(fā)明的目的是通過以下技術(shù)方案來實現(xiàn)的：基于空間克里金插值的高精度多波 匹配方法，包括初始匹配和精細匹配兩個步驟，通過基于空間克里金插值的算法完成對原 始PS數(shù)據(jù)進行初始匹配，包括以下步驟：
[0007]S1、對最初的gamma數(shù)據(jù)體進行三角剖分，構(gòu)建三維Delaunay三角網(wǎng)；
[0008] S2、對于構(gòu)建的Delaunay三角網(wǎng)，通過空間索引的方法對待插值點進行空間定 位；
[0009] S3、對每一個待插值點進行克里金插值運算，計算出每個點的8曰_3值，根據(jù)數(shù)據(jù) 體內(nèi)所有點的ga_a值，采用重采樣算法對PS數(shù)據(jù)進行初始匹配；
[0010] 所述的精細匹配包括以下步驟：
[0011] S4、對初始匹配后的PS數(shù)據(jù)和原始PP數(shù)據(jù)，通過基于層位約束的精細匹配算法， 通過控制層位點以及劃定時窗的方式完成精細匹配。
[0012] 進一步地，所述的步驟S1中，以人機交互設(shè)置的數(shù)據(jù)點和自動追蹤算法計算出來 的數(shù)據(jù)點作為最初的gamma數(shù)據(jù)體；其實現(xiàn)方法具體包括以下子步驟：
[0013] SI1、根據(jù)給定的點集的區(qū)域構(gòu)建一個初始三角形，使該三角形能夠包含所有的數(shù) 據(jù)點，作為初始Delaunay三角網(wǎng)；
[0014]S12、從點集中選擇一個點進行插入操作，在已有的Delaunay三角網(wǎng)中找到包含 該點的三角形，將該點與包含它的三角形的頂點連接形成三個新的三角形；
[0015] S13、根據(jù)Deluanay三角網(wǎng)的外接空圓特性，利用Lawson提出的L0P算法對當前 Delaunay三角網(wǎng)進行優(yōu)化，使之滿足Delaunay特性；
[0016] S14、重復(fù)步驟S12和S13,當點集中所有點都包含在Delaunay三角網(wǎng)中時結(jié)束；
[0017] S15、將與初始三角形連接的所有三角形刪除。
[0018] 進一步地，所述的步驟S2包括以下子步驟：
[0019] S21、對于當前四面體，計算它的四個頂點的最大和最小Z值，進而計算出Z方向的 索引號范圍；
[0020] S22、按照Z值從小到大的順序，將Z固定，依次對該四面體進行切割，計算該Z平 面與四面體的交點；
[0021] S23、若求得交點個數(shù)為1，則該點定位完成，若交點個數(shù)為3或4,則按照Y值大小 將交點排序，交點構(gòu)成三角形或者四邊形；
[0022] S24、計算出Y方向的索引號范圍，按照Y值從小到大的順序，將Y值固定，依次對 三角形或四邊形進行切割，計算該平行于X軸的直線與三角形或四邊形的交點Xi，X2,則落 在Xi，X2范圍內(nèi)的點即為下一個待處理的四面體；
[0023]S25、對于下一個四面體重復(fù)步驟S21~S24,直到遍歷完所有四面體。
[0024] 進一步地，所述的步驟S3包括以下子步驟：
[0025]S31、進行克里金插值運算：
[0027] 其中區(qū)域化變量為Z(x)，待插值點為x。，樣本點記為Xi(i= 1，2,......，η)，η為 樣本點的個數(shù)，在點^處的屬性值記為Ζ(χJ，則點χ。處的屬性值Ζ(χ。）是由各個樣本點屬 性值加權(quán)求和得到的，其中，λ 待定權(quán)系數(shù)；
[0028]S32、在無偏估計和區(qū)域化變量滿足二階平穩(wěn)的條件下，為了使得到的估計方差值 達到最小，得到加權(quán)系數(shù)的方程組，即克里金方程組：
[0030] 在實際插值過程中，為了對方程組求解得到權(quán)重系數(shù)，將方程組轉(zhuǎn)換為矩陣的形 式進行表達，之后利用矩陣運算進行求解，便于計算，矩陣形式如式（3)所示：
[0031] [Κ]·[λ] = [Μ] (3)
[0032] 在式（3)中，Κ矩陣的表達式如下：
[0034] 其中，γ(νη，νη)表示第η個已知點與第η個已知點之間的變差函數(shù)值，在Κ矩陣 中，前η行η列的值是兩兩已知點之間的變差函數(shù)值；
[0035] Μ矩陣的表達式如下：
[0037] γ(νη，ν)表示第η個已知點與未知點之間的變差函數(shù)值，Μ矩陣的前η行表示待 插值點與每個已知點之間的變差函數(shù)值；
[0038]λ矩陣的表達式如下：
[0040]λ矩陣即權(quán)重系數(shù)矩陣，由λ= [Κ] + [Μ]計算得到，[Κ]+表示的是Κ矩陣的逆矩 陣；
[0041]S33、將λ矩陣帶入公式（1)中，計算點χ。處的屬性值Z(x。)，即為該點的gamma 值；
[0042]S34、重復(fù)步驟S31~S33,計算所有待插值點的gamma值；
[0043]S35、根據(jù)所有待插值點的gamma值，利用重采樣算法重新計算該道數(shù)據(jù)上每個點 的數(shù)據(jù)值，完成初始匹配。
[0044] 進一步地，所述的步驟S4具體包括以下子步驟：
[0045]S41、記錄每道PS數(shù)據(jù)的層位點坐標；
[0046]S42、抽取一道PP數(shù)據(jù)和對應(yīng)的一道初始匹配后的PS數(shù)據(jù)構(gòu)成一個N行2列的數(shù) 據(jù)面，其中N為每道PS數(shù)據(jù)上采樣點的個數(shù)；
[0047]S43、設(shè)定時窗參數(shù)：凡、乂、隊，其中，Nw為波形對比窗口的大小，N"為波形對比窗口 的移動量，隊為波形對比窗口的搜索半徑，Nw，N"，Ns均為整數(shù)；以Nw為窗口大小，以N"為窗 口移動量，從上到下，滑動窗口，以窗口中心點作為種子點，計算每個種子點移動量；具體計 算方法為：
[0048]設(shè)is= (k-1) *Nm+Nw/2,ie=iS+Nw，定義集合f= …,義丨丨，k為波形對比窗口 個數(shù)；Z={-Λζ，-Μ+1，··.，0,…，-1，Λζ}，
[0049] 定義長度為（Nw+1)的向量f，f的第i個元素/(/) =ΖΧ(十
[0050] 定義2NS+1個長度為（Nw+1)的向量gl，gl的第i個元素 g,(i) =D(i^+/+z-1, /2), /e?,Ι
[0051] 最優(yōu)移動量為Sa,j2)，其中jJPj2代表第jJPj2道；
[0052] 則S(j\，j2)的計算方法為：
[0054] 對于第k個波形對比窗口，計算拉平種子點移動量，定義長度為J的向量 m用來存儲每道PS數(shù)據(jù)拉平種子點移動量，所述的拉平種子點移動量為：m(j2)= sa，j2)-sa，j);
[0055]S44、增加對層位點