本發(fā)明涉及一種無人機序列影像批處理三維重建方法，特別是融合低精度GPS/IMU信息的無人機序列影像批處理三維重建方法。

背景技術：

無人機能夠連續(xù)獲取重疊度大的高精度序列影像，但獲取的影像會丟失深度信息。基于圖像的三維重建，是指利用多幅數(shù)碼相機圖像全自動恢復出場景三維結構的方法與技術。近年來三維重建技術在視頻、圖像三維重建處理領域獲得了巨大的成功，將其應用到無人機圖像處理領域，對無人機圖像進行全自動重建相關的應用，可以拓展無人機的應用范圍，提高無人機的應用水平。但目前對于無人機序列影像三維重建的研究尚處于起步階段，主要存在以下問題：(1)相對于地面影像，基于無人機序列影像的三維重建一般是大數(shù)據(jù)量大場景的三維重建；(2)大多直接將計算機視覺中成熟的算法應用于無人機序列影像三維重建中；(3)沒有充分利用精度不高的輔助信息。

現(xiàn)在越來越多的成像系統(tǒng)都帶有全球定位系統(tǒng)和慣性測量裝置，可以獲得包含三維世界地理坐標系(如WGS84)下地理坐標信息和相機姿態(tài)信息的序列圖像。但是，這些系統(tǒng)賴高精度的地理定位設備，通過這些設備進行的標定以及獲得的姿態(tài)和位置數(shù)據(jù)的精度一般都比圖像的方式(例如，亞像素級的圖像配準)要高。另一方面，目前的各種地理定位定向系統(tǒng)一般可以提供連續(xù)的但常常精度不高有時甚至不準確的位置和姿態(tài)信息，就像無人機所搭載的地理定位定姿系統(tǒng)。不幸的是，從這些設備得到的GPS/IMU數(shù)據(jù)達不到直接用于三維目標重建及導航等部分計算機視覺工作所要求的像素級的圖像匹配精度要求。這些設備可以提供概略的相機姿態(tài)和位置信息，充分利用這些精度不高的輔助信息有望提高基于無人機序列影像三維重建的效率和精度。

技術實現(xiàn)要素：

本發(fā)明所要解決的技術問題是，克服現(xiàn)有技術的缺點，提供一種無人機序列影像批處理三維重建方法，本發(fā)明通過融合低精度GPS/IMU信息，利用無人機序列影像全自動地恢復出場景三維結構，得到重構點云的相對坐標。

本發(fā)明提供了一種無人機序列影像批處理三維重建方法，包括以下步驟：

步驟一、融合低精度GPS/INS信息的影像匹配：

借助于GPS/IMU系統(tǒng)提供的全局的位置信息，這些容易混淆的影響匹配效果的影像可以被過濾掉；假設一組影像I＝I_i，...，I_n以及對應的概略的位置姿態(tài)信息G＝G_i，...，G_n，以及一個可能匹配的視圖對集的子集V_i；

㈠位置姿態(tài)信息

通過GPS和IMU獲得的位置姿態(tài)信息G_i＝[R_i|t_i]，其中R_i是一個3×3的旋轉矩陣，t_i是一個3維空間向量分別代表相機的位置和姿態(tài)角；標準GPS接收機得到的全局位置坐標是以經(jīng)緯度和高程表示地球某一位置的WGS84坐標系；下一步，就是將GPS數(shù)據(jù)轉換為地心地固坐標系，地心地固坐標系是一種可以在全球表示場景重建的笛卡爾坐標系統(tǒng)；相機的定向可以用偏航角、俯仰角和滾轉角來表示，其中，偏航角是從地磁北起算的；這樣，外部姿態(tài)G_i就包括了轉換到了地心地固坐標系的GPS坐標和三個旋轉角；再加上已知的相機的內(nèi)方位參數(shù)，就可以得到每張影像I_i的完整的投影矩陣了：

該投影矩陣給出了相機的位置和姿態(tài)信息的概略值，這些概略值將用于后續(xù)的處理過程中；無人機所搭載的GPS/IMU系統(tǒng)由于系統(tǒng)跳變等原因有時會獲取一些明顯錯誤的位置姿態(tài)數(shù)據(jù)，所以這些概略值在用于后續(xù)的處理過程之前，必須首先進行預處理；預處理的主要內(nèi)容是剔除明顯錯誤的數(shù)據(jù)，用前后兩張影像的位置姿態(tài)數(shù)據(jù)的平均值來近似代替；

需要注意的是，當錯誤數(shù)據(jù)出現(xiàn)航帶的起始點或終點時，則用前后兩條航帶起始點或終點影像的位置姿態(tài)數(shù)據(jù)的平均值來代替；

㈡視圖選擇

為了識別可能存在共同對應特征點的影像，本發(fā)明為每張影像I_i選擇相應的具有足夠相似度的候選匹配影像集T_i＝T₁...T_k。接下來，影像集將根據(jù)影像對應的GPS/IMU信息得到的概略的重疊區(qū)域標準剔除影像；如果場景的精細三維模型是可見的，影像間的重疊區(qū)域可以很容易地通過視圖I_i和I_j之間相互投影的視錐體來得到；如果場景的精細三維模型是不可顯示的，只能通過估計最大景深S_i來限制影像I_i可見的區(qū)域；例如，給定了數(shù)字高程模型(DEM)，估算的地面高程可以限制相機拍攝地面的最大景深范圍，地面高程可以由公開的全球30米分辨率的DEM數(shù)據(jù)得到；而且，最大景深值S_ij可以通過影像對＜I_i，I_j＞的基線恢復出來。定義S_ij：

S_ij＝t.d(G_i，G_j)，

其中，d(.，.)表示歐氏距離，t是一個決定重建所需精度的參數(shù)。給定了這些約束條件，可以通過影像對＜I_i，I_j＞的重建計算最大景深值S：

S＝min(S_ij，S_i，S_j)，

而且，這些影像必須要有重疊，為了計算一個粗略的重疊標準定義一個平行于影像I_i且到相機中心點G_i的距離為S的平面πⁱ；R_i和R_j表示視圖I_i和I_j投影在平面πⁱ上的圖像范圍；影像的重疊度可以通過下式計算：

其中，a(.)表示投影矩形的面積；

因為SIFT等特征描述子只能適用于旋轉角度小于30度的情況，這就需要視圖的投影矩陣和標準平面πⁱ之間的夾角小于最大旋轉角度α，否則就設置為0；對于每組影像對＜I_i，I_j＞計算其重疊區(qū)域其中I_j∈T_i；如果重疊區(qū)域值大于設定的閾值，I_j就被加入到子集V_i以用于后續(xù)的精細匹配；也就是說每一幅圖像I_i只與同時滿足以下兩個條件的圖像I_j進行匹配：

其中，表示圖像I_i和I_j拍攝時的方向角，threshold為設定的閾值；

步驟二、建立極線圖：

首先，對每幅圖像提取尺度不變特征點；本發(fā)明的方法是采用高效的SIFT提取算子和描述算子，這種方法對于寬基線的圖像匹配具有很好的適應性，尤其本發(fā)明用的是公開的SiftGPU軟件；圖像對之間的特征點匹配應用基于CUBLAS矩陣運算庫的GPU加速的圖像匹配方法；

在對每個候選的視圖I_i匹配上相關的影像集V_i之后，利用五點算法進行幾何驗證；因為根據(jù)特征點描述子的匹配經(jīng)常會出現(xiàn)錯誤的外點，采用RANSAC算法進行剔除；匹配的輸出結果是用極線圖表示的結構圖，極線圖由對應影像的頂點集v＝{I₁...I_N}和邊界集ε＝{e_ij|i，j∈v}組成并是成對重建的，也就是由視圖i和j之間的相對定向e_ij＝＜P_i，P_j＞和各自的影像三角化后的點集組成；其中，

P_i＝K_i[I|O]，，P_j＝K_j[R|t]；

步驟三、計算全局一致性的旋轉集：

給定了極線圖接下來就是確定相機的初始位置和定向信息；根據(jù)兩幅圖像之間的約束，兩個相機的絕對位置姿態(tài)(R_i，t_i)和(R_j，t_j)需要滿足旋轉一致性和平移方向一致性

首先，視圖對i和j之間的相對旋轉集{R_ij}通過解超定方程組可以升級為全局一致性的旋轉集{R_i}，

R_ijR_i＝R_j，

上述是以R_i必須為標準正交的為限制條件的；然后利用SVD(奇異值分解方法)分解使R_i滿足正交約束，得到最終解可以通過解出系統(tǒng)初始的近似旋轉矩陣(不受R_i必須為標準正交的條件限制)并用Frobenius范數(shù)將近似旋轉矩陣投影到最接近的旋轉矩陣的方法解決這個限制問題；

并不是所有的對極幾何都是一樣重要的，也就是通過外極幾何得到的R_ij不是同等重要的；本文通過下面的公式來計算權重：

其中，N＝|F_ij|是視圖i和j之間的內(nèi)點數(shù)量，c_i，c_j是特征覆蓋范圍值，

其中，α是所需內(nèi)點的最小數(shù)量，表示整個影像的面積，A(F_ij，r)是特征點F_ij覆蓋范圍經(jīng)過以為圓半徑進行擴張操作后的面積。除了原始的內(nèi)點數(shù)量可以決定相對定向的可靠性外，覆蓋范圍標準更多地考慮對應特征點的空間分布；結果，具有恰當分布對應點的聚合的視圖對會比那些有同樣多的對應點但隨機分布的視圖對所占權重要高；因此，全局一致性的旋轉集可以拓展為加了權重的形式，

其中，為R_i的列(k＝1，2，3)；上式可以通過稀疏的最小二乘算子，例如應用ARPACK庫，解算出來；

步驟四、初始化相機中心點：

要實現(xiàn)相機的中心點在地心地固坐標系下的初始化，需要進行一個轉換，就是旋轉矩陣R_i必須轉換為適應GPS的方式，這可以通過將相對定向v_ij調(diào)整為相應的GPS定向

其中，v_ij是全局坐標系中圖像I_i和圖像I_j之間的相對平移，是GPS坐標系中圖像I_i和圖像I_j之間的相對平移；這是一個典型的正交的Procrustes問題，可以通過奇異值分解法解算R；

步驟五、生成對應特征點軌跡：

極線圖存儲了一個相對定向集和視圖對＜I_i，I_j＞之間的對應特征點；每張影像I_i都和鄰近一定數(shù)量的影像匹配，匹配的信息被存儲在本地的節(jié)點；需要注意的是，是一個單向圖，I_i→I_j的匹配并不一定包含I_j←I_i的匹配；然后，對于極線圖中的每張圖節(jié)點I_i，節(jié)點被聚合成軌跡(track)其中f＝＜x^k，y^k＞表示特征點在影像I_k中的坐標位置；也就是說，根據(jù)圖像匹配關系，尋找每幅圖像中的每個特征點在其他匹配圖像中對應的特征點，所有這些特征點，構成一個點軌跡，對應現(xiàn)實世界中的一個3D點；因為點軌跡是為每張影像而建，并存儲在本地，起初，點軌跡集是冗余的，例如，影像I_k的一個特征點f可以屬于幾個不同的軌跡；然后，點軌跡應用光束法平差進行整體優(yōu)化；從實用化的角度來看，由于在進行整體優(yōu)化時會涉及更多的參數(shù)，冗余的估算并不可取，因此，本發(fā)明盡可能使用最小的軌跡集來表達；為了這個目的，本發(fā)明確定的一個軌跡的子集，它只包括極線圖上每一個匹配的對應特征點一次；這是一個集合覆蓋問題，是眾所周知的最早的非完全多項式問題之一；本發(fā)明運用貪婪搜索算法來確定軌跡集的最小子集，這個最小子集隨后被用來初始化稀疏的3D結構；

步驟六、初始化3D結構：

通過前面的處理步驟，可以得到相機的方向信息集合(即標定和姿態(tài))和點軌跡還需要根據(jù)每個點軌跡確定3D點的坐標；假如相機的方向信息集總體上達不到像素級的精度且中還有外點，基于的線性三角化將導致隨機的大的重建錯誤，即3D結構初始化錯誤；事實上，可以發(fā)現(xiàn)直接的三角化方法不能保證足夠的結構初始化的精度，甚至經(jīng)常連cheirality約束都不能滿足；但是，極線圖能提供視圖對之間像素級或亞像素級精度的相機方向信息，可以滿足兩視圖三角化的精度；因此，在每個點軌跡中基于相對定向對具有最長基線的視圖對(長基線可以保證相對較低的GPS坐標誤差)進行兩視圖三角化，也就是選取每個點軌跡中特征點所在圖像的GPS坐標相差最大的兩個特征點做三角化，得到初始3D點；

步驟七、光束法平差：

給定一個量測值集合，光束法平差可以通過最小化重投影誤差來優(yōu)化相機的定向和結構恢復；

其中，x_ij是在未知相機P_i中未知3D點X_j的觀測值，v_ij是一個二進制的變量(當點X_j在影像P_i中可見時為1，否則為0)；事實上，光束法平差涉及到通過最小化非線性實值函數(shù)的平方和來調(diào)整每個3D點和相機中心點之間的射線束；光束法平差是一個龐大的，但稀疏的幾何參數(shù)值估算問題，它適用于缺失數(shù)據(jù)的情況(即并不要求每一個3D點必須在每個相機中可見)；在存在高斯噪音的情況下，非線性的最小二乘法可以實現(xiàn)最大似然估計(必要的限制條件是初始值必須足夠接近于整體最小值)；但是圖像量測值高斯噪音的假設只是設定的理想狀況，對于真實世界中基于自然特征的匹配技術的運動恢復結構問題是不確切的；在本發(fā)明中，特征點軌跡是基于極線約束幾何調(diào)整過的，但是誤匹配仍有可能出現(xiàn)；因此，最小平方和并不是一個合適的標準，需要一個魯棒的目標函數(shù)來處理外點。在光束法平差的基本實現(xiàn)形式是最小化矢量的平方和Σ_i||∈||²，其中因此，魯棒的目標函數(shù)可以以重設誤差矢量∈′_i＝w_i∈_i的權重值來實現(xiàn)，如下式所示，

由此得出，上式滿足Σ_iC(||∈_i||)＝Σ_i||∈_i||²的要求，其中，

權值w_i常被稱為衰減因子，因為它用于降低外點的影響；

步驟八、稠密點云重建：

采用PMVS算法進行稠密點云的重建；

步驟9、紋理映射：

對重建的稠密點云基于Poisson算法構建場景或物體的網(wǎng)格表面，并通過遮擋分析，實現(xiàn)紋理自動映射。

本發(fā)明的應用基礎是帶有低精度GPS/IMU信息的無人機序列影像，這里低精度是指定位精度一般在10米左右，姿態(tài)角精度一般在5度左右，這樣的定姿定位數(shù)據(jù)精度不能滿足三維重建中圖像匹配和相機定向的精度要求，無法直接依賴這些低精度的定姿定位數(shù)據(jù)進行三維重建；但是，這些低精度的GPS/IMU數(shù)據(jù)可以提供序列影像的概略地理位置和相互關系，可以提供相機的概略定向信息，這些信息都可以作為輔助信息提高基于無人機序列影像三維重建的效率和精度。

附圖說明

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

圖2為實例中無人機影像數(shù)據(jù)的待匹配集合圖；

圖3為實例中生成的極線圖；

圖4為實例中平差后的相機位置和姿態(tài)信息圖；

圖5為實例中稠密點云重建的結果圖；

圖6為實例中三維重建結果圖。

具體實施方式

實施例1

本實施例提供一種無人機序列影像批處理三維重建方法，以下將結合附圖1-6對本發(fā)明的技術方案進行詳細描述。

如圖1所示，首先要融合低精度GPS/INS信息對序列影像進行圖像特征匹配，然后通過建立極線圖、計算全局一致性的旋轉集、初始化相機中心點、生成對應特征點軌跡、初始化3D結構和光束法平差等步驟完成運動恢復結構重建過程，最后通過稠密點云重建和自動紋理映射，得到了三維重建模型。

以下將通過具體的算例對本發(fā)明的技術方案的實現(xiàn)過程予以說明。

以下結合一組由搭載低精度的GPS/IMU組合傳感器的Quickeye II型無人機獲取的某校園的圖像數(shù)據(jù)集，對本發(fā)明的技術方案進行詳細描述。

拍攝使用的相機為Canon EOS 5D Mark II，無人機照片影像的分辨率為5616×3744。無人機飛行高度為700米，平臺上掛載了定位精度在10米左右的動態(tài)單點定位的GPS，和精度在5度左右的陀螺儀，共獲得153幅無人機照片影像。

本發(fā)明的技術方案通過如下步驟實現(xiàn)：

(1)融合低精度GPS/INS信息的影像匹配

利用低精度的GPS信息和IMU信息，獲取每幅圖像的待匹配的圖像集合，將每幅圖像的對應連接繪制成矩陣的形式，得到如圖2所示的待匹配集合圖。根據(jù)飛行實驗前的無人機航線規(guī)劃數(shù)據(jù)，圖像航向重疊度最大約為70％，旁向重疊度最大約為30％，利用本文提出的視圖選擇方法，并將影像的重疊度閾值設為50％，可以確定最大鄰近視圖數(shù)d＝18，即每幅圖像最多取18幅圖像進行圖像匹配，就可以將圖像匹配時間復雜度從O(153×153)降到了O(153×18)。

(2)建立極線圖

匹配的輸出結果是用極線圖表示的結構圖，如圖3所示，極線圖由對應影像的頂點集v＝{I₁...I_N}和邊界集ε＝{e_ij|i，j∈v}組成并是成對重建的，也就是由視圖i和j之間的相對定向e_ij＝＜P_i，P_j＞和各自的影像三角化后的點集組成。

(3)計算全局一致性的旋轉集

視圖對i和j之間的相對旋轉集{R_ij}通過解超定方程組可以升級為全局一致性的旋轉集{R_i}

R_ijR_i＝R_j

上述是以R_i必須為標準正交的為限制條件的，利用SVD(奇異值分解方法)分解使R_i滿足正交約束，得到最終解

根據(jù)計算，這153張影像全局一致性的旋轉集為：

……

(4)初始化相機中心點

要實現(xiàn)相機的中心點在地心地固坐標系下的初始化，需要進行一個轉換，就是旋轉矩陣R_i必須轉換為適應GPS的方式，這可以通過將相對定向v_ij調(diào)整為相應的GPS定向

根據(jù)計算，這153張影像的初始化相機中心點坐標為：

ImageIndex，CameraPos_x，CameraPos_y，CameraPos_z

0，-3.33200979232788，2.05530166625977，-2.15661430358887

1，-1.72906887531281，2.11266350746155，-1.87004625797272

2，-0.306412011384964，2.10093069076538，-1.73969519138336

3，-0.497344762086868，1.97058546543121，-2.0764467716217

4，1.33266019821167，1.89519202709198，-2.15731501579285

5，1.38418865203857，1.90700125694275，-2.0815737247467

6，3.32213163375854，2.05892443656921，-1.49860370159149

7，3.29912209510803，2.09036254882813，-1.41904950141907

8，3.39662861824036，2.08454847335815，-1.40254628658295

9，3.45938682556152，2.06369590759277，-1.46558439731598

10，6.04636240005493，2.01808667182922，-1.26724004745483

11，8.23372077941895，2.00806713104248，-1.07459223270416

12，11.5055103302002，1.93985688686371，-0.930897831916809

13，13.8602714538574，1.91319286823273，-0.736517250537872

……

146，16.170618057251，1.88167405128479，-0.515337228775024

147，17.8826656341553，1.84787964820862，-0.456986546516418

148，20.083517074585，1.80206346511841，-0.362621396780014

149，21.4840564727783，1.80437672138214，-0.227693825960159

150，22.4221382141113，1.85853540897369，0.0575104840099812

151，23.7377300262451，1.75472116470337，-0.0763002559542656

152，25.1186275482178，1.78453361988068，0.136511072516441

(5)生成對應特征點軌跡

根據(jù)圖像匹配關系，尋找每幅圖像中的每個特征點在其他匹配圖像中對應的特征點，所有這些特征點，構成一個點軌跡，對應現(xiàn)實世界中的一個3D點，并運用貪婪搜索算法來確定軌跡集的最小子集。

根據(jù)計算，軌跡集的最小子集為：

ImageIndex，Point_id，Pos_x，Pos_y，

C000000，0，-0.885，2.106

C000000，1，-0.339，2.076

C000000，2，2.472，1.961

……

C000001，0，-0.883，0.503

C000001，1，-0.352，0.577

C000001，2，2.449，0.810

……

C000152，14876，-3.258，-2.092

C000152，14976，-2.754，-1.976

C000152，15629，-2.643，-1.972

(6)初始化3D結構

在每個點軌跡中基于相對定向對具有最長基線的視圖對(長基線可以保證相對較低的GPS坐標誤差)進行兩視圖三角化，也就是選取每個點軌跡中特征點所在圖像的GPS坐標相差最大的兩個特征點做三角化，得到初始3D點。

根據(jù)計算，得到7862個3D點初始化坐標為：

(7)光束法平差

對于上述步驟得到解算結果，光束法平差可以通過最小化重投影誤差來優(yōu)化相機的定向和結構恢復，平差后的相機位置和姿態(tài)解算結果如圖4所示，具體數(shù)值如下：

PhotoID，X，Y，Z，Omega，Phi，Kappa

0，0.0254479340093484，0.3415414621610514，0.0250366517649538，175.7103213867737800，2.0291928596802138，2.0673573899635556

1，0.0410902629089479，0.1477144927140597，0.0294442811860536，178.3745342696717300，4.0489057627394587，-3.1093510454025033

2，0.0684955405233049，-0.0479907859590106，0.0364048464871272，-176.2261276627340600，2.9858672357529628，-12.5320937664505990

……

152，2.1318594371907205，0.2284795097077054，0.3511295664437739，-179.6044286840305200，13.2343419289542870，158.7878915203033300

(8)稠密點云重建

采用PMVS算法進行稠密點云的重建，重建結果如圖5所示。

根據(jù)計算，共得到1640119個稠密點云，坐標為：

(9)紋理映射

對重建的稠密點云基于Poisson算法構建場景(或物體)的網(wǎng)格表面，并通過遮擋分析，實現(xiàn)紋理自動映射，結果如圖6所示。

以上實施例僅為說明本發(fā)明的技術思想，不能以此限定本發(fā)明的保護范圍，凡是按照本發(fā)明提出的技術思想，在技術方案基礎上所做的任何改動，均落入本發(fā)明保護范圍之內(nèi)。