賈 雨，鄧世武，姚興苗，蔡元菲

（1.成都理工大學(xué) 核技術(shù)與自動化工程學(xué)院，成都610059；2.電子科技大學(xué) 通信與信息工程學(xué)院，成都611731）

克里金（Kriging）插值［1］是基于地質(zhì)統(tǒng)計學(xué)變差函數(shù)模型發(fā)展起來的空間插值方法，是利用區(qū)域化變量的原始數(shù)據(jù)和變差函數(shù)的結(jié)構(gòu)特點，對未采樣點的區(qū)域化變量的值進行最優(yōu)、線性、無偏估計的一種方法，廣泛應(yīng)用于地下水模擬、油氣儲層建模預(yù)測、煤層分布估計等領(lǐng)域。自1951年由南非采礦工程師D.G.Krige提出至今，克里金方法的發(fā)展已形成了一套完整的理論體系，并產(chǎn)生了一些實際有效的程序和軟件。

為了有效地提高插值精度，許多學(xué)者對克里金插值算法進行了改進。嚴(yán)華雯［2］等通過利用加權(quán)最小二乘法優(yōu)化遺傳算法中的適度函數(shù)，改進普通基于遺傳算法優(yōu)化的克里金插值方法；邵才瑞［3］等針對克里金具有平滑性，神經(jīng)網(wǎng)絡(luò)難以反映變量的空間相關(guān)性等缺點，用變差函數(shù)修正神經(jīng)網(wǎng)絡(luò)的目標(biāo)函數(shù)，并利用遺傳算法對神經(jīng)網(wǎng)絡(luò)進行全局優(yōu)化形成了一種遺傳神經(jīng)克里金混合插值方法?？死锝鸩逯档南禂?shù)是以變差函數(shù)的計算為基礎(chǔ)的，早期變差函數(shù)的擬合方法主要依賴于專業(yè)人員的地質(zhì)經(jīng)驗和對地質(zhì)數(shù)據(jù)特征的了解，直接給出理論變差函數(shù)的參數(shù)，沒有客觀的標(biāo)準(zhǔn)。在實際應(yīng)用過程中，為提高變差函數(shù)的擬合精度，學(xué)者們提出了多種變差函數(shù)的擬合方法：加權(quán)最小二乘法、極大似然法、遺傳算法等等。

變差函數(shù)擬合是一個最優(yōu)化問題，加權(quán)最小二乘法和極大似然法常常陷入局部最優(yōu)而無法得到全局最優(yōu)解，同時解的穩(wěn)定性問題也是長期困擾我們的難點問題。遺傳算法有較好的全局尋優(yōu)能力，且對目標(biāo)函數(shù)是否可導(dǎo)沒有限制；但是該方法需要調(diào)整的參數(shù)較多，結(jié)構(gòu)復(fù)雜，影響算法的執(zhí)行效率。為此，本文通過改變粒子群算法中粒子多樣性，結(jié)合地質(zhì)變量的特征和數(shù)據(jù)特征，提出了一種改進的插值方法——基于約束粒子群優(yōu)化的克里金插值算法，通過改進變差函數(shù)的擬合精度，提高克里金插值的精度。

1 克里金插值的基本原理

克里金插值是建立在變差函數(shù)空間分析的基礎(chǔ)上，對有限區(qū)域內(nèi)的區(qū)域化變量取值進行無偏最優(yōu)化估計的一種方法。假設(shè)區(qū)域化變量為Z（x），待插值點為x0，樣本點記為xi（i＝1，2，…，n），在點xi處的屬性值記為Z（xi），則待插值點x0處的屬性值是各個樣本點屬性值的加權(quán)和，記為

其中：λi為采樣點xi的權(quán)重系數(shù)。相距為h的2個空間點x、x＋h處的函數(shù)值Z（x）和Z（x＋h）之間的方差稱為變差函數(shù)，其數(shù)學(xué)表達式為［4］

γ＊（h）為實驗變差函數(shù)，h為滯后距［5］。對于不同的滯后距h，計算出不同的實驗變差函數(shù)值γ＊（h）后需用一個合理的理論變差函數(shù)模型對其進行擬合［6，7］。常用的模型有線形模型、球狀模型、指數(shù)模型和高斯模型［8］等。以高斯模型為例，基本公式如下

其中：C0為塊金常數(shù)；C0＋C為基臺值；C為拱高。當(dāng)時，γ（h）近似達到最大值C0＋C，所以高斯模型的變程約為當(dāng)C0＝0，C＝1時，稱為標(biāo)準(zhǔn)高斯模型。在區(qū)域化變量滿足二階平穩(wěn)的條件下，可推導(dǎo)出克里金方程組

其中：γ（xi，xj）是觀測點xi與xj之間的變差函數(shù)值；γ（xi，x0）是采樣點xi與內(nèi)插點x0之間的變差函數(shù)值；μ是與方差最小化有關(guān)的拉格朗日乘數(shù)。由此方程計算出權(quán)重λi的值，代入公式（1）中即可求出待估點x0處的內(nèi)插值Z（x0）。

2 基于約束粒子群優(yōu)化的克里金算法

2.1 粒子群算法原理

粒子群算法（PSO）是根據(jù)適應(yīng)度函數(shù)使得粒子達到最優(yōu)位置的一種優(yōu)化算法。群體中的每個粒子都有其自身的位置和速度，且粒子的位置和速度是變化的。粒子移動的速度與兩個因素有關(guān)：粒子到目前為止的最好位置和群體中所有粒子到目前為止的最好位置。

假設(shè)粒子在Q維空間中移動，群體中粒子的個數(shù)為N0，則第i個粒子的位置表示為：Xi＝（xi1，xi2，…，xi0），根據(jù)適應(yīng)度函數(shù)，可以得到每個粒子的歷史最優(yōu)位置和群體中所有粒子的最優(yōu)位置，第i個粒子的歷史最優(yōu)位置pbest為：Pi＝（pi1，pi2，…，pi0），群體最優(yōu)位置gbest表示為Pg，是所有Pi（i＝1，2，…，N0）中的最優(yōu)，表達式為：Pg＝（pg1，pg2，…，pg0）。第i個粒子的移動速度為：Vi＝（vi1，vi2，…，vi0）。

在每次迭代時，粒子根據(jù)以下公式更新自己的速度和位置

其中：k為當(dāng)前迭代數(shù)；ω為慣性權(quán)重因子，表示粒子維持原來速度的程度；c1和c2為加速常數(shù)，分別代表將粒子推向個體最優(yōu)位置pbest和群體最優(yōu)位置gbest的權(quán)重；ξ和η是［0，1］區(qū)間內(nèi)均勻分布的隨機數(shù)；γ為收斂因子。

2.2 基于約束粒子群克里金插值算法

粒子群優(yōu)化算法直接應(yīng)用于克里金插值會出現(xiàn)早熟、變差函數(shù)的擬合不符合實際地質(zhì)規(guī)律等問題，為此本文采用高斯變異、地質(zhì)規(guī)律約束等方法進行改進，具體方法如下。

2.2.1 高斯變異

大量的研究表明，由于PSO算法按照追隨種群最優(yōu)粒子的策略進行迭代更新，算法易陷入局部最優(yōu)和早熟收斂等缺陷。鑒于此，本文采用一種基于高斯變異的方法來提高種群的多樣性。在算法出現(xiàn)過早收斂時，能夠使粒子在解空間中的其他區(qū)域進行搜索，跳出局部最優(yōu)，尋找更優(yōu)的解。改變粒子多樣性的方法如下。

在迭代到一半次數(shù)后，開始對粒子進行變異。對每個粒子以概率P進行高斯變異［9］。其中P一般根據(jù)函數(shù)的復(fù)雜性和經(jīng)驗來決定。粒子變異的公式如下

其中：gbestd為全局最優(yōu)在d維的值；σ為高斯白噪聲。

2.2.2 權(quán)重系數(shù)的設(shè)定

一般情況下，用最小二乘法擬合變差函數(shù)的時候，適應(yīng)度函數(shù)由公式（8）給定

其中：F（j）為第j個粒子的適應(yīng)度函數(shù)值；hi，j為第j個粒子的第i個滯后距；γ（hi，j）為第j個粒子在第i個滯后距處的變差函數(shù)值；γ＊（hi，j）表示第j個粒子在第i個滯后距處的實驗變差函數(shù)值。

通過對適應(yīng)度函數(shù)增加權(quán)重系數(shù)來描述對某些實驗變差函數(shù)值的加重，以強化地質(zhì)因素和擬合要求。根據(jù)研究，權(quán)重系數(shù)λi的選擇與以下3個因素有關(guān)。

a.滯后距：實驗變差函數(shù)中滯后距較小的幾個點比較大地反映了區(qū)域化變量的變異程度，在變差函數(shù)擬合時，需要著重考慮滯后距較小的頭幾個實驗變差函數(shù)值，使該處的誤差盡量小，即適應(yīng)度函數(shù)值盡量小。

b.樣本點的密度：在實際的地質(zhì)問題中，樣本點在平面上的分布是極不均勻的。例如在石油勘探的應(yīng)用中，通常鉆井或測井?dāng)?shù)據(jù)是我們的樣本數(shù)據(jù)，而井位在平面上的分布極不均勻，有的區(qū)域井很多，有的區(qū)域井很少。在克里金插值過程中，樣本點少的區(qū)域，其樣本點的權(quán)重相對較大；樣本點多的區(qū)域，其樣本點的權(quán)重相對較小。

c.樣本點的絕對值：樣本點的絕對值可能相差很大，直接使用公式（8）構(gòu)建適應(yīng)度函數(shù)，可能出現(xiàn)絕對值大的樣本點對適應(yīng)度函數(shù)貢獻過大，絕對值小的樣本點對適應(yīng)度函數(shù)貢獻過小的問題。

根據(jù)上述問題，給出新的適應(yīng)度函數(shù)

其中

λi表示每個樣本點的權(quán)重系數(shù)；hi為滯后距；Ni為對應(yīng)滯后距處的樣本對數(shù)；為實驗變差函數(shù)值的平均值。

2.2.3 參數(shù)搜索范圍的約束

由于理論變差函數(shù)模型的每個參數(shù)都具有實際的物理含義，對每個參數(shù)給出搜索范圍，以符合實際規(guī)律。給定待擬合參數(shù)的搜索范圍，每進行一次迭代后，判斷參數(shù)是否在設(shè)定的搜索范圍內(nèi)。若參數(shù)已經(jīng)超過限定的搜索范圍，則采取如下方式進行處理

3 基于約束粒子群優(yōu)化的克里金插值算法的實現(xiàn)

本文選擇高斯模型，根據(jù)粒子群算法，將理論變差函數(shù)中的未知參數(shù)（a，C，C0）看作一個粒子，每個粒子包含a、C和C0這3個分量，基于約束的粒子群優(yōu)化算法步驟如下。

第一步，初始化：設(shè)定粒子個數(shù)m，生成一個粒子群X＝｛x1，x2，…，xm｝，根據(jù)待擬合參數(shù)的物理意義，設(shè)定參數(shù)的取值范圍：0＜a＜兩點之間距離的最大值，0＜C＜實驗變差函數(shù)值的最大值，C0≥0；在各參數(shù)的取值范圍內(nèi)，隨機取m個值，作為每個粒子的初始位置，并將其設(shè)置為當(dāng)前個體最優(yōu)位置pi。設(shè)置粒子各個分量的最大速度值（d表示粒子的第d個分量），在中隨機取值，作為粒子第d個分量的初速度；設(shè)置最大迭代次數(shù)n。

第二步，根據(jù)滯后距hi及對應(yīng)滯后距處的樣本對數(shù)Ni，計算權(quán)重系數(shù)λi。

第三步，根據(jù)適應(yīng)度函數(shù)計算各個粒子的適應(yīng)度函數(shù)值F（j）。

第四步，根據(jù)下面2個公式，確定粒子i的當(dāng)前最優(yōu)位置以及整個粒子群當(dāng)前的全局最優(yōu)位置

其中：表示第i個粒子經(jīng)過k次迭代后的當(dāng)前最優(yōu)位置表示經(jīng)過k次迭代后的全局最優(yōu)位置。

第五步，根據(jù)下面的粒子速度和位置更新公式，更新粒子的位置

第六步，判斷粒子當(dāng)前速度和位置是否超過設(shè)定的范圍，如果是，則在搜索范圍內(nèi)為該粒子重新隨機取值。

第七步，迭代次數(shù)達到最大迭代次數(shù)的一半后，對每個粒子按50%的概率進行高斯變異，防止其陷入局部最優(yōu)。

第八步，返回第三步，重復(fù)該計算過程，直到滿足終止條件（達到最大迭代次數(shù)或預(yù)定的最小適應(yīng)度函數(shù)值），此時獲得了一個理想的最優(yōu)解。

4 實際數(shù)據(jù)測試結(jié)果分析

本文所用的數(shù)據(jù)是采樣得到的實測數(shù)據(jù)，該數(shù)據(jù)包含坐標(biāo)、深度、速度、層位名稱等。對于擬合算法的仿真只是針對二維情況下的插值，從這些數(shù)據(jù)中，只需要提取出某個層位的實測數(shù)據(jù)坐標(biāo)、深度值和速度值。

4.1 擬合效果比較

分別利用最小二乘法、基于約束的PSO算法對實驗變差函數(shù)進行擬合得到的參數(shù)，如表1所示。

實驗變差函數(shù)的散點以及擬合后的理論變差函數(shù)曲線如圖1所示。圖1中，藍色代表計算出來的實驗變差函數(shù)散點；紅色曲線為最小二乘法擬合出的理論變差函數(shù)曲線；綠色曲線代表利用約束的PSO算法得到的理論變差函數(shù)曲線。由圖可知，本文提出的約束粒子群優(yōu)化算法在滯后距較小的地方取得的了較好的效果；而在滯后距較大的地方，變差函數(shù)值更接近于基臺值C0＋C。結(jié)合表1中的數(shù)據(jù)也可得知：相較于最小二乘法，本文提出的方法具有較小的基臺值，更接近實驗變差函數(shù)值的穩(wěn)定值（由后幾個實驗變差函數(shù)值反映出來）。

4.2 插值效果比較及誤差分析

圖1 最小二乘、約束PSO擬合變差函數(shù)對比曲線圖Fig.1 Contrast figure of the variation functions fitted by the least squares method and the constraint PSO

表1 擬合后的理論變差函數(shù)模型參數(shù)值Table 1 Parameters of the theoretical variation functions after fitting

利用上文中得到的理論變差函數(shù)模型對待插值點進行屬性值估計。本文對1 000×1 000個待插值點進行屬性值估計。待插值點的坐標(biāo)表示為

其中：xmin表示已知點x坐標(biāo)的最小值；ymin表示已知點y坐標(biāo)的最小值；xstep代表x方向的步長；ystep代表y方向的步長；i，j分別是從0到999的整數(shù)。運用本文4.1節(jié)中獲得的2種理論變差函數(shù)曲線對待插值點做插值處理得到的插值效果如圖2所示。從圖中可看出，采用約束的PSO算法擬合后變差函數(shù)插值結(jié)果高點更清晰。實際鉆井和地質(zhì)分析也證實了本文方法的有效性。

為進一步比較兩種擬合算法得到變差函數(shù)的精度，本文采用交叉驗證的方法。對已知點進行插值處理，將計算得到的屬性值與真實的屬性值之差作為誤差，如圖3所示。對已知的94個樣本點插值后，在大部分已知點處，利用本文提出的基于約束粒子群優(yōu)化的克里金插值算法得到的誤差更小。圖中最后一條數(shù)據(jù)為平均誤差，從中也可以看出本文提出的方法優(yōu)越于常規(guī)克里金插值算法。

圖2 使用最小二乘法擬合變差函數(shù)的插值（左圖）和使用約束的PSO算法擬合的插值結(jié)果Fig.2 Interpolation renderings with variation functions fitted by the least squares method and the constraint PSO

圖3 最小二乘法和約束的PSO變差函數(shù)插值后誤差分析Fig.3 Error analysis of the least squares method and the constraint PSO

5 結(jié)論

變差函數(shù)擬合結(jié)果直接影響克里金插值效果，選擇合適的方法擬合變差函數(shù)對于改進克里金插值效果具有較大的作用。本文在考慮地質(zhì)變量特征和數(shù)據(jù)特征的基礎(chǔ)上，將基于約束粒子群優(yōu)化算法應(yīng)用于變差函數(shù)擬合中，應(yīng)用實測數(shù)據(jù)計算出實驗變差函數(shù)值，分別利用最小二乘法和基于約束粒子群算法對其進行變差函數(shù)計算和擬合，并運用擬合結(jié)果進行二維插值，實驗結(jié)果表明：基于約束粒子群優(yōu)化的克里金插值算法獲得的插值效果具有較高插值精度。

［1］孫洪泉.地質(zhì)統(tǒng)計學(xué)及其應(yīng)用［M］.徐州：中國礦業(yè)大學(xué)出版社，1990.Sun H Q.Geological Statistics and Its Application［M］.Xuzhou：China University of Mining and Technology Press，1990.（In Chinese）

［2］嚴(yán)華雯，吳健平.加權(quán)最小二乘法改進遺傳克里金插值方法研究［J］.計算機技術(shù)與發(fā)展，2012，22（3）：92－95.Yan H W，Wu J P.Research on genetic algorithm Kriging optimized by weight least square［J］.Computer Technology and Development，2012，22（3）：92－95.（In Chinese）

［3］邵才瑞，印興耀，李洪奇，等.儲層屬性的遺傳神經(jīng)克里金插值方法及其應(yīng)用［J］.中國石油大學(xué)學(xué)報：自然版，2007，31（5）：35－40.Shao C Y，Yin X Y，Li H Q，etal.An integrated genetic-neural-Kriging interpolation method for reservoir property and its application［J］.Journal of China University of Petroleum，2007，31（5）：35－40.（In Chinese）

［4］Gringarten E，Deutsch C V.Teacher's aide variogram interpretation and modeling［J］.Mathematical Geology，2001，33（4）：507－534.

［5］喬金海，潘懋，金毅，等.基于Kriging方法的天然地基承載力三維模擬及分析［J］.北京大學(xué)學(xué)報：自然科學(xué)版，2011，47（5）：812－818.Qiao J H，Pan M，Jin Y，etal.3Dmodeling and analysis of natural foundation bearing capacity based on Kriging method［J］.Acta Scientiarum Naturalium Universitatis Pekinensis，2011，47（5）：812－818.（In Chinese）

［6］Park J J，Shin K I，Lee J H，etal.Detecting and cleaning outliers for robust estimation of variogram models in insect count data［J］.Ecological Research，2012，27（1）：1－13.

［7］Clark R G，Allingham S.Robust resampling confidence intervals for empirical variograms［J］.Mathematical Geosciences，2011，43（2）：243－259.

［8］顏輝武，祝國瑞，徐智勇，等.基于Kriging水文地質(zhì)層的三維建模與體視化［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2004，29（7）：611－614.Yan H W，Zhu G R，Xu Z Y，etal.Volume rendering and 3Dmodeling of hydrogeologic layer based on Kriging algorithm［J］.Geomatics and Information Science of Wuhan University，2004，29（7）：611－614.（In Chinese）

［9］Krohling A K.Gaussian swarm：A novel particle swarm optimization algorithm［C］//Proceedings of the 2004IEEE，Singapore，Conference on Cybernetics and Intelligent Systems，2004：372－376.