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

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

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

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

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

1 克里金插值的基本原理

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

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

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

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

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

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

2.1 粒子群算法原理

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

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

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

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

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

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

2.2.1 高斯變異

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

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

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

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

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

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

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

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

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

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

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

其中

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.1 擬合效果比較

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

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

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

圖1 最小二乘、約束PSO擬合變差函數(shù)對(duì)比曲線圖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ù)模型對(duì)待插值點(diǎn)進(jìn)行屬性值估計(jì)。本文對(duì)1 000×1 000個(gè)待插值點(diǎn)進(jìn)行屬性值估計(jì)。待插值點(diǎn)的坐標(biāo)表示為

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

為進(jìn)一步比較兩種擬合算法得到變差函數(shù)的精度，本文采用交叉驗(yàn)證的方法。對(duì)已知點(diǎn)進(jìn)行插值處理，將計(jì)算得到的屬性值與真實(shí)的屬性值之差作為誤差，如圖3所示。對(duì)已知的94個(gè)樣本點(diǎn)插值后，在大部分已知點(diǎn)處，利用本文提出的基于約束粒子群優(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ù)對(duì)于改進(jìn)克里金插值效果具有較大的作用。本文在考慮地質(zhì)變量特征和數(shù)據(jù)特征的基礎(chǔ)上，將基于約束粒子群優(yōu)化算法應(yīng)用于變差函數(shù)擬合中，應(yīng)用實(shí)測(cè)數(shù)據(jù)計(jì)算出實(shí)驗(yàn)變差函數(shù)值，分別利用最小二乘法和基于約束粒子群算法對(duì)其進(jìn)行變差函數(shù)計(jì)算和擬合，并運(yùn)用擬合結(jié)果進(jìn)行二維插值，實(shí)驗(yàn)結(jié)果表明：基于約束粒子群優(yōu)化的克里金插值算法獲得的插值效果具有較高插值精度。

［1］孫洪泉.地質(zhì)統(tǒng)計(jì)學(xué)及其應(yīng)用［M］.徐州：中國(guó)礦業(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ìn)遺傳克里金插值方法研究［J］.計(jì)算機(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］邵才瑞，印興耀，李洪奇，等.儲(chǔ)層屬性的遺傳神經(jīng)克里金插值方法及其應(yīng)用［J］.中國(guó)石油大學(xué)學(xué)報(bào)：自然版，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é)報(bào)：自然科學(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］顏輝武，祝國(guó)瑞，徐智勇，等.基于Kriging水文地質(zhì)層的三維建模與體視化［J］.武漢大學(xué)學(xué)報(bào)：信息科學(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.