優(yōu)化后的克里格空間插值法在土方平衡中的應(yīng)用

鄧?yán)誓荩?黃靜怡， 劉曉鳳， 馬晉超

(廣西科技大學(xué) 土木建筑工程學(xué)院， 廣西 柳州 545006)

隨著現(xiàn)代城鎮(zhèn)化的迅速發(fā)展和城市化水平的不斷加快，可利用的原有平緩地資源非常有限，土地供給與保護(hù)耕地之間的矛盾日益凸顯。為保護(hù)有限的耕地資源，城市建設(shè)用地逐漸向高處復(fù)雜地形轉(zhuǎn)移。因此，在工程項(xiàng)目中尤其是復(fù)雜地形項(xiàng)目中合理的土方平衡設(shè)計(jì)方案，科學(xué)、精確的計(jì)算方法是非常必要的。

傳統(tǒng)的土方量理論計(jì)算有很多算法，一般是方格網(wǎng)法[1]、斷面法[2]、等高線法[3]等，其中方格網(wǎng)法主要適用于大面積平坦，地勢(shì)起伏變化小的地形；斷面法一般適合道路、隧道、管道等地面高程變化明顯的帶狀地貌；等高線法主要適用于坡度轉(zhuǎn)折較多、地面高程變化較大的較完整規(guī)則地貌。幾種土方計(jì)算方法均有各自的適用性及局限性，但在數(shù)據(jù)采集上均不能全面地反映出地形表面數(shù)據(jù)特征點(diǎn)，計(jì)算精度不高[4,5]。

針對(duì)傳統(tǒng)土方量的計(jì)算結(jié)果精度差、計(jì)算繁雜、效率低下等問(wèn)題，本文在已成熟的地統(tǒng)計(jì)學(xué)空間插值法的基礎(chǔ)上，選定較其他插值法計(jì)算精度更高的普通克里格空間插值法為理論核心計(jì)算方法，就標(biāo)高和克里格權(quán)重兩個(gè)關(guān)鍵因子對(duì)克里格插值法進(jìn)行理論優(yōu)化，將優(yōu)化后的克里格空間插值法與未優(yōu)化的克里格空間插值法進(jìn)行誤差對(duì)比，研究結(jié)果表明，較未優(yōu)化的克里格空間插值法，優(yōu)化后的克里格空間插值法理論計(jì)算精度更高[6]。

1 普通克里格空間插值法

普通克里格空間插值法[7]是地統(tǒng)計(jì)學(xué)空間插值法的一種，其本質(zhì)是將某一區(qū)域范圍內(nèi)的變量作為理論基礎(chǔ)，將變異函數(shù)作為計(jì)算載體，研究空間分布上既有隨機(jī)性又有依賴性或結(jié)構(gòu)性的自然科學(xué)。克里格空間插值法是結(jié)合某一范圍內(nèi)變量的原始數(shù)據(jù)和半方差函數(shù)的結(jié)構(gòu)性，對(duì)未知采樣點(diǎn)區(qū)域化變量進(jìn)行線性無(wú)偏最優(yōu)估計(jì)手段，是在樣點(diǎn)分割間距范圍內(nèi)，考慮到樣點(diǎn)與待估樣點(diǎn)在空間上的相互位置關(guān)系和結(jié)構(gòu)上的變異函數(shù)，再根據(jù)已經(jīng)測(cè)量出來(lái)的樣點(diǎn)數(shù)據(jù)，對(duì)計(jì)劃估算樣點(diǎn)進(jìn)行線性無(wú)偏最優(yōu)估計(jì)的一種手段。數(shù)學(xué)領(lǐng)域的解釋就是對(duì)空間上分布的數(shù)據(jù)進(jìn)行內(nèi)插法線性最優(yōu)無(wú)偏估計(jì)的一種方法。其具體計(jì)算方法如下[8]：

假設(shè)規(guī)定分割間距范圍內(nèi)變量Z(x)同時(shí)滿足本征假設(shè)和二階平穩(wěn)假設(shè)，數(shù)學(xué)期望為m，且協(xié)方差函數(shù)c(h)及變異函數(shù)γ(h)存在，即

(1)

設(shè)Z(x)是一個(gè)二階平穩(wěn)隨機(jī)函數(shù)，對(duì)n個(gè)位置對(duì)應(yīng)有n個(gè)不同取樣：Z(x1)，Z(x2)，…，Z(xn)，則x0點(diǎn)的估量為：

(2)

式中：λi為克里格權(quán)重系數(shù)，表示各空間樣本點(diǎn)實(shí)際測(cè)量值Z(x)對(duì)待估值Z(x0)的貢獻(xiàn)程度。

1.1 構(gòu)造不規(guī)則三角網(wǎng)

設(shè)碎布點(diǎn)Aij(i=1,2,3,…,n，j=1,2,3,…,n)的原始高程數(shù)據(jù)為ZAij，其在平面上的投影位置為XAij。

構(gòu)造不規(guī)則三角網(wǎng)[]首先需要根據(jù)原地形確定凸包，基于原地形數(shù)據(jù)的方格網(wǎng)規(guī)則網(wǎng)絡(luò)，對(duì)原地形進(jìn)行粗略的單元式規(guī)劃，再分割建立成三角網(wǎng)。

從方格網(wǎng)點(diǎn)中取具有最小緯度碎布點(diǎn)的投影視為起始頂點(diǎn)X0，其他投影(X0除外)連接到起始頂點(diǎn)。從X0投影到其他投影點(diǎn)的線之間具有最小角度的投影被認(rèn)為是第一頂點(diǎn)X1；當(dāng)存在多個(gè)這樣的投影點(diǎn)時(shí)，選擇距起始頂點(diǎn)X0最遠(yuǎn)的投影作為第一頂點(diǎn)X1。然后將其他投影(X0，X1除外)連接到第一頂點(diǎn)X1。從Xk投影到線Xk-1Xk之間具有最小角度的投影被認(rèn)為是第k+1個(gè)頂點(diǎn)，其中k=2,3,…,n-2。當(dāng)?shù)趎個(gè)頂點(diǎn)是初始頂點(diǎn)X0，其中n為頂點(diǎn)的數(shù)量，所有頂點(diǎn)根據(jù)數(shù)字序列連接以形成凸多邊形，其被稱為凸包，如圖1。所有頂點(diǎn)都稱為凸包點(diǎn)[9]。

圖1 凸多邊形的建立

其次是基于凸包的平面不規(guī)則三角網(wǎng)的建立。根據(jù)初始頂點(diǎn)X0連接到其他凸包點(diǎn)以形成凸包三角形的步驟操作。隨機(jī)選擇一個(gè)投影并連接到它相鄰的三角形的頂點(diǎn)?？梢酝ㄟ^(guò)對(duì)凸包中的所有投影進(jìn)行相同操作來(lái)建立三角形不規(guī)則網(wǎng)絡(luò)。圖2表示出平面不規(guī)則三角網(wǎng)的建立，其中洋紅色點(diǎn)表示投影點(diǎn)。

圖2 平面不規(guī)則三角網(wǎng)的建立

最后建立空間三角形不規(guī)則網(wǎng)絡(luò)地形，空間三角形不規(guī)則網(wǎng)絡(luò)地形可以根據(jù)平面三角形不規(guī)則網(wǎng)絡(luò)地形和初始高程數(shù)據(jù)構(gòu)造。圖3表示空間三角形不規(guī)則網(wǎng)絡(luò)地形。

圖3 空間三角形不規(guī)則網(wǎng)絡(luò)地形

1.2 克里格權(quán)重

克里格算法中最關(guān)鍵的步驟是對(duì)權(quán)重系數(shù)的計(jì)算，所以，克里格權(quán)重的計(jì)算必須滿足以下兩個(gè)條件：

(1)無(wú)偏性

若使Z*(x)為Z(xi)的無(wú)偏估計(jì)量，即

E[Z*(x)]=E[Z(x)]

(3)

當(dāng)E[Z*(x)]=m時(shí)，即

(4)

時(shí)，則有∑λ=1。

(2)最優(yōu)性

在滿足上述條件下，使待估值Z*(x)和實(shí)際測(cè)量值Z(xi)之差的平方和最小，即

σ2=E[Z(x)-Z*(x)]2=

(5)

用協(xié)方差函數(shù)可以表達(dá)為：

(6)

可根據(jù)拉格朗日乘數(shù)原理得出最小估計(jì)協(xié)方差值，令

(7)

令偏導(dǎo)數(shù)為0，F(xiàn)對(duì)λi和μ求偏導(dǎo)數(shù)，得到克里格方程組：

(8)

解線性方程組(8)，求出權(quán)重系數(shù)λi和拉格朗日乘數(shù)μ，代入式(2)，(5)，分別得估計(jì)值和估計(jì)方差。

在變異函數(shù)存在的條件下，協(xié)方差函數(shù)c(h)及變異函數(shù)γ(h)的關(guān)系式為：c(h)=c(0)-γ(h)。

用變異函數(shù)表示普通克里格方程組和克里格估計(jì)方差，即

(9)

用矩陣形式表示有：

kλ=D

(10)

由式(10)可得：

λ=k-1D

(11)

估計(jì)方差為：

σ2=γ(x,x)-λTD

(12)

2 基于普通克里格空間插值法的優(yōu)化分析

2.1 優(yōu)化不規(guī)則三角網(wǎng)

1.1中所述地形是根據(jù)固定的拓?fù)溥B接關(guān)系建立的，但實(shí)際上不匹配實(shí)際地形特性，因此需要全局優(yōu)化。

平面優(yōu)化是指只優(yōu)化地形數(shù)據(jù)的經(jīng)度和緯度，而空間優(yōu)化除了經(jīng)度和緯度還考慮到了高程。平面優(yōu)化算法即Delaunay優(yōu)化規(guī)則[10,11]，其優(yōu)化思想是基于最小平面內(nèi)角的最大化規(guī)則。平面優(yōu)化規(guī)則僅根據(jù)平面距離而不是實(shí)際地形波動(dòng)建立點(diǎn)的關(guān)系，并且不考慮高程信息?？臻g優(yōu)化規(guī)則根據(jù)地形波動(dòng)考慮高程數(shù)據(jù)。當(dāng)?shù)匦纹教箷r(shí)，優(yōu)化效果對(duì)于平面和空間優(yōu)化都是類似的。然而，當(dāng)?shù)匦蚊黠@波動(dòng)時(shí)，空間優(yōu)化規(guī)則可以更準(zhǔn)確地近似于實(shí)際地形。以下以內(nèi)角的標(biāo)準(zhǔn)偏差優(yōu)化[12]規(guī)則進(jìn)行空間優(yōu)化。

首先，用兩個(gè)毗鄰的空間三角網(wǎng)的四個(gè)頂點(diǎn)構(gòu)造四面體，如圖4，四面體由四個(gè)三角形，兩個(gè)原始三角形(圖中的黃色和綠色三角形)和兩個(gè)另外的三角形(圖中的黃色和綠色三角形)組成。

然后，如果兩個(gè)原始三角網(wǎng)格的內(nèi)角標(biāo)準(zhǔn)差小于兩個(gè)附加三角網(wǎng)格的內(nèi)角，則兩個(gè)原始三角網(wǎng)格組合成的空間三角網(wǎng)被認(rèn)為是真實(shí)的地形。否則，兩個(gè)附加的三角形被認(rèn)為是真實(shí)的地形。地形的內(nèi)角標(biāo)準(zhǔn)偏差按式(13)，(14)計(jì)算。

Sa=

(13)

Sb=

(14)

當(dāng)Sa

圖4 基于內(nèi)部標(biāo)準(zhǔn)偏差的優(yōu)化方法

最后，對(duì)空間三角形不規(guī)則網(wǎng)絡(luò)地形中的所有相鄰三角形進(jìn)行上述步驟操作，建立優(yōu)化的空間三角形不規(guī)則網(wǎng)絡(luò)地形。

2.2 建立變異函數(shù)(結(jié)構(gòu)分析理論)

克里格插值中插值數(shù)據(jù)的選擇和克里格方程中系數(shù)矩陣的計(jì)算都基于變異函數(shù)。因此，必須根據(jù)優(yōu)化后的空間三角形不規(guī)則網(wǎng)絡(luò)地形，首先生成變異函數(shù)，再通過(guò)變異函數(shù)建立變異函數(shù)理論擬合模型。

(15)

接下來(lái)，建立變異函數(shù)理論擬合模型：

(1)球形模型變異函數(shù)

(16)

(2)指數(shù)模型變異函數(shù)

(17)

式中：a，C0，C為球狀模型、指數(shù)模型的主要參數(shù)，其中，變程a=200，塊金常數(shù)C0=2，拱高C=20；h′=lh。

圖5，6所示分別為球形變異函數(shù)和指數(shù)變異函數(shù)擬合模型。

圖5 球形變異函數(shù)擬合模型

圖6 指數(shù)變異函數(shù)擬合模型

2.3 克里格權(quán)重系數(shù)校正

通常采用克里格插值法中的插值方差來(lái)評(píng)價(jià)插值結(jié)果的精度。實(shí)際模擬結(jié)果表明，插值誤差與克里格方差成正比，因此需要校正插值結(jié)果。通過(guò)交叉驗(yàn)證方法對(duì)克里格權(quán)重進(jìn)行系數(shù)校正，從而校正插值結(jié)果[13]，插值位置局部范圍中的插值誤差也用于校正插值結(jié)果[14,15]。然而，不同插值方法會(huì)導(dǎo)致不同的插值誤差，所以下面采用克里格內(nèi)插處理方法處理。

(1)設(shè)置規(guī)則網(wǎng)格的間距。在初始地形原始數(shù)據(jù)的所有投影中，網(wǎng)格劃分首先設(shè)置為方格網(wǎng)形式，對(duì)接下來(lái)建立不規(guī)則三角網(wǎng)有一定的簡(jiǎn)化作用。

(2)克里格空間插值法的誤差校正。在數(shù)字地面模型(Digital Terrain Model，DTM)[16]研究領(lǐng)域中已經(jīng)提出了許多插值方法，而克里格空間插值法因其計(jì)算精度高是優(yōu)選的理論計(jì)算方法。但是當(dāng)?shù)匦蚊黠@變化時(shí)，插值誤差通常也會(huì)隨之變大，所以克里格空間插值法必須進(jìn)行插值校正。

(18)

(4)克里格權(quán)重修正。在校正方法中，負(fù)權(quán)重被設(shè)置為零，并且基于標(biāo)準(zhǔn)化方法校正其他權(quán)重如下：

(19)

(20)

(21)

(22)

2.4 優(yōu)化后的曲線與原始曲線比較分析

運(yùn)用ANSYS有限元軟件[17]建立場(chǎng)地模型[18,19]并進(jìn)行高精度網(wǎng)格劃分，提取出各節(jié)點(diǎn)坐標(biāo)(高程點(diǎn))，進(jìn)行網(wǎng)格優(yōu)化，對(duì)其進(jìn)行內(nèi)角的標(biāo)準(zhǔn)偏差優(yōu)化和空間形狀優(yōu)化，再揉合優(yōu)化后的克里格空間插值法為核心理論算法，最后運(yùn)用Matlab軟件[20]編程先對(duì)比分析出適合的變異函數(shù)理論模擬模型。

設(shè)部分已知高程點(diǎn)未知，根據(jù)Matlab編程軟件模擬，首先將克里格插值球形變異函數(shù)擬合值、指數(shù)變異函數(shù)擬合值和原變異值進(jìn)行對(duì)比，如圖7所示。確定變異函數(shù)類型，然后通過(guò)變異函數(shù)類型的擇優(yōu)選取，將具有和不具有校正的克里格插值曲線模擬數(shù)值與實(shí)際數(shù)值進(jìn)行對(duì)比，如圖8。

圖7顯示了實(shí)際曲線(黑色點(diǎn))，球狀模型(黑色曲線)擬合曲線和指數(shù)模型(洋紅色曲線)擬合曲線的離散點(diǎn)的三種變差函數(shù)曲線，結(jié)果顯示球形模型具有比指數(shù)模型更高的擬合精度。所以選擇球狀曲線擬合模型為克里格插值法理論計(jì)算模型。

根據(jù)圖8，未優(yōu)化校正克里格插值誤差范圍基本集中在-5～7 m范圍內(nèi)，其值普遍大于范圍值集中于-3.8～4 m的已優(yōu)化校正克里格插值誤差，可知較于未優(yōu)化校正克里格插值法，具有校正的克里格插值法更有效地減少了插值算法誤差。

3 土方平衡算例

案例來(lái)自廣西百色某學(xué)院總體規(guī)劃，總用地面積為41.93 ha。該學(xué)院場(chǎng)地位于百色盆地內(nèi)，原場(chǎng)地為森林用地，地貌單元屬于右江左岸Ⅰ級(jí)階地，地貌為緩丘坡地，地形起伏較大，地面高程為103.40～180.90 m，地形坡度10°～30°。場(chǎng)地總體為北東較高，西南地勢(shì)較低，兩條呈北東-南西走向的沖溝，溝底地面高程115.40 ～117.20 m，西南側(cè)坡腳為魚塘，魚塘水位約116.10 m，塘低地面高程約103.20～116.10 m，水深一般為5.0～13.0 m。

3.1 基于Matlab的開(kāi)挖土方量計(jì)算

根據(jù)2.4節(jié)結(jié)論，采用具有校正的克里格插值法，運(yùn)用Matlab 編程運(yùn)算得各樓塊開(kāi)挖土方量，如表1。

表1 1#～20#樓土方挖填量

3.2 土方調(diào)配

土方調(diào)配是指在工程施工前，對(duì)路基土石方的挖填方數(shù)量進(jìn)行統(tǒng)計(jì)并合理分配路線的過(guò)程。設(shè)有m個(gè)挖方區(qū)Wi(i=1,2,…,m)，設(shè)計(jì)的待挖土方量分別為Si；有n個(gè)填方區(qū)Tj(j=1,2,…,n)，設(shè)計(jì)需要填方量分別為dj；從i到j(luò)平均運(yùn)輸距離為Cij。若用Xij表示從挖方區(qū)i運(yùn)往填方區(qū)j的土方量，則在挖填平衡的條件下，土方的總運(yùn)輸量最小(若Cij為從i到j(luò)運(yùn)輸單位土方的單價(jià)，則為總運(yùn)費(fèi)最小) 的調(diào)運(yùn)方案的數(shù)學(xué)模型為：

挖填平衡時(shí)，線性規(guī)劃目標(biāo)函數(shù)：

(23)

(24)

(25)

xij≥0

(26)

式(23)為目標(biāo)函數(shù)，即土方總運(yùn)輸量最小約束函數(shù)；式(24)為挖方區(qū)的供給約束；式(25)為填方區(qū)的需求約束；式(26)為變量非負(fù)約束[21]。

當(dāng)挖填不平衡時(shí)，土方調(diào)配的線性規(guī)劃目標(biāo)模型分別表示如下：

(1)挖方量大于填方量：

(27)

(28)

(29)

xij≥0

(30)

(2)挖方量小于填方量：

(31)

(32)

(33)

xij≥0

(34)

線性規(guī)劃屬運(yùn)籌學(xué)范疇，從數(shù)學(xué)角度看，線性規(guī)劃可以用來(lái)求解變量滿足線性約束條件時(shí)，多變量線性函數(shù)的最優(yōu)解，同時(shí)，線性規(guī)劃也可以利用計(jì)算機(jī)技術(shù)來(lái)實(shí)現(xiàn)。將線性規(guī)劃模型和Matlab的編程運(yùn)行功能結(jié)合，可以快捷簡(jiǎn)便地得出土方調(diào)配結(jié)果，所以這里土方調(diào)配問(wèn)題通過(guò)線性規(guī)劃來(lái)解決。

為了求解方便，需要根據(jù)已知各挖填土方可調(diào)配量和挖填方之間的平均運(yùn)距，繪制挖、填方土方量-運(yùn)距表(表2)，將各挖填土方可調(diào)配量和各平均運(yùn)距填入其中，作為土方調(diào)配數(shù)據(jù)依據(jù)。

表2 填挖土方量-運(yùn)距表

通過(guò)挖、填土方量-運(yùn)距表，進(jìn)行數(shù)學(xué)表達(dá)，建立合適的線性規(guī)劃模型，為Matlab編程運(yùn)行提供理論支持。依據(jù)表2，則線性規(guī)劃數(shù)學(xué)模型應(yīng)如下：

目標(biāo)函數(shù)：

MinZ=189.4X11+343.5X12+349.3X13+…+278.7X14,5+538X14,6+98X14,7

(35)

約束條件為：

(36)

變量非負(fù)約束為：

Xij≥0，i=1,2,3,4，j=1,2,3

(37)

根據(jù)上述目標(biāo)函數(shù)和變量約束，可知有14×7個(gè)未知量，有14×7-1個(gè)獨(dú)立方程，方程的解有無(wú)窮多組。而通過(guò)Matlab規(guī)劃運(yùn)行求解可求出一組目標(biāo)函數(shù)值最小的土方調(diào)運(yùn)最優(yōu)解，最終，調(diào)運(yùn)結(jié)果繪成土方調(diào)運(yùn)路線圖(圖9)。

圖9 土方調(diào)運(yùn)路線

4 結(jié) 語(yǔ)

本文基于普通克里格空間插值法，就標(biāo)高和克里格權(quán)重兩個(gè)關(guān)鍵因子對(duì)克里格插值法進(jìn)行理論優(yōu)化，通過(guò)Matlab數(shù)值模擬，將優(yōu)化后的克里格空間插值法與未優(yōu)化的克里格空間插值法進(jìn)行誤差對(duì)比，研究結(jié)果表明：優(yōu)化后的克里格空間插值法誤差小于未優(yōu)化的克里格空間插值法，即優(yōu)化后的克里格空間插值法計(jì)算精度更高。通過(guò)以上理論優(yōu)化后的 Matlab 數(shù)值模擬分析可知，空間克里格空間插值法較于傳統(tǒng)土方計(jì)算理論具有一定的適用性、科學(xué)性和高精度性，為土方工程的計(jì)算提供了一定的參考依據(jù)，其土方平衡優(yōu)化的結(jié)果也可作為豎向規(guī)劃合理可行的有利依據(jù)

另外，本文建立克里格空間插值法變異函數(shù)Matlab數(shù)值模型時(shí)，只是根據(jù)傳統(tǒng)變異函數(shù)進(jìn)行了簡(jiǎn)單優(yōu)化，并沒(méi)有具體地將整個(gè)廣西百色某學(xué)院規(guī)劃區(qū)域的東—西、東南—西北、東北—西南三個(gè)方向，從而分析研究不同方向情況下的不同變異函數(shù)擬合精度問(wèn)題，這些都有待進(jìn)一步研究。