王桂林，向林川，孫帆

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

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

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

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

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

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

1 方法原理

1.1 普通克里金插值法

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

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

(1)

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

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

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

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

(2)

(3)

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

1.2 協(xié)同克里金法

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

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

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

(4)

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

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

(5)

(6)

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

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

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

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

(7)

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

(8)

式中：i為第i個(gè)粒子，j為第j個(gè)維度；t為更新次數(shù)；v為粒子速度，x為粒子位置；r1、r2為0到1均勻分布的隨機(jī)數(shù)；xpbest粒子自身最好位置，xgbest全局最好位置；ω為慣性權(quán)重值。

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

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

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

xd=gbest(d)×(0.5+σ)

(9)

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

1.5 預(yù)測(cè)精度評(píng)價(jià)方法

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

(10)

(11)

式中：Z(xi)與Z*(xi)分別為xi處的勘察值與估計(jì)值；N為樣本數(shù)量。

2 工程實(shí)例

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

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

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

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

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

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

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

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

2.2 空間插值模型試驗(yàn)

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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