王桂林，向林川，孫帆

(重慶大學(xué) a.山地城鎮(zhèn)建設(shè)與新技術(shù)教育部重點實驗室；b.土木工程學(xué)院, 重慶 400045)

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

基于上述方法在估計空間屬性值時對變量考慮的單一性，學(xué)者提出了一種能考慮多個相關(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ù)學(xué)模型,由其確定擬合模型參數(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ù)

(1)

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

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

圖1 半變異函數(shù)圖

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

(2)

(3)

式中：λi為待定的權(quán)重；u為拉格朗日乘子；r(hij)為xi與xj兩點的半變異函數(shù)；N為樣本點數(shù)量；x0為待插值點，Z為其屬性值。

1.2 協(xié)同克里金法

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

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

[Z2(xi)-Z2(xi+h)]

(4)

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

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

(5)

(6)

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

粒子群優(yōu)化算法(PSO)是一種進化演變技術(shù)，最早是由美國的心理學(xué)研究者 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。

圖2 粒子群優(yōu)化算法流程圖Fig.2 The flow chart of particle swarm

(7)

xij(t+1)=xij(t)+vij(t+1)

(8)

式中：i為第i個粒子，j為第j個維度；t為更新次數(shù)；v為粒子速度，x為粒子位置；r1、r2為0到1均勻分布的隨機數(shù)；xpbest粒子自身最好位置，xgbest全局最好位置；ω為慣性權(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]。公式為

xd=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ù)測精度的提高程度。

(10)

(11)

式中：Z(xi)與Z*(xi)分別為xi處的勘察值與估計值；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。

圖3 試驗區(qū)鉆孔分布圖Fig.3 The distribution of

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

圖4 土層厚度和高程關(guān)系圖Fig.4 The correlation between soil thickness and

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

表1 SPSS統(tǒng)計軟件計算結(jié)果Table 1 The result calculated by SPSS statistical software

注：**水平(雙側(cè))上顯著相關(guān)。

2.2 空間插值模型試驗

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

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

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

圖5 最小二乘法的O-Kriging插值結(jié)果與真實土層厚度Fig.5 O-Kriging based on the least squares method interpolation results and the true soil

圖6 粒子群優(yōu)化算法的O-Kriging插值結(jié)果與真實土層厚度Fig.6 O-Kriging based on PSO method interpolation results and the true soil

圖7 最小二乘法的CO-Kriging插值結(jié)果與真實土層厚度Fig.7 CO-Kriging based on the least squares method interpolation results and the true soil

圖8 粒子群優(yōu)化算法的Co-Kriging插值結(jié)果與真實土層厚度Fig.8 Co-Kriging based on PSO method interpolation results and the true soil

從圖5～圖8中可以看出，在上述的4類方法中，基于最小二乘法的普通克里金插值法與勘察值之間的誤差是最大的；而考慮高程的協(xié)同克里金法的插值的精度要高于基于粒子群優(yōu)化算法的普通克里金法；基于粒子群優(yōu)化算法的協(xié)同克里金插值法的精度在4類方法中是最高的，除了極個別的土層厚度較厚的鉆孔點外，其他鉆孔點均與實際情況吻合得較好。

用數(shù)理統(tǒng)計的方法來比較4種插值結(jié)果的精度的提高程度。其結(jié)果見表2。

表2 O-Kriging與Co-Kriging預(yù)測精度的比較Table 2 Error analysis of O-Kriging and Co-Kriging

從表2可以看出，當(dāng)采用普通克里金插值法對邊坡的土層厚度進行空間插值時，基于粒子群優(yōu)化算法的普通克里金插值法較基于最小二乘法的普通克里金插值法的均方根誤差減少了28.79%；當(dāng)均采用最小二乘法去擬合半變異函數(shù)時；協(xié)同克里金法的均方根誤差較普通克里金法減少了39.32%；當(dāng)都采用PSO法去擬合模型參數(shù)時；協(xié)同克里金法的均方根誤差較普通克里金法減少了19.66%；同時還可以從表中得出結(jié)論：當(dāng)考率單一因素時，考慮高程的協(xié)同克里金法較基于粒子群優(yōu)化算法的普通克里金插值法的插值精度提高了10.53%。

3 結(jié)論

以重慶某坡地為研究對象，利用該區(qū)域勘察報告中的124個鉆孔點的勘察數(shù)據(jù)為試驗樣本，采用交叉驗證的方法分別用普通克里金插值法和考慮高程的協(xié)同克里金插值法對土層厚度進行虛擬差值。還對在擬合半變異函數(shù)模型參數(shù)時存在的誤差進行了研究，用粒子群優(yōu)化算法擬合半變異函數(shù)模型，并對4種不同方法的插值結(jié)果進行了對比分析，得到的結(jié)論如下：

1)將最小二乘法擬合的半變異函數(shù)用于普通克里金法插值時會產(chǎn)生較大誤差。

2)對土層厚度進行克里金插值時，考慮高程的協(xié)同克里金法較普通的克里金法的插值精度高。

3)通過交叉驗證表明，粒子群優(yōu)化算法能在一定程度上提高普通克里金插值法和協(xié)同克里金插值法的插值精度。