融合無人機LiDAR和高分影像的路面平整度監測方法與流程
本發明涉及公路路面平整度遙感監測技術，是利用無人機激光雷達(lightdetectionandranging,lidar)和高分辨率遙感技術針對公路路面進行大范圍、高效率、低成本的路面平整度監測與評估。
背景技術：
機載激光雷達(airbornelidar)是一種主動式對地觀測系統，集全球導航衛星系統、慣性導航、激光測距及計算機處理技術于一身，能夠直接獲取高精度、高分辨率的數字地面模型以及地物的三維空間信息，具有傳統攝影測量方法無法取代的優越性，并在地形測繪、道路勘測設計、植被參數反演、災害監測、城市規劃及三維建模等領域得到了較為廣泛的應用。公路作為一種現代化的運輸通道在社會經濟中發揮著越來越重要的作用，對沿線的物流、資源開發、招商引資、產業結構調整、橫向經濟聯合等起到積極推動作用，因而其質量好壞直接影響到地區和國家的經濟發展。文獻[1](sunl,zhangzm,ruthj.modelingindirectstatisticsofsurfaceroughness[j].journaloftransportationengineering-asce,2001,127(2):105-111)中，1960年美國國家公路與交通運輸行政聯合會開展的道路試驗表明：約95％的路面服務性能取決于道路表面的平整度。路面平整度作為評價路面質量的主要指標之一，也是執行公路路面養護管理的重要決策依據。據《2015年交通運輸行業發展統計公報》統計，2015年中國公路總里程達到457.73萬公里，公路養護里程446.56萬公里，占公路總里程97.6％，且呈逐年快速增長趨勢。隨著我國公路建設的飛速發展，如何實現大范圍、高效率、低成本的公路路面健康狀況監測就成為一個亟待解決的問題。
平整度監測是常用的公路路面健康狀況監測與評價的重要指標之一。平整度評價指標的測量有3米直尺測定最大間隙(h)、國際平整度指數(internationalroughnessindex,iri)、行使質量指數(runningqualityindex,rqi)、顛簸累積值(bumpcumulativevalue,vbi)、平整度標準差(σ)、功率譜密度(powerspectraldensity,psd)等。現有傳統的路面平整度監測方法有：3m直尺、連續式平整度儀、手推式斷面儀、車載式顛簸累積儀以及近年出現的激光平整度檢測儀(包括集成激光平整度儀的道路檢測車)，其各自的特點概括于表1。這些方法的共性之一就是只能依靠地面測量獲得為數很少的道路縱剖面信息，對道路平整度的描述也較粗略，而且無法快速獲取道路表面的三維幾何信息，使得管理部門對道路整體狀況的了解非常有限。
表1.常用的路面平整度監測方法
低空無人機lidar技術的出現不僅為獲取道路表面三維信息提供了可能，而且該技術較少受地形限制、高精度、高效率、低成本(尤其是大范圍監測)、可同時獲取地面遙感影像的優勢更是為路面平整度監測提供了一種新的有效解決方案。目前利用lidar技術進行路面平整度監測的相關案例和研究少見。文獻[2](china,olsenmj.evaluationoftechnologiesforroadprofilecapture,analysis,andevaluation[j].journalofsurveyingengineering,2014,141(1):4014011)記載了利用地面激光雷達技術對高速公路進行了平整度監測，但是，該方法并未對路面平整度的空間可視化方法進行嘗試，也未對點云的分類方法和路面表面模型的構建方法進行探索，該技術的另外一個致命點就是只能在某個固定點位置附近處獲取小范圍(半徑<50m)的點云數據。如果要獲取大范圍的路面信息則需大量的點云拼接、坐標糾正、點云配準、坐標轉換等一系列復雜、誤差不斷積累的后處理過程，不僅流程復雜而且非常低效。文獻[3](郭姣.基于車載激光的道路平整度檢測系統研究[d].首都師范大學,2013)記載了利用車載激光雷達技術分析道路剖面平整度特征并定位道路損壞超標位置，但是，該方法同樣沒有對點云的分類方法進行探索，而是借助傳統的點云處理平臺來完成分類工作。此外，該方法僅僅依靠單個數值進行表達，未能對路面質量的空間分布特征進行探索，并要求檢測車行駛速度穩定、gps信號穩定，而這在某些建筑物遮擋、崎嶇不平的復雜路段難以應用。
綜上所述，現有傳統的公路路面平整度監測方法大多存在效率低下、自動化程度低，無法進行大范圍全路段監測，并且只能獲得為數很少的道路縱剖面信息，不能得到道路表面的三維結構信息。現有的激光雷達技術的平整度監測方法則存在測量范圍小、作業效率低、缺乏空間可視化表達、數據源單一的缺陷，并且會受到地形條件限制。然而公路路面平整度監測迫切需要探索新的技術途徑，提高大范圍路面平整度監測效率與準確性。
技術實現要素：
針對現有公路路面平整度監測方法中的不足，本發明提出一種融合無人機lidar和高分辨率影像的路面平整度監測方法及其實施方案。首先，對無人機lidar系統獲取的點云數據和高分影像數據進行預處理，然后從高分影像數據提取光譜特征，對點云數據提取多尺度幾何特征。將所得到的幾何特征和光譜特征作為點云分類的特征變量，并對這些特征變量進行去冗余和降維操作，在特征級別上對兩類數據進行融合，利用機器學習中的隨機森林分類算法進行點云分類；再進一步通過濾波操作和粗差糾正得到比較精準的道路點云數據集。然后利用三維插值法構建高精度的精細路面數字表面模型(digitalsurfacemodel,dsm)。基于iri指數原理計算路面平整度，再利用gis進行空間可視化表達與平整度質量評價。相比于以往的路面平整度監測技術只能用單個iri值來表達整個路段信息而言，該方法則可給出某條路不同路段不同車道iri值的空間分布，提高了路面養護管理工作決策的科學性，有效提高了無人機遙感技術在公路養護管理應用中的實用性，可廣泛應用于城市道路及其他低等級道路的路面平整度監測需要。
本發明的原理是：無人機lidar遙感系統獲取的點云數據中包含公路及周邊地物的三維坐標信息，通過選擇合適的點云去噪手段可有效去除噪聲干擾，而同期獲得的高分辨率影像數據富含地物的光譜信息且非常易于目視判讀。點云數據則因為點集分散而不易直觀判讀，為此將高分影像與點云數據進行融合就能同時將光譜信息與空間幾何信息進行集成，為點云分類提供更多的可利用信息。目前，針對高分影像的地物分類存在一個挑戰，就是對“異物同譜”的地物目標無法區分，而點云數據的空間幾何信息則可以對此起到很好的補充。同樣的，點云數據自身由于點集非常離散而導致在一些邊緣交界處以及幾何結構特征非常相似的地物之間出現混淆，而高分影像則能對此進行很好的輔助作用。基于此，本發明一方面從高分影像提取典型地物類別的光譜特征，另一方面從點云上提取典型地物類別的幾何特征，通過構建多尺度幾何特征來強化地物類型的空間形態差異，然后采用高效且抗噪聲強的機器學習分類算法對點云進行分類，繼而對分類后的點云進行粗差糾正和類別濾波即可得到較精準的道路點云數據，再結合點云內插算法構建精細的道路dsm。在此基礎上，從道路表面選取一系列縱剖面曲線，分段計算國際平整度值并進行空間可視化表達，最后對平整度值進行分等定級，從而實現對道路平整度質量的評價。本發明方法經過精度評價證明可用來進行大范圍公路路面的快速自動化平整度質量狀況監測。
本發明提供的技術方案：
一種融合無人機lidar與高分影像的路面平整度監測方法，包括：數據獲取與預處理、特征提取與融合、點云分類、路面dsm模型的構建及路面平整度評價等技術環節。
其中，數據獲取與預處理的目的是對低空無人機lidar系統獲得的點云數據進行去噪，對同期獲取的高分影像進行配準與拼接，得到實驗區的正射影像數據，為數據融合做準備。數據融合的目的是將點云數據的幾何空間信息與高分影像數據的光譜信息進行集成，從而使得點云數據同時具備空間幾何信息和光譜信息。點云分類的目的是為了得到道路點云數據，數據融合的過程為點云分類提供了可利用信息，對實驗區進行典型地物類別劃分，構建一系列幾何特征參數來表達地物類型的差異，并結合光譜特征參數就構成了分類器的特征變量，接著使用隨機森林算法實現對點云數據的分類，得到路面點云數據。路面dsm模型的構建可以為路面三維空間結構的表達提供基礎，而且dsm模型的每個柵格單元存儲的都是高程信息，便于平整度計算，也為平整度的空間可視化表達提供了條件。路面平整度評價的目的就是通過選取一定的平整度指標對路面的平整度狀況進行量化描述，并建立分等定級體系對實驗區的道路狀況進行客觀評價，從而對路面當前的整體狀況給出前瞻性分析和決策。
本發明提供的路面平整度監測方法具體包括以下步驟：
1)利用低空無人機lidar系統(包括：無人機平臺、lidar掃描儀、多光譜相機、姿態定位單元、自動駕駛儀等單元的集成)，獲取實驗區的lidar點云數據和高分辨率影像數據；
2)對獲取的lidar點云數據轉換為當前通用的點云數據文件格式(*.las)，然后對點云數據和無人機拍攝的高分辨率影像數據進行預處理，包括點云的去噪、影像的拼接與配準處理及正射影像數據的生成；
3)對步驟2)中得到的正射影像的波段進行分析并提取反映地物差異的光譜特征參數，連同高分辨率影像原來的波段，構成包含光譜信息的多個光譜特征集；
4)對得到的光譜特征數據與點云數據進行特征融合，并以點云為空間參考基準進行地理空間配準；
5)對實驗區的主要地物類型進行劃分，并在步驟4)得到的點云數據基礎上經過抽稀得到抽稀的點云數據，然后對抽稀的點云數據選取一系列的幾何特征參數進行計算，并通過改變空間分析的尺度得到包含多尺度幾何特征參數和光譜特征參數的點云數據，接著從中選取部分相應地物類型的點云數據作為后續點云分類的訓練樣本，將剩余的抽稀點云數據作為待分數據。點云的抽稀過程減少了幾何特征參數計算過程中的運算量，同時也降低了分類所需處理的數據量，大大提高了數據處理效率；
6)步驟5)中得到的特征變量通常會很多，不可避免的會出現某些變量之間存在較強的相關性，導致特征變量冗余。為了找到關鍵特征變量，從而對數據起到降維作用，進而簡化分類器構建過程的復雜度，在特征級別上，融合多光譜影像的光譜信息與激光點云數據的空間幾何信息，構成新的點云分類特征集；
7)利用隨機森林模型基于步驟6)的特征集對步驟5)的訓練樣本進行模型訓練，通過不斷的調整模型參數取得滿意的分類效果后，應用于待分類樣本點云數據進行分類，然后將樣本點云的類別屬性按照空間最近鄰插值原理賦給經去噪后的點云數據，從而完成全部點云的分類。分類結果可能還包含不屬于路面的點云，此時要對分類后的點云數據進行人工檢查和粗差糾正，并進一步通過類別屬性過濾即得到道路點云數據；
8)對步驟7)得到的道路點云數據的高程信息，采用不規則三角網模擬路面起伏變化，然后通過自然鄰域插值法構建規則格網的高精度精細道路路面的dsm；
9)從步驟8)中得到的dsm模型數據中選取一系列的縱剖面線，運行iri指數模型計算各剖面線對應的iri值并經過柵格化處理，得到路面的iri值的空間分布圖；
10)對步驟9)中得到的iri值，依據道路平整度的分等定級標準進行分類，從而得到待測路段對應的平整度質量評價結果。
針對上述公路路面平整度監測方法，進一步地，步驟2)中涉及的lidar點云數據，具體為滿足以下4個條件的點云數據：
a)坐標系：點云數據的坐標系為wgs84坐標系(經度、緯度)或者與當地相匹配的地方平面投影坐標系(x，y)；
b)屬性信息：點云數據屬性信息除經度、緯度外，還應包括高程、掃描角、反射強度；
c)點云密度：點云密度在400點/平方米以上(即激光腳點半徑在50mm以內)；
針對上述平整度監測方法，進一步地，對步驟2)中涉及的數據預處理過程，主要是指以下4種：lidar點云的噪聲點過濾、lidar點云的高程異常點過濾、lidar點云的掃描角過濾與高分辨率影像的拼接與配準。針對上述平整度監測方法，進一步地，對步驟2)中涉及的高分辨率影像數據，是指市面上性能良好的多光譜數碼相機所拍攝的高清數碼照片，并且無人機能提供相匹配的姿態與定位(pos)信息文件。針對上述平整度監測方法，進一步地，對步驟3)中涉及的光譜特征參數，為以下3類：
a)高分影像的紅、綠、藍波段的灰度值。
b)高分影像的紅、綠、藍波段的灰度值之標準差；
c)激光點云的反射強度值(lidar掃描儀發射的激光脈沖也屬于電磁波的一種，同樣具有波長信息)。
針對上述平整度監測方法，進一步地，對步驟4)中涉及的特征融合，具體是指將點云數據與高分辨率影像數據的光譜特征參數進行融合，并且融合的過程是以點云數據的空間參考為基準的，采用仿射逆變換模型將影像上每個像元的光譜值匹配到點云數據坐標系；
針對上述平整度監測方法，進一步地，對步驟5)中涉及的幾何特征參數，為以下3類：
a)局部粗糙度(localdegreeofroughness,ldr)。該特征變量的具體含義是指點云中的某個參考點到該參考點所在的某個空間尺度下的鄰域點所形成的最佳擬合平面的距離；
b)局部維度特征(localdimensionfeature,ldf)。文獻[4](brodun,lagued.3dterrestriallidardataclassificationofcomplexnaturalscenesusingamulti-scaledimensionalitycriterion:applicationsingeomorphology[j].isprsjournalofphotogrammetryandremotesensing,2012,68:121-134)記載了局部維度特征ldf，該特征變量的具體含義是指點云中的某個參考點與該參考點的鄰域點之間在空間上是以一維、二維還是三維分布的特征度量值，用一個數值來描述它分屬于這三種維度的大小程度，該值就是局部維度特征參數；
c)局部高程差(localheightdifference,lhd)。該特征變量的具體含義是指點云中某個參考點與該參考點在某個空間尺度下的鄰域點所構成的全體點集中，到達該點集所形成的最佳擬合面的最大與最小距離之差，即為局部高程差，起伏越大的平面該值也越大。
針對上述平整度監測方法，進一步地，對步驟5)中涉及的多尺度幾何特征，其具體含義是指，在計算某個幾何特征值時其空間分析尺度有多個(固定某個空間分析尺度會得到一個對應的空間幾何特征參數)，不同地物類型在不同的空間分析尺度下其幾何特征值會有不同。
針對上述平整度監測方法，進一步地，對步驟5)中涉及的點云數據類別，主要包括道路點云(含交通標線)和低矮植被(農作物和草地)、土壤、樹木、移動目標、電力線(包括電線桿)這6大類。
針對上述平整度監測方法，進一步地，對步驟5)中涉及的點云抽稀過程，是在空間中按照一定的空間間隔進行采樣，該過程要保證點云在空間的分布相對于原始點云而言是均勻的(而不是零亂的分布)。
針對上述平整度監測方法，進一步地，對步驟6)中涉及的特征選擇，主要包括特征子集的搜索和特征子集的評價兩個過程。特征子集的搜索采用的是貪心策略，即：初始的最優子集為空，所有待評估的特征變量則組成備選特征集，從備選特征集中篩選解釋能力最強的特征變量加入到最優特征子集中，一次循環只添加一個變量；某一特征變量添加結束后需要從備選特征集中剔除該特征變量，重新下一次循環。最優特征子集的評價標準：新加入的特征變量與已選中的特征變量之間的相關性較低。
針對上述平整度監測方法，進一步地，對步驟8)中涉及的dsm模型的分辨率要盡可能高，并且其高程精度應控制在15mm以內。具體的水平和高程精度評估準則為：從現場的路面中選擇一些具有明顯特征的目標物(如坑槽)并記錄其對應的gps坐標，通過測量其長寬高屬性，然后在高分辨率影像中找到對應的目標物，利用影像與lidar點云的融合，易于在dsm模型中找到該目標的對應位置，進而可利用長度測量工具交互式地測定該目標的長寬高，再統計分析完成對dsm模型的水平和高程的精度評估。
針對上述平整度監測方法，進一步地，對步驟9)中涉及的道路縱剖面線，其間距越小則對應的路面平整度空間分布越精細。iri指數模型的理論依據是四分之一車模型，在計算過程中需要嚴格控制采樣間隔(50mm以內)和測量單元的長度(10m或20m)。最佳的采樣間隔和測量單元長度可通過蒙特卡羅模擬方法定量分析在高程誤差一定時iri的誤差隨采樣間隔和測量單元長度的變化規律來確定。
針對上述平整度監測方法，進一步地，步驟9)中涉及的柵格化是指將路面縱剖面線的iri值按照最鄰近原則賦給周圍的柵格單元，從而得到待測路段的平整度空間分布圖。
針對上述平整度監測方法，進一步地，通過在待測路段選擇一系列包含不同平整度質量等級的路面樣本作為參考數據，可對本方法的評價結果進行精度評定，具體步驟如下：
1)根據設定的平整度質量分級標準，同期在待測路段選擇一系列(包含平整度質量較好與較差，覆蓋健康與病害類型，每種路面選擇一定數量的地面實測樣本)的路面樣本作為參考數據；
2)對參考數據選取外觀、表面粗糙度、車轍深度、病害面積、損壞程度這5項指標使用專家評分法對其進行平整度質量評價，并作為參考值；
3)運用本方法同樣的對此路面實測樣本進行平整度值進行計算，依據平整度質量分級標準，得到本方法的路面平整度質量的評價結果；
4)對平整度質量結果的評價可以類比于遙感影像分類結果的評價，而混淆矩陣是用來評價遙感影像分類精度的最常用形式，因此本方法也采用混淆矩陣對路面平整度質量結果進行評價。引進精度評價指標總體精度(overallaccuracy,oa)和kappa系數對結果進行分析，得到本方法的精度。
與現有技術相比，本發明的有益效果是：
本發明基于低空無人機平臺獲取的lidar數據和高分辨率遙感影像，首次提出融合影像的光譜特征與lidar點云的多尺度幾何特征構建點云分類特征集的方法，解決當前路面點云與非路面點云分類的難題；再利用路面點云構建精細三維路面數字表面模型(dsm)，通過引入國際平整度指數的計算原理，快速計算獲得公路路面的平整度參數，實現路面質量的快速監測與評估。本發明首次提出了一種融合低空無人機lidar和高分影像的公路路面快速高效的平整度監測方法，有效解決了現有路面平整度監測方法中效率低下、自動化程度較低，不適合大范圍全路段監測，且只能獲得為數很少的道路縱剖面平整度信息，不能得到道路三維結構信息等不足。本方法可快速掌握大范圍的路面平整度狀況，從而為路面養護維修決策提供科學依據。同時，本發明設計的點云數據處理流程簡單易實現，效果顯著。另外，本發明提供的平整度質量空間分布圖能全面反映公路路面的平整度狀況，相比于以往的平整度監測技術只能依靠單個數值來表達整個路段的平整度信息而言，有力地促進了道路養護管理的信息化與精細化發展，提高了路面養護管理工作決策的科學性。
本發明融合影像和lidar進行路面平整度監測，與現有技術相比，一是采用了低空無人機lidar數據源，并首次提出融合光譜特征信息與多尺度幾何特征信息來優化道路點云分類，在此基礎上構建精細的dsm，快速提取具有空間分布的路面平整度信息；二是數據處理過程簡潔易實現，使用的數學模型易于理解，效果顯著，所有操作均可在個人普通電腦上完成，對硬件的要求不高；三是在提供可行性的操作流程上，進一步給出了優化數據處理的方法，并通過實驗完成驗證，具有很好的可重復性。
附圖說明
圖1為本發明提供的路面平整度監測方法的流程圖
圖2為局部粗糙度計算方法示意圖。
圖3為局部維度特征的簡化表達示意圖；
其中，(a)為一維空間分布特征，(b)為二維空間分布特征，(c)為三維空間分布特征。
圖4為局部高程差計算示意圖。
圖5為本方法中進行點云插值的自然鄰域插值算法示意圖。
圖6為本方法中使用的國際平整度指數(iri)的算法示意圖；
其中，zs表示簧載質量的垂直位移，zu表示非簧載質量的垂直位移，ms表示簧載質量的大小，mu表示非簧載質量的大小，ks表示連接簧載質量和非簧載質量的彈簧的勁度系數，cs表示連接簧載質量和非簧載質量的線性阻尼系數，ks表示連接輪胎的彈簧的精度系數，可模擬輪胎所產生的減震效果，y(x)表示路面高程，b代表的是輪胎對地面的包容長度(輪胎與地面接觸的部分)。
圖7為本發明實施例中地面采集的不同平整度特征的路面樣本照片
其中，(a1)坑槽，特征為石料暴露且呈大塊狀分布，形狀近似為圓形或者橢圓形，平整度質量評價為極差(failed)；(a2)坑槽，特征為石料暴露但面積相對較小，形狀近似圓形或橢圓形，平整度質量評價為極差(failed)；(b1)塌陷，特征為條帶狀分布且深度較大(石料沒有暴露)、覆蓋面積大，平整度質量評價為極差(failed)；(b2)塌陷，特征為條帶狀分布且深度較小(石料沒有暴露)、覆蓋面積小，平整度質量評價為差(poor)或極差(failed)，具體取決于沉陷深度以及覆蓋面積；(c)裂縫(包括細裂縫和粗裂縫)，特征為狹長狀分布、裂縫間隙在10mm以內，平整度質量評價為中等(fair)或差(poor)，具體取決于裂縫寬度和面積；(d)麻面，特征為路面非常粗糙或者有許多細微的凹坑，平整度質量評價為中等(fair)或差(poor)，具體取決于表面破壞的程度和面積；(e)沙石覆蓋，特征為路面出現塊狀破損且被沙石覆蓋，平整度質量評價為差(poor)；(f)平整度一般的舊路面，特征為路面紋理較粗糙但整體表面流暢，平整度質量評價為中等(fair)；(g)平整度較好的舊路面，特征為路面色澤較為光亮、紋理致密、表面非常流暢，平整度質量評價為好(good)。
圖8為本發明實施例中實驗區內獲取并經過格式轉換lidar點云數據。
圖9為本發明實施例經影像拼接與配準處理的實驗區多光譜正射影像。
圖10為本發明實施例中融合影像波段信息的lidar點云數據(圖中顯示的是藍波段)。
圖11為對本發明實施例中使用5個光譜信息和36個多尺度幾何特征信息進行特征選擇后所保留的23個關鍵特征變量及訓練樣本在這些特征變量下的屬性值分布特點；
其中，(a)不同地物的反射強度、r-g-b波段灰度標準差、r-g-b波段灰度值等5個光譜特征；(b)不同地物的局部粗糙度(ldr)與局部維度特征(ldf)；(c)不同地物的局部維度特征(ldf)和局部高程差(lhd)；(d)不同地物的局部高程差(lhd)。括弧里的長度數字代表空間分析尺度大小，ldf_1表示一維的局部維度特征，ldf_1表示二維的局部維度特征
圖12為本發明實施例中對路面點云插值生成tin表達的dsm。
圖13為本發明實施例中利用自然鄰域插值算法得到的道路dsm。
圖14為本發明實施例利用iri指數模型計算得到的路面平整度空間分布圖。
圖15為本發明實施例對平整度值進行重分類后得到的路面平整度質量空間分布圖。
圖16為本發明實施例中不同路面質量類型的空間分布統計圖。
圖17為評估本發明實施案例所建立的平整度質量評價結果的混淆矩陣示意圖。
具體實施方式
下面結合附圖，并通過具體實例進一步描述本發明的具體實施過程，但不以任何方式限制本發明的適用范圍。
圖1為本發明提出的融合無人機lidar和高分影像的路面平整度監測方法的流程框圖，包括以下步驟：
步驟1：利用無人機lidar系統同時獲取待測路段lidar點云數據和高分辨率影像數據；
步驟2：由于不同的lidar系統生產廠商提供的點云數據格式不同(生產廠商大都有一套自定義的存儲格式)，為了能在一些通用軟件上進行處理，因此需要對lidar點云數據進行格式轉換，在輸出時注意投影坐標系的選擇以及保留的字段屬性；
步驟3：對步驟2中得到的通用點云文件(*.las)進行去噪處理，初步去除干擾噪聲信息，從而得到去燥后的點云數據。具體的去噪處理方式包括：
a)噪聲點過濾。主要是對lidar點云中由于激光脈沖的折射或者多路徑效應產生的異常點進行剔除；
b)高程異常點過濾。主要是剔除局部高程異常點(粗差)，具體的判斷方法為：設置一定的搜索半徑閾值r，如果某點與該點所在半徑r內的鄰域點的高程平均值相差超過3倍的標準差，則可視為高程異常點并加以剔除；
c)掃描角過濾。由于在無人機航跡規劃時，為了能全面獲取路段的表面信息，無人機的飛行路線應該在路段正上方，因此利用點云的掃描角屬性可以對數據進行過濾從而減少數據處理量，節省不必要的數據處理時間。假設某一點到飛機正下方的水平距離為d，飛機的飛行高度為h，則該點的掃描角θ定義為式1：
θ＝arctan(d/h)(式1)
根據道路的邊界范圍可以確定道路的最大掃描角θmax，因此設置合適的掃描角閾值就可將道路區域分離。考慮到航線有一定的偏差，實驗中需要盡可能將掃描角閾值設置稍大一些；
步驟4：利用無人機影像拼接軟件(如本實例中使用的由瑞士pix4d公司開發的商用軟件pix4dmapper)，對待測路段內無人機拍攝的高分辨率影像數據進行影像拼接和配準處理，可得到待測路段的正射遙感影像；
步驟5：對高分辨率影像進行光譜特征提取，本發明用到的光譜特征參數包括以下4個：紅、綠、藍波段的灰度值以及這三個波段的標準差，標準差的數學模型為式2：
其中，std.代表標準差，gi(i＝1,2,3)是紅、綠、藍三波段的灰度值。
步驟6：對包含多波段高分影像與點云數據進行特征融合，以點云數據的地理空間參考為基準，通過仿射變換將影像的像素坐標映射到該像素中心點的地理坐標，從而將影像上對應的光譜信息匹配到點云數據。仿射變換的數學模型為式3：
但在實際操作中，由于影像的像素數量遠遠大于點云的數據量，而且兩者也不完全重合，因此為了減少計算量，本發明采用的是從地理坐標轉換到圖像坐標的過程，即仿射反變換，表示為式4：
式3、式4中，g即是仿射變換矩陣，ai,bi,ci(i＝1,2)是仿射變換系數，mapx、mapy分別是點云的x坐標、y坐標，row、col分別是影像的行、列值。geotiff或其他帶有地理坐標信息的影像文件的元數據中一般就有仿射系數信息，另外也可通過影像上的4個邊界點坐標及其對應的行列數(像素坐標)信息通過構建方程來對仿射系數進行求解而得到。
步驟7：對點云數據構建八叉樹進行存儲，以八叉樹格網單元的長度為最小空間間隔對點云進行抽稀，得到抽稀的點云數據。然后對抽稀的點云數據進行幾何特征參數的計算，本發明中用到的幾何特征參數包括以下三個：
a)局部粗糙度(localdegreeofroughness,ldr)，如圖2所示，假設整個三維點云點集為c＝{pi|pi＝(xi,yi,zi),i＝1,2,3,...,n}，當前計算點為p＝(x,y,z)∈c，三維空間中其半徑為r的鄰域點集為p＝{pj|||pj-p||≤r,pj≠p,j＝1,2,3,...,n}，設平面t:ax+by+cz+d＝0是點集q＝{p∪p}的最優擬合平面，那么ldr的數學形式可以表述為式5：
式5中的a,b,c,d分別代表平面方程的基本參數，x,y,z分別代表當前搜索點的x,y,z坐標。本發明中的最優擬合平面t是指點集中所有點投射到x-y平面上時，使得z值誤差平方和s為最小的平面，具體的數學形式為式6：
式6中xi,yi,zi分別代表第i個點的x,y,z坐標，其他參數含義同式5；
b)局部維度特征(localdimensionfeature,ldf)，如圖3所示，三維點云點集和鄰域點集的構造類似于局部粗糙度，假設當前搜索點為p＝(x,y,z)∈c，令點集q＝{p∪p}＝[xyz]，首先對q進行主成分變換，得到矩陣q的三個主成分系數μ1,μ2,μ3(μ1≥μ2≥μ3)，通過式7進一步對這三個主成分系數進行歸一化：
其中，λ1、λ2、λ3分別對應于當前搜索點在該鄰域內服從一維、二維、三維空間分布的程度，也即一維、二維、三維特征值，由于λ1+λ2+λ3＝1，因此只需要前兩個特征參數即可表達當前搜索點的局部維度特征，這樣可以提高點云數據的幾何特征計算效率；
c)局部高程差(localheightdifference,lhd)，如圖4所示，其中最優擬合平面的定義與局部粗糙度計算中的最優擬合平面定義一致，假設當前搜索點為p＝(x,y,z)∈c，三維空間中其半徑為r的鄰域點集為p＝{pj|||pj-p||≤r,pj≠p,j＝1,2,3,...,n}，點集q＝{p∪p}是包含當前搜索點的鄰域點集，那么lhd的數學形式可以表述為式8：
式9中
實驗過程中為了加快計算搜索點的局部高程差的速度，對q進行svd分解可以快速得到平面t的4個關鍵參數a,b,c,d，然后尋找距離平面t最大的點和距離平面t最小的點進行距離求差，從而得到搜索點的局部高程差。
通過不斷改變搜索鄰域的半徑r的大小，并依次計算各個幾何特征(ldr、ldf和lhd)在不同搜索空間半徑r下的對應值，可以得到關于這三個幾何特征的一組列向量，這組列向量就是對地物目標(搜索點所在的物體)的多尺度幾何特征參數。
步驟8：對高分辨率影像與點云數據進行融合后，此時的點云數據包含了大量的特征變量(幾何特征參數和光譜特征參數)，形成一個高維的矩陣(每個點本身就是一個3維向量)，本身這些特征變量之間也可能存在相關性，導致信息會產生一定的冗余，去除冗余信息的過程就是特征選擇的過程。本發明使用的特征選擇策略分為兩個過程：特征子集搜索和特征子集評價。特征子集搜索的策略為前向搜索，而特征子集的評價則主要是借助相關系數來評估子集內部的相關性程度。具體的實現過程如下：
1)首先建立兩個集合形式的數據結構，假設分別為p和q，p代表候選子集，q代表最優子集，類別屬性變量為y，除類別屬性外的其他屬性為x＝{x1,x2,x3,...,xn}。將除y外的其他特征變量添加到候選子集中，設置子集內部的相關性閾值為δ；
2)如果
式10中，xi、y分別代表第i個特征、分類屬性，
3)將待選特征變量xc與q內的每個特征變量進行相關性評價，如果xc與q內的任意一個特征變量的相關性均小于δ，則將xc添加到q內，同時從p內刪除該特征變量，否則直接從p內刪除該特征變量，繼續步驟2)；
4)輸出q，結束特征子集搜索。
特別說明的是，如果特征變量個數不是很多(比如30個以內)該步驟也可省略。
步驟9：獲得抽稀后的點云數據后，交互選擇其中的一部分作為訓練樣本，訓練樣本要求覆蓋所有的地物類型，并且訓練樣本不能過于集中，適當分散在點云所在空間內，每種地物類型的訓練樣本數量不宜太少，否則分類器無法構建分類模型。隨機森林模型具有對數據集的適應能力強(能同時處理離散型、數值型數據)、訓練速度快、抗噪聲能力強、不容易陷入過擬合、模型參數少、分類效果好等優點而成為機器學習領域內最受歡迎的算法之一，本發明借鑒其優點來對同時融合多尺度幾何特征和光譜特征的點云數據進行地物類別分類。
該算法的基本思想為：假設輸入樣本數據的個數為n(每個樣本對應一行)，特征變量的個數為m(每個特征變量對應一列，即原始樣本數據大小為n*m的矩陣)，初始時隨機有放回的從n個樣本中選擇n個樣本作為一棵樹(隨機森林中每一棵樹都是二叉樹)的訓練集，該樹的根節點就是訓練集；隨機地從m個特征中選擇m個特征變量，然后從這m個特征變量中按照純度最小的原則選擇特征變量進行分裂，從而可將該樹分為兩個節點，這兩個節點各自包含訓練子集的一部分，如此遞歸的進行劃分，直至該樹無法繼續分裂(到達二叉樹的最大遞歸深度)或者這m個特征使用完畢為止。重復剛才構建樹的過程可得到含有任意棵樹的森林(隨機森林)，并且每棵樹的訓練樣本都不完全一樣，使用的特征變量也不完全一樣。在預測類別階段，當輸入一個待分類樣本時，森林中的每棵樹對該樣本都會得到一個決策類別，對決策類別進行統計便可得到待分類樣本在各個預測類別下的分布頻率，輸出其中頻率最大的類別作為其預測類別，從而完成分類。
本發明中使用的純度度量指標為基尼系數，采用文獻[5](陳華舟,陳福,石凱,等.基于隨機森林的魚粉蛋白近紅外定量分析[j].農業機械學報,2015(05):233-238)記載的計算方法，其數學模型為：
1)假設訓練樣本集包含p個類別，第i個類別的樣本數比例為pi，則其基尼系數gini(p)值為：
對該訓練樣本集按照某個特征變量劃分為m個子集，第i個子集的樣本個數為ni，該訓練樣本集的樣本總數為n，則劃分后子集的基尼系數ginisplit(p)為：
步驟9：對抽稀后的點云進行分類，依據類別屬性對其進行過濾可獲得路面點云數據，進一步通過相關的點云內插算法就可以得到待測路段的dsm數據。為減少高程噪聲對表面模型效果及后續平整度計算的影響，本方法先利用不規則三角網(tin)生成三維地形數據，然后采用自然鄰域插值法構建道路的dsm模型。文獻[6](sibsonr.abriefdescriptionofnaturalneighbourinterpolation[j].interpretingmultivariatedata,1981:21-36)記載的自然鄰域插值法通過式13計算得到待插值點的高程值：
式13中的wi代表鄰域內第i個高程點的權重系數，zi則是鄰域內第i個高程點的高程，i代表鄰域內第i個高程點，z為待插值點的高程值。該方法確定鄰近點的原則和距離無關，而是與待插值點的泰森(voronoi)多邊形是否有交集來確定鄰近關系，而且鄰近點的權重則是待插值點的泰森多邊形與鄰近點的泰森多邊形的公共面積所占比例(如圖5)。需要注意的是：待插值點的泰森多邊形構建過程中考慮的是所有點，而鄰近點的泰森多邊形構建過程中則將待插值點排除在外。該插值方法的基本特點是它具有局部性，僅使用待插值點周圍的樣本點子集，權重系數保證了插值得到的結果在所使用的樣本值范圍之內。其優點在于不會出現外推趨勢且不會生成輸入樣本點尚未表示的山峰、凹地、山脊或山谷，而且得到的表面在除輸入樣本位置之外的其他所有位置均是平滑的；
步驟10：利用得到的道路dsm選取一系列一定間距(如0.2m)的縱剖面線，然后利用國際平整度指數iri模型來表達實驗路段的平整度值。本發明具體實施中，采用文獻[7](sayersmw,gillespietd,queirozca.theinternationalroadroughnessexperiment[j].1986)記載的方法計算得到路段平整度值。iri指數的核心原理為四分之一車模型(如圖6)，當四分之一車輛以一定的速度沿道路縱剖面行駛時，在路面坡度的激勵作用下，系統將產生振動，計算其以規定速度(如80km/h)行駛一定距離后(如1km)懸掛系統的累積位移值即為iri，單位為m/km。為求解該懸掛系統的相對位移，建立二階振動微分方程如式14～式15：
其中，y是路面高程，zs和zu分別代表簧載質量位移和非簧載質量位移，μ、c、k1、k2是系數，且依據文獻[7]記載，有：
其中，ms表示簧載質量的大小，mu表示非簧載質量的大小，ks表示連接簧載質量和非簧載質量的彈簧的勁度系數，cs表示連接簧載質量和非簧載質量的線性阻尼系數，ks表示連接輪胎的彈簧的精度系數。
構造狀態變量
其中：
式16是一個非齊次線性微分方程，其對應的狀態方程為：
其中，
pr＝a-1(st-i)b,(式19)
式18、式19中的i是單位矩陣，結合
文獻[8](王建鋒,宋宏勛,馬榮貴.路面平整度評價指標iri的影響因素[j].重慶交通大學學報(自然科學版),2012,31(06):1145-1148)記載了上述計算iri的方法，但使用的數據為傳統地面測量數據。本實例中為實現iri指數的計算過程，首先需要將dsm數據進行等距離間隔采樣處理得到路面高程點序列y＝{y0,y1,...,yn}，并將其作為模型的輸入數據。然后依據矩陣a和b計算狀態矩陣st和pr，然后對y進行一階差分運算得到
步驟11：步驟6的計算結果僅僅是得到了單個剖面線的iri值序列，為了能模擬整個路段表面的平整度值，可對道路dsm上的每個柵格單元按照最鄰近原則分配對應的iri值，這種處理方式的誤差大小主要取決于縱剖面線的間距大小，間距越小則誤差也越小。經過上述過程便可得到整個待測路段的平整度空間分布。
步驟12：參照如下標準(如表2)對得到的路面iri值進行分類，即完成對路面平整度質量的評價。該分級設置標準需要結合我國《公路瀝青路面養護技術規范》的路面養護質量的相關要求來設置，本方法是結合《公路瀝青路面養護技術規范》(jtj073.2-2001)規范來確定的。通過統計待測路段各個質量等級的像元數目可計算其對應的分布比例，并結合道路養護維修的相關規范，從而對該路段需要何種維護措施做出科學決策。
表2路面平整度質量分級設置標準
為了對本方法的可靠性進行評價，采用下述步驟對該方法路面平整度質量監測結果進行精度評價。
1)根據設定的平整度質量分級標準，在待測路段選擇一系列(包含平整度質量較好與較差，包含健康與病害類型，每種路面選擇一定數量的地面實測樣本數據)的路面樣本作為參考數據；
2)對參考數據選取外觀、表面粗糙度、車轍深度、病害面積、損壞程度這5項指標使用專家評分法對其進行平整度質量評價，并作為參考值；
3)運用本方法同樣的對此路面樣本進行平整度值進行計算，依據平整度質量分級標準，得到本方法的對應評價結果；
4)對平整度質量結果的精度評價可以類比于遙感影像分類結果的精度評價，而混淆矩陣是用來評價遙感影像分類精度的最常用形式。因此本方法也采用混淆矩陣對路面平整度質量結果進行評價。借助分類精度評價指標總體精度(oa)和kappa系數對結果進行分析，得到本方法的精度，具體的精度評價指標計算方法為：
式21、式22中，n是樣本總數，nii代表本方法的評價結果與參考評價結果均為i等級的樣本數量，ni.則表示參考評價結果為i等級的樣本總數，n.i則代表本方法的評價結果為i等級的樣本總數。
以下以某市一條縣級公路為例說明本方法具體實施及精度評價過程。
(1)數據獲取
該縣級公路所在地理位置為北緯44°24′47″、東經85°53′47″附近，道路寬度約為8m。研究人員于2016年6月23日借助瑞士生產的scoutb1-100無人直升機搭載的激光掃描儀系統(rigelvux-1lr)在實驗區正上方進行實驗數據采集，該激光掃描儀的基本參數信息如表3所示(來源：http://www.riegl.com/products/newriegl-vux-1-series/newriegl-vux-1lr/)。為了使得獲取的點云數據密度足夠高，實驗中無人機飛行的距離路面高度設置為30m、飛行速度為5m/s、掃描角為110°、掃描頻率為550khz，最后得到的激光腳點密度為300-600pts，掃描幅寬為85.7m，數據獲取當天晴朗無(少)云。
表3激光掃描儀系統的關鍵參數信息
在獲取lidar點云數據的同一天，研究人員在飛行區域內進行地物特征點與路面樣本數據的采集。地物特征點數據主要是沿道路周邊采集，通過選取一些關鍵的地物目標(如車身、標志線、坑槽、大塊裂縫)進行幾何尺寸度量，利用實時差分gps設備進行定位并進行筆錄存檔；然后在飛行區域區道路表面上選取不同類型的路面樣本(包含平整度質量較好與較差，覆蓋健康與病害類型)進行取樣拍照并作特征描述與分類統計，所有的數據采集工作都進行位置記錄與拍照保存。本研究總共采集了10個關鍵地物特征點，9種不同類型的路面樣本(樣本點總數52，每種路面樣本大約為5-7個)。路面樣本類型包括：坑槽(分大塊狀和小塊狀)、塌陷(分石料完全暴露和石料未暴露)、裂縫(大裂縫和小裂縫)、麻面、沙石覆蓋、路面平整度一般及平整度良好，如圖7所示。
(2)數據處理
數據的具體處理步驟如圖1所示，具體的步驟執行過程如下。
第一步：由rigelvux-1lr激光掃描儀所獲取的lidar點云數據包含的屬性有：三維坐標、掃描角、回波數量、回波次數、激光強度。原始點云文件并非通用的點云文件格式，借助rigel公司提供的配套軟件對點云數據進行格式轉換得到通用的.las格式，如圖8所示；
第二步：對上述lidar點云數據進行去噪處理，實驗過程中先通過目視識別初步將噪聲點進行剔除，然后設置搜索半徑為0.2m對局部高程異常點進行剔除，接著根據道路寬度和無人機飛行高度計算得到的路面點云的掃描角范圍在±8°內，考慮到航線有一定的偏差，實驗過程中設置掃描角閾值為±24°對點云數據進行過濾；
第三步：利用瑞士pix4d公司開發的商用無人機數據處理軟件pix4dmapper對待測路段內無人機拍攝的高分辨率影像數據進行影像拼接和配準處理后可得到待測路段的正射影像圖(圖9)。然后對該正射影像進行特征提取，計算紅綠藍三波段的灰度值標準差作為新的特征波段，并與紅綠藍三波段組合，得到含有光譜信息的4個特征波段；
第四步：采用仿射變換模型將點云中的點與影像中的像素進行匹配，然后將影像上的波段值融合到點云數據中，結果如圖10所示(以藍波段信息為例)；
第五步：對點云數據與高分辨率影像數據進行特征融合后，需要對該融合數據進行多尺度幾何特征參數計算。由于點云數據量非常大，為節省時間，同時又不影響結果的可靠性，本例先對點云數據的存儲結構進行優化，即構建八叉樹存儲結構，設置一個比較適宜的數值(本例中使用的是0.3m)作為空間存儲單元長度進行重采樣，得到抽稀后的點云數據；
第六步：以上述點云數據的每個點為中心，設置搜索半徑在[0.2m,1.0m]內且以0.1m為間隔單位進行搜索，搜索的鄰域點則是去噪后的點云數據中與當前搜索點距離小于等于搜索半徑的點集，然后分別計算當前尺度下的局部粗糙度(ldr)、局部維度特征(ldf)、局部高程差(lhd)。隨著空間搜索尺度的調整，即可得到不同尺度幾何特征信息的點云數據；
第七步：對步驟六得到的點云數據通過目視選擇出覆蓋所有地物類型的部分點云作為訓練樣本，然后通過特征選擇對特征變量進行降維(如圖11)，容易看出特征選擇后多尺度幾何特征參數lhd相對于其他特征變量而言對各類地物的區分能力最強。進一步借助隨機森林模型進行分類器構建。本例在構建過程中設置樹的個數為200，樹的深度為40。其次利用所構建的隨機森林模型對將抽稀后的點云數據進行分類，從而完成抽稀點云的分類過程；
第八步：由于到目前為止只完成了對抽稀的點云進行了分類，需要采取一定的方法對去噪后的全部點云進行分類。由于空間距離越近的兩個物體其屬性越相似，因此本例按照空間最近鄰原則給去噪后的點云數據的類別屬性賦以距離該點最近的已分類點的類別屬性值；
第九步：對分類后的點云數據進行類別篩選，得到道路點云數據，再通過不規則三角網構建道路表面tin模型(圖12)，進一步利用自然鄰域插值得到道路表面的dsm模型，為了滿足iri指數計算需要同時兼顧數據處理需要，柵格分辨率設為50mm，結果如圖13；
第十步：在arcgis10中通過交互矢量化方式，在道路內部沿道路方向數字化生成間隔為0.125m的一系列縱剖面線，然后利用iri指數模型計算這些剖面線上的iri值。為了反映各路段縱向上的差別，在計算iri時需設置合適的計算單元長度，本實例設置為10m。為了可視化效果再將計算的iri值進行柵格化得到整個路段的平整度空間分布結果，如圖14；
第十一步：依據本方法設定的平整度質量分級設置標準，對得到的路面平整度空間分布圖進行重分類，得到路面的平整度質量分布圖，如圖15所示，通過對不同質量等級的柵格單元進行統計分析，可確定該路段的各質量等級分布比例統計圖，如圖16所示；
(3)結果分析
總體上來看，縱向方面的特點為道路南段的質量評價結果主要以“fair”為主，少部分區域是“failed”；道路中部的質量評價最差，評價結果主要為“poor”或“failed”；而道路北段的質量評價結果則以“fair”或“failed”居多，這與平整度的空間分布趨勢一致。橫向方面的特點則是中心線附近的質量等級要比兩側車行道低，以“failed”為主，而行車道附近的評價結果則主要以“fair”或“poor”為主，道路邊界地帶的路面質量等級則相對較低，為“poor”或“failed”，這部分區域由于路面沉陷、砂土混合等原因導致表面粗糙不平。結合各質量等級分布圖來看(圖15)，可知該路段整體的平整度質量以“fair”為主，其分布比例達41.1％，質量為“poor”或“failed”的路段分布比例居其次，各為27.4％和24.9％，而質量為“good”的路段所占比例最小，僅6.6％。而且質量為“poor”或者“failed”的分布比例總和占整個路段的一半以上，這就表明該路段整體的質量偏差，需要做出重大維修。
(4)精度評價
以實地采集的52個路面樣本為分析數據，通過選取外觀、表面粗糙度、車轍深度、病害面積、損壞程度這5項指標使用專家評分法對其進行平整度質量評價，并將其作為參考值。然后對這些路面樣本依據實地測量的gps點坐標并輔助以高分辨率影像數據在dsm模型上進行地圖定位，接著根據路面平整度空間質量空間分布圖對這些路面樣本的平整度質量等級進行識別，再對所有路面樣本的識別結果進行統計，并由此得到用于進行精度評定的混淆矩陣(如圖17)。最后利用式21和式22分別計算得到了本方法的總體分類精度為75％，kappa系數為0.65，表明該方法的整體評價結果與實際結果大體相同，且一致性程度較高。
需要注意的是，公布實施例的目的在于幫助進一步理解本發明，但是本領域的技術人員可理解：在不脫離本發明及所附權利要求的精神和范圍內，各種替換和修改都是可能的。因此，本發明不應局限于實施例所公開的內容，本發明要求保護的范圍以權利要求書界定的范圍為準。
- 還沒有人留言評論。精彩留言會獲得點贊！