劉尚蔚, 楊陽, 趙繼偉
(華北水利水電大學(xué),河南 鄭州 450046)
隨著傾斜攝影技術(shù)的不斷發(fā)展,無人機(jī)攝影測量技術(shù)在水利工程中得到了廣泛應(yīng)用[1]。無人機(jī)攝影測量具有低成本、高效率的特點(diǎn),彌補(bǔ)了傳統(tǒng)近景攝影測量布控困難的問題,在逆向建模過程中具有獨(dú)特的優(yōu)勢。面對水利工程體積龐大、基巖表面不平整的特點(diǎn),無人機(jī)傾斜攝影很好地解決了效率的問題,但由于傾斜角度的限制,造成了部分地面視角數(shù)據(jù)的丟失。對此,文獻(xiàn)[2]研究了LiDAR點(diǎn)云輔助的航空影像空中三角測量,提高了航空影像數(shù)據(jù)的定位精度;文獻(xiàn)[3]在傾斜攝影的基礎(chǔ)上加入了雙手機(jī)成像攝影系統(tǒng),實(shí)現(xiàn)了空地一體的三維重建;文獻(xiàn)[4]利用無人機(jī)和激光掃描儀對建筑物進(jìn)行了聯(lián)合三維建模。這些傾斜攝影與近景攝影結(jié)合的逆向建模三維重建方法,克服了單一數(shù)據(jù)源的局限,進(jìn)一步完善了逆向建模的數(shù)據(jù)來源。
然而,無人機(jī)在采集圖像的過程中,易受到外界環(huán)境因素和圖像采集設(shè)備本身的影響,尤其是水利工程通常所處測量環(huán)境差,導(dǎo)致圖像幾何信息略微偏差,對三維重建產(chǎn)生了較大的影響。文獻(xiàn)[5]基于經(jīng)驗(yàn)小波變換算法思想對無人機(jī)采集的圖像進(jìn)行了直接降噪,實(shí)現(xiàn)了較好的圖像降噪效果,但該方法僅停留在平面,未就三維重建后的模型進(jìn)行處理分析。文獻(xiàn)[6]根據(jù)點(diǎn)云濾波非地面數(shù)據(jù)獲取建筑物高度等信息,結(jié)合無人機(jī)采樣紋理進(jìn)行三維重建,實(shí)現(xiàn)參數(shù)化建模,但這種方法嚴(yán)重依賴人工干預(yù),工作量巨大,無法滿足水利工程快速建模的需求。文獻(xiàn)[7]采用二次濾波的方法減輕了三維匹配塊在平滑區(qū)產(chǎn)生的抓痕現(xiàn)象,但該方法無法處理邊緣噪點(diǎn)產(chǎn)生的鋸齒狀紋理。文獻(xiàn)[8-10]對邊緣檢測算子進(jìn)行了討論,分別提出了邊緣保持的方法,但面向信息量較大的水利工程模型時(shí)效率很低,不適合推廣。
針對以上問題,本文以引導(dǎo)濾波算法[11]為基本思想,從近景攝影和無人機(jī)傾斜攝影三維重建的模型出發(fā),對該模型進(jìn)行簡化等預(yù)處理,提出了一種適用于水利工程依據(jù)邊緣信息進(jìn)行引導(dǎo)濾波的自適應(yīng)算法,以提高逆向建模的精度。
對于水利工程,本文采用傾斜攝影測量和近景攝影測量結(jié)合的方法獲取原始模型信息,對傾斜攝影測量進(jìn)行運(yùn)動(dòng)結(jié)構(gòu)恢復(fù),并輔以地面近景攝影測量進(jìn)行矯正,實(shí)現(xiàn)對實(shí)體水利工程的逆向建模。
本文利用無人機(jī)對水利工程進(jìn)行傾斜攝影測量,獲取測區(qū)目標(biāo)物幾何信息,同時(shí)記錄無人機(jī)定位信息和相機(jī)焦距等屬性;使用廣義成像儀實(shí)現(xiàn)地面上的近景攝影測量,在地面視角實(shí)現(xiàn)無人機(jī)無法獲取的凹陷部信息獲取,完善無人機(jī)傾斜攝影的不足。將獲取的傾斜攝影和近景攝影數(shù)據(jù)導(dǎo)入ContextCapture軟件,通過空中三角測量計(jì)算獲取空中影像外方位元素,結(jié)合地面近景攝影信息構(gòu)件不規(guī)則三角網(wǎng)(Triangulated Irregular Network,TIN)模型,最后在空間中進(jìn)行紋理映射生成三維重建模型,以obj格式文件導(dǎo)出。上述三維重建的技術(shù)流程如圖1所示,以某拱壩樞紐為測區(qū)所建立的逆向模型如圖2所示。
圖1 傾斜攝影協(xié)同近景攝影三維重建技術(shù)路線
圖2 傾斜攝影協(xié)同近景攝影逆向建模生成的某拱壩樞紐模型
由采集誤差導(dǎo)致部分未經(jīng)上述過程成功配準(zhǔn)的像素點(diǎn)在模型中生成了噪點(diǎn),形成了三維重建后的誤差。無人機(jī)采集信息過程中的誤差來源于采集設(shè)備和外界環(huán)境的變化,該誤差服從均值為0的正態(tài)分布[12]。由此,為進(jìn)一步提高三維重建的精度,本文采用濾波的原理對上述過程生成的三維重建模型進(jìn)行進(jìn)一步處理。
曲率下采樣,即在點(diǎn)云曲率越大的地方采樣點(diǎn)越多。取某點(diǎn)為中心的鄰域進(jìn)行研究,計(jì)算點(diǎn)與鄰域法線的夾角值,夾角值越大說明曲率越大。設(shè)定一角度閾值,若某點(diǎn)鄰域夾角值大于該閾值,則認(rèn)為該點(diǎn)處為特征明顯的區(qū)域;反之,則為特征不明顯區(qū)域。對特征明顯區(qū)域和特征不明顯區(qū)域進(jìn)行獨(dú)立的均勻采樣,采樣數(shù)分別為S(1-U)與SU,其中S為目標(biāo)采樣數(shù),U為采樣均勻性系數(shù)。
點(diǎn)T鄰域窗口大小的控制決定了局部紋理的性質(zhì),其中包含的采樣點(diǎn)個(gè)數(shù)t為重要參數(shù)。對此,本文在點(diǎn)K處設(shè)一半徑為R的球面,球面以內(nèi)的點(diǎn)為所求鄰域點(diǎn),該球半徑R為自適應(yīng)取值。
這種采樣由于在水工逆向模型中特征明顯區(qū)域的采樣密度大,故計(jì)算效率高。通過幾何特征區(qū)域的劃分,使得采樣結(jié)果具有強(qiáng)的抗噪性、高的穩(wěn)定性。
2.2 離群統(tǒng)計(jì)消除
假設(shè)模型服從高斯分布,其結(jié)構(gòu)由均值和標(biāo)準(zhǔn)差決定,平均距離在指定范圍(視為全局平均距離和方差的函數(shù))之外的點(diǎn)可視為離群點(diǎn)剔除。為提高檢索速度,可應(yīng)用KD樹對點(diǎn)云進(jìn)行索引[13]。
某點(diǎn)到所有K鄰域點(diǎn)的平均距離Da按式(1)計(jì)算:
(1)
式中:K為鄰域中所有點(diǎn)的數(shù)量;i為鄰域中點(diǎn)的索引值;Di為該點(diǎn)到鄰域中第i個(gè)點(diǎn)的空間距離。
統(tǒng)計(jì)分析所有點(diǎn)平均距離Da的平均值D以及方差S2,分別按下式計(jì)算:
(2)
(3)
式中:n為模型中點(diǎn)的數(shù)量;j為模型中點(diǎn)的索引值。
離群點(diǎn)判定閾值DT定義為:
DT=D+αS2。
(4)
式中α為松弛系數(shù),按照試驗(yàn)與訓(xùn)練的需求進(jìn)行設(shè)定。
若某點(diǎn)的Da>DT,則剔除該點(diǎn);反之,則保留該點(diǎn)。
2.3 引導(dǎo)濾波原理
引導(dǎo)濾波是一種局部線性濾波器,對圖像平滑濾波的同時(shí)還保持了邊緣細(xì)節(jié),通常用于二維平面圖像,稍加修改可用于圖2中三維重建的水工模型。首先定義一個(gè)通用的線性平移變量濾波的過程,包括一個(gè)引導(dǎo)圖像I、一個(gè)濾波輸入圖像p和一個(gè)輸出圖像q。引導(dǎo)濾波的關(guān)鍵是引導(dǎo)圖像I與輸出圖像q之間的局部線性模型,假設(shè)q在I中以點(diǎn)k為中心的區(qū)間wk進(jìn)行線性變換:
qi=Akpi+Bk,?i∈ωk。
(5)
式中:i、k為點(diǎn)的索引值;ωk為中心位于k時(shí)的窗口,半徑為Ra,Ra為輸入?yún)?shù);qi為濾波后輸出的三維點(diǎn);pi為輸入三維點(diǎn)(假設(shè)與引導(dǎo)圖像存在對應(yīng)關(guān)系,替代引導(dǎo)圖像對應(yīng)點(diǎn));Ak為一個(gè)3×3的方陣;Bk為一個(gè)3×1的矩陣。
從濾波器的約束中輸入p,在模型輸出q的過程中,除去輸入p中不需要的成分n(如噪聲和多余的紋理):
qi=pi-ni。
(6)
上述基本原理如圖3所示。
圖3 引導(dǎo)濾波原理
引導(dǎo)濾波模型中回歸模型是假設(shè)的解析表達(dá)式,此時(shí)圖像為二維函數(shù),可用系數(shù)ak和bk分別替換式(5)中的Ak和Bk。對于該回歸模型,要求有合理的系數(shù)ak和bk使輸入值p與輸出值q之間的差距最小,轉(zhuǎn)為最優(yōu)化問題,即:
(7)
式中:e為向量(1,1,1);ε為ak的規(guī)整因子,為調(diào)節(jié)引導(dǎo)濾波效果的重要參數(shù),用于防止所求的ak過大。
(8)
(9)
(10)
2.4 Marr-Hildreth算子改進(jìn)的引導(dǎo)濾波
原生的引導(dǎo)濾波,在不同的窗口采用了相同的規(guī)整因子ε,顯然是未考慮到不同窗口中點(diǎn)云的紋理存在的差異。根據(jù)2.3節(jié)的結(jié)論,對于紋理多變、邊緣梯度信息豐富的窗口,線性模型系數(shù)ak的值偏大,需要較小的規(guī)整因子ε來修正;而對于邊緣梯度平緩的窗口,需要調(diào)整更小的規(guī)整因子ε來減小逼近誤差。因此,有必要根據(jù)不同的窗口對規(guī)整因子ε進(jìn)行自適應(yīng)調(diào)節(jié)。LI Z等[14]對引導(dǎo)濾波進(jìn)行了改進(jìn),該方法將以窗口內(nèi)的方差定義加權(quán)因子,對規(guī)整因子ε進(jìn)行自適應(yīng)調(diào)整。實(shí)際上,方差定義的加權(quán)因子并不能很好地反映圖形的邊緣信息。
本文將加權(quán)引導(dǎo)濾波的思路引入至水工模型的修正中,同時(shí)考慮到水利工程模型信息量大的特點(diǎn),利用Marr-Hildreth算子重新定義加權(quán)因子。這種算子將高斯濾波和拉普拉斯算子結(jié)合[15],可以給出封閉的邊緣邊界,并且能夠避開遞歸滯后的計(jì)算過程,從而減少了堆棧要求。定義邊緣權(quán)重Marr-Hildreth算子如下:
(11)
式中:Mi為當(dāng)前窗口的Marr-Hildreth邊緣檢測算子;Mi′為引導(dǎo)圖像的Marr-Hildreth邊緣檢測算子;N為引導(dǎo)圖像中像素點(diǎn)數(shù)量;i′為對應(yīng)索引值;i為當(dāng)前窗口中點(diǎn)的索引值;η為固定值。
由于Marr-Hildreth算子為二階Laplacian算子的一種,在邊緣處具有較大的幅值,因此在邊緣處點(diǎn)的ri>1,在平滑處點(diǎn)的ri<1。經(jīng)過試驗(yàn)與訓(xùn)練發(fā)現(xiàn),自適應(yīng)于Mi值變化的η有利于提高算法的魯棒性,對此,本文中η取為|Mi|最大值的0.1倍,效果較好。
將式(11)稱為Marr-Hildreth加權(quán)引導(dǎo)濾波算子,替換2.3節(jié)中的ε為ε/ri,式(7)替換為本文求解引導(dǎo)濾波線性參數(shù)方案,即:
(12)
對于初步建立的三維模型,為減少后續(xù)工作的計(jì)算量,并盡可能保留局部線性變化效果,對該模型進(jìn)行曲率下采樣工作。由于離散的點(diǎn)云之間不存在拓?fù)浣Y(jié)構(gòu),為了得到其幾何屬性,通過各點(diǎn)的鄰域結(jié)構(gòu)來識(shí)別點(diǎn)云的幾何特征,進(jìn)行低尺度離群點(diǎn)的剔除,得到引導(dǎo)模型?;谝龑?dǎo)濾波原理利用引導(dǎo)模型對原模型進(jìn)行平滑處理,引入邊緣檢測算子實(shí)現(xiàn)該原理對模型細(xì)節(jié)紋理的自適應(yīng)調(diào)整。具體降噪方法步驟如下:
1)對于水工逆向建模obj格式的三維模型,使用pcl_mesh_sampling程序進(jìn)行整體點(diǎn)云采樣,獲得模型原始點(diǎn)云數(shù)據(jù)。由于水利工程模型曲率變化不均勻不連續(xù),這種整體采樣的均勻采樣方法會(huì)在曲率較小處獲取同樣密度的點(diǎn),造成后續(xù)工作需求計(jì)算量變大。為減少這種曲率不同產(chǎn)生的計(jì)算量浪費(fèi),本文對點(diǎn)云原始模型依據(jù)曲率大小按照2.1節(jié)中的原理進(jìn)行二次采樣,以降低后續(xù)工作計(jì)算量。
2)為消除設(shè)備收集數(shù)據(jù)時(shí)產(chǎn)生的部分噪聲和偽影造成過大的點(diǎn)云偏差,引入低尺度統(tǒng)計(jì)離群消除濾波算法,按照2.2節(jié)中的原理,依據(jù)各點(diǎn)到其鄰域內(nèi)所有點(diǎn)的平均距離,將超出閾值距離過大的離群點(diǎn)剔除。按照該方法對曲率下采樣后的模型進(jìn)行離群點(diǎn)的剔除,減小水工逆向模型中偏差過大的離群噪點(diǎn),獲得引導(dǎo)濾波中所需要的引導(dǎo)模型。
3)利用統(tǒng)計(jì)離群點(diǎn)處理后的引導(dǎo)模型對采樣后的待濾波模型進(jìn)行引導(dǎo),依據(jù)2.3節(jié)中原理,按照式(1)的方法進(jìn)行引導(dǎo)濾波。在求解引導(dǎo)濾波的過程中,依據(jù)2.4節(jié)中改進(jìn)原理,以式(12)為依據(jù)求線性變換參數(shù)ak與bk,通過局部線性變換與約束條件完成對待濾波水工模型的平滑處理。
上述方法技術(shù)路線如圖4所示。
圖4 本文處理方法技術(shù)路線
Python提供了高效的數(shù)據(jù)結(jié)構(gòu),可以簡單地面向?qū)ο缶幊?能夠在多數(shù)平臺(tái)上運(yùn)行,同時(shí)方便擴(kuò)展新功能,本文調(diào)用了Numpy和Open3D擴(kuò)展程序庫。Numpy代表了“Numeric Python”,是由多維數(shù)組對象及用于處理數(shù)組的例程集合組成的程序庫,相對于Python list不需要做類型檢查,可減少代碼的冗余,提高程序運(yùn)行速度。Open3D是用于處理三維數(shù)據(jù)的程序庫,支持快速開發(fā)和處理3D數(shù)據(jù),可對點(diǎn)云進(jìn)行操作及可視化。本文基于引導(dǎo)濾波的水利工程逆向模型降噪方法的部分步驟中關(guān)鍵Python程序如下:
1)將ply格式的點(diǎn)云文件與程序放入同一目錄中,具體程序代碼如下:
import numpy as np
import open3d as o3d
pcd = o3d.io.read_point_cloud(′ex.ply′)
2)根據(jù)曲率判定點(diǎn)的特征是否明顯,分別對特征明顯點(diǎn)和特征不明顯點(diǎn)進(jìn)行均勻采樣,拼接點(diǎn)云,主要程序代碼如圖5所示。
圖5 曲率下采樣主要程序代碼
3)通過KD樹索引計(jì)算點(diǎn)云密度,主要程序代碼如圖6所示。
圖6 計(jì)算點(diǎn)云密度主要程序代碼
4)通過定義的閾值剔除離群的點(diǎn)云。主要代碼如下:
c1,ind = pcd.remove_radius_outlier(K,D)
radius_cloud = pcd.select_by_index(ind)
5)定義引導(dǎo)濾波函數(shù),在每個(gè)窗口內(nèi)求解線性變換參數(shù),對點(diǎn)進(jìn)行線性變換,循環(huán)處理每個(gè)窗口,并輸出點(diǎn)云,部分程序代碼如圖7所示。
圖7 部分引導(dǎo)濾波函數(shù)代碼
為驗(yàn)證本文方法的有效性,利用Python語言對上文的方法編程實(shí)現(xiàn)。試驗(yàn)使用的計(jì)算機(jī)硬件配置為AMD Ryzen 7 4800H的CPU,32G內(nèi)存,使用NVIDIA GeForce RTX 2060的GPU進(jìn)行加速,操作系統(tǒng)Win10 64位,運(yùn)行環(huán)境Python 3.7。將obj模型和pcl_mesh_sampling.exe文件放入同一文件夾下,在cmd中打開該文件夾目錄,運(yùn)行pcl_mesh_sampling.exe程序,輸入obj模型并輸出pcd格式點(diǎn)云文件,在MATLAB中利用pcwrite命令將pcd格式轉(zhuǎn)換為ply格式的點(diǎn)云文件,如圖8(a)所示。按照3.1節(jié)中的方法,將py文件與ply點(diǎn)云文件放入同一目錄下運(yùn)行Python程序:
1)按照曲率大小對點(diǎn)云采樣,效果如圖8(b)所示,模型坐標(biāo)點(diǎn)由36 004個(gè)降至28 533個(gè);
2)對采樣后的點(diǎn)云模型統(tǒng)計(jì)離群點(diǎn)并剔除,效果如圖8(c)所示;
3)利用剔除離群點(diǎn)后的模型對曲率下采樣后的模型做引導(dǎo)濾波,效果如圖8(d)所示。
圖8 點(diǎn)云文件降噪過程示意圖
本文使用無人機(jī)與廣義成像儀對某拱壩樞紐進(jìn)行采樣,采用ContextCapture軟件逆向建立該工程的obj格式三維模型(如圖2),按照上述方法及操作對該逆向模型進(jìn)行處理,試驗(yàn)中涉及的參數(shù)見表1。
表1 試驗(yàn)參數(shù)設(shè)定
將處理后的模型與原始的模型紋理部分邊緣細(xì)節(jié)進(jìn)行對比,觀察模型噪點(diǎn)變化以及紋理效果,并按《三維地理信息模型數(shù)據(jù)產(chǎn)品規(guī)范》[16]中的規(guī)定核實(shí)精度要求,與本文方法處理前的模型進(jìn)行對比,驗(yàn)證本文方法的有效性。在該水利工程中選取30個(gè)檢查點(diǎn)進(jìn)行實(shí)際坐標(biāo)測量,作為精度檢測的依據(jù)。
水利工程三維重建模型中的地形和水工建筑物表面幾何結(jié)構(gòu)、曝光度和紋理清晰程度需要達(dá)到較高水平。為此,將本文所提出的降噪方法處理前后的三維重建模型進(jìn)行對比分析??紤]到實(shí)驗(yàn)中水利樞紐工程結(jié)構(gòu)較大,為保證對比效果,本文選取4個(gè)代表部位結(jié)構(gòu)進(jìn)行對比,結(jié)果見表2。
表2 本文方法對某拱壩逆向模型處理前后紋理對比
由表2中試驗(yàn)結(jié)果對比圖可以看出:由于運(yùn)動(dòng)結(jié)構(gòu)恢復(fù)過程中形成了一些離群點(diǎn)和噪點(diǎn),導(dǎo)致原始建立的模型邊緣出現(xiàn)一些鋸齒狀紋理,甚至部分曲率較小的曲面上出現(xiàn)了不光滑平整的現(xiàn)象。通過本文所述方法處理后,大尺度的離群噪點(diǎn)被統(tǒng)計(jì)離群點(diǎn)的方法剔除,小尺度的噪點(diǎn)在本文的自適應(yīng)引導(dǎo)濾波的算法下得以修正,曲面處理得更加平順,有效改善了邊緣鋸齒的現(xiàn)象,在主觀上達(dá)成了更優(yōu)秀的視覺體驗(yàn)。
在客觀評價(jià)上,引入結(jié)構(gòu)相似度指數(shù)(Structure Similarity Index Measure,SSIM)與均方誤差(Mean Squared Error,MSE)進(jìn)行定量評價(jià)[17]。SSIM反映了兩個(gè)模型的相似程度,取值在0到1之間,數(shù)值越大表示模型三維重建越完整,質(zhì)量越好。MSE反映了處理前后的模型變化大小,該值越小說明方法效果愈顯著。本文將處理前后的三維模型進(jìn)行比對,得到的SSIM參數(shù)為0.86,即處理后的模型較為完整地保留了原始模型的信息;得到的MSE參數(shù)為0.15,說明降噪效果比較明顯。
根據(jù)《三維地理信息模型數(shù)據(jù)產(chǎn)品規(guī)范》中三維模型精度規(guī)定的Ⅰ級(jí)要求,平面精度中誤差(σXY)不高于0.3 m,按式(13)計(jì)算;高度精度中誤差(σZ)不高于0.5 m,按式(14)計(jì)算。
(13)
(14)
式中:X實(shí)、Y實(shí)、Z實(shí)為點(diǎn)的實(shí)際坐標(biāo),由實(shí)測獲??;X模、Y模、Z模為對應(yīng)的點(diǎn)在模型中的坐標(biāo);n為檢查點(diǎn)的數(shù)量,本試驗(yàn)中為30。
由式(13)和式(14)計(jì)算可得,未經(jīng)處理的三維模型在水平精度上的中誤差為0.037 m,高度精度中誤差為0.031 m;經(jīng)本文方法處理后的三維模型在水平精度上的中誤差為0.032 m,高度精度中誤差為0.028 m,水平精度中誤差降低了13.5%,高度精度中誤差降低了9.7%,合計(jì)精度中誤差降低了11.9%。上述結(jié)果均滿足規(guī)范中的要求,但顯然本文方法處理后的模型更加精準(zhǔn),處理前后兩模型各檢查點(diǎn)三維誤差向量統(tǒng)計(jì)如圖9所示,并以平均大小為半徑做球面作為參考,經(jīng)本文方法處理后的誤差均值由0.046 m降低至0.041 m,降低了11.3%。各檢查點(diǎn)累計(jì)誤差歸一化后的方向向量分別為(0.224,-0.588,-0.187)和(0.217,-0.604,-0.189),兩向量夾角弧度為0.018 5,基本保持一致,即經(jīng)本文處理后的模型既降低了原有模型的誤差,也保留了原有誤差的方向?qū)傩浴?/p>
圖9 本文方法處理模型前后誤差統(tǒng)計(jì)
本文在多源采集數(shù)據(jù)建立的模型基礎(chǔ)上,采用曲率下采樣減小不必要的工作量,并以低尺度的統(tǒng)計(jì)離群點(diǎn)剔除獲取引導(dǎo)圖像,引入Marr-Hildreth算子采用引導(dǎo)濾波對模型進(jìn)行處理。試驗(yàn)結(jié)果表明,本文提出的方法在滿足國家相關(guān)規(guī)范的基礎(chǔ)上進(jìn)一步提高了三維重建的精度(平均誤差降低11.3%),模型細(xì)節(jié)更平順,改善了邊緣鋸齒情況,基本保留了采集數(shù)據(jù)的完整性,同時(shí)滿足了水利工程快速逆向建模的需求。更高精度的三維重建,可為水利工程施工期信息提供載體,為施工可視化仿真提供更多可能。另外,受當(dāng)前計(jì)算機(jī)計(jì)算能力的限制,難以快速做到設(shè)計(jì)模型對攝影采集模型的自動(dòng)引導(dǎo)匹配。未來,隨著計(jì)算機(jī)計(jì)算能力的提高,匹配分塊相關(guān)算法還需進(jìn)一步研究,水利工程的逆向建模精度還有提升空間。