基于三維激光掃描數據的不規則實體表面積和體積計算方法

郭景仁， 王艷林， 于 蕾

(中國人民解放軍61175部隊， 山東淄博 255020)

地面激光掃描，簡稱TLS(Terrestrial laser scanning),是一種允許數字采集并通過點云在三維空間中以數字形式再現真實物體的技術.激光掃描操作的主要原理是基于作為目標本身的表面反射的激光脈沖傳輸到目標和分析返回脈沖[1-2].地面三維激光掃描儀的發展和應用，為人們獲取豐富的局部地面空間信息提供了一種全新的技術手段，它是一種非接觸式主動測量系統，可進行大面積高密度空間三維數據的采集，具有測量點位精度高、采集空間點的密度大、速度快等特點[3].地面三維激光掃描技術的應用為解決土方量計算的空間三維數據采集，提高土方量計算精度提供了有效的方法，土方工程是土地整理的主要項目之一, 其主要目的是通過填挖土方使平整地塊的平面位置和高程滿足設計要求[4].

目前土方量計算的常用方法有DTM法，方格網法，等高線法，平均高程法和斷面法.DTM模型計算的土方量精度相對較高，它是根據實地測定的地面點坐標和設計高程計算得到土方量，測量點的數量及其精度決定了土方量計算的精度，但是數據點越多，占用內存越大，不適合大范圍土方的計算[5].方格網法是一種基于柵格的解決方法, 其特點是數據結構和數學模型簡單, 易于程序實現, 但格網間距的大小直接影響著計算精度, 而對于復雜地形很難確定合適的格網間距.因此, 計算結果往往與現實情況出入較大, 不能滿足實際應用的需要[6].平均高程法操作比較簡便，但是該方法誤差較大.等高線法的工作內容與步驟和方格法大致相同，不同之處在于計算場地平均高程的方法，它適合地面起伏較大、坡度變化較多時采用.斷面法只需要知道兩端橫斷面的面積，因而計算很簡單．但如果施工場地范圍很大，那么求面積工作量將會非常大，而且容易出錯，同時大量的數據容易在計算時帶來較大的誤差，特別是地形較復雜時，算出的面積值誤差會很大[7].

以上介紹的這四種土方量計算的方法，都有其缺點，為了進一步提高土方量解算精度，簡化土方量計算的工作程序，滿足土方量計算的實際應用要求，現提出有限元的思想.

1 有限元體積及表面積計算原理

任何立方結構都可以通過數學解剖模型分解為有限個體積元.對于土方工程量計算來說，可以先對該土方填挖處進行三維激光地面掃描，利用IDL軟件構建三維立體模型.三維實體重構的首要任務是將測量數據按實物原型的幾何特征進行分割成不同的數據塊, 使得位于同一數據塊內的數據點可以用一張特定的曲面來表示, 然后針對不同數據塊采用不同的曲面建構方案進行曲面重建, 最后將這些曲面塊拼接成實體[8].將重建后的曲面實體進行數學有限元積分，求出其精確的表面面積[9-10].再將此立體模型根據實際地形分布劃分為有限個體積元，分析每個體積元中的土方分布，從而計算其土方量.

1.1 有限元體積計算

1.1.1 有限元體積計算原理

將復雜實體三維空間劃分為有限個體積大小相同的立方元(圖1).利用積分原理，先擬合出所求立方元外部曲面并根據曲面上已知點云數據計算該曲面法線向量，再利用柱面積分原理計算體積.在一定范圍內，分割的立方元體積越小，得出的體積精度越高.

圖1 三維立體結構有限元分割圖

用有限元思想解算掃描實體模型土石方量的步驟為：

(1)根據點云數據計算空間直角坐標系X軸方向的最大值Xmax和最小值Xmin，Y軸方向的最大值Ymax和最小值Ymin.

(2)如圖2所示，將模擬大點云分割為形狀、大小相等的有限個虛擬的小立方元，對單個立方元建立空間直角坐標系o-xyz，x′軸、y′軸、z′軸分別與x軸、y軸、z軸平行.

圖2 立方元分割示意圖

圖3 立方元法向量示意圖

(3)對該立方元內曲面下的區域進行雙重積分求體積.

1.1.2 曲面擬合

曲線曲面擬合是一種古老而常用的技術，在工程、實驗、統計和計算機圖形等方面有著廣泛的應用[11].建立最佳的擬合參數, 運用最小二乘法對其進行曲面擬合, 通過線性變換精確解出控制點的坐標值[12].擬合的基本思路為先選擇基本反映曲面形狀而且分布均勻的數據點,利用曲面蒙皮技術[13]生成一個基本反映型面形狀而且性態良好的曲面,稱之為初始曲面.然后以某種映射關系找到數據點在初始曲面上的對應點,并把該點的參數值稱為數據點在初始曲面上的參數值,同時建立超定線性方程組.最后利用最小二乘法求解超定線性方程組,得到數據點的最佳擬合曲面[14].本文中我們給出三種基本曲面擬合原理公式，其中法向矢量均以z軸方向為例.

十參數曲面擬合原理[10]的二次曲面方程為

Ax2+By2+Cz2+D+

Exy+Fxz+Gx+Hyz+Jy+Kz=0

(1)

三維空間R3法向矢量為

(2)

其中：

根據法向量與三個坐標軸之間的夾角，選擇式(3)或式(4)中的簡化方程.雙曲面S上的點表示為S={Pk}?R3,k=1,2,…，n,k點滿足方程：

Sx?R3:xk-A-Byk-Czk-Dykzk=0

(3a)

Sy?R3:yk-A-Bxk-Czk-Dxkzk=0

(3b)

Sz?R3:zk-A-Bxk-Cyk-Dxkzk=0

(3c)

地面平坦時用(3c)，陡峭地方用(3a)或(3b)擬合.當為二次拋物面擬合時，用如下方程：

x-A-By-Cz-Dy2-Ez2-Fyz=0

(4a)

y-A-Bx-Cz-Dx2-Ez2-Fxz=0

(4b)

z-A-Bx-Cy-Dx2-Ey2-Fxy=0

(4c)

上式使用原理同式(3)所述，首先判斷法向量與各軸之間的夾角.如當法向量與x軸之間的夾角較小時，在下一步積分時就應該使用關于x軸的方向矢量.

1.1.3 幾種體積計算方法

以沿z方向上積分為例，給出幾種體積積分公式

(1)十參數擬合曲面體積積分公式：

Fxz+Gx+Hyz+Jy+Kz)dxdy=

(Cz2+D+Kz)(x2-x1)(y2-y1)+

(5)

(2)四參數擬合曲面體積積分公式：

A(x2-x1)(y2-y1)+

(6)

(3)六參數擬合曲面體積積分公式：

Dx2+Ey2+Fxy)dxdy=

A(x2-x1)(y2-y1)+

(7)

其中Ω代表實體體積，公式中x2、x1分別為空間直角坐標系x軸方向的最大值和最小值，y2、y1分別為空間直角坐標系y軸方向的最大值和最小值.

1.2 有限元表面積計算

曲面積分可以看做是無限個子區域的平面積分，即為有限元數學積分原理.

有限元表面積計算與Bézier曲面面積計算原理相似.令Bézier曲面面積為S,分割后得到多個小的曲面片,設第i個小曲面片的面積為Si.由曲面分割過程知,原Bézier曲面經過1次分割后得到4個曲面片,經過2次分割得到16個曲面片,一直進行下去,經過k次分割后可得到4k個曲面片[9]，如圖4所示.

圖4 曲面有限元分割

有限元表面面積計算公式：

(8)

SΔij為k次分割后得到的第i個曲面片的控制多面體中第j個三角形的面積.

2 實驗及分析

2.1 實驗模型及數據介紹

本實驗數據為不同的分辨率(相鄰兩個掃描點之間的距離)點云數據組成的半球，半徑為10m.表1是不同的分辨率下模擬半球對應的點個數.

表1 不同的分辨率下模擬半球對應的點個數

模擬的半球如圖5所示，半球點個數約130萬、分辨率為0.038m,數據處理采用IDL語言開發的基于激光點云的EEXLT(地圖要素提取)軟件.該實驗默認坐標系為空間直角坐標系，各軸之間的位置要求符合右手法則.

本實驗中設定的有限元單位立方元長0.2m、寬0.2m、高0.2m.

圖5 半球模型

2.2 結果分析

以分辨率為0.038m的半球模型為例，說明本方法的可行性.表2為分辨率0.038m的半球模型的部分點云數據，其高程起算點為零點，建立坐標系原點為(0,0,0)，既以X軸、Y軸所在的平面為起算平面.將模擬半球大點云按長0.2m、寬0.2m、高0.2m的有限元單位立方元分割，得到約2萬個單位立方元，表3為部分小立方元內的點個數及體積.

表3 分辨率為0.038m的模擬半球部分小立方元內的點個數及體積

表4 有限元思想解算不同分辨率下半徑為10m的模擬半球的體積、表面積及相對誤差

相對誤差計算公式如下：

(9)

2.3 有限元體積計算精度評定

有限元法體積和表面積計算的精度用相對誤差來評定，見圖6.

由圖6可以看出，點云掃描分辨率越高，體積和面積的計算結果越精確.圖7可以看出，模型點云總數越多，體積相對誤差越小，結果越精確.表5為分別用等高線法、DTM法、格網法、斷面法及有限元法計算分辨率為0.038m的模型體積的結果與精度對比.

圖6 不同分辨率下體積和表面積相對誤差

圖7 不同分辨率下點云總數的模型體積相對誤差

計算方法所求體積/m3體積相對誤差/%有限元法2094.3511/240000等高線法1471.9151/3DTM法2105.31/40格網法2075.01/100斷面法2106.9511/40

顯然，相對于其它已知方法，有限元法大大提高了實體模型體積的求解精度，為工程預算提供了可靠的數據.

3 結論與展望

(1)本文提出的有限元法計算不規則實體的表面積和體積是通過有限元積分將其劃分為足夠小的區域，根據無限逼近的方法進行解求實體的表面積和體積.

(2)有限元法所分得不規則實體的體積元與點云獲取的分辨率大小有密切的關系，分辨率越大，其體積元就越小，計算的體積就越精確，反之亦然.

(3)有限元法計算的體積比等高線法、DTM法、格網法、斷面法的精度要高幾個數量級.

(4)影響有限元積分法計算體積精度的因素為：擬合面上點的個數及分割立方元的大小.所以應用有限元積分原理計算體積時，應首先利用三維激光掃描儀正確并充分掃描得到實體表面無噪聲點云數據，并根據實體實際情況劃分合適的有限元.

(5)三維激光掃描得到的點云數據必須能夠覆蓋整個實體表面，點云數據在使用之前必須去除噪音點.在進行曲面擬合，獲取曲面參數的時候必須保證有足夠多的點云落在該擬合曲面上.

(6)在有限元分割時依據實際情況將實體分割為盡可能小的立方元，以便于進一步提高體積計算的精度.

(7)有限元積分法用于計算工程土石方量，不同實體參數下的曲面擬合，有限元分割有較大差異，具體的分類方法還需進一步研究.

[1] Casula G,Mora P,Bianchi M G. Detection of terrain morphologic features using GPS, TLS, and land surveys:“ttana della volpe” blind valley case study[J]. Journal of Surveying Engineering， 2009,136(3):132-138.

[2] Pirotti F,Guarnieri A,Vettore A.Ground filtering and vegetation mapping using multi-return terrestrial laser scanning[J].ISPRS Journal of Photogrammetry and Remote Sensing,2013,76:56-63.

[3] 朱清海,黃承亮,李凱.基于EPS 2008及地面三維激光掃描點云數據進行斷面線提取[J].城市勘測,2013(2):89-91.

[4] 王曉莉，王琳琳.工程土方量計算中土方量精度與DEM分辨率的關系[J].科技情報開發與經濟，2007,17(7):191-192.

[5] 胡奎.三維激光掃描在土方計算中的應用[J].礦山測量，2013,(1):71-72.

[6] 楊友長，景妮.構TIN法填挖方計算方法研究[J].金屬礦山.2008(7):88-91.

[7] 劉鴻劍，肖偉紅.基于DEM的工程土方量算法研究[J].江西理工大學報，2009,30(4):18-21.

[8] 胡敏捷,閻利,張毅.三維激光掃描技術在大型艙容測量中的應用研究[J].船舶設計通訊訊,2011(1):60-64.

[9] 王勝兵,戴明強,黃登斌.基于雙線性插值擬合的山形曲面面積計算[J].兵工自動化, 2012, 3(3): 42-43.

[10] 曹俊,陳普春,徐瑩,等.基于斷層掃描思想的曲面面積計算方法研究[J].現代電子技術, 2012, 35(6): 131-133.

[11] 何芳.移動曲面擬合法在復雜曲面造型中的研究與應用[D].武漢:武漢理工大學, 2008.

[12] 李二濤,張國煊,曾虹.基于最小二乘的曲面擬合算法研究 [J][J].杭州電子科技大學學報, 2009, 29(2): 48-51.

[13] 李玉梅.基于三次 B樣條的曲線、曲面逼近算法的研究[D].南京:南京信息工程大學, 2013.

[14] 吳小剛.幾何造型技術中逆向柔性曲面重構技術研究[D].成都：電子科技大學, 2012.