面點(diǎn)，確定所有數(shù)據(jù)點(diǎn)的三維軸 向外包圍盒（Axis Aligned Bounding Box，AABB)大小，即三維點(diǎn)云坐標(biāo)的最大最小值， (Xmin, Xmax, Ymin, Ymax, Zmin, Zmax)。計(jì)算方法是，迭代處理所有地面點(diǎn)，對(duì)每個(gè)點(diǎn)，分別比 較該點(diǎn)三維坐標(biāo)（X，Y，Z)與當(dāng)前AABB范圍，若該點(diǎn)落入AABB之外，則對(duì)AABB進(jìn)行更新；
[0086] (312)根據(jù)金字塔層級(jí)信息，DEM格網(wǎng)原點(diǎn)以及行列數(shù)，DEM窗口大小w，正則化約 束權(quán)重λ等信息，并且初始化DEM和彎曲能量柵格結(jié)構(gòu)。正則化約束權(quán)重λ需要與地面 點(diǎn)包含粗差比例成正比，隨著由粗到精的金字塔濾波的逐級(jí)遞進(jìn)，地面點(diǎn)包含噪聲比例也 隨之增大。在最粗的金字塔層級(jí)，假設(shè)無(wú)粗差點(diǎn)，窗口大小為《_，此時(shí)正則化約束權(quán)重λ =0 ;而在最精細(xì)的金字塔層級(jí)，此時(shí)粗差點(diǎn)最多，窗口大小為《_，此時(shí)正則化約束權(quán)重達(dá) 到最大值λ = λ_;中間層級(jí)，則采用線性內(nèi)插方式獲取，g卩
初始 化DEM和彎曲能量柵格圖時(shí)，需要利用AABB平面最小值（Xmin，Ymin)，窗口大小w，以及行 列數(shù)目；
[0087] (313)對(duì)地面點(diǎn)的二維平面坐標(biāo)，構(gòu)建KD-樹的空間索弓丨，用于迭代處理DEM格 網(wǎng)的近鄰索引，內(nèi)插對(duì)應(yīng)的格網(wǎng)信息；具體方法為，根據(jù)獲取所有地面點(diǎn)的二維平面坐標(biāo) (X，Y)，采用開源算法FLANN建立二維地面點(diǎn)數(shù)據(jù)的KD-樹索引。其中KD-樹中的索引序號(hào) 與原始三維點(diǎn)序號(hào)相同；
[0088] (314)對(duì)每一 DEM格網(wǎng)，求取該格網(wǎng)處的二維平面坐標(biāo)（X，Y)，在KD-樹索引中按 一定規(guī)則搜索一定數(shù)目的近鄰點(diǎn)，并且獲取其對(duì)應(yīng)的三維坐標(biāo)（X，Υ，Ζ)，作為樣條函數(shù)的控 制點(diǎn)；在搜索過程中需要滿足一下條件：1)樣條函數(shù)控制點(diǎn)應(yīng)分布均勻；2)控制點(diǎn)數(shù)目應(yīng) 滿足最小要求；3)控制點(diǎn)應(yīng)選和該格網(wǎng)最相關(guān)的三維點(diǎn)，即和格網(wǎng)點(diǎn)最靠近的點(diǎn)。為滿足 上述要求，具體實(shí)時(shí)方法如下：
[0089] (3141)計(jì)算DEM格網(wǎng)點(diǎn)處的二維平面坐標(biāo)（X，Υ);
[0090] (3142)獲取當(dāng)前窗口大小w，并且以此初始化一個(gè)數(shù)組r = [w，2w，3w，4w，5w]，對(duì) 應(yīng)半徑搜索窗口大??；
[0091] (3143)從最小的搜索半徑r = W開始，利用KD-樹索引，搜索半徑r中包含的點(diǎn)的 索引號(hào)，并且獲取對(duì)應(yīng)的三維坐標(biāo)（X，Y，Z);
[0092] (3144)將所有三維點(diǎn)分為4個(gè)象限，每個(gè)象限中的三維點(diǎn)，按與DEM格網(wǎng)點(diǎn)距離大 小由小到大排序；
[0093] (3145)若每象限點(diǎn)數(shù)均大于等于4,則每個(gè)象限選取離格網(wǎng)點(diǎn)最近的三維點(diǎn)，作 為樣條函數(shù)控制點(diǎn)；
[0094] 若其中一個(gè)象限點(diǎn)數(shù)小于4,則增大搜索半徑r = Γι+1，并且重復(fù)步驟（3143)~ (3145)，直至 r = 5w ;
[0095] 若搜索半徑r = 5w，近鄰點(diǎn)仍小于9,則將該格網(wǎng)點(diǎn)設(shè)置為無(wú)效之-9999,并跳過后 續(xù)步驟（315)和（316);
[0096] (3146)返回搜索到的近鄰點(diǎn)的索引號(hào)集合，并且獲取對(duì)應(yīng)的三維坐標(biāo)，用于計(jì)算 樣條函數(shù)的控制點(diǎn)。
[0097] (315)傳統(tǒng)的DEM內(nèi)插，通常是通過擬合地面控制點(diǎn)P = {pi = (Xi，Yi，Zi) I i = 1，2, 3,…，η}，得到一張參數(shù)化的曲面Z = f (X，Y)，在擬合過程中，由于地面控制點(diǎn)中不包 含粗差，因此，僅需要保持對(duì)地面控制點(diǎn)有較好的擬合度，就可取得較好的效果。即使得數(shù) 據(jù)擬合度edata(data term)具有較小的值，或完全經(jīng)過所有控制點(diǎn)，如下式所示
[0099] 然而，在點(diǎn)云濾波過程中，地面點(diǎn)不可避免的會(huì)引入非地面粗差點(diǎn)，因此此時(shí)即使 ε ^3為〇,也并不意味著對(duì)地面有較好的擬合程度。采用地表一致性約束，可以一定程度 上，降低噪聲敏感性，提高內(nèi)插算法的魯棒性。彎曲能量es_th，如下式所示，
[0101] 可以表達(dá)一張參數(shù)化曲面的二階連續(xù)程度，因此若將ε s_th作為一個(gè)正則化約束 引入地面內(nèi)插過程中，則可隱式的蘊(yùn)含著地表一致性約束，抵抗噪聲點(diǎn)的干擾。即通過優(yōu)化 具有ε = edata+X es_th形式的能量函數(shù)，進(jìn)行DEM內(nèi)插。
[0102] 給定控制點(diǎn)集P，該優(yōu)化模型為的解為薄板樣條函數(shù)，即下式所示，
[0104] 為獲取該樣條函數(shù)的參數(shù)羅和歹，，可通過以下形式的線性方程組求解：
[0106] 其中，上述矩陣中的各元素具有以下形式，并且可直接采用QR分解的方式獲取最 終的參數(shù)聲和α &
CN 105118090 A I兄明 8/9 頁(yè)
[0112] (316)利用樣條函數(shù)參數(shù)，計(jì)算DEM格網(wǎng)對(duì)應(yīng)二維平面坐標(biāo)（X，Y)處的高程值z(mì)并 保存入DEM的柵格結(jié)構(gòu)Z中；其具體步驟如下：
[0113] (3161)計(jì)算DEM高程值：將DEM格網(wǎng)點(diǎn)坐標(biāo)（X，Υ)和樣條函數(shù)參數(shù)彥:和足，帶入 函數(shù)Z (Χ，Υ)之中，即求得EffiM對(duì)應(yīng)格網(wǎng)處的高程值，并存入DHM對(duì)應(yīng)的行列號(hào)之中；
[0114] (3162)計(jì)算DEM格網(wǎng)處的彎曲能量：根據(jù)樣條函數(shù)參數(shù)#以及矩陣Κ，計(jì)算樣條函 數(shù)的彎曲能量二戶'/φ 5并存入彎曲能量對(duì)應(yīng)的行列號(hào)之中；
[0115] (3163)直至計(jì)算完所有格網(wǎng)點(diǎn)的信息，并將計(jì)算完成的DEM與彎曲能量傳出至下 一階段處理。
[0116] (317)重復(fù)步驟（314)至（316)直至迭代處理所有DEM格網(wǎng)，獲取最終的DEM與彎 曲能量柵格圖。
[0117] 如圖5所示，步驟（32)中自適應(yīng)復(fù)雜地形結(jié)構(gòu)的參數(shù)優(yōu)選方法包括以下步驟：
[0118] (321)依據(jù)金字塔層級(jí)信息s，計(jì)算當(dāng)前層級(jí)對(duì)應(yīng)的尺度補(bǔ)償閾值Ts ;考慮到 地形坡度影響，在較粗的金字塔層級(jí)，由于窗口 w較大，即使平坦區(qū)域也可能會(huì)存在微小 的地形起伏，因此需要根據(jù)金字塔層級(jí)S和窗口大小W補(bǔ)償?shù)匦纹露绕鸱跋瘛＞唧w 而言，在最粗的金字塔層級(jí)，補(bǔ)償閾值最大為Ts = Ts_max ;在最精細(xì)的金字塔層級(jí)，補(bǔ) 償閾值為〇;而在中間層級(jí)的數(shù)據(jù)，補(bǔ)償閾值根據(jù)窗口大小w采用線性內(nèi)插方法獲取，BP
[0119] (322)將彎曲能量柵格圖采用分段線性內(nèi)插轉(zhuǎn)化為對(duì)應(yīng)的彎曲能量補(bǔ)償閾值Tb。 顧及彎曲能量的參數(shù)自適應(yīng)優(yōu)選方法的基本思想是，在地形起伏較大的地方，如山脊線、地 面斷裂處等區(qū)域，濾波閾值應(yīng)更大；而在地形起伏較小的區(qū)域，如平坦地表等區(qū)域，濾波閾 值應(yīng)更小。由于在地形欺負(fù)較大的地方，地形變化彎曲能量也相應(yīng)的會(huì)變大，因此濾波閾值 應(yīng)與彎曲能量成正比例變化。因此，其具體獲取彎曲能量補(bǔ)償閾值Tb的具體步驟為；
[0120] (3221)對(duì)所有DEM格網(wǎng)點(diǎn)的彎曲能量從小至大排序；
[0121] (3222)將所有彎曲能量均分為10段，并記錄分割處的彎曲能量大小
[b0, bu....., b10]；
[0122] (3223)設(shè)置一最大補(bǔ)償閾值Tbmax，對(duì)所有彎曲能量值b，進(jìn)行分段內(nèi)插，若其落入 第i段之內(nèi)，貝1J有
[0123] (323)內(nèi)插待分類點(diǎn)所處格網(wǎng)的高程均值Zni與彎曲能量補(bǔ)償值T b。由于DEM在濾 波過程中可能會(huì)存在一定噪聲，因此不能僅用該格網(wǎng)處的DEM高程值，而考慮采用3X3鄰 域DEM的高程值z(mì)和彎曲能量補(bǔ)償值Tb。首先，根據(jù)待分類點(diǎn)坐標(biāo)（X，Y，Z)，計(jì)算其所對(duì)應(yīng) 的格網(wǎng)行列號(hào)，r= (Y-Ynin)/w，C= (X-Xnin)/￥。隨后，在DEM和補(bǔ)償閾值中獲取對(duì)應(yīng)3X3 鄰域內(nèi)的高程均值z(mì)n和補(bǔ)償閾值T b;
[0124] (324)計(jì)算自適應(yīng)濾波閾值，該閾值綜合考慮金字塔層級(jí)的尺度信息s，彎曲能量 信息b，以及初始濾波閾值t，故最終的濾波閾值應(yīng)是上述三者的函數(shù)T = f (s，b，t)，即T = Tb+Ts+t ;
[0125] (325)依據(jù)待分類點(diǎn)高程，DEM高程與濾波閾值之間關(guān)系判斷該點(diǎn)為地面或非地 面點(diǎn)；在3x3鄰域內(nèi)，判斷該待分類點(diǎn)和鄰域內(nèi)的9個(gè)格網(wǎng)點(diǎn)的相互關(guān)系。記N為待分類點(diǎn) 高程值Z〈z+T的數(shù)目，若N多5,則該待分類點(diǎn)為地面點(diǎn)；否則仍然保留為待分類點(diǎn)，在后續(xù) 處理中繼續(xù)