劉尚蔚, 楊陽, 趙繼偉

(華北水利水電大學,河南 鄭州 450046)

隨著傾斜攝影技術的不斷發(fā)展,無人機攝影測量技術在水利工程中得到了廣泛應用[1]。無人機攝影測量具有低成本、高效率的特點,彌補了傳統(tǒng)近景攝影測量布控困難的問題,在逆向建模過程中具有獨特的優(yōu)勢。面對水利工程體積龐大、基巖表面不平整的特點,無人機傾斜攝影很好地解決了效率的問題,但由于傾斜角度的限制,造成了部分地面視角數(shù)據(jù)的丟失。對此,文獻[2]研究了LiDAR點云輔助的航空影像空中三角測量,提高了航空影像數(shù)據(jù)的定位精度；文獻[3]在傾斜攝影的基礎上加入了雙手機成像攝影系統(tǒng),實現(xiàn)了空地一體的三維重建；文獻[4]利用無人機和激光掃描儀對建筑物進行了聯(lián)合三維建模。這些傾斜攝影與近景攝影結合的逆向建模三維重建方法,克服了單一數(shù)據(jù)源的局限,進一步完善了逆向建模的數(shù)據(jù)來源。

然而,無人機在采集圖像的過程中,易受到外界環(huán)境因素和圖像采集設備本身的影響,尤其是水利工程通常所處測量環(huán)境差,導致圖像幾何信息略微偏差,對三維重建產(chǎn)生了較大的影響。文獻[5]基于經(jīng)驗小波變換算法思想對無人機采集的圖像進行了直接降噪,實現(xiàn)了較好的圖像降噪效果,但該方法僅停留在平面,未就三維重建后的模型進行處理分析。文獻[6]根據(jù)點云濾波非地面數(shù)據(jù)獲取建筑物高度等信息,結合無人機采樣紋理進行三維重建,實現(xiàn)參數(shù)化建模,但這種方法嚴重依賴人工干預,工作量巨大,無法滿足水利工程快速建模的需求。文獻[7]采用二次濾波的方法減輕了三維匹配塊在平滑區(qū)產(chǎn)生的抓痕現(xiàn)象,但該方法無法處理邊緣噪點產(chǎn)生的鋸齒狀紋理。文獻[8-10]對邊緣檢測算子進行了討論,分別提出了邊緣保持的方法,但面向信息量較大的水利工程模型時效率很低,不適合推廣。

針對以上問題,本文以引導濾波算法[11]為基本思想,從近景攝影和無人機傾斜攝影三維重建的模型出發(fā),對該模型進行簡化等預處理,提出了一種適用于水利工程依據(jù)邊緣信息進行引導濾波的自適應算法,以提高逆向建模的精度。

1 水利工程逆向建模技術

對于水利工程,本文采用傾斜攝影測量和近景攝影測量結合的方法獲取原始模型信息,對傾斜攝影測量進行運動結構恢復,并輔以地面近景攝影測量進行矯正,實現(xiàn)對實體水利工程的逆向建模。

本文利用無人機對水利工程進行傾斜攝影測量,獲取測區(qū)目標物幾何信息,同時記錄無人機定位信息和相機焦距等屬性；使用廣義成像儀實現(xiàn)地面上的近景攝影測量,在地面視角實現(xiàn)無人機無法獲取的凹陷部信息獲取,完善無人機傾斜攝影的不足。將獲取的傾斜攝影和近景攝影數(shù)據(jù)導入ContextCapture軟件,通過空中三角測量計算獲取空中影像外方位元素,結合地面近景攝影信息構件不規(guī)則三角網(wǎng)(Triangulated Irregular Network,TIN)模型,最后在空間中進行紋理映射生成三維重建模型,以obj格式文件導出。上述三維重建的技術流程如圖1所示,以某拱壩樞紐為測區(qū)所建立的逆向模型如圖2所示。

圖1 傾斜攝影協(xié)同近景攝影三維重建技術路線

圖2 傾斜攝影協(xié)同近景攝影逆向建模生成的某拱壩樞紐模型

由采集誤差導致部分未經(jīng)上述過程成功配準的像素點在模型中生成了噪點,形成了三維重建后的誤差。無人機采集信息過程中的誤差來源于采集設備和外界環(huán)境的變化,該誤差服從均值為0的正態(tài)分布[12]。由此,為進一步提高三維重建的精度,本文采用濾波的原理對上述過程生成的三維重建模型進行進一步處理。

2 引導濾波及相關原理

2.1 曲率下采樣

曲率下采樣,即在點云曲率越大的地方采樣點越多。取某點為中心的鄰域進行研究,計算點與鄰域法線的夾角值,夾角值越大說明曲率越大。設定一角度閾值,若某點鄰域夾角值大于該閾值,則認為該點處為特征明顯的區(qū)域；反之,則為特征不明顯區(qū)域。對特征明顯區(qū)域和特征不明顯區(qū)域進行獨立的均勻采樣,采樣數(shù)分別為S(1-U)與SU,其中S為目標采樣數(shù),U為采樣均勻性系數(shù)。

點T鄰域窗口大小的控制決定了局部紋理的性質,其中包含的采樣點個數(shù)t為重要參數(shù)。對此,本文在點K處設一半徑為R的球面,球面以內的點為所求鄰域點,該球半徑R為自適應取值。

這種采樣由于在水工逆向模型中特征明顯區(qū)域的采樣密度大,故計算效率高。通過幾何特征區(qū)域的劃分,使得采樣結果具有強的抗噪性、高的穩(wěn)定性。

2.2 離群統(tǒng)計消除

假設模型服從高斯分布,其結構由均值和標準差決定,平均距離在指定范圍(視為全局平均距離和方差的函數(shù))之外的點可視為離群點剔除。為提高檢索速度,可應用KD樹對點云進行索引[13]。

某點到所有K鄰域點的平均距離Da按式(1)計算：

(1)

式中：K為鄰域中所有點的數(shù)量；i為鄰域中點的索引值；Di為該點到鄰域中第i個點的空間距離。

統(tǒng)計分析所有點平均距離Da的平均值D以及方差S2,分別按下式計算：

(2)

(3)

式中：n為模型中點的數(shù)量；j為模型中點的索引值。

離群點判定閾值DT定義為：

DT=D+αS2。

(4)

式中α為松弛系數(shù),按照試驗與訓練的需求進行設定。

若某點的Da>DT,則剔除該點；反之,則保留該點。

2.3 引導濾波原理

引導濾波是一種局部線性濾波器,對圖像平滑濾波的同時還保持了邊緣細節(jié),通常用于二維平面圖像,稍加修改可用于圖2中三維重建的水工模型。首先定義一個通用的線性平移變量濾波的過程,包括一個引導圖像I、一個濾波輸入圖像p和一個輸出圖像q。引導濾波的關鍵是引導圖像I與輸出圖像q之間的局部線性模型,假設q在I中以點k為中心的區(qū)間wk進行線性變換：

qi=Akpi+Bk,?i∈ωk。

(5)

式中：i、k為點的索引值；ωk為中心位于k時的窗口,半徑為Ra,Ra為輸入?yún)?shù)；qi為濾波后輸出的三維點；pi為輸入三維點(假設與引導圖像存在對應關系，替代引導圖像對應點)；Ak為一個3×3的方陣；Bk為一個3×1的矩陣。

從濾波器的約束中輸入p,在模型輸出q的過程中,除去輸入p中不需要的成分n(如噪聲和多余的紋理)：

qi=pi-ni。

(6)

上述基本原理如圖3所示。

圖3 引導濾波原理

引導濾波模型中回歸模型是假設的解析表達式,此時圖像為二維函數(shù),可用系數(shù)ak和bk分別替換式(5)中的Ak和Bk。對于該回歸模型,要求有合理的系數(shù)ak和bk使輸入值p與輸出值q之間的差距最小,轉為最優(yōu)化問題,即：

(7)

式中：e為向量(1,1,1)；ε為ak的規(guī)整因子,為調節(jié)引導濾波效果的重要參數(shù),用于防止所求的ak過大。

(8)

(9)

(10)

2.4 Marr-Hildreth算子改進的引導濾波

原生的引導濾波,在不同的窗口采用了相同的規(guī)整因子ε,顯然是未考慮到不同窗口中點云的紋理存在的差異。根據(jù)2.3節(jié)的結論,對于紋理多變、邊緣梯度信息豐富的窗口,線性模型系數(shù)ak的值偏大,需要較小的規(guī)整因子ε來修正；而對于邊緣梯度平緩的窗口,需要調整更小的規(guī)整因子ε來減小逼近誤差。因此,有必要根據(jù)不同的窗口對規(guī)整因子ε進行自適應調節(jié)。LI Z等[14]對引導濾波進行了改進,該方法將以窗口內的方差定義加權因子,對規(guī)整因子ε進行自適應調整。實際上,方差定義的加權因子并不能很好地反映圖形的邊緣信息。

本文將加權引導濾波的思路引入至水工模型的修正中,同時考慮到水利工程模型信息量大的特點,利用Marr-Hildreth算子重新定義加權因子。這種算子將高斯濾波和拉普拉斯算子結合[15],可以給出封閉的邊緣邊界,并且能夠避開遞歸滯后的計算過程,從而減少了堆棧要求。定義邊緣權重Marr-Hildreth算子如下：

(11)

式中：Mi為當前窗口的Marr-Hildreth邊緣檢測算子;Mi′為引導圖像的Marr-Hildreth邊緣檢測算子;N為引導圖像中像素點數(shù)量;i′為對應索引值;i為當前窗口中點的索引值;η為固定值。

由于Marr-Hildreth算子為二階Laplacian算子的一種,在邊緣處具有較大的幅值,因此在邊緣處點的ri>1,在平滑處點的ri<1。經(jīng)過試驗與訓練發(fā)現(xiàn),自適應于Mi值變化的η有利于提高算法的魯棒性,對此,本文中η取為|Mi|最大值的0.1倍,效果較好。

將式(11)稱為Marr-Hildreth加權引導濾波算子,替換2.3節(jié)中的ε為ε/ri,式(7)替換為本文求解引導濾波線性參數(shù)方案,即：

(12)

3 基于引導濾波的水工逆向模型降噪方法及實現(xiàn)

3.1 基于引導濾波的逆向建模降噪方法

對于初步建立的三維模型,為減少后續(xù)工作的計算量,并盡可能保留局部線性變化效果,對該模型進行曲率下采樣工作。由于離散的點云之間不存在拓撲結構,為了得到其幾何屬性,通過各點的鄰域結構來識別點云的幾何特征,進行低尺度離群點的剔除,得到引導模型?；谝龑V波原理利用引導模型對原模型進行平滑處理,引入邊緣檢測算子實現(xiàn)該原理對模型細節(jié)紋理的自適應調整。具體降噪方法步驟如下：

1)對于水工逆向建模obj格式的三維模型,使用pcl_mesh_sampling程序進行整體點云采樣,獲得模型原始點云數(shù)據(jù)。由于水利工程模型曲率變化不均勻不連續(xù),這種整體采樣的均勻采樣方法會在曲率較小處獲取同樣密度的點,造成后續(xù)工作需求計算量變大。為減少這種曲率不同產(chǎn)生的計算量浪費,本文對點云原始模型依據(jù)曲率大小按照2.1節(jié)中的原理進行二次采樣,以降低后續(xù)工作計算量。

2)為消除設備收集數(shù)據(jù)時產(chǎn)生的部分噪聲和偽影造成過大的點云偏差,引入低尺度統(tǒng)計離群消除濾波算法,按照2.2節(jié)中的原理,依據(jù)各點到其鄰域內所有點的平均距離,將超出閾值距離過大的離群點剔除。按照該方法對曲率下采樣后的模型進行離群點的剔除,減小水工逆向模型中偏差過大的離群噪點,獲得引導濾波中所需要的引導模型。

3)利用統(tǒng)計離群點處理后的引導模型對采樣后的待濾波模型進行引導,依據(jù)2.3節(jié)中原理,按照式(1)的方法進行引導濾波。在求解引導濾波的過程中,依據(jù)2.4節(jié)中改進原理,以式(12)為依據(jù)求線性變換參數(shù)ak與bk,通過局部線性變換與約束條件完成對待濾波水工模型的平滑處理。

上述方法技術路線如圖4所示。

圖4 本文處理方法技術路線

3.2 運用Python語言的程序設計

Python提供了高效的數(shù)據(jù)結構,可以簡單地面向對象編程,能夠在多數(shù)平臺上運行,同時方便擴展新功能,本文調用了Numpy和Open3D擴展程序庫。Numpy代表了“Numeric Python”,是由多維數(shù)組對象及用于處理數(shù)組的例程集合組成的程序庫,相對于Python list不需要做類型檢查,可減少代碼的冗余,提高程序運行速度。Open3D是用于處理三維數(shù)據(jù)的程序庫,支持快速開發(fā)和處理3D數(shù)據(jù),可對點云進行操作及可視化。本文基于引導濾波的水利工程逆向模型降噪方法的部分步驟中關鍵Python程序如下：

1)將ply格式的點云文件與程序放入同一目錄中，具體程序代碼如下：

import numpy as np

import open3d as o3d

pcd = o3d.io.read_point_cloud(′ex.ply′)

2)根據(jù)曲率判定點的特征是否明顯,分別對特征明顯點和特征不明顯點進行均勻采樣,拼接點云,主要程序代碼如圖5所示。

圖5 曲率下采樣主要程序代碼

3)通過KD樹索引計算點云密度,主要程序代碼如圖6所示。

圖6 計算點云密度主要程序代碼

4)通過定義的閾值剔除離群的點云。主要代碼如下:

c1,ind = pcd.remove_radius_outlier(K,D)

radius_cloud = pcd.select_by_index(ind)

5)定義引導濾波函數(shù),在每個窗口內求解線性變換參數(shù),對點進行線性變換,循環(huán)處理每個窗口,并輸出點云,部分程序代碼如圖7所示。

圖7 部分引導濾波函數(shù)代碼

4 試驗與分析

4.1 試驗方法

為驗證本文方法的有效性,利用Python語言對上文的方法編程實現(xiàn)。試驗使用的計算機硬件配置為AMD Ryzen 7 4800H的CPU,32G內存,使用NVIDIA GeForce RTX 2060的GPU進行加速,操作系統(tǒng)Win10 64位,運行環(huán)境Python 3.7。將obj模型和pcl_mesh_sampling.exe文件放入同一文件夾下,在cmd中打開該文件夾目錄,運行pcl_mesh_sampling.exe程序,輸入obj模型并輸出pcd格式點云文件,在MATLAB中利用pcwrite命令將pcd格式轉換為ply格式的點云文件,如圖8(a)所示。按照3.1節(jié)中的方法,將py文件與ply點云文件放入同一目錄下運行Python程序：

1)按照曲率大小對點云采樣,效果如圖8(b)所示,模型坐標點由36 004個降至28 533個；

2)對采樣后的點云模型統(tǒng)計離群點并剔除,效果如圖8(c)所示；

3)利用剔除離群點后的模型對曲率下采樣后的模型做引導濾波,效果如圖8(d)所示。

圖8 點云文件降噪過程示意圖

本文使用無人機與廣義成像儀對某拱壩樞紐進行采樣,采用ContextCapture軟件逆向建立該工程的obj格式三維模型(如圖2),按照上述方法及操作對該逆向模型進行處理,試驗中涉及的參數(shù)見表1。

表1 試驗參數(shù)設定

將處理后的模型與原始的模型紋理部分邊緣細節(jié)進行對比,觀察模型噪點變化以及紋理效果,并按《三維地理信息模型數(shù)據(jù)產(chǎn)品規(guī)范》[16]中的規(guī)定核實精度要求,與本文方法處理前的模型進行對比,驗證本文方法的有效性。在該水利工程中選取30個檢查點進行實際坐標測量,作為精度檢測的依據(jù)。

4.2 試驗結果

水利工程三維重建模型中的地形和水工建筑物表面幾何結構、曝光度和紋理清晰程度需要達到較高水平。為此,將本文所提出的降噪方法處理前后的三維重建模型進行對比分析?？紤]到實驗中水利樞紐工程結構較大,為保證對比效果,本文選取4個代表部位結構進行對比,結果見表2。

表2 本文方法對某拱壩逆向模型處理前后紋理對比

由表2中試驗結果對比圖可以看出：由于運動結構恢復過程中形成了一些離群點和噪點,導致原始建立的模型邊緣出現(xiàn)一些鋸齒狀紋理,甚至部分曲率較小的曲面上出現(xiàn)了不光滑平整的現(xiàn)象。通過本文所述方法處理后,大尺度的離群噪點被統(tǒng)計離群點的方法剔除,小尺度的噪點在本文的自適應引導濾波的算法下得以修正,曲面處理得更加平順,有效改善了邊緣鋸齒的現(xiàn)象,在主觀上達成了更優(yōu)秀的視覺體驗。

在客觀評價上,引入結構相似度指數(shù)(Structure Similarity Index Measure,SSIM)與均方誤差(Mean Squared Error,MSE)進行定量評價[17]。SSIM反映了兩個模型的相似程度,取值在0到1之間,數(shù)值越大表示模型三維重建越完整,質量越好。MSE反映了處理前后的模型變化大小,該值越小說明方法效果愈顯著。本文將處理前后的三維模型進行比對,得到的SSIM參數(shù)為0.86,即處理后的模型較為完整地保留了原始模型的信息；得到的MSE參數(shù)為0.15,說明降噪效果比較明顯。

根據(jù)《三維地理信息模型數(shù)據(jù)產(chǎn)品規(guī)范》中三維模型精度規(guī)定的Ⅰ級要求,平面精度中誤差(σXY)不高于0.3 m,按式(13)計算；高度精度中誤差(σZ)不高于0.5 m,按式(14)計算。

(13)

(14)

式中：X實、Y實、Z實為點的實際坐標,由實測獲??；X模、Y模、Z模為對應的點在模型中的坐標；n為檢查點的數(shù)量,本試驗中為30。

由式(13)和式(14)計算可得,未經(jīng)處理的三維模型在水平精度上的中誤差為0.037 m,高度精度中誤差為0.031 m；經(jīng)本文方法處理后的三維模型在水平精度上的中誤差為0.032 m,高度精度中誤差為0.028 m,水平精度中誤差降低了13.5%,高度精度中誤差降低了9.7%,合計精度中誤差降低了11.9%。上述結果均滿足規(guī)范中的要求,但顯然本文方法處理后的模型更加精準,處理前后兩模型各檢查點三維誤差向量統(tǒng)計如圖9所示,并以平均大小為半徑做球面作為參考,經(jīng)本文方法處理后的誤差均值由0.046 m降低至0.041 m,降低了11.3%。各檢查點累計誤差歸一化后的方向向量分別為(0.224,-0.588,-0.187)和(0.217,-0.604,-0.189),兩向量夾角弧度為0.018 5,基本保持一致,即經(jīng)本文處理后的模型既降低了原有模型的誤差，也保留了原有誤差的方向屬性。

圖9 本文方法處理模型前后誤差統(tǒng)計

5 結語

本文在多源采集數(shù)據(jù)建立的模型基礎上,采用曲率下采樣減小不必要的工作量,并以低尺度的統(tǒng)計離群點剔除獲取引導圖像,引入Marr-Hildreth算子采用引導濾波對模型進行處理。試驗結果表明,本文提出的方法在滿足國家相關規(guī)范的基礎上進一步提高了三維重建的精度(平均誤差降低11.3%),模型細節(jié)更平順,改善了邊緣鋸齒情況,基本保留了采集數(shù)據(jù)的完整性,同時滿足了水利工程快速逆向建模的需求。更高精度的三維重建,可為水利工程施工期信息提供載體,為施工可視化仿真提供更多可能。另外,受當前計算機計算能力的限制,難以快速做到設計模型對攝影采集模型的自動引導匹配。未來，隨著計算機計算能力的提高,匹配分塊相關算法還需進一步研究,水利工程的逆向建模精度還有提升空間。