邓朗妮， 黄静怡， 刘晓凤， 马晋超

(广西科技大学 土木建筑工程学院， 广西 柳州 545006)

随着现代城镇化的迅速发展和城市化水平的不断加快，可利用的原有平缓地资源非常有限，土地供给与保护耕地之间的矛盾日益凸显。为保护有限的耕地资源，城市建设用地逐渐向高处复杂地形转移。因此，在工程项目中尤其是复杂地形项目中合理的土方平衡设计方案，科学、精确的计算方法是非常必要的。

传统的土方量理论计算有很多算法，一般是方格网法[1]、断面法[2]、等高线法[3]等，其中方格网法主要适用于大面积平坦，地势起伏变化小的地形；断面法一般适合道路、隧道、管道等地面高程变化明显的带状地貌；等高线法主要适用于坡度转折较多、地面高程变化较大的较完整规则地貌。几种土方计算方法均有各自的适用性及局限性，但在数据采集上均不能全面地反映出地形表面数据特征点，计算精度不高[4,5]。

针对传统土方量的计算结果精度差、计算繁杂、效率低下等问题，本文在已成熟的地统计学空间插值法的基础上，选定较其他插值法计算精度更高的普通克里格空间插值法为理论核心计算方法，就标高和克里格权重两个关键因子对克里格插值法进行理论优化，将优化后的克里格空间插值法与未优化的克里格空间插值法进行误差对比，研究结果表明，较未优化的克里格空间插值法，优化后的克里格空间插值法理论计算精度更高[6]。

1 普通克里格空间插值法

普通克里格空间插值法[7]是地统计学空间插值法的一种，其本质是将某一区域范围内的变量作为理论基础，将变异函数作为计算载体，研究空间分布上既有随机性又有依赖性或结构性的自然科学。克里格空间插值法是结合某一范围内变量的原始数据和半方差函数的结构性，对未知采样点区域化变量进行线性无偏最优估计手段，是在样点分割间距范围内，考虑到样点与待估样点在空间上的相互位置关系和结构上的变异函数，再根据已经测量出来的样点数据，对计划估算样点进行线性无偏最优估计的一种手段。数学领域的解释就是对空间上分布的数据进行内插法线性最优无偏估计的一种方法。其具体计算方法如下[8]：

假设规定分割间距范围内变量Z(x)同时满足本征假设和二阶平稳假设，数学期望为m，且协方差函数c(h)及变异函数γ(h)存在，即

(1)

设Z(x)是一个二阶平稳随机函数，对n个位置对应有n个不同取样：Z(x1)，Z(x2)，…，Z(xn)，则x0点的估量为：

(2)

式中：λi为克里格权重系数，表示各空间样本点实际测量值Z(x)对待估值Z(x0)的贡献程度。

1.1 构造不规则三角网

设碎布点Aij(i=1,2,3,…,n，j=1,2,3,…,n)的原始高程数据为ZAij，其在平面上的投影位置为XAij。

构造不规则三角网[]首先需要根据原地形确定凸包，基于原地形数据的方格网规则网络，对原地形进行粗略的单元式规划，再分割建立成三角网。

从方格网点中取具有最小纬度碎布点的投影视为起始顶点X0，其他投影(X0除外)连接到起始顶点。从X0投影到其他投影点的线之间具有最小角度的投影被认为是第一顶点X1；当存在多个这样的投影点时，选择距起始顶点X0最远的投影作为第一顶点X1。然后将其他投影(X0，X1除外)连接到第一顶点X1。从Xk投影到线Xk-1Xk之间具有最小角度的投影被认为是第k+1个顶点，其中k=2,3,…,n-2。当第n个顶点是初始顶点X0，其中n为顶点的数量，所有顶点根据数字序列连接以形成凸多边形，其被称为凸包，如图1。所有顶点都称为凸包点[9]。

图1 凸多边形的建立

其次是基于凸包的平面不规则三角网的建立。根据初始顶点X0连接到其他凸包点以形成凸包三角形的步骤操作。随机选择一个投影并连接到它相邻的三角形的顶点。可以通过对凸包中的所有投影进行相同操作来建立三角形不规则网络。图2表示出平面不规则三角网的建立，其中洋红色点表示投影点。

图2 平面不规则三角网的建立

最后建立空间三角形不规则网络地形，空间三角形不规则网络地形可以根据平面三角形不规则网络地形和初始高程数据构造。图3表示空间三角形不规则网络地形。

图3 空间三角形不规则网络地形

1.2 克里格权重

克里格算法中最关键的步骤是对权重系数的计算，所以，克里格权重的计算必须满足以下两个条件：

(1)无偏性

若使Z*(x)为Z(xi)的无偏估计量，即

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

(3)

当E[Z*(x)]=m时，即

(4)

时，则有∑λ=1。

(2)最优性

在满足上述条件下，使待估值Z*(x)和实际测量值Z(xi)之差的平方和最小，即

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

(5)

用协方差函数可以表达为：

(6)

可根据拉格朗日乘数原理得出最小估计协方差值，令

(7)

令偏导数为0，F对λi和μ求偏导数，得到克里格方程组：

(8)

解线性方程组(8)，求出权重系数λi和拉格朗日乘数μ，代入式(2)，(5)，分别得估计值和估计方差。

在变异函数存在的条件下，协方差函数c(h)及变异函数γ(h)的关系式为：c(h)=c(0)-γ(h)。

用变异函数表示普通克里格方程组和克里格估计方差，即

(9)

用矩阵形式表示有：

kλ=D

(10)

由式(10)可得：

λ=k-1D

(11)

估计方差为：

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

(12)

2 基于普通克里格空间插值法的优化分析

2.1 优化不规则三角网

1.1中所述地形是根据固定的拓扑连接关系建立的，但实际上不匹配实际地形特性，因此需要全局优化。

平面优化是指只优化地形数据的经度和纬度，而空间优化除了经度和纬度还考虑到了高程。平面优化算法即Delaunay优化规则[10,11]，其优化思想是基于最小平面内角的最大化规则。平面优化规则仅根据平面距离而不是实际地形波动建立点的关系，并且不考虑高程信息。空间优化规则根据地形波动考虑高程数据。当地形平坦时，优化效果对于平面和空间优化都是类似的。然而，当地形明显波动时，空间优化规则可以更准确地近似于实际地形。以下以内角的标准偏差优化[12]规则进行空间优化。

首先，用两个毗邻的空间三角网的四个顶点构造四面体，如图4，四面体由四个三角形，两个原始三角形(图中的黄色和绿色三角形)和两个另外的三角形(图中的黄色和绿色三角形)组成。

然后，如果两个原始三角网格的内角标准差小于两个附加三角网格的内角，则两个原始三角网格组合成的空间三角网被认为是真实的地形。否则，两个附加的三角形被认为是真实的地形。地形的内角标准偏差按式(13)，(14)计算。

Sa=

(13)

Sb=

(14)

当Sa

图4 基于内部标准偏差的优化方法

最后，对空间三角形不规则网络地形中的所有相邻三角形进行上述步骤操作，建立优化的空间三角形不规则网络地形。

2.2 建立变异函数(结构分析理论)

克里格插值中插值数据的选择和克里格方程中系数矩阵的计算都基于变异函数。因此，必须根据优化后的空间三角形不规则网络地形，首先生成变异函数，再通过变异函数建立变异函数理论拟合模型。

(15)

接下来，建立变异函数理论拟合模型：

(1)球形模型变异函数

(16)

(2)指数模型变异函数

(17)

式中：a，C0，C为球状模型、指数模型的主要参数，其中，变程a=200，块金常数C0=2，拱高C=20；h′=lh。

图5，6所示分别为球形变异函数和指数变异函数拟合模型。

图5 球形变异函数拟合模型

图6 指数变异函数拟合模型

2.3 克里格权重系数校正

通常采用克里格插值法中的插值方差来评价插值结果的精度。实际模拟结果表明，插值误差与克里格方差成正比，因此需要校正插值结果。通过交叉验证方法对克里格权重进行系数校正，从而校正插值结果[13]，插值位置局部范围中的插值误差也用于校正插值结果[14,15]。然而，不同插值方法会导致不同的插值误差，所以下面采用克里格内插处理方法处理。

(1)设置规则网格的间距。在初始地形原始数据的所有投影中，网格划分首先设置为方格网形式，对接下来建立不规则三角网有一定的简化作用。

(2)克里格空间插值法的误差校正。在数字地面模型(Digital Terrain Model，DTM)[16]研究领域中已经提出了许多插值方法，而克里格空间插值法因其计算精度高是优选的理论计算方法。但是当地形明显变化时，插值误差通常也会随之变大，所以克里格空间插值法必须进行插值校正。

(18)

(4)克里格权重修正。在校正方法中，负权重被设置为零，并且基于标准化方法校正其他权重如下：

(19)

(20)

(21)

(22)

2.4 优化后的曲线与原始曲线比较分析

运用ANSYS有限元软件[17]建立场地模型[18,19]并进行高精度网格划分，提取出各节点坐标(高程点)，进行网格优化，对其进行内角的标准偏差优化和空间形状优化，再揉合优化后的克里格空间插值法为核心理论算法，最后运用Matlab软件[20]编程先对比分析出适合的变异函数理论模拟模型。

设部分已知高程点未知，根据Matlab编程软件模拟，首先将克里格插值球形变异函数拟合值、指数变异函数拟合值和原变异值进行对比，如图7所示。确定变异函数类型，然后通过变异函数类型的择优选取，将具有和不具有校正的克里格插值曲线模拟数值与实际数值进行对比，如图8。

图7显示了实际曲线(黑色点)，球状模型(黑色曲线)拟合曲线和指数模型(洋红色曲线)拟合曲线的离散点的三种变差函数曲线，结果显示球形模型具有比指数模型更高的拟合精度。所以选择球状曲线拟合模型为克里格插值法理论计算模型。

根据图8，未优化校正克里格插值误差范围基本集中在-5～7 m范围内，其值普遍大于范围值集中于-3.8～4 m的已优化校正克里格插值误差，可知较于未优化校正克里格插值法，具有校正的克里格插值法更有效地减少了插值算法误差。

3 土方平衡算例

案例来自广西百色某学院总体规划，总用地面积为41.93 ha。该学院场地位于百色盆地内，原场地为森林用地，地貌单元属于右江左岸Ⅰ级阶地，地貌为缓丘坡地，地形起伏较大，地面高程为103.40～180.90 m，地形坡度10°～30°。场地总体为北东较高，西南地势较低，两条呈北东-南西走向的冲沟，沟底地面高程115.40 ～117.20 m，西南侧坡脚为鱼塘，鱼塘水位约116.10 m，塘低地面高程约103.20～116.10 m，水深一般为5.0～13.0 m。

3.1 基于Matlab的开挖土方量计算

根据2.4节结论，采用具有校正的克里格插值法，运用Matlab 编程运算得各楼块开挖土方量，如表1。

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

3.2 土方调配

土方调配是指在工程施工前，对路基土石方的挖填方数量进行统计并合理分配路线的过程。设有m个挖方区Wi(i=1,2,…,m)，设计的待挖土方量分别为Si；有n个填方区Tj(j=1,2,…,n)，设计需要填方量分别为dj；从i到j平均运输距离为Cij。若用Xij表示从挖方区i运往填方区j的土方量，则在挖填平衡的条件下，土方的总运输量最小(若Cij为从i到j运输单位土方的单价，则为总运费最小) 的调运方案的数学模型为：

挖填平衡时，线性规划目标函数：

(23)

(24)

(25)

xij≥0

(26)

式(23)为目标函数，即土方总运输量最小约束函数；式(24)为挖方区的供给约束；式(25)为填方区的需求约束；式(26)为变量非负约束[21]。

当挖填不平衡时，土方调配的线性规划目标模型分别表示如下：

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

(27)

(28)

(29)

xij≥0

(30)

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

(31)

(32)

(33)

xij≥0

(34)

线性规划属运筹学范畴，从数学角度看，线性规划可以用来求解变量满足线性约束条件时，多变量线性函数的最优解，同时，线性规划也可以利用计算机技术来实现。将线性规划模型和Matlab的编程运行功能结合，可以快捷简便地得出土方调配结果，所以这里土方调配问题通过线性规划来解决。

为了求解方便，需要根据已知各挖填土方可调配量和挖填方之间的平均运距，绘制挖、填方土方量-运距表(表2)，将各挖填土方可调配量和各平均运距填入其中，作为土方调配数据依据。

表2 填挖土方量-运距表

通过挖、填土方量-运距表，进行数学表达，建立合适的线性规划模型，为Matlab编程运行提供理论支持。依据表2，则线性规划数学模型应如下：

目标函数：

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

(35)

约束条件为：

(36)

变量非负约束为：

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

(37)

根据上述目标函数和变量约束，可知有14×7个未知量，有14×7-1个独立方程，方程的解有无穷多组。而通过Matlab规划运行求解可求出一组目标函数值最小的土方调运最优解，最终，调运结果绘成土方调运路线图(图9)。

图9 土方调运路线

4 结 语

本文基于普通克里格空间插值法，就标高和克里格权重两个关键因子对克里格插值法进行理论优化，通过Matlab数值模拟，将优化后的克里格空间插值法与未优化的克里格空间插值法进行误差对比，研究结果表明：优化后的克里格空间插值法误差小于未优化的克里格空间插值法，即优化后的克里格空间插值法计算精度更高。通过以上理论优化后的 Matlab 数值模拟分析可知，空间克里格空间插值法较于传统土方计算理论具有一定的适用性、科学性和高精度性，为土方工程的计算提供了一定的参考依据，其土方平衡优化的结果也可作为竖向规划合理可行的有利依据

另外，本文建立克里格空间插值法变异函数Matlab数值模型时，只是根据传统变异函数进行了简单优化，并没有具体地将整个广西百色某学院规划区域的东—西、东南—西北、东北—西南三个方向，从而分析研究不同方向情况下的不同变异函数拟合精度问题，这些都有待进一步研究。