基于Delaunay三角網(wǎng)法的土石方量計(jì)算精度探討

候漢霜

(深圳市長(zhǎng)勘勘察設(shè)計(jì)有限公司，廣東 深圳 518003)

1 引 言

土石方量的計(jì)算結(jié)果關(guān)系到工程項(xiàng)目的各方利益，目前我們常采用的土石方量計(jì)算方法有斷面法、方格網(wǎng)法、等高線法、平均高程法和DTM(數(shù)字地面模型)法等。其中，每一種計(jì)算方法都有其適用條件和適用范圍，斷面法計(jì)算土石方量適合在線狀地形使用，比如道路、河流渠道等；方格網(wǎng)法在比較平坦的地區(qū)或者地形起伏不大的場(chǎng)地上更經(jīng)濟(jì)實(shí)用；而在地形起伏較大、精度要求高的區(qū)域則需使用DTM法來(lái)計(jì)算土石方量較為合理。DTM法由于考慮了地形特征(采用不規(guī)則三角網(wǎng)TIN法建模)，并以原始測(cè)量采集的數(shù)據(jù)作結(jié)點(diǎn)的三角形為最小計(jì)算單元，從理論上來(lái)說(shuō)是計(jì)算土石方量最為精確的方法。而通過(guò)建立Delaunay三角網(wǎng)保證得到的是最優(yōu)的TIN模型，一般情況下是計(jì)算土石方量的最佳選擇。目前，主流的Delaunay三角網(wǎng)生成算法有三種：分割合并法、逐點(diǎn)插入法和三角網(wǎng)生長(zhǎng)法。本文采用三角網(wǎng)生長(zhǎng)法編寫基于Delaunay三角網(wǎng)的土石方計(jì)算程序，探究Delaunay三角網(wǎng)法土石方計(jì)算的方案與南方CASS軟件中傳統(tǒng)三角網(wǎng)法方案對(duì)比在土方量計(jì)算精度的偏差。

2 基于Delaunay三角網(wǎng)法的土石方量計(jì)算方法

土石方量的計(jì)算是通過(guò)地形表面模型和設(shè)計(jì)表面模型進(jìn)行疊加，求出疊加后所夾空間區(qū)域的體積從而得到土石方量。這里需要分別建立地形表面TIN模型和設(shè)計(jì)表面TIN模型。

2.1 地形表面TIN模型的構(gòu)建

要建立地形表面TIN模型必須先確定TIN模型中三角形的形成法則，即三角剖分準(zhǔn)則。常用的三角形剖分準(zhǔn)則為以下6種：①?gòu)埥亲畲鬁?zhǔn)則：一點(diǎn)到基線邊的張角是最大的；②空外接圓準(zhǔn)則：就是TIN中過(guò)每個(gè)三角形的外接圓均不包含點(diǎn)集除了該三角形的頂點(diǎn)之外的任何其他點(diǎn)；③最大最小角準(zhǔn)則：每?jī)上噜徣切嗡纬傻耐顾倪呅沃?，這兩個(gè)三角形的最小內(nèi)角一定大于交換該凸四邊形對(duì)角線后新形成的兩三角形的最小內(nèi)角；④最短距離和準(zhǔn)則：即一點(diǎn)到基線邊兩個(gè)端點(diǎn)的距離之和是最小的；⑤對(duì)角線準(zhǔn)則：若是兩三角形組成凸四邊形的兩條對(duì)角線之比超過(guò)限定閾值時(shí)，則需要對(duì)三角形進(jìn)行優(yōu)化；⑥面積比準(zhǔn)則：就是三角形的面積與周長(zhǎng)平方之比最小(或者三角形內(nèi)切圓的面積與三角形的面積之比最小)。在理論上已經(jīng)證明了張角最大準(zhǔn)則、空外接圓準(zhǔn)則和最大最小角準(zhǔn)則是等價(jià)的，但其他的則不是。其中對(duì)角線準(zhǔn)則受主觀因素影響，已經(jīng)很少用到。

通常將在空外接圓準(zhǔn)則、最大最小角準(zhǔn)則下進(jìn)行的三角剖分稱為Delaunay三角剖分(簡(jiǎn)稱DT)。并且空外接圓準(zhǔn)則與最小角最大準(zhǔn)則是Delaunay三角網(wǎng)的兩個(gè)基本性質(zhì)。實(shí)際上，任何三角剖分準(zhǔn)則下得到的TIN，只要利用LOP法則(局部?jī)?yōu)化過(guò)程，Local optimal procedure，Lawson于1977年提出)來(lái)進(jìn)行優(yōu)化處理，都能得到一個(gè)唯一的DT三角網(wǎng)。

2.2 設(shè)計(jì)表面TIN模型的建立

在各種土方工程中，設(shè)計(jì)表面可以圖形的方式給出，也可以數(shù)學(xué)函數(shù)的形式給出，如場(chǎng)地平整平面，整個(gè)區(qū)域可有不同的表達(dá)式，例如道路設(shè)計(jì)表面中的左右邊坡、路面等。

若設(shè)計(jì)表面以數(shù)學(xué)表達(dá)式給出，則地形表面TIN模型和設(shè)計(jì)表面模型TIN在平面位置上具有相同的三角網(wǎng)形，而在同一平面位置上地面模型與設(shè)計(jì)模型的高程值不同。但是設(shè)計(jì)表面TIN模型的高程可由模型的三角形頂點(diǎn)坐標(biāo)代入到設(shè)計(jì)表面函數(shù)計(jì)算得到。設(shè)地形表面TIN模型三角形頂點(diǎn)P坐標(biāo)為(xp，yp)，設(shè)計(jì)表面函數(shù)為HD=g(x，y)，則P點(diǎn)的設(shè)計(jì)高程為Hd(P)=g(xp，yp)。根據(jù)設(shè)計(jì)表面函數(shù)求出每個(gè)地形表面三角形頂點(diǎn)所對(duì)應(yīng)的設(shè)計(jì)高程，則可建立相應(yīng)的設(shè)計(jì)表面TIN模型。

若設(shè)計(jì)表面以圖件形式或散點(diǎn)形式給出，則可先對(duì)設(shè)計(jì)表面建立TIN。但由于地形表面模型和設(shè)計(jì)表面模型具有不同的三角網(wǎng)形結(jié)構(gòu)。這時(shí)要進(jìn)行數(shù)據(jù)加密處理，即把設(shè)計(jì)表面中的點(diǎn)內(nèi)插到地形表面上，并求出設(shè)計(jì)表面散點(diǎn)對(duì)應(yīng)的地面高程；同時(shí)，將地形表面上的點(diǎn)內(nèi)插到設(shè)計(jì)表面中，并求出地形點(diǎn)對(duì)應(yīng)的設(shè)計(jì)高程。另外，不同的TIN模型上，三角形邊相交處的坐標(biāo)和高程也相應(yīng)求出。這樣就建立了具有相同結(jié)構(gòu)的兩個(gè)TIN模型。

2.3 土石方量的計(jì)算

得到地形表面TIN模型和設(shè)計(jì)表面模型TIN后，通過(guò)疊加產(chǎn)生地形表面與設(shè)計(jì)表面高程差異，從而在三維空間上形成不同結(jié)構(gòu)的棱柱體。土石方量計(jì)算其實(shí)就是計(jì)算三角棱的體積，通過(guò)計(jì)算三角網(wǎng)中每個(gè)三棱錐柱的填挖方量，再求和得到整個(gè)三角網(wǎng)的總填挖土石方量。

根據(jù)三角形三個(gè)頂點(diǎn)地形高程和設(shè)計(jì)高程情況，每個(gè)三角形區(qū)域分為全填全挖或有填有挖兩種情況考慮。當(dāng)三角形三個(gè)頂點(diǎn)的高程都大于(或小于)設(shè)計(jì)高程。如圖1(a)所示。

(1)

式(1)中：S為三角形投影到水平面的面積；H1，H2，H3分別為三角形三個(gè)頂點(diǎn)與設(shè)計(jì)面的高差(取絕對(duì)值)。

當(dāng)三角形三個(gè)頂點(diǎn)的地形高程值既有大于也有小于設(shè)計(jì)高程。這時(shí)設(shè)計(jì)面將三角形分成兩部分，一部分為底面是三角形的錐體，另一部分是底面為四邊形的楔體。如圖1(b)所示。錐體部分和楔體部分的體積計(jì)算公式分別為：

圖1 三角棱柱的體積計(jì)算

(2)

(3)

式(2)和(3)中：S為三角形投影到水平面的面積；H1，H2，H3為三角形三個(gè)頂點(diǎn)與設(shè)計(jì)面的高差(取絕對(duì)值)，其中H1恒指錐體頂點(diǎn)與設(shè)計(jì)面的高差(取絕對(duì)值)。另外，有填方也有挖方在實(shí)際計(jì)算過(guò)程中主要分為6種情況考慮，即：①H1<0，H2>0，H3>0；②H1>0，H2<0，H3>0；③H1>0，H2>0，H3<0；④H1>0，H2<0，H3<0；⑤H1<0，H2>0，H3<0；⑥H1<0，H2<0，H3>0。

3 程序?qū)崿F(xiàn)與實(shí)驗(yàn)分析

3.1 程序?qū)崿F(xiàn)

為了分析Delaunay三角網(wǎng)法計(jì)算土石方的效果，筆者采用三角網(wǎng)生長(zhǎng)算法編寫基于Delaunay三角網(wǎng)法計(jì)算土石方的程序?？紤]到便于對(duì)所剖分的離散點(diǎn)集和形成的三角形進(jìn)行管理，程序設(shè)計(jì)3個(gè)單向鏈表類。它們分別是：①點(diǎn)表類(CPointList)，用于管理的離散點(diǎn)集。參數(shù)有點(diǎn)名(PointName)、點(diǎn)號(hào)(indexP)、X坐標(biāo)(x)、Y坐標(biāo)(y)和Z坐標(biāo)(z)，還有該點(diǎn)是否為封閉點(diǎn)(closePoint)以及一點(diǎn)所在的圓環(huán)號(hào)(RingNum)。②邊表類(CEdgeList)，用于管理三角剖分后生成的每一條邊。類的參數(shù)有邊號(hào)(indexE)、起始點(diǎn)(startPt)、終止點(diǎn)(endPt)和該邊的已被幾個(gè)三角形共用(useCount)。其中起始點(diǎn)(startPt)和終止點(diǎn)(endPt)均屬于點(diǎn)表類(CPointList)。③三角形表類(CTriList)，用于管理三角剖分后生成的每一個(gè)三角形。類的參數(shù)有三角形號(hào)(indexT)、三角形的L0邊(L0)、三角形的L1邊(L1)和三角形的L2邊(L2)。其中三角形的L0邊(L0)、三角形的L1邊(L1)和三角形的L2邊(L2)均屬于邊表類(CEdgeList)。其程序?qū)崿F(xiàn)過(guò)程如下：

(1)創(chuàng)建初始三角形

在數(shù)據(jù)集中找到距幾何中心最近的點(diǎn)A，再找到距A點(diǎn)最近的點(diǎn)B，由點(diǎn)A、B組成第一條基線AB。然后對(duì)數(shù)據(jù)集進(jìn)行搜索，找到與AB邊有最大夾角的點(diǎn)C，將A、B、C三點(diǎn)按逆時(shí)針鏈接形成初始三角形，并且給邊AB(L0邊)、BC(L1邊)、CA(L2邊)的擴(kuò)展次數(shù)(useCount)賦為1。并規(guī)定以后每生成一個(gè)三角形，都按照逆時(shí)針?lè)较蜴溄尤齻€(gè)頂點(diǎn)。

(2)擴(kuò)展邊以形成整個(gè)區(qū)域的無(wú)約束Delaunay三角網(wǎng)

有了初始三角形，后面的工作就是開(kāi)始循環(huán)逐步擴(kuò)展生成三角網(wǎng)，設(shè)從第triNum個(gè)三角形開(kāi)始，分別對(duì)三角形triNum的三邊L0、L1、L2做如下處理(以L0邊為例)：

①查找可用的頂點(diǎn)。要擴(kuò)展生成下一個(gè)三角形，下一個(gè)點(diǎn)必定不是三角形triNum上的點(diǎn)，而且也不是與C點(diǎn)在同一側(cè)。

②找出可用點(diǎn)中與邊L0形成夾角最大的那個(gè)點(diǎn)。然后與L0兩端點(diǎn)連接完成三角形的一次擴(kuò)展，并給L0邊的useCount值+1。

分別處理L1、L2邊，并以此循環(huán)，直到三角形的個(gè)數(shù)不再增長(zhǎng)，則完成整個(gè)數(shù)據(jù)域的無(wú)約束Delaunay三角網(wǎng)生成。

(3)將生成的點(diǎn)表、邊表、三角形表傳遞到繪圖類、土方計(jì)算類中。

(4)繪圖類中將每一條邊繪制，連接成網(wǎng)。

(5)土方計(jì)算類中，對(duì)每個(gè)三角形進(jìn)行計(jì)算，最后計(jì)算所有三角形的填挖方總和。這里主要有三種情況：僅有填方、僅有挖方和既有填方也有挖方。其中既有填方也有挖方又細(xì)分六種情況考慮。

3.2 土石方量計(jì)算結(jié)果比較分析

為了說(shuō)明Delaunay三角網(wǎng)法土石方量計(jì)算本文方法精度，進(jìn)行了仿真測(cè)試精度分析和實(shí)際測(cè)量精度分析。

(1)仿真測(cè)試精度分析

利用曲面函數(shù)構(gòu)造模型進(jìn)行測(cè)試(圖2)，首先對(duì)該函數(shù)進(jìn)行積分計(jì)算(本處指定z=0，a=16，b=16，-50≤x≤50，-50≤y≤50)。因?yàn)檫@里指定a=b，所以該曲面是對(duì)稱的，當(dāng)|x|>|y|時(shí)z<0，即此時(shí)為填方，反之為挖方，且填挖方量相等。而|x|=|y|將區(qū)域分成4個(gè)部分，4個(gè)部分的計(jì)算結(jié)果絕對(duì)值相等：

圖2 雙曲拋物面示意圖

V總=4×V=4×8138.021=32552.084

圖3 劃定雙曲拋物面的測(cè)試區(qū)域范圍

圖4 函數(shù)的三維模型視圖

在CASS7.0中生成三角網(wǎng)，并用DTM法計(jì)算土方，如圖5所示。由圖5可以看到，土石方總量=15632.268+15659.817=31292.085，與積分計(jì)算得到的結(jié)果(32552.084)有3.87%(1259.999方)的偏差。原因在于CASS7.0中三角網(wǎng)不是嚴(yán)格的Delaunay三角網(wǎng)，另外CASS生成的三角網(wǎng)的最外層邊連成的多邊形不是嚴(yán)格凸的。

圖5 在CASS7.0中用DTM法計(jì)算1 009個(gè)點(diǎn)的土石方量

同樣采用上面生成的1 009個(gè)點(diǎn)數(shù)據(jù)文件，在基于Delaunay三角網(wǎng)法程序中計(jì)算土石方量，如圖6所示。由圖6可以看到，這里土石方總量=16245.10+16274.10=32519.20，與積分計(jì)算得到的結(jié)果(32552.084)僅有0.10%(32.884方)的偏差。

圖6 Delaunay三角網(wǎng)算法下計(jì)算1009個(gè)點(diǎn)的土石方量

由于該函數(shù)對(duì)應(yīng)的曲面是一個(gè)光滑曲面，參與計(jì)算的隨機(jī)點(diǎn)越多，計(jì)算得到的土方量與積分結(jié)果越接近。為了證明點(diǎn)數(shù)與計(jì)算結(jié)果的精確性的關(guān)系，筆者還分別測(cè)試隨機(jī)生成109個(gè)點(diǎn)、509個(gè)點(diǎn)和 2 009個(gè)點(diǎn)計(jì)算的土石方量作為比較的依據(jù)(圖7)：

圖7 測(cè)點(diǎn)數(shù)量與誤差關(guān)系圖

①109個(gè)隨機(jī)點(diǎn)的測(cè)試結(jié)果為：土石方總量=15223.76+15623.66=30847.42，與積分計(jì)算得到的結(jié)果(32552.084)偏差5.24%(1704.664方)。

②509個(gè)隨機(jī)點(diǎn)的測(cè)試結(jié)果為：土石方總量=16067.79+16060.36=32128.15，與積分計(jì)算得到的結(jié)果(32552.084)偏差1.30%(423.934方)。

③2 009個(gè)隨機(jī)點(diǎn)的測(cè)試結(jié)果為：土石方總量=16265.75+16299.50=32565.25，與積分計(jì)算得到的結(jié)果(32552.084)僅偏差了0.04%(13.166方)。

通過(guò)這幾次測(cè)試可以看到，同一個(gè)測(cè)試區(qū)域內(nèi)，隨著離散點(diǎn)數(shù)量的增加，計(jì)算得到的土石方量與積分得到的土石方量誤差也越來(lái)越小。

(2)實(shí)際數(shù)據(jù)精度分析

為了說(shuō)明基于Delaunay三角網(wǎng)法計(jì)算土石方量本文方法精度，筆者外業(yè)采用測(cè)繪儀器在野外選取一片地形稍有起伏的區(qū)域?qū)嵉夭杉?1個(gè)點(diǎn)位坐標(biāo)數(shù)據(jù)，如表1所示。

某區(qū)域?qū)崪y(cè)數(shù)據(jù) 表1

內(nèi)業(yè)利用南方CASS7.0生成DTM后計(jì)算的土石方量，其結(jié)果如圖8所示，本文方法計(jì)算的土石方量如圖9所示。

圖8 CASS7.0的DTM法生成的三角網(wǎng)

圖9 Delaunay三角網(wǎng)法生成的三角網(wǎng)

由圖8、圖9比較看出：基于Delaunay三角網(wǎng)法計(jì)算的挖方量與南方CASS7.0的不同。原因則是由于三角網(wǎng)的網(wǎng)形不相同，CASS7.0生成的三角網(wǎng)不符合Delaunay三角剖分準(zhǔn)則(圖8中畫圈的四邊形不滿足空外接圓準(zhǔn)則，需要交換對(duì)角線)。又根據(jù)前文(圖5、圖6對(duì)比)證明，本次實(shí)際數(shù)據(jù)測(cè)試基于Delaunay三角網(wǎng)法計(jì)算的精度更可靠。

4 結(jié) 語(yǔ)

本文通過(guò)數(shù)據(jù)測(cè)試表明，基于Delaunay三角網(wǎng)法計(jì)算得到的土石方量結(jié)果相比傳統(tǒng)三角網(wǎng)計(jì)算的土石方量有較大精度上的提高，尤其在采樣的數(shù)據(jù)點(diǎn)不是很密集的時(shí)候，表現(xiàn)更為明顯。因此基于Delaunay三角網(wǎng)的TIN模型更接近實(shí)際地面模型，計(jì)算出的土石方量準(zhǔn)確性更高。