蔡來良 康洪躍 杜 莊 張志斌

(1.河南理工大學(xué)測繪與國土信息工程學(xué)院,河南 焦作 454003;2.中國南水北調(diào)集團(tuán)中線有限公司渠首分公司,河南 南陽 473000)

近年來,三維激光掃描技術(shù)逐漸在高效率獲取空間立體信息應(yīng)用中占據(jù)主要地位,成為地理空間數(shù)據(jù)獲取的一種重要的技術(shù)手段[1]。三維激光掃描技術(shù)在我國礦山測量領(lǐng)域得到了大量研究和應(yīng)用,特別是在露天開采邊坡監(jiān)測、露天礦儲(chǔ)量管理、井工開采邊坡監(jiān)測、采空區(qū)安全監(jiān)測、地表沉陷監(jiān)測等方面應(yīng)用較為廣泛[2]。相較于傳統(tǒng)的全站儀、RTK等方式,三維激光掃描技術(shù)具有掃描速度快、數(shù)據(jù)獲取全面等優(yōu)點(diǎn),尤其是在懸崖、邊坡等危險(xiǎn)地點(diǎn),該技術(shù)可以實(shí)現(xiàn)無接觸測量,有助于確保測量人員安全。國內(nèi)許多學(xué)者在基于點(diǎn)云的地表變形分析方法方面取得了豐富的研究成果,大致有點(diǎn)云直接對(duì)比求取變形法、點(diǎn)云模型對(duì)比求取變形法、點(diǎn)云特征點(diǎn)線提取變形法這3類方法。

在點(diǎn)云直接對(duì)比求取變形的研究方面,劉昌軍等[3]通過對(duì)Hausdorff距離算法進(jìn)行改進(jìn),提出了基于八叉樹結(jié)構(gòu)的點(diǎn)云與點(diǎn)云的精細(xì)直接比較算法,實(shí)現(xiàn)了變形區(qū)域多個(gè)點(diǎn)云數(shù)據(jù)的自動(dòng)比較。王海城等[4]提出了一種基于兩期點(diǎn)云直接計(jì)算變形的新方法,即通過對(duì)掃描點(diǎn)云進(jìn)行特征分割、坐標(biāo)變換與投影、格網(wǎng)劃分、最短坐標(biāo)距離計(jì)算與排序、中位值段取平均以及逆變換等一系列處理過程后,最終得到整個(gè)變形體高密度監(jiān)測點(diǎn)的變形值。陳弘奕等[5]針對(duì)點(diǎn)云數(shù)據(jù)的特點(diǎn),提出了點(diǎn)云和點(diǎn)云直接比較變形分析方法,采用最佳擬合面法來提取變化信息。魏振全[6]采用三維激光掃描儀與全站儀相組合的測量方法,對(duì)田陳煤礦某工作面上方地表移動(dòng)和建筑物變形情況進(jìn)行了監(jiān)測和分析,獲取了控制點(diǎn)三維坐標(biāo),進(jìn)而求取了變形量。

在點(diǎn)云模型對(duì)比求取變形算法方面,柏雯娟[7]基于礦區(qū)的兩期三維激光掃描數(shù)據(jù),利用Kriging算法分別建立了礦區(qū)地表下沉盆地?cái)?shù)字高程模型,對(duì)礦區(qū)地表沉陷進(jìn)行分析,最后選取礦區(qū)若干有代表性的監(jiān)測點(diǎn)的開采沉陷監(jiān)測值與對(duì)應(yīng)的水準(zhǔn)測量結(jié)果進(jìn)行了對(duì)比分析。于海洋等[8]在LiDAR DEM精度與不確定性分析的基礎(chǔ)上,利用模糊推理方法建立了基于坡度、點(diǎn)云密度和地表粗糙度的誤差相關(guān)表面,用于差值DEM不確定性的量化與DEM最小變化閾值探測,采用基于權(quán)重濾波窗口的貝葉斯估計(jì)判定與修正,較精確地獲取了研究區(qū)地表形變信息。呂志鵬等[9]獲取了變形體變形前后的點(diǎn)云數(shù)據(jù),并利用移動(dòng)最小二乘法進(jìn)行了均勻格網(wǎng)擬合內(nèi)插來提取變形信息。孟萬利等[10]以相同網(wǎng)格劃分開采前后點(diǎn)云數(shù)據(jù),再利用改進(jìn)的反距離權(quán)重法插值求取網(wǎng)格節(jié)點(diǎn)坐標(biāo),匹配相同網(wǎng)格節(jié)點(diǎn)作為同名點(diǎn)并求取下沉量,最后獲取了下沉曲線與下沉DEM。梁周雁等[11]對(duì)DEM模型求差法和Hausdorff距離法計(jì)算的監(jiān)測結(jié)果進(jìn)行了比較分析,驗(yàn)證了兩種方法在地表變形監(jiān)測中的優(yōu)勢。王堃宇等[12]采用三維激光掃描技術(shù)實(shí)現(xiàn)了對(duì)某邊坡位移的監(jiān)測,獲取了整個(gè)邊坡的三維模型及位移云圖。何榮等[13]采用S-G濾波對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行去噪處理,建立了沉陷區(qū)三維數(shù)字高程模型,利用局部表面擬合方法求取了特征點(diǎn)法向量,計(jì)算了地表變形前后法向量的變化量,并與特征點(diǎn)真實(shí)傾斜角進(jìn)行了對(duì)比分析。楊帆等[14]通過對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行處理,形成了相應(yīng)時(shí)段的尾礦邊坡三維模型,采用等高線比較、三維質(zhì)檢軟件比較等方法,提取了邊坡變形數(shù)據(jù)。劉衛(wèi)南等[15]將離散的三維點(diǎn)云轉(zhuǎn)化為二維的密度圖像,再利用粒子圖像測速技術(shù)分析了位移前后兩幅點(diǎn)云密度圖像的相關(guān)性,從而計(jì)算了柵格圖像中各子集的相對(duì)位移值,得到了目標(biāo)區(qū)域的平面位移場。盧遙等[16]通過對(duì)礦區(qū)兩期點(diǎn)云數(shù)據(jù)提取的DEM作差運(yùn)算,得出該段時(shí)期內(nèi)地表的垂直位移量,而后通過設(shè)置高差閾值和面積閾值,并采用基于高差分析方法,提取了地表沉陷信息。何倩等[17]利用地面三維激光掃描點(diǎn)云數(shù)據(jù)構(gòu)建高分辨率、高精度DEM,再由低分辨率SRTM DEM內(nèi)插補(bǔ)充數(shù)據(jù)缺失區(qū)域以獲取整個(gè)研究區(qū)域的DEM,最后將其與高分辨率SAR影像進(jìn)行時(shí)序差分處理,獲取礦區(qū)地表動(dòng)態(tài)沉降監(jiān)測數(shù)據(jù)。

在基于點(diǎn)云特征點(diǎn)線提取變形的研究方面,張小青[18]通過提取變形監(jiān)測點(diǎn)的三維坐標(biāo)信息,進(jìn)行多期數(shù)據(jù)的監(jiān)測點(diǎn)坐標(biāo)信息比較,獲取了局部的變形信息。徐進(jìn)軍等[19]根據(jù)確定的邊界,對(duì)點(diǎn)云進(jìn)行網(wǎng)格分割并通過基于最短距離取“中位值”的變形計(jì)算方法,大致確定了同名變形部位,最后利用該部位的兩次掃描點(diǎn)云來提取同名變形點(diǎn),并計(jì)算其變形。司夢元等[20]設(shè)計(jì)了一套基于邊坡的永久性三棱錐作為控制點(diǎn),通過各期點(diǎn)云數(shù)據(jù)中三棱錐特征提取相應(yīng)控制點(diǎn),運(yùn)用算法程序處理從而提取了邊坡變形信息。趙不釩等[21]針對(duì)地面三維激光掃描技術(shù)在變形監(jiān)測中的應(yīng)用需求,結(jié)合點(diǎn)云序列特征,提出了利用點(diǎn)云特征要素提取變形信息的方法。李強(qiáng)等[22]選取了研究區(qū)變形特征明顯的墻體作為特征點(diǎn)提取對(duì)象,通過提取兩期墻體特征點(diǎn)坐標(biāo)值,獲得了研究區(qū)的典型特征點(diǎn)變形結(jié)果,并提取了研究區(qū)Z軸方向剖面線,可直觀地反映出地面沉降線性變形特征。呂寶雄等[23]利用地面三維激光掃描技術(shù)提取了崩塌滑坡災(zāi)害體特征點(diǎn)三維坐標(biāo)、虛擬斷面線,通過點(diǎn)、線、面的綜合分析,可獲得災(zāi)害體的局部細(xì)節(jié)和整體位移量信息,實(shí)現(xiàn)對(duì)災(zāi)害體的變形監(jiān)測。

上述分析表明,近年來利用三維激光掃描技術(shù)監(jiān)測地表移動(dòng)的研究取得了顯著進(jìn)展,但還存在一定的不足。例如:在利用三維激光掃描儀開展變形監(jiān)測時(shí),在缺乏人工標(biāo)靶的情況下,無法準(zhǔn)確采集到同名點(diǎn)坐標(biāo)。本研究對(duì)此進(jìn)行了相關(guān)探索,提出了同名特征區(qū)域概念,建立了特征區(qū)域最近點(diǎn)迭代地表移動(dòng)分析算法(Surface Movement Analysis Algorithm Based on Feature Region Iteration Closest Points,FRICP),并通過試驗(yàn)分析驗(yàn)證了該算法的可行性。

1 特征區(qū)域最近點(diǎn)迭代地表移動(dòng)分析算法(FRICP)原理

為保證觀測人員安全并提高工作效率,采用三維激光掃描儀測量山地陡坡的地表地形,并通過對(duì)比多期觀測數(shù)據(jù),分析地表移動(dòng)。由于每期掃描激光發(fā)射方向不同,無法在各期觀測值中找到相同觀測點(diǎn),進(jìn)而也無法使用常規(guī)的定點(diǎn)觀測數(shù)據(jù)處理方法計(jì)算地表移動(dòng)。為了克服多期觀測無同名點(diǎn)這一技術(shù)難點(diǎn),本研究將山區(qū)地表分布的巖石、樹木等地物視為特征區(qū)域,構(gòu)建同名特征區(qū)域提取算法;并在此基礎(chǔ)上建立多期同名特征區(qū)域的匹配算法,最后利用ICP算法分析同名特征區(qū)域的地表移動(dòng)量,算法流程如圖1所示。

圖1 FRICP算法流程Fig.1 Flow of FRICP algorithm

1.1 點(diǎn)云數(shù)據(jù)預(yù)處理

使用地面三維激光掃描儀獲得的數(shù)據(jù)無法直接進(jìn)行變形分析,需進(jìn)行點(diǎn)云的預(yù)處理,主要步驟如下:

(1)點(diǎn)云去噪。將點(diǎn)云中的噪聲點(diǎn)去除,降低噪聲點(diǎn)對(duì)后續(xù)分析影響。

(2)點(diǎn)云配準(zhǔn)。對(duì)采集的第一期和第二期的多站點(diǎn)云數(shù)據(jù)進(jìn)行配準(zhǔn),建立完整的測區(qū)點(diǎn)云。

(3)兩期數(shù)據(jù)坐標(biāo)轉(zhuǎn)換。以第一期配準(zhǔn)后的點(diǎn)云數(shù)據(jù)坐標(biāo)系統(tǒng)為基準(zhǔn)坐標(biāo)系統(tǒng),將第二期配準(zhǔn)結(jié)束的點(diǎn)云數(shù)據(jù)坐標(biāo)轉(zhuǎn)換到基準(zhǔn)坐標(biāo)系統(tǒng)。

(4)點(diǎn)云濾波。將第一期點(diǎn)云數(shù)據(jù)和第二期的點(diǎn)云數(shù)據(jù)進(jìn)行濾波,分離地面點(diǎn)和非地面點(diǎn),為提取同名地面特征區(qū)域做準(zhǔn)備。

1.2 地面特征區(qū)域點(diǎn)云提取

山區(qū)地表起伏較大,與平原地區(qū)相比,地表特征十分明顯。在點(diǎn)云中提取特征區(qū)域可為形變分析奠定基礎(chǔ)。

1.2.1 使用地面點(diǎn)云提取地面特征區(qū)域

地表點(diǎn)云法向量是描述地表凹凸起伏的關(guān)鍵幾何參數(shù)。平坦地區(qū)的點(diǎn)云法向量變化較小,非平坦區(qū)域(特征區(qū)域)的法向量變化較大,如圖2所示。可以利用該點(diǎn)的法向量與周圍點(diǎn)法向量的夾角和的平均值來判斷該點(diǎn)是否為特征點(diǎn),如果該點(diǎn)法向量與鄰域內(nèi)其他點(diǎn)的法向量夾角和的平均值較大(設(shè)置角度閾值為20°),認(rèn)為該點(diǎn)是特征點(diǎn);反之,該點(diǎn)不是特征點(diǎn)。

圖2 不同地形的點(diǎn)云法向量Fig.2 Point cloud normal vectors of different terrains

經(jīng)過法向量關(guān)系提取出的特征點(diǎn)離散性強(qiáng),沒有建立聚集關(guān)系,無法形成特征區(qū)域。本研究采用一種基于密度對(duì)噪聲魯棒的空間聚類算法(Density-based Spatial Clustering of Applications with Noise,DBSCAN)來進(jìn)行點(diǎn)云聚類。該算法流程如下:

(1)讀入點(diǎn)云數(shù)據(jù),設(shè)定初始參數(shù)鄰域半徑為0.1,最小鄰域點(diǎn)個(gè)數(shù)為6。

(2)遍歷所有點(diǎn)云數(shù)據(jù)對(duì)點(diǎn)云進(jìn)行鄰域計(jì)算,若某一點(diǎn)鄰域的個(gè)數(shù)大于最小鄰域點(diǎn)個(gè)數(shù),則認(rèn)為該點(diǎn)為核心點(diǎn)。

(3)如果點(diǎn)q在核心點(diǎn)p的鄰域內(nèi),則稱點(diǎn)q和點(diǎn)p是直接密度可達(dá)的。如果存在數(shù)據(jù)點(diǎn)集p1,p2,…,pn是從點(diǎn)p1出發(fā)直接密度可達(dá)的,則稱點(diǎn)pn和點(diǎn)p1是密度可達(dá)的。任意選取一個(gè)核心點(diǎn)找出由其密度可達(dá)的所有點(diǎn)云數(shù)據(jù),生成點(diǎn)云簇。

(4)迭代執(zhí)行上述步驟,直到所有的核心點(diǎn)都計(jì)算完為止。

DBSCAN算法流程如圖3所示。

圖3 DBSCAN聚類流程Fig.3 Flow of DBSCAN clustering

聚類效果如圖4所示。地面特征點(diǎn)聚類完成之后就形成了帶有聚類屬性的特征區(qū)域,不再是離散的特征點(diǎn)。

圖4 DBSCAN算法聚類效果示意Fig.4 Schematic of the clustering effects by DBSCAN algorithm

1.2.2 使用植被點(diǎn)云提取地面特征區(qū)域

山地測區(qū)除了地形特征明顯之外,植被信息也十分豐富。本研究將樹木視為天然地面標(biāo)靶,對(duì)其鄰近范圍地面進(jìn)行標(biāo)定。詳細(xì)步驟為:首先采用濾波算法獲得植被點(diǎn)云;其次設(shè)置高程閾值分離低矮植被;然后對(duì)高大樹木進(jìn)行單木分割,獲得高大樹木根部附近的點(diǎn)云;最后對(duì)點(diǎn)云進(jìn)行聚類,樹木聚類效果如圖5所示。以樹根點(diǎn)云上一點(diǎn)為球心,r為半徑畫球,球內(nèi)所包含的地面點(diǎn)云即為樹木所標(biāo)識(shí)的地面特征區(qū)域。

圖5 植被點(diǎn)云聚類示意Fig.5 Schematic of vegetation point cloud clustering

1.3 同名特征區(qū)域匹配

獲得各期監(jiān)測的特征區(qū)域之后,對(duì)多期特征區(qū)域進(jìn)行同名匹配,找到不同觀測時(shí)期數(shù)據(jù)中的相同觀測目標(biāo)。下文分別探討地面特征區(qū)域和樹木特征區(qū)域的同名匹配算法。

1.3.1 地面點(diǎn)云提取的同名特征區(qū)域匹配

首先求出第一期地面特征區(qū)域的重心,然后以第一期地面特征區(qū)域的重心為球心,以1.2～1.5倍最大變形預(yù)計(jì)量為半徑,搜索第二期地面特征區(qū)域重心。搜索到第二期區(qū)域內(nèi)重心點(diǎn)后,對(duì)所有出現(xiàn)在該搜索區(qū)域內(nèi)的地面特征區(qū)域求其特征度(點(diǎn)云數(shù)量、包圍盒大小等),分別與第一期地面特征區(qū)域的特征度進(jìn)行對(duì)比分析。若存在符合閾值的第二期地面特征區(qū)域,則認(rèn)為該特征區(qū)域?yàn)橥麉^(qū)域,并將該區(qū)域點(diǎn)云記錄下來。

1.3.2 使用樹木點(diǎn)云進(jìn)行同名地面特征區(qū)域匹配

樹冠隨季節(jié)變化而變化,樹干在橫向上生長緩慢,樹干在水平方向上比較穩(wěn)定且接近圓形,可以作為匹配的依據(jù)。使用距地面0.5倍樹高的水平截面截取樹木,對(duì)該截面內(nèi)的樹干點(diǎn)云利用最小二乘法擬合圓,利用擬合的圓心和半徑匹配同名植被。首先以第一期植被點(diǎn)云擬合得到的圓心為搜索圓心,以r為搜索半徑,建立圓型搜索區(qū)域;再將第二期的植被點(diǎn)云數(shù)據(jù)擬合求得的圓心放到第一期數(shù)據(jù)中,利用圓形搜索區(qū)域進(jìn)行同名點(diǎn)匹配;最后將第一期和第二期的半徑作差,并設(shè)定閾值,如果符合閾值,則認(rèn)為是同名植被。

以擬合的截面圓心為球心,以大于截面高度的長度為半徑畫球,球內(nèi)包含就是植被點(diǎn)云所標(biāo)注的地面特征區(qū)域。對(duì)兩期同名植被所標(biāo)注的地面特征區(qū)域進(jìn)行特征度比較,例如點(diǎn)云數(shù)量、包圍盒大小等,如果符合閾值,則保留。

1.4 利用ICP算法進(jìn)行變形分析

1.4.1 ICP算法原理

ICP算法是在源點(diǎn)云P中選取點(diǎn)pi集,在目標(biāo)點(diǎn)云Q中找到對(duì)應(yīng)的鄰近點(diǎn)qi集,找到并由此解算旋轉(zhuǎn)矩陣和平移矩陣,使得誤差函數(shù)值最小,并將旋轉(zhuǎn)矩陣和平移矩陣用于源點(diǎn)云中的點(diǎn)集,得到新的對(duì)應(yīng)點(diǎn)集,計(jì)算新的對(duì)應(yīng)點(diǎn)集和目標(biāo)點(diǎn)云中點(diǎn)集的平均距離,直到平均距離小于某一給定的閾值,或者迭代到一定的次數(shù),即可完成配準(zhǔn)。ICP算法配準(zhǔn)流程如圖6所示。

圖6所示的旋轉(zhuǎn)矩陣R中,α、β、γ分別為沿X,Y,Z軸的旋轉(zhuǎn)角度。

圖6 ICP配準(zhǔn)流程Fig.6 Flow of ICP registration

平移矩陣T為

式中,tx,ty,tz分別表示沿Z,Y,Z軸的平移量。

誤差函數(shù)E(R,T)為

式中,n為最鄰近點(diǎn)集個(gè)數(shù);pi為源點(diǎn)云點(diǎn)集pi中一點(diǎn);qi為目標(biāo)點(diǎn)云點(diǎn)集qi中與pi對(duì)應(yīng)的最近點(diǎn);R為旋轉(zhuǎn)矩陣;T為平移矩陣。

式中,pi′為源點(diǎn)云點(diǎn)集中經(jīng)過旋轉(zhuǎn)和平移矩陣計(jì)算后的新源點(diǎn)云點(diǎn)集;pi′為新源點(diǎn)云點(diǎn)集中的一點(diǎn);d為新源點(diǎn)云點(diǎn)集pi′與對(duì)應(yīng)點(diǎn)集qi的平均距離。

1.4.2 變形量計(jì)算

ICP算法共有5個(gè)參數(shù),分別為源點(diǎn)云、目標(biāo)點(diǎn)云、距離閾值、粗配準(zhǔn)的初始變換矩陣(旋轉(zhuǎn)矩陣和平移矩陣)以及迭代次數(shù)。如圖7所示,源點(diǎn)云設(shè)置為第一期點(diǎn)云特征區(qū)域C1,目標(biāo)點(diǎn)云為第二期對(duì)應(yīng)的同名點(diǎn)云特征區(qū)域C2,距離閾值設(shè)置為1.5。由于兩期點(diǎn)云已經(jīng)統(tǒng)一了坐標(biāo)系,不存在旋轉(zhuǎn)變換,所以初始變換矩陣中的旋轉(zhuǎn)矩陣為單位矩陣,平移矩陣為兩期點(diǎn)云特征區(qū)域的重心值相減組成的3×1列矩陣,迭代次數(shù)默認(rèn)為30次。運(yùn)用ICP算法進(jìn)行配準(zhǔn)時(shí),將第一期的點(diǎn)云特征區(qū)域C1配準(zhǔn)到第二期的同名點(diǎn)云特征區(qū)域C2上,用配準(zhǔn)后的第一期數(shù)據(jù)的地面特征區(qū)域C1′的重心(x′,y′,z′)減去C1區(qū)域的重心(x,y,z),該差值(Δx,Δy,Δz)即為C1區(qū)域的位移值。特征區(qū)域位移前后的重心坐標(biāo)差值既可以作為ICP算法的初始平移參數(shù),又可作為下沉和水平位移量的計(jì)算依據(jù),是本研究算法中重要的幾何參數(shù)。

圖7 特征區(qū)域ICP配準(zhǔn)前后對(duì)比Fig.7 Comparison of the feature area before and after ICP registration

求出地表多個(gè)特征區(qū)域的位移值后,可繪制地表移動(dòng)等值線圖,求取地表形變場。本研究算法基于Python語言,引入Open3D庫編程實(shí)現(xiàn),完成點(diǎn)云數(shù)據(jù)處理系統(tǒng)開發(fā)。

2 實(shí)例分析

2.1 試驗(yàn)數(shù)據(jù)采集

本研究以沁水煤田2301工作面采空區(qū)上方的邊坡為例進(jìn)行分析,邊坡長約29 m,寬約28 m,最大高差約16 m,邊坡較為陡峭,坡度約35°。邊坡表面植被茂盛,大多數(shù)為雜草和樹木,地面起伏高低不平,石頭等雜物比較多,邊坡中下部坡度變化較大,且有類似臺(tái)階狀的陡坎,上部坡度變化平緩。地下煤炭資源開采引起該邊坡發(fā)生沉陷和水平移動(dòng)。采用RIEGL VZ-1000型地面激光掃描儀,對(duì)區(qū)域進(jìn)行了變形前后兩期掃描,如圖8所示。

圖8 兩期觀測結(jié)果Fig.8 Observing results of two periods

2.2 地面點(diǎn)云特征區(qū)域提取及匹配

首先,對(duì)兩期點(diǎn)云進(jìn)行去噪、配準(zhǔn)以及坐標(biāo)轉(zhuǎn)換和濾波,為后續(xù)變形分析做準(zhǔn)備。然后,基于法向量對(duì)地面點(diǎn)云進(jìn)行特征點(diǎn)提取,以某一點(diǎn)為圓心,以半徑0.8 m以內(nèi)的點(diǎn)作為該點(diǎn)的鄰域,進(jìn)行平面擬合,求取該點(diǎn)的法向量與鄰域內(nèi)點(diǎn)法向量的夾角和的平均值,并設(shè)置夾角和閾值為20°。第一期地面點(diǎn)云特征點(diǎn)提取出25 130個(gè),第二期地面點(diǎn)云特征點(diǎn)提取出19 889個(gè)?；贒BSCAN算法對(duì)地面點(diǎn)云進(jìn)行了聚類分析,結(jié)果如圖9所示。

圖9 地面特征區(qū)域聚類結(jié)果Fig.9 Region clustering results of ground feature region

以第一期地面特征區(qū)域的重心為球心,建立球型搜索區(qū)域搜索第二期的地面特征區(qū)域,其中搜索半徑設(shè)為2 m,點(diǎn)云數(shù)量差設(shè)為30,包圍盒大小為0.2 m?；诘孛纥c(diǎn)云共搜索到8對(duì)同名特征區(qū)域。

2.3 植被點(diǎn)云特征區(qū)域提取及匹配

將植被點(diǎn)云進(jìn)行低矮植被濾除之后,保留的植被點(diǎn)云基本都是高大的樹木點(diǎn)云。在此基礎(chǔ)上,基于八鄰域的搜索方法進(jìn)行了植被點(diǎn)云聚類,結(jié)果如圖10所示。

圖10 植被點(diǎn)云聚類結(jié)果Fig.10 Clustering results for vegetation point clouds

以0.8 m高度、0.03 m厚度的圓柱去截取樹干,將得到的植被點(diǎn)云進(jìn)行最小二乘擬合,求取半徑和圓心。以同樣的方法建立圓形搜索區(qū)域,匹配第二期的擬合植被圓心,搜索半徑為1,半徑閾值為0.3,共匹配到12對(duì)同名植被點(diǎn)云。

將擬合的圓心置于地面數(shù)據(jù)之中,以擬合圓心為球心、以2 m為半徑建立球形區(qū)域。同名植被點(diǎn)云搜索到的地面點(diǎn)云中,符合點(diǎn)云數(shù)量差值為50、點(diǎn)云包圍盒為0.2的即為同名特征區(qū)域,共找到8對(duì)同名地面特征區(qū)域。

2.4 變形分析

對(duì)從點(diǎn)云中提取的同名特征區(qū)域采用ICP算法進(jìn)行變形分析。用配準(zhǔn)后第一期同名地面特征區(qū)域C1′的重心坐標(biāo)減去配準(zhǔn)前的第一期同名地面特征區(qū)域C1的重心坐標(biāo)得到地表移動(dòng)值。如表1所示,使用點(diǎn)云的特征區(qū)域求取了山地邊坡局部位移,包括下沉和水平移動(dòng)信息,在此基礎(chǔ)上繪制了水平移動(dòng)矢量圖,如圖11所示。

圖11 邊坡水平位移矢量圖Fig.11 Vector diagram of horizontal displacement of slope

表1 特征區(qū)域重心點(diǎn)坐標(biāo)配準(zhǔn)前后變化值Table 1 Variation values of barycentric point coordinates before and after registration in feature region

由表1和圖11可知:水平位移最大值0.46 m位于dm1點(diǎn),最小值0.075 m位于zb6點(diǎn),水平位移矢量方向指向西南方向。西南方向?yàn)楸O(jiān)測對(duì)象的下坡方向,也是采空區(qū)所在方向,求取結(jié)果符合礦區(qū)開采沉陷規(guī)律。

利用表1中16個(gè)點(diǎn)的下沉值生成了地表下沉等值線圖,如圖12所示。由表1和圖12可知:地表最大下沉值-1 m位于dm8點(diǎn),地表最小下沉值-0.65 m位于zb6點(diǎn);西側(cè)地表下沉量多,東側(cè)地表下沉量少,該結(jié)果與西側(cè)距離采空區(qū)較近、東側(cè)距離采空區(qū)較遠(yuǎn)有密切關(guān)系,監(jiān)測到的下沉等值線符合開采沉陷規(guī)律。

圖12 地表下沉等值線圖(單位:m)Fig.12 Contour of surface subsidence

為了分析FRICP算法計(jì)算變形的精度,采用人工精確判讀方法進(jìn)行對(duì)比,以16個(gè)重心點(diǎn)中的8個(gè)地面特征區(qū)域的重心點(diǎn)為參考,用Cloudcompare軟件在兩期原始點(diǎn)云中手工選出對(duì)應(yīng)點(diǎn)云同名特征區(qū)域,使用該軟件附帶的ICP算法進(jìn)行第一期和第二期同名點(diǎn)云特征區(qū)域配準(zhǔn),用配準(zhǔn)后的第一期點(diǎn)云區(qū)域的重心減去原始的第一期點(diǎn)云特征區(qū)域重心,求出重心的差值。將其與FRICP算法計(jì)算的重心變形值進(jìn)行了對(duì)比,結(jié)果見表2。

表2 FRICP算法與人工精確判讀法精度對(duì)比Table 2 Comparison between the precise of FRICP algorithm and manual accurate interpretation method

由表2可知:FRICP算法自動(dòng)計(jì)算和人工精確判讀的地表位移值最大偏差為12 mm,最小偏差為3 mm,三維坐標(biāo)位移中誤差為7.882 mm,X方向中誤差為14.966 mm,Y方向中誤差為10.862 mm,Z方向中誤差為7.818 mm。由此可見:FRICP算法對(duì)于地表移動(dòng)值的計(jì)算精度接近固定測點(diǎn)的監(jiān)測水平,對(duì)于分析煤礦地表變形有一定的實(shí)用價(jià)值。

3 結(jié) 論

(1)針對(duì)變形監(jiān)測工作中三維激光掃描技術(shù)存在的多期測量無法直接獲得同名點(diǎn)的不足,提出了特征區(qū)域概念,并建立了特征區(qū)域最近點(diǎn)迭代地表移動(dòng)分析算法(FRICP)。該算法通過特征區(qū)域自動(dòng)提取、同名匹配、ICP轉(zhuǎn)換等步驟,實(shí)現(xiàn)了地表移動(dòng)量計(jì)算。試驗(yàn)表明:該方法對(duì)于地表移動(dòng)的監(jiān)測精度能達(dá)到定點(diǎn)測量的技術(shù)水平。

(2)FRICP算法擴(kuò)展了傳統(tǒng)ICP算法的應(yīng)用場景,突破了三維激光掃描技術(shù)應(yīng)用于變形監(jiān)測領(lǐng)域的技術(shù)瓶頸,具有一定的理論意義,在山區(qū)開采引起的地表移動(dòng)監(jiān)測等方面有較好的應(yīng)用前景。

(3)本研究算法的較好應(yīng)用場景是被監(jiān)測范圍內(nèi)有明顯的特征區(qū)域。在地表平坦、起伏不大且標(biāo)志物較少的工況下,不明顯特征區(qū)域的提取算法有待進(jìn)一步研究。