本發(fā)明涉及三維點云獲取技術領域，具體涉及一種聯合激光掃描和影像匹配的三維點云生成方法。

背景技術：

21世紀是一個信息化的時代，突出表現在將人類世界以信息化的方式進行顯示、統(tǒng)計、分析和應用。美國前副總統(tǒng)戈爾在1998年的演講中，提出了“數字地球”的概念，掀起了一股數字地球研究的浪潮。2005年，谷歌公司推出了GoogleEarth產品，將航片、衛(wèi)片數據映射到虛擬的數字地球上，便于用戶瀏覽和使用世界各地的地理信息產品。在全世界正從數字地球向智慧地球發(fā)展的過程中，我國也在推動“數字中國”、“智慧城市”等的快速發(fā)展，在為政府部門、企事業(yè)單位和社會公眾等提供公共管理、突發(fā)事件應急、科學決策等服務方面具有重要意義。

數字表面模型(Digital Surface Model,DSM)是以數字化的形式，表達地球自然表面和人工地物的三維模型，是“智慧城市”重要的顯示和分析工具。近十年來，隨著“數字中國”、“智慧城市”的快速發(fā)展，對DSM生成的精度、速度和分辨率提出了越來越高的要求。激光掃描和立體影像密集匹配技術是目前兩種主流的地形地物三維信息獲取手段，通過將相機、LiDAR等傳感器搭載于車輛、無人機、航空飛機乃至衛(wèi)星等不同高度的遙感平臺，間接或直接地快速獲取測區(qū)地表的三維測繪信息。

LiDAR(Light Detection And Ranging，激光探測與測量點云)系統(tǒng)由全球定位系統(tǒng)(GPS)、慣性導航系統(tǒng)(IMU)以及激光器(Laser Scanner)共同組成，以主動方式快速獲取地面模型及空間點云信息，具有速度快、時效性強、作業(yè)范圍大等優(yōu)點。但是LiDAR系統(tǒng)成本昂貴，且點云密度較低(與對應的影像分辨率比較而言)，無法很好地描述線狀地物；并且由于自然地表和建筑物的材質問題，在某些區(qū)域無法形成反射點，無法很好地為智慧城市等三維重建應用服務。

影像密集匹配通過覆蓋一個測區(qū)的多張立體影像，根據同名光線對對相交的原理，從影像的二維信息恢復整個測區(qū)的三維空間信息。與LiDAR激光點云相比，立體影像具有密集匹配點云密度大、平面精度高、粗差剔除技術完善、影像上的地物幾何特征更明顯、影像數據獲取成本低等突出優(yōu)勢。但是影像匹配的結果取決于地物紋理情況以及影像輻射質量。影像密集匹配技術在紋理貧乏區(qū)域、重復紋理區(qū)域表現較弱。較差的影像輻射質量以及陰影區(qū)域也會嚴重影響匹配結果。影像密集匹配往往速度較慢，很難達到實時的要求，嚴重制約數字城市的重建速度。

由于上述兩種技術各自的不足，單純地采用激光掃描，或者單純地采用影像密集匹配，均無法很好地描述測區(qū)地表的三維信息。為了滿足高精度、高分辨率、高實時性的三維重建的需求，需要將兩種技術相結合，充分發(fā)揮各自的優(yōu)勢，得到完整可靠的點云產品，能夠滿足快速、大范圍、高精度、高分辨率三維重建的需求，為數字地球、智慧城市的發(fā)展提供技術支撐，為社會的可持續(xù)發(fā)展服務。

技術實現要素：

本發(fā)明的目的在于提供一種聯合激光掃描和影像匹配的三維點云生成方法，該方法充分發(fā)揮激光掃描技術和影像密集匹配技術各自的優(yōu)勢，利用激光掃描技術直接快速獲取測區(qū)的三維點云，對于線狀地物、無反射點區(qū)域(如水面)等掃描困難區(qū)域，采用影像密集匹配技術獲取掃描困難區(qū)域的三維點云，最終快速生成稠密的高精度的三維點云，為智慧城市等應用服務。

為解決上述技術問題，本發(fā)明公開的一種聯合激光掃描和影像匹配的三維點云生成方法，其特征在于，它包括如下步驟：

步驟1：根據原始激光探測與測量點云數據，生成數字表面模型，提取數字表面模型中的無效區(qū)域，將該無效區(qū)域中最小外接矩形的四個角點，分別投影到三維重建測區(qū)的影像上，生成一個三維重建測區(qū)影像上的四邊形區(qū)域，該四邊形區(qū)域即可認為是激光探測與測量點云數據中的無效區(qū)域在三維重建測區(qū)影像上對應的范圍；

步驟2：將激光探測與測量點云數據中的無效區(qū)域在三維重建測區(qū)影像上對應的范圍作為三維重建測區(qū)影像的密集匹配區(qū)域，根據該三維重建測區(qū)影像密集匹配區(qū)域構建影像金字塔，從金字塔頂層開始，采用半全局密集匹配方法(Semi-global Matching,SGM)進行由粗到精的分級匹配，獲得密集匹配區(qū)域的二維視差圖，并采用前方交會的方式，生成三維重建測區(qū)影像的密集匹配區(qū)域所對應的三維點云數據；

步驟3：采用Canny邊緣檢測算子，在三維重建測區(qū)影像上提取線特征，將每個線特征作為一個密集匹配區(qū)域，分別進行影像密集匹配，由于線特征只有一維，因此匹配方式采用沿著線特征方向的一維動態(tài)規(guī)劃的優(yōu)化方式，獲取密集匹配區(qū)域的視差圖，并采用前方交會的方式，生成線特征對應的三維點云數據；

步驟4：在原始激光探測與測量點云數據的基礎上，采用點云融合技術將步驟2生成的三維重建測區(qū)影像密集匹配區(qū)域所對應的三維點云數據，以及步驟3生成的線特征對應的三維點云數據進行三維點云融合，生成融合激光探測與測量點云和密集匹配點云的三維重建點云。

本發(fā)明的有益效果：

本發(fā)明公開的一種聯合激光掃描和影像匹配的三維點云生成方法，能夠充分利用影像分辨率高、灰度特征明顯的優(yōu)勢，解決LiDAR點云中無反射點區(qū)域、線狀地物描述等重建問題，能夠極大地減少影像密集匹配的計算復雜度，快速獲取更加完整精確的點云產品，能夠為智慧城市、數字城市等應用服務。

附圖說明

圖1為本發(fā)明的流程圖；

圖2為激光探測與測量點云中的無效區(qū)域；

圖3為半全局密集匹配的掃描線方向；

圖4為無效區(qū)域的密集匹配點云；

圖5為激光探測與測量點云對線狀地物的定位精度；

圖6為線特征匹配的代價積聚方向。

具體實施方式

以下結合附圖和具體實施例對本發(fā)明作進一步的詳細說明：

本發(fā)明針對現有的主流三維點云獲取技術：激光掃描和影像匹配各自的特點，充分發(fā)揮兩者各自的優(yōu)勢，提出一種聯合激光掃描和影像匹配的三維點云生成方法，能夠解決激光掃描在弱/無反射區(qū)域、線狀地物區(qū)域的三維點云獲取困難問題，解決影像密集匹配時間復雜度高的問題，生成的三維點云具有完整度高、精度高、點云密度大、獲取速度快等優(yōu)勢，如圖1所示，它包括如下步驟：

步驟1：激光探測與測量點云完整度評估；

根據原始激光探測與測量點云數據，生成數字表面模型(Digital Surface Model，DSM)，提取數字表面模型中的無效區(qū)域，將該無效區(qū)域中最小外接矩形的四個角點，分別投影到三維重建測區(qū)的影像上，生成一個三維重建測區(qū)影像上的四邊形區(qū)域，該四邊形區(qū)域即可認為是激光探測與測量點云數據中的無效區(qū)域在三維重建測區(qū)影像上對應的范圍；

步驟2：弱/無反射區(qū)域三維點云生成；

將激光探測與測量點云數據中的無效區(qū)域在三維重建測區(qū)影像上對應的范圍作為三維重建測區(qū)影像的密集匹配區(qū)域，根據該三維重建測區(qū)影像的密集匹配區(qū)域構建影像金字塔，從金字塔頂層開始，采用半全局密集匹配方法進行由粗到精的分級匹配，獲得密集匹配區(qū)域的二維視差圖，并采用前方交會的方式，生成三維重建測區(qū)影像密集匹配區(qū)域對應的三維點云數據；

步驟3：線特征地物的三維點云生成；

采用Canny邊緣檢測算子，在三維重建測區(qū)影像上提取線特征，將每個線特征作為一個密集匹配區(qū)域，分別進行影像密集匹配，由于線特征只有一維，因此匹配方式采用沿著線特征方向的一維動態(tài)規(guī)劃的優(yōu)化方式，獲取密集匹配區(qū)域的視差圖，并采用前方交會的方式，生成線特征對應的三維點云數據；

步驟4：激光探測與測量點云和密集匹配點云融合；

在原始激光探測與測量點云數據的基礎上，采用點云融合技術將步驟2生成的三維重建測區(qū)影像密集匹配區(qū)域對應的三維點云數據，以及步驟3生成的線特征對應的三維點云數據進行三維點云融合，生成融合激光探測與測量點云和密集匹配點云的三維重建點云。

上述技術方案的步驟1中，本發(fā)明的輸入是激光探測與測量點云和測區(qū)的光學影像集合。為了保證影像密集匹配順利實施，需要具有重疊度的影像至少兩張，即能夠組成至少一個立體像對。激光探測與測量是主流的獲取地形地貌三維信息的手段之一，得到了廣泛的應用。但是在水面等無法反射激光的區(qū)域，無法獲得激光探測與測量點云，導致激光探測與測量點云中存在無效區(qū)域，如圖2(a)中的河流區(qū)域所示。本文定義僅僅有極少量的，甚至不存在任何三維點的區(qū)域為無效區(qū)域。因此，需要對激光探測與測量點云進行完整度評估，找出無效區(qū)域，在后續(xù)處理中，采用影像密集匹配技術，對這些無效區(qū)域進行處理，生成三維點云。

上述技術方案的步驟1中，為了方便統(tǒng)計激光探測與測量點云中的無效區(qū)域，首先根據激光探測與測量點云，生成數字表面模型(DSM)，定義激光探測與測量點云的點間距為s，一般點間距s由工程文件提供。計算無效區(qū)域中最小外接矩形的四個角點的X、Y坐標的方式如下：

X_lb＝min{Xl_i|i＝1…t} Y_lb＝min{Yl_i|i＝1…t}

X_rt＝max{Xl_i|i＝1…t} Y_rt＝max{Yl_i|i＝1…t}

其中，t表示原始激光探測與測量點云數據中三維點的數目；(Xl_i,Yl_i)表示第i個激光探測與測量點的X、Y坐標；(X_lb,Y_lb)表示無效區(qū)域中最小外接矩形左下角的角點坐標；(X_rt,Y_rt)表示無效區(qū)域中最小外接矩形右上角的角點坐標，(X_lb,Y_rt)表示無效區(qū)域中最小外接矩形左上角的角點坐標，(X_rt,Y_lb)表示無效區(qū)域中最小外接矩形右下角的角點坐標；

將激光探測與測量點云數據的最小外接矩形范圍定義為數字表面模型范圍，數字表面模型的起點為無效區(qū)域中最小外接矩形左下角的角點，定義數字表面模型的正方形網格的大小為s_D x s_D，s_D＝σxs，其中，s表示激光探測與測量點云數據中點和點之間的平均間距，σ表示數字表面模型的網格大小與上述激光探測與測量點云數據中點和點之間平均間距之間的比值，σ要取大于1的值，一般取3，目的就是為了保證每個網格內都存在激光探測與測量點，如果不存在，則該網格可能處于激光探測與測量點云無效區(qū)域；

則數字表面模型的寬W_D和高H_D分別為：

W_D＝(int)(X_rt-X_lb)/s_D；H_D＝(int)(Y_rt-Y_lb)/s_D

其中，int表示取整操作，從而將數字表面模型定義為一個H_D x W_D大小的規(guī)則格網，該規(guī)則格網的起點為(X_lb,Y_lb)；如果數字表面模型的網格內不存在任何激光探測與測量點，則該數字表面模型的網格定義為無效區(qū)域；否則，該數字表面模型的網格定義為有效區(qū)域，在有效區(qū)域中，每個數字表面模型網格都存在至少一個激光探測與測量點，每個數字表面模型網格的高程為對應激光探測與測量點的最大高程，如下式所示：

Z(m,n)＝max(Z_i|(int)(X_i-X_lb)/s_D＝n；(int)(Y_i-Y_lb)/s_D＝m)

其中，(int)(X_i-X_lb)/s_D＝n；(int)(Y_i-Y_lb)/s_D＝m表示坐標為(X_i,Y_i)的激光探測與測量點落入數字表面模型第m行第n列的網格中；Z_i表示第i個激光探測與測量點的Z坐標；Z(m,n)表示數字表面模型網格(m,n)的高程，取落入網格中所有激光探測與測量點的最大高程；

在有效區(qū)域，每個DSM網格都存在有效的高程值；但是由于無效區(qū)域不存在LiDAR點，無效區(qū)域的網格同樣也不存在有效高程。因此，可以根據DSM網格內是否存在有效高程，定義該區(qū)域是否為無效區(qū)域，如圖2(b)中河流區(qū)域所示。根據所述無效區(qū)域內每個數字表面模型網格的X、Y坐標，定義無效區(qū)域的最小外接矩形，如下式所示：

N_lb＝min{n_i|i∈Ω} M_lb＝min{m_i|i∈Ω}

N_rt＝max{n_i|i∈Ω} M_rt＝max{m_i|i∈Ω}

其中，N_lb表示數字表面模型網格所有n列中的最小值，M_lb表示數字表面模型網格所有n列中的最大值，N_rt表示數字表面模型網格所有m行中的最小值，M_rt表示數字表面模型網格所有m行中的最大值，N_lb、M_lb、N_rt和M_rt共同定義了無效區(qū)域最小外接矩形的范圍；Ω表示無效區(qū)域中數字表面模型網格的集合；n表示數字表面模型的網格的列數；m表示數字表面模型的網格的行數；具體外接矩形如圖2(b)虛線框所示。

根據無效區(qū)域外接矩形四個角點的數字表面模型行列坐標，獲得四個角點的物方空間三維坐標，如下式所示：

X_c＝X_lb+n x s_D Y_c＝Y_lb+m x s_D Z_c＝Z(m,n)

其中，(m,n)表示無效區(qū)域外接矩形角點的數字表面模型行列坐標；X_c、Y_c、Z_c分別表示無效區(qū)域外接矩形角點的物方三維坐標；

最后，通過共線方程，將無效區(qū)域外接矩形角點反投影到三維重建測區(qū)影像上得到無效區(qū)域在三維重建測區(qū)影像上的范圍，生成無效區(qū)域的影像，即三維重建測區(qū)影像密集匹配區(qū)域，如下式所示：

其中，x_c、y_c表示無效區(qū)域外接矩形角點在三維重建測區(qū)影像上的坐標；x₀、y₀表示相機內像主點(相機攝影中心與像平面的垂線與像平面的交點)坐標參數，f表示相機的焦距參數；Xs、Ys、Zs表示相機的外方位線元素；a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃表示攝影坐標系與大地坐標系之間的旋轉矩陣的九個元素，最終無效區(qū)域在影像上的顯示范圍如圖2(c)所示。上述旋轉矩陣由外方位元素計算得來的旋轉矩陣，一般通過測量平差手段獲得精確值。

上述技術方案的步驟2中，影像劃分為有效區(qū)域和無效區(qū)域。有效區(qū)域已經存在對應的激光探測與測量點云，不需要影像匹配的方法重新生成三維點云；而無效區(qū)域不存在對應的激光探測與測量點云，需要借助影像密集匹配的方法來生成無效區(qū)域的三維點云。無效區(qū)域往往是水面等紋理貧乏區(qū)域，為了減少匹配歧義，提高匹配魯棒性，需要對無效區(qū)域影像構建金字塔，本發(fā)明采用均值濾波的方式構建影像金字塔，將三維重建測區(qū)影像的密集匹配區(qū)域構建影像金字塔的方法為：定義S為影像金字塔重采樣的尺度，上一級影像金字塔的一個像素，會對應下一級影像金字塔S×S的區(qū)域，影像金字塔構建方式如下式所示；

式中，gⁿ表示上一級第n層影像金字塔的像素灰度；表示下一級第n-1層影像金字塔第i個像素的灰度，S²＝S×S；

在構建影像金字塔后，采用由粗到精的分級匹配策略，獲得密集匹配區(qū)域的二維視差圖，在影像金字塔頂層開始，采用半全局的密集匹配策略，每個像素受到八個方向掃描線的代價累積，如圖3所示，在每條掃描線上，代價積聚方式為：

式中，L_r(p,d)表示像素p在當前掃描線對應視差d的累積代價；L_r(p-1,d)表示像素p-1在當前掃描線對應視差d的累積代價；L_r(p-1,d-1)表示像素p-1在當前掃描線對應視差d-1的累積代價；L_r(p-1,d+1)表示像素p-1在當前掃描線對應視差d+1的累積代價；m_kinL_r(p-1,k)表示像素p-1在當前掃描線對應視差k的累積代價中，最小的累積代價，其中k的絕對值大于1；表示像素p-1在當前掃描線對應視差i的累積代價中，最小的累積代價；r表示掃描線的方向；C(p,d)表示當前像素p對應視差d的代價；p–1表示在當前路徑上，像素p的前一個像素；P₁表示視差平滑項懲罰因子；P₂表示視差階躍項懲罰因子；

單條掃描線的累積結果是不穩(wěn)定的，容易產生“條紋”問題。為了使匹配結果更加穩(wěn)定可靠，對無效區(qū)域影像進行8個方向的代價累積，最后對各個方向的累積結果相加，得到最終的代價積聚結果，如下式所示：

式中，S(p,d)表示對各個方向代價積聚結果相加后，得到的總體代價積聚結果，采用Winner Takes All(WTA)策略獲取頂層影像的初始視差圖，將初始視差圖傳遞到下一級影像金字塔，約束下一級影像金字塔的視差搜索范圍，然后再次采用半全局的密集匹配方法，獲取視差圖，并向下一級影像金字塔傳遞，直至計算到影像金字塔底層為止，根據視差圖，可以快速獲得立體影像的同名點坐標，如下式所示：

x_r＝x_l-d y_r＝y(tǒng)_l

式中，(x_l,y_l)、(x_r,y_r)分別表示左右影像的同名點坐標；d表示視差值，根據左右影像的同名點坐標，確定對應的三維點坐標：

B_u＝X_S2-X_S1B_v＝Y_S2-Y_S1B_w＝Z_S2-Z_S1

X＝X_S1+U₁＝X_S2+U₂

Y＝Y_S1+V₁＝Y_S2+V₂

Z＝Z_S1+W₁＝Z_S2+W₂

式中，X、Y、Z表示物方點的三維坐標；f表示相機的焦距參數；R_i(i＝1,2)表示旋轉矩陣；X_s1、Y_s1、Z_s1、X_s2、Y_s2、Z_s2表示外方位線元素，U₁、V₁、W₁表示地面點在左影像像空間輔助坐標系中的坐標；U₂、V₂、W₂表示地面點在右影像像空間輔助坐標系中的坐標；u₁、v₁、w₁表示像點在左影像像空間輔助坐標系中的坐標；u₂、v₂、w₂表示像點在右影像像空間輔助坐標系中的坐標；B_u、B_v、B_w表示相機之間的基線分量；N₁、N₂表示點投影系數；只要無效區(qū)域在影像上是可見的，即可根據影像密集匹配技術恢復無效區(qū)域的三維點云，如圖4所示。從圖4可以看出，河流區(qū)域已經充滿了密集匹配點云。

上述技術方案的步驟3，由于激光探測與測量點云本身密度較低，直接從LiDAR點云提取的建筑物邊緣并不精確，產生邊緣“偏移”問題，如圖5所示。圖5中，上面一條線條表示激光探測與測量點云提取建筑物邊緣的投影結果；而下面一條直線表示影像上建筑物的邊緣。可以看出，激光探測與測量點云中的建筑物邊緣并不精確，存在一定的偏移。為了進一步提高三維點云的精度，有必要借助影像匹配技術，生成建筑物精確的“真實邊緣”；

采用canny算子，在三維重建測區(qū)影像上提取線特征，將每個線特征作為一個密集匹配區(qū)域，分別進行影像密集匹配，具體方法為，沿著線特征的方向，采用一維動態(tài)規(guī)劃方法，進行代價積聚，如下式所示：

式中，L_r(p,d)表示像素p在當前線特征方向上對應視差d的累積代價；r表示線特征的方向；C(p,d)表示當前像素p對應視差d的代價；p–1表示在當前路徑上，像素p的前一個像素；

由于線特征一般存在正反兩個方向，因此上式所示的代價積聚也存在兩個方向，如圖6所示，將兩個方向的代價積聚結果累加，得到最終的代價積聚結果，采用Winner Takes All策略獲取每條線特征的視差值，計算對應的同名點對，并生成線特征對應的三維點云數據。

上述技術方案的步驟4中，通過影像密集匹配技術，可以生成無效區(qū)域和線狀地物區(qū)域的三維點云。將密集匹配點云和LiDAR點云融合，可以生成更加完整、精度更高的三維點云產品。在無效區(qū)域匹配過程中，不可避免地會將一部分有效區(qū)域進行重復匹配。由于有效區(qū)域已經存在LiDAR點云，因此重復匹配的三維點是不需要融合的，在步驟2中的三維重建測區(qū)影像密集匹配區(qū)域對應的三維點云數據，以及步驟3中的線特征對應的三維點云數據進行三維點云融合前，需要判斷哪些匹配點位于激光探測與測量點云有效區(qū)域，哪些匹配點位于激光探測與測量點云無效區(qū)域，如下式所示：

m＝(X_i-X_lb)/s_D n＝(Y_i-Y_lb)/s_D

其中，P(X_i,Y_i,Z_i)表示坐標為(X_i,Y_i,Z_i)的三維點；None表示該點不參與融合；Exist表示該點參與融合；Z_DSM表示該點所在數字表面模型網格的高程；(m,n)表示該點所在網格的行列號，valid表示有效值，即數字表面模型網格內存在激光點；invalid表示無效值，即數字表面模型網格內不存在任何激光點；(X_i,Y_i)表示三維點的平面坐標；(X_lb,Y_lb)表示數字表面模型的起點；S_D表示數字表面模型網格的大小。上式表明，當該匹配點落入的數字表面模型網格存在有效高程，則說明該匹配點是有效區(qū)域的匹配點，不參與融合過程；否則，說明該匹配點是無效區(qū)域的匹配點，參與融合。

本發(fā)明充分利用激光掃描技術和影像密集匹配技術各自的優(yōu)勢，生成完整、高精度、稠密的三維點云，在弱\無反射區(qū)域，線狀地物區(qū)域有著可靠精細的三維描述，能夠為數字城市、智慧城市等應用提供技術支撐。

本說明書未作詳細描述的內容屬于本領域專業(yè)技術人員公知的現有技術。