一種適用于泛克里金趨勢(shì)函數(shù)的優(yōu)化方法?

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

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

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

1 引言

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

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

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

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

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

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

2.1 泛克里金方法概述

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

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

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

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

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

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

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

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

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

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

式中，c0是塊金值，c是基臺(tái)值，a是變程。

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

2.2 殘差分析基本原理

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

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

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

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

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

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

圖1 插值建模流程

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

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

3.1 擬合趨勢(shì)函數(shù)

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

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

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

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

3.3 插值準(zhǔn)備

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

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

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

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

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

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

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

4 實(shí)例驗(yàn)證

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

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

圖2 實(shí)驗(yàn)數(shù)據(jù)平面分布圖

4.1 實(shí)驗(yàn)數(shù)據(jù)檢驗(yàn)

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

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

4.2 擬合趨勢(shì)函數(shù)

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

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

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

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

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

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

圖5 去奇異點(diǎn)殘差分析圖

4.3 變異函數(shù)擬合

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

4.4 插值計(jì)算

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

圖6 變異函數(shù)圖

4.5 結(jié)果與分析

4.5.1 統(tǒng)計(jì)量分析

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

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

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

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

4.5.2 交叉驗(yàn)證

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

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

5 結(jié)語(yǔ)

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

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

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

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

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