粒子群優化 協同克里金法在確定山地斜坡土層厚度中的應用

王桂林 向林川 孫帆

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

關鍵詞：克里金;協同克里金;土層厚度;空間插值;粒子群優化算法

中圖分類號：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

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

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

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

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

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

1 方法原理

1.1 普通克里金插值法

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

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

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

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

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

通過式（1）確定半變異函數云圖，然后再擬合出相應的模型表達式，應用式（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兩點的半變異函數;N為樣本點數量;x 0為待插值點，Z為其屬性值。

1.2 協同克里金法

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

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

γ 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分開的成對樣本點的數量;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樣本點數量;N 2為變量2樣本點數量;γ 11 為變量1的半方差函數;γ 22 為變量2的半方差函數;γ 12 為變量1與變量2的協半方差函數，γ 12 = γ 21 ;Z 1（x 1i ）為變量1的屬性值;Z 2（x 2j ）為變量2的屬性值;x 0為待插值點，Z* 2（x 0） 為未知點的屬性值

1.3 粒子群優化算法

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

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

1.4 基于粒子群優化算法的克里金插值法

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

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

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

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

1.5 預測精度評價方法

為了能夠對上述方法的空間插值精度進行比較，使各方法間的結果更具有可比性，可通過計算鉆孔點的勘察數據與預測數據的誤差來評估各種方法的優劣[19] 。均方根誤差（RMSE）可以用來評價預測值和真實值之間的接近程度。利用協同克里金法與普通克里金的均方根誤差減少的百分數（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 為樣本數量。

2 工程實例

2.1 研究區域概況與數據處理

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

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

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

2.2 空間插值模型試驗

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

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

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