王桂林 向林川 孫帆

摘 要：克里金法是廣泛應(yīng)用的空間插值方法，但僅考慮單一因素的普通克里金法在確定山地斜坡土層厚度中存在較大誤差。針對普通克里金法中的不足之處，提出了一種確定土層厚度的基于粒子群優(yōu)化的協(xié)同克里金法。該方法首先用粒子群優(yōu)化算法擬合半變異函數(shù)，然后將該函數(shù)用于以高程值作為輔助變量的協(xié)同克里金法中，并根據(jù)均方根誤差來評價土層厚度的不確定性。將該方法應(yīng)用于重慶萬盛某邊坡土層厚度的確定，通過交叉驗證，結(jié)果表明：與普通克里金插值法相比較，考慮高程的協(xié)同克里金法插值的均方根誤差降低了39.32%;基于粒子群優(yōu)化的普通克里金法和協(xié)同克里金法的均方根誤差分別降低了28.79%和48.45%?；诹Ｗ尤簝?yōu)化的協(xié)同克里金插值法對提高土層厚度的插值精度有較大作用。

關(guān)鍵詞：克里金;協(xié)同克里金;土層厚度;空間插值;粒子群優(yōu)化算法

中圖分類號：TU191.1

文獻標志碼：A? 文章編號：1674-4764（2018）06-0060-07

Application of cooperative Kriging method based on particle swarm

optimization in estimation slope soil thickness

Wang Guilin1，2 ， Xiang Linchuan2，Sun Fan2

（a.Key Laboratory of New Technology for Construction of Cities in Mountain Area， Ministry of Education;

b.School of Civil Engineering， Chongqing University， Chongqing 400045， P.R. China）

Abstract：?The Kriging method is widely used for the spatial interpolation. However， the conventional Kriging method which only considers a singer factor generally leads to a considerable inaccuracy. In this paper， a cooperative kriging method based on particle swarm optimization is proposed to estimate the soil thickness distribution. Estimation is divided into two steps. Firstly， the particle swarm optimization is used to fit the semi-variance function. Secondly， the cooperative Kriging method which uses the altitude as an auxiliary variable is employed for estimation. In addition， a root mean square error is obtained to evaluate the estimation uncertainty of soil thickness. The proposed method is applied to estimate the soil thickness of a slope in Wansheng， Chongqing. It shows that compared with the conventional method， the cooperative Kriging approach improves the estimation accuracy by reducing the standard deviation by 39.32%， indicating that the proposed method is advantageous in improving the accuracy of spatial interpolation.

Keywords：Kriging; Co-Kriging; soil thickness; spatial interpolation; particle swarm optimization

隨著中國城鎮(zhèn)化進程的不斷推進，為了更有效地保護優(yōu)質(zhì)耕地，國土資源部提出了開展低丘緩坡等未利用土地開發(fā)試點的工作。殘坡積土是低丘緩坡重要的組成部分，高質(zhì)量的殘坡積土層厚度信息在土地資源管理、工程建設(shè)、地質(zhì)災(zāi)害預(yù)警預(yù)報等方面具都有重要意義。獲取比較精確的土層厚度數(shù)據(jù)的方法之一是布設(shè)高密度的鉆孔，但由于受到經(jīng)濟水平、技術(shù)手段和地形條件的限制，鉆孔點的數(shù)量是一定的，這為確定土層厚度帶來了較大的誤差。隨著地統(tǒng)計學理論在工程界不斷應(yīng)用，空間插值的方法也日益成為解決上述問題的經(jīng)典方法。將利用統(tǒng)計學與地理信息系統(tǒng)相結(jié)合的方法，基于已知鉆孔點的觀測數(shù)據(jù)進行空間插值來獲得估計點的土層厚度。

空間插值是利用已知的部分空間樣本信息，對未知地理空間的特征屬性值進行估計，是地理信息系統(tǒng)的重要功能模塊之一[1] ，并在礦業(yè)、水文、氣候預(yù)報、農(nóng)業(yè)等領(lǐng)域有著廣泛的應(yīng)用。目前已發(fā)展了較多的空間插值方法：如泰森多邊形法[2] 、克里金（Kriging）插值法[3] 、反距離加權(quán)平均法[4] 、趨勢面分析法[5] 、多項式回歸法[6] 等。但這些方法只是局限于已知點的單一屬性值，沒有考慮到其他因素對待估計點屬性值的影響（如高程對土層厚度的影響等）。

基于上述方法在估計空間屬性值時對變量考慮的單一性，學者提出了一種能考慮多個相關(guān)變量互相影響的協(xié)同克里金法。Yates等[7] 將裸土表面溫度和土壤沙粒含量作為協(xié)變量，利用協(xié)同克里金法（Co-Kriging）估計質(zhì)量含水率;Ghadermazi 等[8] 利用pH值作為輔助變量，估計了飲用水中硝酸鹽的含量。胡丹桂等[9] 用考慮降雨量的協(xié)克里金法來研究空氣的濕度的空間變化;許晶玉等[10] 用考慮土壤粗砂含量和全氮含量的協(xié)克里金法來研究山東省種植區(qū)地下水硝態(tài)氮污染空間變異及分布規(guī)律;杜文鳳等[11] 將地震數(shù)據(jù)引入?yún)f(xié)克里金法中估計煤層厚度的分布規(guī)律;黃大年等[12] 首次將重力梯度引入?yún)f(xié)克里金法中來研究巖脈傾向的問題。

半變異函數(shù)是克里金法中反映區(qū)域化變量空間變化特征的有效數(shù)學模型，由其確定擬合模型參數(shù)直接影響插值精度。嚴華雯等[13] 通過利用加權(quán)最小二乘法優(yōu)化遺傳算法中的適度函數(shù)，改進普通基于遺傳算法優(yōu)化的克里金插值方法;張強等[14] 采用基于線性遞減權(quán)值的粒子群算法估計半變異函數(shù)參數(shù)的方法，提高了插值的精度;賈雨等[15] 將約束粒子群優(yōu)化算法用于克里金的插值研究，結(jié)果表明能獲得較好的插值精度。

本文采用粒子群優(yōu)化算法與協(xié)同克里金法相結(jié)合的方法對土層厚度的空間插值進行研究，并對插值結(jié)果進行對比，以期提高土層厚度插值分布的精度。

1 方法原理

1.1 普通克里金插值法

克里金插值法又稱空間自協(xié)方差最佳插值法，是以南非礦業(yè)工程師D.G.Krige的名字命名的一種最優(yōu)內(nèi)插法。該方法首先考慮的是空間屬性在空間位置上的變異分布，確定對一個待插點值有影響的距離范圍，然后用此范圍內(nèi)的采樣點來估計待插點的屬性值。該方法考慮了已知信息與待估計點相互間的空間結(jié)構(gòu)特性。為達到線性、無偏和最小估計方差的估計，而對每一個樣品賦予一定的權(quán)重，最后進行加權(quán)平均來估計待預(yù)測點屬性值的一種方法。

在克里金法中，用來衡量各個樣本點之間空間相關(guān)程度的是半變異函數(shù)

r（h）= 1 2N ·∑ N i=1 [Z（x i）-Z（x i+h）]2 （1）

式中： h 為兩點之間距離; N是由h 分開的成對樣本點的數(shù)量; Z 是點的屬性值。

一種典型的半變異函數(shù)圖像如圖1所示，半變異值隨距離增大而增大，其中有兩個非常重要的點：間距為0時的點以及函數(shù)趨于平穩(wěn)時的拐點，前者表示的是兩個非常接近的樣本點之間的誤差及空間變異，后者表示兩樣本點超過此間離后將不存在空間相關(guān)性。

通過式（1）確定半變異函數(shù)云圖，然后再擬合出相應(yīng)的模型表達式，應(yīng)用式（2）確定內(nèi)插所需要的權(quán)重，并通過式（3）進行未知樣本點屬性值的估計。 ?∑ N i=1 λ i·r（h ij ）+u=r（h 0j ）∑ N i=1 λ i=1??? （2）

Z*（x 0）=∑ N i=1 λ i…Z i（x） （3）

式中：λ i為待定的權(quán)重;u為拉格朗日乘子;r（h ij ）為x i與x j兩點的半變異函數(shù);N為樣本點數(shù)量;x 0為待插值點，Z為其屬性值。

1.2 協(xié)同克里金法

協(xié)同克里金插值法可以利用同一變量在不同時空或不同變量在同一時空上的協(xié)同區(qū)域化性質(zhì)，用易于測定的變量來對那些難以測定的屬性或變量進行估值，或者用樣品多的變量對樣品少的變量進行估值;如果兩種以及多種屬性具有顯著的空間相關(guān)性，則可以優(yōu)先考慮協(xié)同克里金插值法。

在協(xié)同克里金法中，γ 11 、γ 22 分別表示變量1與變量2的半變異模型函數(shù)，其計算方法與普通克里金插值法中的半方差計算方法相同;γ 12 為協(xié)半變異函數(shù)，是用來衡量兩個變量各個樣本點之間空間相關(guān)程度的表達式

γ 12 = 1 2N ·∑ N i=1? Z 1（x i）-Z 1（x i+h） ·

[Z 2（x i）-Z 2（x i+h）] （4）

式中： h為兩點之間距離;N是由h分開的成對樣本點的數(shù)量;Z 1是點關(guān)于變量1屬性值;Z? 2是點關(guān)于變量2屬性值。

擬合后用式（5）即可確定空間插值所需要的權(quán)重，并通過式（6）進行未知樣本點屬性值的估計。

∑ N 1 i=1 λ 1i γ 11 （x 1i -x i）+∑ N 2 j=1 λ 2j γ 21 （x 2j -

x i）+u 2=γ 21 （x 0-x i）∑ N 1 i=1 λ 1i γ 21 （x 1i -x j）+∑ N 2 j=1 λ 2j γ 22 （x 2j -

x j）+u 2=γ 22 （x 0-x j）∑ N 1 i=1 λ 1i =0，∑ N 1 i=1 λ 1i =0?? （5）

Z* 2（x 0）=∑ N 1 i=1 λ 1i Z 1（x 1i ）+∑ N 2 j=1 λ 2j Z 2（x 2j ） （6）

式中： λ 1i 為變量1待定的權(quán)重;λ 2j 為變量2待定的權(quán)重;N 1為變量1樣本點數(shù)量;N 2為變量2樣本點數(shù)量;γ 11 為變量1的半方差函數(shù);γ 22 為變量2的半方差函數(shù);γ 12 為變量1與變量2的協(xié)半方差函數(shù)，γ 12 = γ 21 ;Z 1（x 1i ）為變量1的屬性值;Z 2（x 2j ）為變量2的屬性值;x 0為待插值點，Z* 2（x 0） 為未知點的屬性值

1.3 粒子群優(yōu)化算法

粒子群優(yōu)化算法（PSO）是一種進化演變技術(shù)，最早是由美國的心理學研究者 Kennedy博士和從事計算智能研究的Eberhart博士受到人工生命的研究結(jié)果啟發(fā)于 1995 年提出的一種基于群智能的優(yōu)化算法[16-17] 。

粒子群算法的實質(zhì)是一種信息共享，粒子根據(jù)自身的個體最優(yōu)信息及群體的最優(yōu)信息不斷更新自己的速度和位置，最后收斂于全局最優(yōu)解。每個粒子都有一個有待優(yōu)化問題所決定的適應(yīng)度值用來評價該粒子的優(yōu)良程度，一個粒子還有一個速度用于決定它的運動距離即速度，粒子根據(jù)自身的運動經(jīng)驗和群體中其它粒子的運動經(jīng)驗來調(diào)整自己的運動。粒子從此時刻運動到下一時刻的速度由式（7）確定，位置由式（8）計算得到。具體流程見圖2。

optimization algorithm

vij （t+1）=ω·vij （t）+A+BA=c 1·rij? 1（t）·（xij?? pbest? （t）-xij （t））B=c 2·rij? 2（t）·（xj? gbest? （t）-xij （t））?? （7）

xij （t+1）=xij （t）+vij （t+1）? （8）

式中： i為第i個粒子，j為第j個維度;t為更新次數(shù);v為粒子速度，x為粒子位置;r 1、r 2為0到1均勻分布的隨機數(shù);x? pbest? 粒子自身最好位置，x? gbest 全局最好位置; ω 為慣性權(quán)重值。

1.4 基于粒子群優(yōu)化算法的克里金插值法

擬合半變異函數(shù)的模型有很多種，但無論是哪種擬合模型，都涉及到3個擬合參數(shù)，分別是圖1中的核（nugget）、變程（range）和梁（sill）。在通常情況下，無論是普通克里金法還是協(xié)同克里金法，在確定上述3個參數(shù)時用的均為最小二乘法，即在滿足方差最小的條件下得到的，這就會使得3個參數(shù)存在較大的誤差。利用粒子群優(yōu)化算法能快速尋找到全局最優(yōu)解的特點，將該優(yōu)化算法應(yīng)用到確定半變異模型3個參數(shù)上。其中該優(yōu)化算法中的適應(yīng)度函數(shù)為插值結(jié)果的均方根誤差。

同時，大量研究表明， PSO算法易陷入局部最優(yōu)和早熟收斂等缺陷。采用一種基于高斯變異的方法來提高種群的多樣性[15] 。在算法出現(xiàn)過早收斂時，能夠使粒子在其他區(qū)域進行搜索，跳出局部最優(yōu)，尋找更優(yōu)的解。即在迭代到一半次數(shù)后，開始對粒子進行變異，對每個粒子進行高斯變異[18] 。公式為

x d= gbest （d）×（0.5+σ）? （9）

式中： gbest （d）為全局最優(yōu)在 d 維的值;σ為高斯白噪聲。

1.5 預(yù)測精度評價方法

為了能夠?qū)ι鲜龇椒ǖ目臻g插值精度進行比較，使各方法間的結(jié)果更具有可比性，可通過計算鉆孔點的勘察數(shù)據(jù)與預(yù)測數(shù)據(jù)的誤差來評估各種方法的優(yōu)劣[19] 。均方根誤差（RMSE）可以用來評價預(yù)測值和真實值之間的接近程度。利用協(xié)同克里金法與普通克里金的均方根誤差減少的百分數(shù)（RRMSE）來表示預(yù)測精度的提高程度。

RMSE =?? ∑ N i=1? Z（x i）-Z*（x i） 2 /N??? （10）

RRMSE = ?RMSE OK -RMSE COK ?RMSE OK? ??（11）

式中： Z（x i）與Z*（x i）分別為x i 處的勘察值與估計值; N 為樣本數(shù)量。

2 工程實例

2.1 研究區(qū)域概況與數(shù)據(jù)處理

2.1.1 研究區(qū)域概況? 試驗區(qū)位于重慶萬盛經(jīng)濟開發(fā)區(qū)，地貌單元為丘陵斜坡地貌;總體趨勢為西側(cè)高，東側(cè)低;地表總體坡度2°～34°。上圖中的紅色部分為鉆孔點的位置，已知數(shù)據(jù)有勘察區(qū)的平面圖、剖面圖及鉆孔柱狀圖。試驗區(qū)最高海拔為418.23 m， 最低海拔為334.77 m;共有124個鉆孔點。具體情況見圖3。

2.1.2 鉆孔點數(shù)據(jù)處理與分析? 對試驗區(qū)124個鉆孔的高程數(shù)據(jù)和土層厚度分析，結(jié)果顯示：土層厚度和鉆孔點的高程有一定的規(guī)律（為了便于與土層厚度比較，高程基準點設(shè)為320 m），具體表現(xiàn)為：在一定的高程范圍內(nèi)，鉆孔的土層厚度大致和高程值成正相關(guān)（見圖4）。

協(xié)同克里金法中，作為輔助變量的前提條件是該變量與待估計變量之間存在著相關(guān)性。用SPSS統(tǒng)計軟件中的數(shù)據(jù)分析功能得到高程值和土層厚度的相關(guān)系數(shù)為0.552，且在0.01的顯著性水平的條件下，相關(guān)性程度為顯著性相關(guān)（表1）。因此，可將高程值作為提高土層厚度插值精度的輔助變量。

2.2 空間插值模型試驗

采用交叉驗證的方法[20] 來估計土層厚度，并用數(shù)理統(tǒng)計的方法來比較不同估計模型的估值精度，即首先將待預(yù)測點的屬性值 Z（x i） 暫時剔除，然后將最小二乘法得到的半變異函數(shù)表達式用于普通克里金法和考慮高程的協(xié)同克里金法來預(yù)測 Z（x i） 的值。結(jié)合相應(yīng)的公式并利用MATLAB編制程序求出拉克朗日系數(shù)和待估計點的影響范圍內(nèi)的各個鉆孔點的權(quán)重系數(shù)，最后利用式（3）與式（6）分別求出待估計點的屬性值 Z*（x i）。 重復(fù)以上過程直到對所有的124個觀測點進行估計。

基于粒子群優(yōu)化算法的克里金插值法與基于最小二乘法的克里金插值法區(qū)別在于：粒子群優(yōu)化算法中的半變異模型中的3個參數(shù)是未知的;隨機的給每個粒子賦值，這樣每一個粒子就對應(yīng)了一個半變異函數(shù)的模型;然后將每個粒子所對應(yīng)的半變異函數(shù)模型，帶入用MATLAB編制的相應(yīng)插值程序中，求出均方根誤差，并且將該均方根誤差作為該粒子的適應(yīng)度函數(shù)的適應(yīng)值。當循環(huán)結(jié)束后，最小的適應(yīng)度值的粒子即為所求。進而求出了半變異函數(shù)的表達式。

無論是哪一種方法，每一個鉆孔點可以得到真實值和估計值兩類數(shù)據(jù)。將各個插值方法得到的插值結(jié)果與真實值做成火柴棍圖，見圖5～圖8。