楊金璇 潘 懋 劉鈺洋 雷征東 吳耕宇

（1.北京大學(xué)地球與空間科學(xué)學(xué)院造山帶與地殼演化教育部重點實驗室 北京 100871）

（2.北京大學(xué)石油與天然氣研究中心 北京 100871）（3.中國石油勘探開發(fā)研究院 北京 100083）

1 引言

近幾年，三維地質(zhì)屬性建模中的確定性建模被廣泛應(yīng)用在油氣勘探、開發(fā)等領(lǐng)域。確定性建模由儲層地震學(xué)方法、儲層沉積學(xué)方法和地質(zhì)統(tǒng)計學(xué)克里金插值方法組成［1］。地質(zhì)統(tǒng)計學(xué)方法根據(jù)待估點周圍的若干已知信息，應(yīng)用變異函數(shù)對待估點未知值做出無偏最優(yōu)估計。目前廣泛應(yīng)用于空間域或時間域自然變量的定量化描述等眾多領(lǐng)域［1］。簡單克里金和普通克里金應(yīng)用最廣泛［2］，但先驗條件嚴(yán)苛，大部分地質(zhì)數(shù)據(jù)都難以滿足，因此這兩種方法的插值效果難以保證。

泛克里金插值是處理隨機變量非平穩(wěn)插值問題的方法，結(jié)合地質(zhì)變量在一定范圍內(nèi)的變化規(guī)律進(jìn)行建模，并且沒有嚴(yán)苛先驗條件，適用于絕大部分非均質(zhì)性比較強的儲層。楊功流［5］等對泛克里金插值方法和普通克里金插值方法進(jìn)行對比研究，一方面在非平穩(wěn)的地磁場數(shù)據(jù)的插值處理中應(yīng)用兩種方法，另一方面對比兩種方法的插值效果。實例驗證結(jié)果顯示，使用泛克里金插值得到的地磁圖更加符合地磁場的特征，證明對非平穩(wěn)地質(zhì)數(shù)據(jù)插值問題，泛克里金方法更適用。

然而，在廣泛的實際應(yīng)用中，泛克里金方法的理論優(yōu)勢并不明顯。相較于其他簡單快速算法，其精度優(yōu)勢比較微弱。目前，對于泛克里金插值精度的研究已成為熱點。李龍［6］結(jié)合回歸分析，應(yīng)用泛克里金插值方法對大氣參數(shù)進(jìn)行插值，獲得更高精度的AOD融合效果。王長虹等［8］在巖土參數(shù)隨機場分析中，將泛克里金插值方法與多重分型理論相結(jié)合，度量局部空間的奇異性。趙愛梅等［11］研究泛克里金插值效果對變異函數(shù)的敏感性。何濤等［12］基于泛克里金插值的時間復(fù)雜度，研究趨勢函數(shù)參數(shù)的估計。

本文考慮，對于泛克里金方法的使用前提是樣品點數(shù)據(jù)滿足正態(tài)分布，但是對于泛克里金方法來說，其插值結(jié)果依賴于反應(yīng)地質(zhì)變化規(guī)律的趨勢函數(shù)，而滿足正態(tài)分布不能保證研究區(qū)數(shù)據(jù)具有一定變化規(guī)律，即無法保證趨勢函數(shù)的精度。本文對于趨勢函數(shù)的模擬建立在殘差分析的基礎(chǔ)上，對樣本數(shù)據(jù)進(jìn)行重采樣，提高插值樣品數(shù)據(jù)的相關(guān)性及結(jié)構(gòu)性，切合泛克里金方法的本質(zhì)。

趨勢函數(shù)是泛克里金插值方法中用于表達(dá)地質(zhì)變量變化趨勢的函數(shù)［4］，趨勢函數(shù)的精度直接影響插值結(jié)果的空間性和隨機性，甚至插值結(jié)果的有效性，因此趨勢函數(shù)的擬合是泛克里金插值過程中最重要的問題之一。傳統(tǒng)的趨勢函數(shù)擬合建立在簡單的多元線性函數(shù)擬合的基礎(chǔ)上，本文提出，在傳統(tǒng)擬合趨勢函數(shù)的基礎(chǔ)上，應(yīng)用殘差分析，將在95%置信度下判別為異常點的實驗點剔除，將提升趨勢函數(shù)的擬合精度。

2 數(shù)學(xué)原理介紹

2.1 泛克里金方法概述

泛克里金方法是一種線性無偏最優(yōu)的空間插值計算方法，研究的非平穩(wěn)區(qū)域化變量，包括趨勢與殘差兩個部分，其數(shù)學(xué)表達(dá)如下：

式中，m(x)是趨勢函數(shù)，R(x)是殘差。趨勢函數(shù)表示區(qū)域化變量在研究區(qū)內(nèi)的某種明確的變化規(guī)律，是隨機變量Z(x)的期望，有：

殘差是表示趨勢函數(shù)波動情況的局部量，是隨機變量在趨勢附近擺動的隨機誤差。在理論上，其期望值為零，且滿足二階平穩(wěn)條件。

在泛克里金方法中，對于估計值的估計，就轉(zhuǎn)化為對m(x)與R(x)的估計。對于趨勢函數(shù)，可以用多元線性函數(shù)對其進(jìn)行擬合。對殘差的R(x)的估計，使用泛克里金方程組進(jìn)行求解，即：

式中，R*(x0)是待插點的最優(yōu)估計值，R(xi)表示已知樣品的殘差計算值，λi表示R(xi)的泛克里金權(quán)重系數(shù)。

泛克里金的目標(biāo)函數(shù)為

式中，γ是變異函數(shù)。

變異函數(shù)γ是度量隨機變量空間相關(guān)性的工具，在泛克里金插值方法中有重要意義［3］。實驗變異函數(shù)是兩個樣品點的函數(shù)，公式如下：

式中，h是滯后距，γ*(h)為實驗變異函數(shù)，N(h)為滯后距為h時的點對數(shù)，R(xi)與R(xi+h)是實驗點xi和 xi+h的殘差值，i=1，…，N(h)。

變異函數(shù)模型刻畫實驗變異函數(shù)值與滯后距的相關(guān)關(guān)系。球形模型是地質(zhì)統(tǒng)計學(xué)中最常用的變異函數(shù)模型，本文選用球形模型進(jìn)行計算。表達(dá)公式如下：

式中，c0是塊金值，c是基臺值，a是變程。

對式（6）分別關(guān)于ul和λi求導(dǎo)，并令其等于0，組成泛克里金算法的方程組。待插點的估計值正是通過解矩陣形式的方程組獲得的。

2.2 殘差分析基本原理

在統(tǒng)計學(xué)中，殘差是指在回歸分析中的實際觀察值與回歸估計值之差?？梢岳斫猓卸嗌賹?shù)據(jù)，就有多少個殘差。通過殘差所提供的信息，分析出數(shù)據(jù)的可靠性、周期性或其他干擾稱為殘差分析。

殘差分析的基本原理：殘差遵從正態(tài)分布N(δ，σ2)；δ與 σ 之比稱為標(biāo)準(zhǔn)化殘差，以 σ*表示；σ*遵從標(biāo)準(zhǔn)正態(tài)分布N(0，1)。若實驗點的標(biāo)準(zhǔn)化殘差落在置信區(qū)間以外的概率小于等于0.05，則稱可在95%置信區(qū)間將其判別為異常實驗點。

置信區(qū)間是指由樣本統(tǒng)計量所構(gòu)造的總體參數(shù)的估計區(qū)間［17］，展現(xiàn)的是這個參數(shù)的真實值有一定概率落在測量結(jié)果的周圍的程度，其給出的是被測量參數(shù)的測量值的可信程度。通俗來講，置信區(qū)間反映的是一種規(guī)律，當(dāng)不斷改變樣本的時候，有95%的機會，真實值落在置信區(qū)間里，而不僅僅局限在特定樣本。

在擬合趨勢函數(shù)部分，利用置信區(qū)間的殘差分析，可以幫助我們?nèi)サ纛A(yù)測值正確率小于0.05的樣本，增加趨勢函數(shù)的結(jié)構(gòu)性和穩(wěn)定性，以此達(dá)到提升插值效果的目的。這與泛克里金方法的適用條件是相一致的，即研究區(qū)域地質(zhì)變量的變化有一定規(guī)律性，且可用簡單函數(shù)進(jìn)行模擬。

2.3 泛克里金插值技術(shù)路線

插值的樣品數(shù)據(jù)是一系列已知具有某種屬性和三維空間坐標(biāo)的點，這些三維屬性點來源于鉆孔上的樣品數(shù)據(jù)。整體建模流程分四步（圖1（左）），步驟二和三是為插值提供已知信息和目標(biāo)信息，步驟三的實現(xiàn)是本文的研究重點，可以進(jìn)一步劃分為5個步驟（圖1（右））。殘差分析應(yīng)用在對趨勢函數(shù)的擬合上，趨勢函數(shù)的對R(x)的計算至關(guān)重要。

圖1 插值建模流程

3 基于殘差分析擬合趨勢函數(shù)

一般情況下，泛克里金插值方法的使用前提為數(shù)據(jù)滿足正態(tài)分布［3］。在地質(zhì)中，數(shù)據(jù)滿足正態(tài)分布并不能代表數(shù)據(jù)在研究區(qū)具有一定的變化趨勢［14］，不代表可以用數(shù)學(xué)函數(shù)或模型對趨勢進(jìn)行表達(dá)。由此，對于進(jìn)滿足正態(tài)分布，但對趨勢函數(shù)沒有進(jìn)行一定處理的泛克里金插值，其效果是不能保證的。

3.1 擬合趨勢函數(shù)

根據(jù)泛克里金方法的應(yīng)用要求，實驗數(shù)據(jù)的空間分布必須滿足正態(tài)分布。因此，使用泛克里金方法的第一步是對數(shù)據(jù)進(jìn)行正態(tài)分布檢驗。對實驗數(shù)據(jù)進(jìn)行正態(tài)分析，在空間數(shù)據(jù)滿足正態(tài)分布的前提下，可基于泛克里金方法基本原理及變異函數(shù)擬合方法，對數(shù)據(jù)進(jìn)行空間變異結(jié)構(gòu)分析，繪制實驗變異函數(shù)曲線圖，并根據(jù)變異函數(shù)的理論模型擬合實驗變異函數(shù)。

趨勢函數(shù)表達(dá)研究變量與空間坐標(biāo)的相關(guān)性的線性程度，趨勢函數(shù)的精度直接關(guān)系到R(x)的準(zhǔn)確度。如果R(x)的計算結(jié)果不準(zhǔn)確，插值結(jié)果將無意義。在多數(shù)情況下，地質(zhì)變量的變化趨勢都可以用一階多元線性函數(shù)進(jìn)行模擬［13］。本文選用一階多元線性函數(shù)對趨勢進(jìn)行模擬。

3.2 對趨勢函數(shù)進(jìn)行殘差分析

基于之前擬合出的趨勢函數(shù)，對實驗數(shù)據(jù)進(jìn)行回歸分析，計算出所有點對數(shù)據(jù)的殘差和其對應(yīng)的95%置信區(qū)間。用Matlab畫出殘差分析圖，奇異點的定義為置信區(qū)間不包含0的點。按照定義，將奇異點剔除。對剔除奇異點的實驗數(shù)據(jù)再次進(jìn)行回歸分析，驗證奇異點已被剔除，實驗數(shù)據(jù)在95%的置信區(qū)間內(nèi)符合趨勢函數(shù)。

3.3 插值準(zhǔn)備

基于已知采樣點屬性值和擬合出的趨勢函數(shù)，按照式（3），計算出各采樣點的R(xi)。按照定義，R(x)為期望為0，協(xié)方差存在且平穩(wěn)，且只依賴于兩點之間的相對距離，與絕對位置無關(guān)的變量。估計待插點的R*(x)，則要按照克里金插值的步驟。首先計算實驗變異函數(shù)值，量化數(shù)據(jù)間的相關(guān)關(guān)系，選擇變異函數(shù)模型，得到最終變異函數(shù)。具體步驟如下：

1）對得到的樣品數(shù)據(jù)進(jìn)行空間性量化，根據(jù)空間位置坐標(biāo)計算數(shù)據(jù)點對的歐氏距離；

2）按照距離對數(shù)據(jù)進(jìn)行排序分組，計算每組的實驗變異函數(shù)值；

3）在h-γ*(h)圖（圖8）上標(biāo)出各點 (h，γ*(h))，得到實驗變異函數(shù)大致走勢圖；

4）選取擬合相關(guān)系數(shù)較高的球行變異函數(shù)理論模型，應(yīng)用線性規(guī)劃法求解變異函數(shù)模型中的各個參數(shù)；

5）推斷出一個統(tǒng)一的、由各向同性的結(jié)果表達(dá)式，得到最終的變異函數(shù)。

得到變異函數(shù)，便可計算協(xié)方差，進(jìn)行泛克里金方程組求解，計算R*(x)的λi權(quán)系數(shù)。計算得到全部待插點的R*(x)，可按照式（3）得到全部待插點的估計值，以完成插值全過程。

4 實例驗證

本節(jié)對某油田解釋的測井?dāng)?shù)據(jù)進(jìn)行孔隙度的泛克里金插值，共152口井，19508個采樣點，層厚95m，南北方向17772.38m，東西方向15230.5m。將采用本文方法進(jìn)行插值的結(jié)果與未去奇異點的泛克里金方法的插值結(jié)果進(jìn)行統(tǒng)計量對比和交叉驗證對比。

在實際建模中，首先建立結(jié)構(gòu)模型，表達(dá)地質(zhì)體表面的情況；生成網(wǎng)格模型；擬合趨勢函數(shù)，計算變異函數(shù)；求解泛克里金方程組，計算全部待插點的估計值（見圖2）。

圖2 實驗數(shù)據(jù)平面分布圖

4.1 實驗數(shù)據(jù)檢驗

通過研究區(qū)孔隙度的統(tǒng)計特征分析結(jié)果（圖3），可以確定研究區(qū)的孔隙度數(shù)據(jù)符合正態(tài)分布，滿足泛克里金方法的應(yīng)用要求。

圖3 數(shù)據(jù)分析直方圖

4.2 擬合趨勢函數(shù)

本文使用多元線性函數(shù)進(jìn)行擬合，得到擬合度為0.42的趨勢函數(shù)：

m(x，y，z)=-0.079x+0.21y+0.147z+9.39 （7）

對去除奇異點之后的實驗數(shù)據(jù)再次進(jìn)行多元線性擬合趨勢函數(shù)，得到擬合度為0.74的趨勢函數(shù)：

m(x，y，z)=-0.081x+0.18y+0.14z+9.09 （8）

為驗證去奇異點效果，再次進(jìn)行回歸分析，得到殘差分析結(jié)果（見圖5）。

圖4 采樣數(shù)據(jù)殘差分析圖

圖5 去奇異點殘差分析圖

4.3 變異函數(shù)擬合

本文選用球形模型來擬合變異函數(shù)，得到塊金值為2831.33，基臺值為9857.28，變程為1768.57。

4.4 插值計算

泛克里金插值方法的最后一步是計算插值結(jié)果，基于前文的各項計算和擬合出的函數(shù)，根據(jù)公式，計算得到待估點周圍相關(guān)樣品數(shù)據(jù)的權(quán)重值，并進(jìn)行加權(quán)求和，計算出待估點的估計值。對全部待估點重復(fù)以上步驟，可以獲得全部待估點的屬性值，最終形成連續(xù)的三維地質(zhì)屬性模型。

圖6 變異函數(shù)圖

4.5 結(jié)果與分析

4.5.1 統(tǒng)計量分析

比較統(tǒng)計量的方法表達(dá)數(shù)據(jù)的穩(wěn)健性和分布情況是評價估計量好壞的基本方法［15］。T1、T2、T3分別是期望、方差和協(xié)方差。統(tǒng)計量如下：

從期望和方差的角度出發(fā)，期望越接近真實值，無偏性越好，方差越小，有效性越好，協(xié)方差用于衡量兩個變量的總體誤差［18］。在此為比較采樣數(shù)據(jù)和本文方法的差異，以及采樣數(shù)據(jù)和未去奇異點的泛克里金方法的差異，分別計算以采樣數(shù)據(jù)為Z(xi)，以本文方法和泛克里金方法為Z(yi)的協(xié)方差。

由表1可知，由本文算法所得到的插值結(jié)果的統(tǒng)計量結(jié)果與采樣數(shù)據(jù)的期望較接近；相較于未去奇異點的克里金，本文方法的方差更小，且更接近于采樣數(shù)據(jù)方差，反應(yīng)本文插值結(jié)果表達(dá)的地質(zhì)變化規(guī)律更加平穩(wěn)，與采樣數(shù)據(jù)反映的地質(zhì)變化規(guī)律具有一致性。由此可證，去奇異點對插值效果的有效性是有意義的。

表1 統(tǒng)計量結(jié)果

4.5.2 交叉驗證

交叉驗證是指在插值過程中，去掉一個或一部分樣品點，對去掉的樣品點處進(jìn)行插值，并將插值結(jié)果與樣品點處原有值進(jìn)行對比，以驗證插值結(jié)果的準(zhǔn)確性［18］。對本文方法和未去奇異點的泛克里金方法進(jìn)行全部樣品數(shù)據(jù)的交叉驗證，結(jié)果為表2。可以看出，本文方法的均方誤差小于未去奇異點的交叉驗證結(jié)果。由此可知，本文方法具有更高的準(zhǔn)確率。

表2 交叉驗證結(jié)果

5 結(jié)語

本文介紹了泛克里金和殘差分析的基本原理，并介紹了泛克里金插值的技術(shù)路線，并指出本文的研究意義在于提升整體泛克里金插值效果。本文利用測井解釋的孔隙度數(shù)據(jù)，對本文方法和未去奇異點的泛克里金插值方法進(jìn)行分析，得出以下結(jié)論。

1）通過對實驗數(shù)據(jù)趨勢函數(shù)進(jìn)行殘差分析，得到對實際地質(zhì)變量變化規(guī)律更高精度的數(shù)學(xué)表達(dá)，為建立有效精確的三維地質(zhì)屬性模型提供理論指導(dǎo)。

2）將本文方法應(yīng)用到實際測井解釋的孔隙度數(shù)據(jù)，獲得合理有效的插值結(jié)果。將本文方法的插值結(jié)果與未去奇異點的泛克里金插值結(jié)果做統(tǒng)計量對比，結(jié)果表示本文方法插值結(jié)果的統(tǒng)計量結(jié)果更優(yōu)，證明去奇異點對插值效果的有效性是有意義的。

3）對本文方法與未去奇異點的泛克里金插值結(jié)果進(jìn)行交叉驗證，本文方法交叉驗證的結(jié)果的平均誤差和均方根誤差均小于未去奇異點的泛克里金，證明本文方法的插值效果的真實可靠性高于未去奇異點的泛克里金的插值效果。

4）研究得到的趨勢函數(shù)模擬方法適用于所有針對帶趨勢的地質(zhì)建模方法，因此可以應(yīng)用于礦產(chǎn)資源預(yù)測、儲層油氣預(yù)測與評價等多個地質(zhì)領(lǐng)域。