一種自適應(yīng)復(fù)雜地形結(jié)構(gòu)的點(diǎn)云濾波方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于地理空間信息系統(tǒng)技術(shù)領(lǐng)域，主要涉及一種自適應(yīng)復(fù)雜地形結(jié)構(gòu)的點(diǎn) 云濾波方法。
【背景技術(shù)】
[0002] 隨著多平臺激光雷達(dá)和多角度影像密集匹配技術(shù)的快速發(fā)展，大量快速可得的點(diǎn) 云數(shù)據(jù)的高分辨率高精度的地形表面三維重建提供了有效可靠的數(shù)據(jù)支撐，點(diǎn)云數(shù)據(jù)處理 已成為國際研究和工業(yè)界的熱點(diǎn)問題。如何從點(diǎn)云數(shù)據(jù)中分離地面與非地面點(diǎn)，即點(diǎn)云數(shù) 據(jù)濾波，是眾多后續(xù)應(yīng)用的首要步驟，因此一直是熱點(diǎn)前沿問題。
[0003] 在實(shí)際生產(chǎn)實(shí)踐中，由于地表結(jié)構(gòu)的多樣性和復(fù)雜性，同一套濾波參數(shù)難以適應(yīng) 錯綜復(fù)雜的地形特征，當(dāng)前國內(nèi)外點(diǎn)云處理的主流商業(yè)軟件，如TerraScan、LasTools等， 其自動處理精度難以達(dá)到實(shí)際生產(chǎn)需求，通常采用半自動濾波加人工交互后處理的方法進(jìn) 行濾波處理，需要耗費(fèi)大量的人力、物力；并且由于人為主觀因素影響，有時還需要諸如真 正射影像等信息輔助處理，并可能出現(xiàn)不同程度的誤差，極大的影響了生產(chǎn)效率。因此，自 動、高效、穩(wěn)健的點(diǎn)云濾波成為學(xué)術(shù)界和工業(yè)界面臨的重大挑戰(zhàn)。
[0004] 通常點(diǎn)云濾波需要假設(shè)在局部鄰域內(nèi)的最低點(diǎn)為地面點(diǎn)，通過建立臨時的地面結(jié) 構(gòu)，通過一定閾值參數(shù)對地面點(diǎn)進(jìn)行加密，并區(qū)分地面-非地面點(diǎn)。該方法主要存在缺陷包 括：一方面，地形起伏變化、城市結(jié)構(gòu)形態(tài)復(fù)雜；另一方面多平臺激光雷達(dá)和多角度影像密 集匹配產(chǎn)生的點(diǎn)云數(shù)據(jù)，無論在點(diǎn)云密度還是數(shù)據(jù)特性上均有較大差異；而現(xiàn)有的濾波算 法通常對同一區(qū)域采用一套參數(shù)，難以適應(yīng)錯綜復(fù)雜的地面-非地面物體結(jié)構(gòu)，并且通常 對參數(shù)變化較為敏感，導(dǎo)致誤分類問題嚴(yán)重。鑒于此，構(gòu)建對地面起伏變化具有自適應(yīng)能力 的濾波參數(shù)選擇方法，成為可能的有效途徑之一。另一方面由于點(diǎn)云中的噪聲、地面點(diǎn)加密 過程中的誤差，在由粗到精的金字塔濾波過程中，難免會受誤分類點(diǎn)的影響，如不有效處理 這些誤分類點(diǎn)將在漸進(jìn)式的濾波過程中被不斷放大，大幅增加人工交互式后處理工作量。 因此，如何進(jìn)行點(diǎn)云內(nèi)插使其具有良好的局部抗噪性是目前亟待解決的問題。
[0005] 目前，現(xiàn)有的主流點(diǎn)云濾波算法可以分為兩類：基于內(nèi)插的濾波方法和基于形態(tài) 學(xué)運(yùn)算的濾波方法：
[0006] (1)基于內(nèi)插的濾波方法
[0007] 基于內(nèi)插的濾波方法的基本思想是，首先在一較大的窗口內(nèi)，提取最低點(diǎn)作為地 面種子點(diǎn)，構(gòu)建不規(guī)則三角網(wǎng)（TIN)或數(shù)字地表模型（DEM)用來表達(dá)地面，然后根據(jù)適當(dāng)閾 值漸進(jìn)的對地面點(diǎn)進(jìn)行加密。該方法可十分簡潔的擴(kuò)展為由粗到精的金字塔濾波策略，點(diǎn) 云金字塔中每一層級存儲的是對應(yīng)窗口內(nèi)的最低值。然而，在不斷內(nèi)插地表模型與迭代濾 波過程中，地面點(diǎn)不可避免包含粗差點(diǎn)，若不有效處理將被不斷放大，嚴(yán)重影響最終濾波結(jié) 果。目前，Terrascan采用一種漸進(jìn)三角網(wǎng)加密的濾波策略，LasTools則采用線性預(yù)測的濾 波策略，兩者均屬于此類方法。
[0008] (2)基于形態(tài)學(xué)運(yùn)算的濾波方法
[0009] 基于形態(tài)學(xué)的濾波方法的基本思想來源于數(shù)學(xué)上的形態(tài)學(xué)計(jì)算一一用形態(tài)學(xué)開 運(yùn)算和閉運(yùn)算從灰度影像上識別目標(biāo)。建立在最低高程假設(shè)之上，非地面點(diǎn)的高程通常高 于周圍的地面點(diǎn)，如果將點(diǎn)云數(shù)據(jù)按高程值轉(zhuǎn)化為灰度影像，由于灰度差異，非地面目標(biāo)將 被識別出來。對點(diǎn)云數(shù)據(jù)進(jìn)行開運(yùn)算，可以將局部窗口內(nèi)高程值最小的點(diǎn)識別出來。通過判 斷點(diǎn)云與窗口內(nèi)的最小高程點(diǎn)之間的坡度差異與角度差異作為濾波準(zhǔn)則，移動濾波窗口， 實(shí)現(xiàn)對全部點(diǎn)云數(shù)據(jù)的濾波?；谛螒B(tài)學(xué)運(yùn)算的濾波方法的關(guān)鍵點(diǎn)在于選取一個合適的運(yùn) 算窗口：若選擇窗口過小，則只有較小的非地面目標(biāo)（比如樹木、汽車）能被濾除，而對于較 大地物則會在膨脹中再次保留；而若選擇窗口過大，則在坡度較大區(qū)域，部分地面點(diǎn)也會被 誤分為非地面點(diǎn)。該方法雖然效率較高，然而由于對參數(shù)選擇十分敏感，難以處理復(fù)雜地形 點(diǎn)云，實(shí)用性不大。

【發(fā)明內(nèi)容】

[0010] 本發(fā)明的目的在于克服上述缺陷，提供一種將能顯式刻畫地形變化特性的彎曲能 量引入到點(diǎn)云數(shù)據(jù)濾波過程中，使濾波參數(shù)隨地形變化在局部范圍內(nèi)自適應(yīng)優(yōu)化調(diào)整；并 將地表一致性約束引入內(nèi)插算法之中，采用彎曲能量的規(guī)則化約束，克服噪聲對濾波精度 的影響，并建立一套科學(xué)合理的點(diǎn)云自動濾波實(shí)施方案。
[0011] 為了實(shí)現(xiàn)上述目的，本發(fā)明采用的技術(shù)方案如下：
[0012] -種自適應(yīng)復(fù)雜地形結(jié)構(gòu)的點(diǎn)云濾波方法，包括以下步驟：
[0013] (1)構(gòu)建點(diǎn)云金字塔；
[0014] (2)采用格網(wǎng)窗口最大的金字塔層級數(shù)據(jù)初始化地面點(diǎn)G，并將下一級數(shù)據(jù)加入 待分類點(diǎn)U，更新待分類點(diǎn)；
[0015] (3)迭代處理點(diǎn)云金字塔的所有層級，處理完所有待分類點(diǎn)，并且保存新增加地面 點(diǎn)
[0016] (31)利用已有地面點(diǎn)，采用具有局部抗噪性的DEM內(nèi)插算法獲取地面DEM，并同時 獲取對應(yīng)區(qū)域的彎曲能量；
[0017] (32)根據(jù)DEM、彎曲能量以及金字塔尺度信息，采用自適應(yīng)復(fù)雜地形結(jié)構(gòu)的參數(shù) 優(yōu)選方法獲取濾波閾值并進(jìn)行濾波，區(qū)分地面點(diǎn)與非地面點(diǎn)；
[0018] (33)重復(fù)步驟（31)~（32)，直至無新增地面點(diǎn)；
[0019] (4)存儲分類后的地面點(diǎn)與非地面點(diǎn)，其中，地面點(diǎn)標(biāo)記為2,非地面點(diǎn)標(biāo)記為1。
[0020] 在構(gòu)建點(diǎn)云金字塔之前，剔除粗差數(shù)據(jù)，其具體方式如下：
[0021] 首先，對點(diǎn)云數(shù)據(jù)構(gòu)建KD-樹的索引結(jié)構(gòu)；然后，對每個點(diǎn)進(jìn)行半徑搜索，計(jì)算在5 米半徑值內(nèi)的點(diǎn)的數(shù)目，統(tǒng)計(jì)其均值μ和標(biāo)準(zhǔn)差σ，將鄰域內(nèi)的點(diǎn)數(shù)目小于μ-3 O的點(diǎn)視 為孤立點(diǎn)，進(jìn)行移除。
[0022] 所述步驟（1)的具體方式如下：
[0023] (11)根據(jù)用戶輸入的最大窗口大小Wmax與最小窗口大小Wmin，均勾獲取一系列金 字塔窗口大小w ;
[0024] (12)從最粗層級開始，對原始數(shù)據(jù)進(jìn)行規(guī)則格網(wǎng)劃分，每個格網(wǎng)中記錄該格網(wǎng)的 最低點(diǎn)，并從原始點(diǎn)云中移除；
[0025] (13)迭代處理所有層級，建立點(diǎn)云金字塔。
[0026] 所述步驟（31)中具有局部抗噪性的DEM內(nèi)插算法包括以下步驟：
[0027] (311)在獲取所有地面點(diǎn)之后，迭代處理所有地面點(diǎn)，確定所有數(shù)據(jù)點(diǎn)的三維軸向 外包盒大小，用于后續(xù)確定DEM范圍、行列數(shù)等信息；
[0028] (312)初始化EffiM和彎曲能量柵格結(jié)構(gòu)，并根據(jù)金字塔層級信息，計(jì)算正則化約束 權(quán)重λ ;DEM格網(wǎng)原點(diǎn)、行列數(shù)、DEM窗口大小w等信息；該正則化約束權(quán)重與地面點(diǎn)中包含 的粗差點(diǎn)比例有關(guān)：若粗差比例越大，則權(quán)重越大；反之權(quán)重越小，直至λ = 〇
[0029] (313)對地面點(diǎn)的二維平面坐標(biāo)，構(gòu)建KD-樹的空間索引，用于迭代處理DEM格網(wǎng) 的近鄰索引，內(nèi)插對應(yīng)的格網(wǎng)信息；
[0030] (314)對每一 DEM格網(wǎng)，求取該格網(wǎng)處的二維平面坐標(biāo)（X，Y)，在KD-樹索引中按 自適應(yīng)半徑增長方法搜索至少9個近鄰點(diǎn)，并且獲取其對應(yīng)的三維坐標(biāo)（X，Υ，Ζ)，作為樣條 函數(shù)的控制點(diǎn)；
[0031] (315)根據(jù)控制點(diǎn)三維坐標(biāo)，更新樣條函數(shù)矩陣，用于計(jì)算樣條函數(shù)參數(shù)，樣條函 數(shù)的解具有線性方程的閉合解形式，采用QR分解求解樣條函數(shù)的參數(shù)信息；
[0032] (316)利用樣條函數(shù)參數(shù)，計(jì)算DEM格網(wǎng)對應(yīng)二維平面坐標(biāo)（Χ，Υ)處的高程值ζ并 保存入DEM的柵格結(jié)構(gòu)Z中；
[0033] (317)重復(fù)步驟（314)至（316)直至迭代處理所有DEM格網(wǎng)，獲取最終的DEM與彎 曲能量柵格圖。
[0034] 所述步驟（313)的具體方式如下：根據(jù)獲取所有地面點(diǎn)的二維平面坐標(biāo)（X，Υ)，采 用開源算法建立二維地面點(diǎn)數(shù)據(jù)的KD-樹索引；其中KD-樹中的索引序號與地面點(diǎn)序號相 同。
[0035] 所述步驟（314)的具體方式如下：
[0036] (3141)計(jì)算DEM格網(wǎng)點(diǎn)處的二維平面坐標(biāo)（X，Υ);
[0037] (3142)獲取當(dāng)前窗口大小w，并且以此初始化一個數(shù)組r = [w, 2w, 3w, 4w