王桂林 向林川 孫帆

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

關鍵詞：克里金;協(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)質耕地，國土資源部提出了開展低丘緩坡等未利用土地開發(fā)試點的工作。殘坡積土是低丘緩坡重要的組成部分，高質量的殘坡積土層厚度信息在土地資源管理、工程建設、地質災害預警預報等方面具都有重要意義。獲取比較精確的土層厚度數(shù)據的方法之一是布設高密度的鉆孔，但由于受到經濟水平、技術手段和地形條件的限制，鉆孔點的數(shù)量是一定的，這為確定土層厚度帶來了較大的誤差。隨著地統(tǒng)計學理論在工程界不斷應用，空間插值的方法也日益成為解決上述問題的經典方法。將利用統(tǒng)計學與地理信息系統(tǒng)相結合的方法，基于已知鉆孔點的觀測數(shù)據進行空間插值來獲得估計點的土層厚度。

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

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

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

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

1 方法原理

1.1 普通克里金插值法

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

在克里金法中，用來衡量各個樣本點之間空間相關程度的是半變異函數(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)時的拐點，前者表示的是兩個非常接近的樣本點之間的誤差及空間變異，后者表示兩樣本點超過此間離后將不存在空間相關性。

通過式（1）確定半變異函數(shù)云圖，然后再擬合出相應的模型表達式，應用式（2）確定內插所需要的權重，并通過式（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為待定的權重;u為拉格朗日乘子;r（h ij ）為x i與x j兩點的半變異函數(shù);N為樣本點數(shù)量;x 0為待插值點，Z為其屬性值。

1.2 協(xié)同克里金法

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

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

γ 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是點關于變量1屬性值;Z? 2是點關于變量2屬性值。

擬合后用式（5）即可確定空間插值所需要的權重，并通過式（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待定的權重;λ 2j 為變量2待定的權重;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）是一種進化演變技術，最早是由美國的心理學研究者 Kennedy博士和從事計算智能研究的Eberhart博士受到人工生命的研究結果啟發(fā)于 1995 年提出的一種基于群智能的優(yōu)化算法[16-17] 。

粒子群算法的實質是一種信息共享，粒子根據自身的個體最優(yōu)信息及群體的最優(yōu)信息不斷更新自己的速度和位置，最后收斂于全局最優(yōu)解。每個粒子都有一個有待優(yōu)化問題所決定的適應度值用來評價該粒子的優(yōu)良程度，一個粒子還有一個速度用于決定它的運動距離即速度，粒子根據自身的運動經驗和群體中其它粒子的運動經驗來調整自己的運動。粒子從此時刻運動到下一時刻的速度由式（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 全局最好位置; ω 為慣性權重值。

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

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

同時，大量研究表明， 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 預測精度評價方法

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

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ù)據處理

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

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

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

2.2 空間插值模型試驗

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

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

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