結(jié)合Forstner與NCC約束的UAV圖像配準研究

賀一楠，耿 娟，秦 軍，劉 晨，楊 輝

(西南交通大學(xué)軌道交通工程遙感聯(lián)合研究中心，成都 610031)

0 引言

近年來，無人機(unmanned aerial vehicle，UAV)正成為航空遙感圖像獲取的重要手段，特別是在人跡罕至、環(huán)境比較惡劣的地區(qū)，UAV的優(yōu)越性更為明顯。但因受限于UAV自身配重較輕、抗風(fēng)能力弱、姿態(tài)不穩(wěn)定等因素影響，UAV圖像存在著復(fù)雜的幾何變形，導(dǎo)致其在圖像配準過程中存在很大的難度。分析現(xiàn)有研究情況，UAV圖像的配準、拼接和糾正處理可以分為2大類:①針對UAV圖像特點、參考傳統(tǒng)的攝影測量流程進行處理。圖像拼接前需要按照嚴格定位模型生成正射圖像，研究主要集中在拼接縫消除方面，這樣處理的結(jié)果精度高;但同時需要有高精度的姿態(tài)參數(shù)和地面控制點，工作周期長，難以滿足重大環(huán)境與自然災(zāi)害快速評估和應(yīng)急響應(yīng)的需求[1];②在不需要地面控制點的情況下，對已有圖像配準后進行拼接和糾正處理，主要包括基于頻率域和空間域的配準與加權(quán)融合[2-3]?；陬l率域的圖像配準算法(如Fourier-Mellin算法)可以很好地解決平移、旋轉(zhuǎn)和縮放的問題，但對圖像的尺度和重合度要求較高[4]?；诳臻g域的配準算法如尺度不變特征轉(zhuǎn)換(scale-invariant feature transform，SIFT)算法，拼接后的圖像質(zhì)量很高，但在特征點搜索時采用固定核尺寸，耗時長，效率不高[5];而結(jié)合Harris算子和局部相關(guān)系數(shù)(local area correlation coefficient，LACC)的利用金字塔多層匹配的多約束條件配準算法，則復(fù)雜度高且效率低[6]。針對上述情況，本文提出一種基于Forstner算子與歸一化互相關(guān)(normalized cross correlation，NCC)約束的配準方法實現(xiàn)UAV圖像的配準，并通過實驗證明該方法的有效性。

1 特征點提取與匹配

在UAV圖像配準和拼接研究中需要重點關(guān)注的是拼接算法的效率，最大可能地在既有參數(shù)條件下實現(xiàn)準確而且快速的圖像配準，才能發(fā)揮UAV的靈活高效的優(yōu)點[7]。

本文采用基于UAV圖像特征信息的特征相關(guān)匹配方法，利用圖像中明顯的特征信息(如邊緣、角點、交叉點等)不變的局部屬性對特征點進行提取，然后在特征點之間進行匹配[8-9];繼而利用最小二乘法對已匹配的特征點進行最優(yōu)化，從而實現(xiàn)UAV圖像的配準和鑲嵌。其技術(shù)流程如圖1所示。

圖1 無人機圖像配準技術(shù)流程Fig.1 Workflow of UAV image registration

UAV圖像配準的主要工作由2部分組成:①特征點提取和圖像分塊;②基于NCC匹配算法的特征點匹配。在特征點匹配之前，需要進行圖像重疊度的估算、重疊區(qū)的特征點提取以及基于信息熵的圖像分塊。

1.1 圖像重疊度估算

與傳統(tǒng)的航空遙感圖像相比，UAV圖像具有更高的重疊度，通常能夠達到航向重疊70% ～80%，旁向重疊50%以上。直接利用UAV的定位定向系統(tǒng)(position orientation system，POS)數(shù)據(jù)、基于傳統(tǒng)的攝影測量方法進行空三測量時，為保證圖像配準的效果，需要通過高精度的數(shù)字地面模型(digital surface model，DSM)生成正射圖像，而高精度的DSM則需要大量的地面控制點，處理工作量大且耗時費力。因此，本文提出基于UAV的POS數(shù)據(jù)的圖像重疊度估算，以便于圖像特征點的搜尋。像對同名點處在大面積的重疊區(qū)域內(nèi)，同時在圖像重疊區(qū)內(nèi)尋找像對的同名點比遍歷整景圖像尋找同名點要節(jié)省時間;因此在進行特征點匹配之前，估算圖像的重疊度非常有利于后續(xù)的運算。假設(shè)依照POS數(shù)據(jù)中基準圖像中心地理平面坐標為(x1，y1)，待匹配圖像坐標為(x2，y2)，根據(jù)已知的飛行航高H與鏡頭焦距f可以計算出2景圖像中心的坐標偏移量S，即

因為2景圖像的像素大小a×b與圖像尺寸p保持一致，設(shè)基準圖像的像片中心點為坐標原點(0，0)，X方向為飛行方向，Y方向為垂直飛行方向，在像平面坐標下可以得到2景圖像的4個角點的坐標，按照4個象限依次排列，即基準圖像4個角點的坐標分別為(a·p/2，b·p/2)，(-a·p/2，b·p/2)，(-a·p/2，-b·p/2)和(a·p/2，-b·p/2);待匹配圖像4 個角點的坐標分別為(a·p/2+Sx，b·p/2+Sy)，(-a·p/2+Sx，b·p/2+Sy)，(-a·p/2+Sx，-b·p/2+Sy)和(a·p/2+Sx，-b·p/2+Sy)。

沿著UAV飛行方向，組成這2景圖像重疊區(qū)域的角點坐標。按照4個象限依次排列，即基準圖像4 個角點的坐標分別為(a·p/2，b·p/2)，(-a·p/2+Sx，b·p/2)，(-a·p/2+Sx，-b·p/2)和(a·p/2，-b·p/2);待匹配圖像4個角點的坐標分別為(a·p/2，b·p/2+Sy)，(-a·p/2+Sx，b·p/2+Sy)，(-a·p/2+Sx，-b·/2+Sy)和(a·p/2，-b·p/2+Sy)。

在實際處理中，為防止特征點提取時出現(xiàn)遺漏，同時提高特征點匹配的成功率，可以在保證飛行方向的重疊區(qū)域不變的情況下，適當增加垂直飛行方向的重疊區(qū)域面積。

1.2 特征點提取

特征點是反映地物類型或區(qū)域地理分布特征的地物點，在地圖上具有準確的地理位置和明確的地理屬性及含義。在圖像處理中，特征點大致可以分為狹義特征點和廣義特征點2種。廣義特征點可以認為是一個抽象的特征區(qū)，而狹義特征點的位置本身則具有常規(guī)的屬性(例如角點、交叉點等)。

在已有的特征點提取方法中，基于模板的最小核值相似區(qū)(small univalue segment assimilating nucleus，SUSAN)算法極其耗時，不適用于大量特征點的提取，在大范圍空間尺度的數(shù)字攝影測量中很少應(yīng)用;尺度不變特征轉(zhuǎn)換(SIFT)算法包含4個步驟(尺度空間極值點探測、極值點精確定位、方向參數(shù)指定和關(guān)鍵點描述符)，自身復(fù)雜度較高，尤其是采用固定的高斯探測核，導(dǎo)致該算法的效率也很低;基于圖像噪聲檢測的Harris算法對圖像尺度很敏感，只能適用于沒有旋轉(zhuǎn)、沒有仿射變換且尺度變化比較小的UAV圖像;而著名的點特征提取算子——Moravec算子提取特征點的速度是幾種算法中最快的，但精度最差，而且受空間尺度變異和噪聲的影響很大。

本文采用攝影測量中著名的點定位算子——Forstner算子。該算子在現(xiàn)有的特征提取方法中復(fù)雜度最低，而且效率和相對精度最高。該算子通過計算各像素的Robert梯度和以像素為中心的一個窗口的灰度協(xié)方差矩陣，在圖像中尋找盡可能小而接近圓的誤差橢圓的點作為特征點[10]。在利用Forstner算子提取特征點的過程中，需要確定2個閾值:①差分算子絕對值均值的閾值;②感興趣值的閾值。閾值的選取將對特征點的提取及后續(xù)的配準精度產(chǎn)生很大的影響。閾值設(shè)置太小，將增加特征點選取及后續(xù)匹配過程的計算量;閾值太大，會增加特征點漏選、錯選的幾率。

1.3 特征點匹配

1.3.1 基于信息熵的圖像分塊

完成像對重疊區(qū)的特征點提取后，在進行特征點匹配之前，先針對這些特征點進行重疊區(qū)圖像分塊，以提高特征點匹配的速度。本文采用基于信息熵的方法對重疊區(qū)圖像進行分塊。

先將重疊區(qū)的圖像均勻地分為4塊，每塊圖像的大小為(m/2)×(n/2)，其中m為重疊區(qū)圖像的寬度，n為重疊區(qū)圖像的長度;然后計算每塊圖像矩陣的信息熵H，其中最大值為Hmax，最小值為Hmin。若最大、最小信息熵的差值△H=(Hmax-Hmin)≥th(“th”為分塊閾值，用于描述分塊的精細程度)，則說明分割有效;再將分割所得的每塊圖像作為新的處理對象進行均勻分割，直至每塊圖像矩陣的大小滿足 min(m/2a，n/2b)≥minsize，{a|a ＞1 且 a∈N}，{b|b＞1且b∈N}(“minsize”為每塊處理對象的熵的下限)，則停止分割。本文經(jīng)過反復(fù)試驗，確定待處理對象的熵的下限minsize=3.47，分塊閾值th=0.51。

1.3.2 基于NCC的特征點匹配

進行歸一化互相關(guān)(NCC)處理是為了解決圖像互相關(guān)中像素灰度隨位置變化導(dǎo)致匹配度降低、匹配系數(shù)依賴于模板大小等問題。NCC的相關(guān)系數(shù)[11]為

式中:δ(u，v)為匹配圖像和模板之間的NCC系數(shù);m 和 n 分別為匹配模板的長和寬;gi+u，j+v和 g′i，j分別為參加匹配的2景圖像中在(i+u，j+v)和(i，j)處的灰度值。2景圖像越相似，其NCC越接近1。

NCC相關(guān)系數(shù)的分子實際上是代表灰度的2個矩陣中對應(yīng)點的值兩兩相乘再求和，在沒有旋轉(zhuǎn)的情況下，可以得到比較好的的結(jié)果[11]。但是UAV圖像存在一定角度的航偏旋轉(zhuǎn)，自帶的云臺所采集的旋轉(zhuǎn)角度的精度不高，該角度數(shù)據(jù)不能直接使用。由于在待匹配圖像中選取的特征點和探測窗口與在基準圖像中的實時窗口已經(jīng)不同，加之窗口中像素的對應(yīng)關(guān)系已經(jīng)不存在，因此，NCC方法不能直接用于UAV圖像的特征點匹配。但如果將探測窗口按照一定的步長進行旋轉(zhuǎn)，則總會找到一種情況與實時窗口基本相同，并且在這種條件下NCC系數(shù)達到最大。探測窗口旋轉(zhuǎn)后，原來探測窗口中的點與旋轉(zhuǎn)后的點之間的對應(yīng)關(guān)系比較復(fù)雜，不能用統(tǒng)一的坐標轉(zhuǎn)換公式來表示，本文方法根據(jù)UAV圖像的特點分2種情況:

1)0°＜θ≤90°

2)-90°＜θ≤0°

式(3)(4)中:(i，j)為原始探測窗口的坐標;(i′，j′)為旋轉(zhuǎn)后探測窗口的坐標;m1和n1分別為原始探測窗口的長和寬;θ為探測窗口逆時針旋轉(zhuǎn)的角度。為確定旋轉(zhuǎn)后的每個像素點的灰度值，選擇雙向線性內(nèi)插法得到旋轉(zhuǎn)后像素點的灰度值。

為了精確地獲取與控制點對應(yīng)的匹配點，經(jīng)過多次試驗，設(shè)定NCC的門限閾值為0.98。采用逐次逼近算法，首先讓滑窗在-90°≤θ≤90°內(nèi)以1°為步長逆時針旋轉(zhuǎn)，同時計算NCC系數(shù);然后截取NCC取得最大值時旋轉(zhuǎn)角θ所在鄰域[θ-1°，θ+1°]作為新的旋轉(zhuǎn)域約束角，繼續(xù)讓滑窗在該區(qū)域內(nèi)改為以1′為步長進行逆時針旋轉(zhuǎn)，同時計算NCC系數(shù);依此類推，不斷地改變步長，通過約束逼近求得使NCC的值滿足閾值時的同名匹配點。實現(xiàn)特征點匹配的技術(shù)流程如圖2所示。

圖2 特征點匹配流程Fig.2 Workflow of feature point matching

1.3.3 誤差點剔除

通過NCC系數(shù)匹配法匹配的特征點對中會存在一些錯誤的匹配點對，必須將這些錯誤的點對剔除，否則會導(dǎo)致解算出的轉(zhuǎn)換模型存在很大的誤差。在正常飛行條件下，由UAV拍攝的2景相鄰圖像的旋轉(zhuǎn)畸變和縮放畸變相對于飛行運動很小，所以連接每對同名點所構(gòu)成的向量之間距離的誤差也很小。如果2對匹配點是正確的匹配點，則它們構(gòu)成的向量之間距離的誤差應(yīng)該很小，反之則較大。因此可以利用該閾值進行誤差點剔除，其計算步驟包括:①計算圖像中所有匹配點對所構(gòu)成的向量的距離，并求其均值E;②定義該均值E作為條件閾值;③剔除匹配點對中與該閾值E相差過大的匹配點對(圖3)。

圖3 誤差匹配點消除示意圖((A1，B1)，(A2，B2)，(A3，B3)及(A4，B4)為匹配點對)Fig.3 Sketch map for elimination of false matching points

圖3 中匹配點對(A3，B3)所構(gòu)成的向量A3B3的距離與閾值E相差過大，故作為誤差點被剔除。

近年來,小兒肺炎發(fā)病率居高不下,嚴重影響患兒身體健康。臨床發(fā)現(xiàn),小兒肺炎發(fā)病的主要原因為病原體通常從呼吸道侵入,引發(fā)肺泡、小支氣管以及肺間質(zhì)炎癥,造成小支氣管管腔變得狹窄甚至阻塞,導(dǎo)致通氣障礙;而炎癥會造成呼吸膜厚度增加,造成很多炎癥滲出物存于肺泡腔內(nèi),導(dǎo)致?lián)Q氣障礙。通氣與換氣障礙會導(dǎo)致機體缺氧與二氧化碳潴留,從而引發(fā)一系列病理與生理性改變,嚴重者甚至?xí)霈F(xiàn)呼吸衰竭癥狀,威脅患兒生命安全[2]。因此,尋找一種安全有效的治療方案具有重要的臨床意義。

2 轉(zhuǎn)換模型

圖像配準就是在2景UAV圖像之間建立映射關(guān)系，確定相應(yīng)的變換參數(shù)，并進行仿射變換?；鶞首鴺讼?xr，yr)和校正坐標系(xc，yc)存在一定的映射關(guān)系[12]，即

式中:β為旋轉(zhuǎn)角度;△x和△y分別為橫向位移和縱向位移;s為縮放因子，在固定飛行高度和固定相機焦距條件下，s=1。

采用最小二乘方法，通過式(6)可以求得各參數(shù)，即

式中:

計算步驟如下:①為了快速地進行匹配和坐標轉(zhuǎn)換計算，選擇多個探測窗口進行匹配，把每個探測窗口左上角的點作為匹配點，利用全局搜索策略和NCC匹配算法在基準圖像中尋找相匹配的同名點;在搜索過程中，因不知道旋轉(zhuǎn)的角度，故探測窗口按照逐步逼近算法進行旋轉(zhuǎn)，并計算NCC系數(shù)，取NCC系數(shù)滿足閾值位置的點作為這一控制點的匹配點;②不斷重復(fù)步驟①，直到待匹配圖像中的特征點在基準圖像中都找到對應(yīng)的匹配點;③根據(jù)原始參考圖像和旋轉(zhuǎn)后參考圖像的位置關(guān)系，建立上述同名點在原始參考圖像中的位置轉(zhuǎn)換模型;④根據(jù)位置轉(zhuǎn)換模型，利用最小二乘法計算轉(zhuǎn)換關(guān)系參數(shù)(包括橫向位移、縱向位移和旋轉(zhuǎn)角)。

3 實驗與分析

3.1 實驗數(shù)據(jù)與設(shè)備

本次實驗采用拍攝于四川省阿壩藏族羌族自治州東南部的UAV圖像(圖4)和同步獲取的POS數(shù)據(jù)(表1)。

圖4 實驗用UAV圖像像對Fig.4 Image pair of UAV images for experiment

表1 實驗圖像的POS數(shù)據(jù)Tab.1 POS data of experimental images

仿真實驗算法基于IDL語言;圖像處理工作站配置為CPU英特爾Xeon志強E3-1 230，3.30 Hz四核處理器，內(nèi)存單元為8 GB，GPU為英偉達Quadro K2000。

3.2 實驗結(jié)果與分析

根據(jù)與UAV圖像同步獲取的POS數(shù)據(jù)，利用前文介紹的UAV圖像重疊度估算方法，結(jié)合式(1)及角點坐標得到UAV圖像的重疊區(qū)域(圖5)。

圖5 圖像重疊度估算結(jié)果(綠色矩形框內(nèi)為圖像重疊區(qū))Fig.5 Estimation result of image overlap area

在利用Forstner算子提取特征點的過程中，經(jīng)過反復(fù)試驗，設(shè)置初選差分閾值為45，計算感興趣值窗口為5像素×5像素，非最大抑制窗口為7像素×7像素，感興趣閾值為1.2。在圖像重疊區(qū)內(nèi)提取特征點的結(jié)果如圖6所示。

圖6 圖像重疊區(qū)內(nèi)特征點提取結(jié)果(綠色矩形框內(nèi)為圖像重疊區(qū))Fig.6 Extraction result of feature points in mage overlap area

利用UAV圖像在不同的探測矩陣尺寸下進行大量實驗，發(fā)現(xiàn)隨著探測矩陣尺寸的逐漸增大，運算的耗時量在增加，匹配點對數(shù)量在減少;當探測矩陣尺寸增大到一定程度時，在耗時量增加的同時，匹配點對數(shù)目趨于穩(wěn)定。因此，在耗時量和提取的匹配點對數(shù)量上必然有最優(yōu)的探測矩陣尺寸。為了找到最優(yōu)探測矩陣尺寸，本文選用多個航空攝影區(qū)域的UAV圖像進行試驗，比較在不同探測矩陣尺寸下基于NCC存在旋轉(zhuǎn)的特征點匹配的耗時以及提取的匹配點對數(shù)量。受篇幅所限，僅列舉部分實驗數(shù)據(jù)(表2)進行說明。其中，當探測矩陣尺寸小于9像素×9像素時，不利于探測矩陣在旋轉(zhuǎn)過程中的坐標轉(zhuǎn)換計算，因此只考慮探測矩陣尺寸大于或等于9像素×9像素的情況。

表2 不同探測矩陣尺寸實驗結(jié)果Tab.2 Experimental result in different detection matrix size

由表2可以看出，在探測矩陣尺寸從9像素×9像素逐漸增到27像素×27像素的過程中，計算出的匹配點對的數(shù)量的臨界點在19像素×19像素，超過該臨界點后耗時增多，匹配點對數(shù)目趨于穩(wěn)定。因此，選擇尺寸為19像素×19像素的探測矩陣作為最優(yōu)尺寸探測矩陣。若尺寸選擇過小，雖然耗時較少，但會增加解算轉(zhuǎn)換模型的運算量;若尺寸選擇過大，提取的匹配點對雖趨于穩(wěn)定，但耗時過多。

在采用基于NCC值進行特征點匹配后，剔除像對中的誤差點對，得到匹配點分布圖(圖7)。

圖7 剔除誤差匹配點后的匹配點分布(綠色矩形框內(nèi)為圖像重疊區(qū))Fig.7 Distribution of matching points after removing false matching points

最終利用本文設(shè)計的方法得到大面積區(qū)域 UAV圖像配準后的鑲嵌結(jié)果(圖8)。

圖8 無人機圖像配準鑲嵌結(jié)果Fig.8 Mosaic result of UAV images after registration

4 結(jié)論

1)利用Forstner算子提取特征點的效果受初選差分閾值的影響很大，而且對某些銳度比較低、連續(xù)大面積相似地物的圖像而言，若閾值設(shè)置不合理，會導(dǎo)致提取出的特征點過少或過多，致使在隨后的特征點匹配中無法找到足夠的同名點對，并為特征點匹配帶來繁重的計算負擔(dān)。

2)采用基于信息熵的自適應(yīng)圖像分塊算法，并結(jié)合逐次逼近的互相關(guān)系數(shù)NCC的特征點匹配，很好地控制了特征點的均勻性問題，并減小了匹配點對的計算量。由于圖像分塊的閾值取決于圖像分塊過程中信息熵的大小，因此需要根據(jù)不同圖像進行圖像分塊閾值的優(yōu)化設(shè)置。

3)通過實驗證明，在UAV圖像特征點匹配的有效性、適應(yīng)性和抗噪性等方面，本文提出的結(jié)合Forstner與NCC約束的UAV圖像配準方法是切實有效的，并保持了較好的魯棒性。

[1]張永軍.無人駕駛飛艇低空遙感影像的幾何處理[J].武漢大學(xué)學(xué)報:信息科學(xué)版，2009，34(3):284-288.Zhang Y J.Geometric processing of low altitude remote sensing images captured by unmanned airship[J].Geometrics and Information Science of Wuhan University，2009，34(3):284-288.

[2]蘇俊英.SIFT特征匹配無人飛艇多光譜影像拼接[J].應(yīng)用科學(xué)學(xué)報，2010，28(6):616-620.Su J Y.Mosaicing of multiple spectrum images acquired from unmanned airship with SIFT feature matching[J].Journal of Applied Sciences-Electronics and Information Engineering，2010，28(6):616-620.

[3]狄穎辰，陳云坪，陳瑩瑩，等.無人機圖像拼接算法綜述[J].計算機應(yīng)用，2011，31(1):170-174.Di Y C，Chen Y P，Chen Y Y，et al.Survey on image mosaic algorithm of unmanned aerial vehicle[J].Journal of Computer Applications，2011，31(1):170-174.

[4]王 斌，王偉鋒.一種基于局部灰度匹配的無人機圖像拼接算法[J].中國石油大學(xué)學(xué)報:自然科學(xué)版，2009，33(2):169-173.Wang B，Wang W F.Unmanned aerial vehicle image mosaic algorithm based on local gray fitting[J].Journal of China University of Petroleum:Edition of Natural Science，2009，33(2):169-173.

[5]白廷柱，侯喜報.基于SIFT算子的圖像匹配算法研究[J].北京理工大學(xué)學(xué)報，2013，33(6):622-627.Bai T Z，Hou X B.An improved image matching algorithm based on SIFT[J].Transactions of Beijing Institute of Technology，2013，33(6):622-627.

[6]唐 敏，李永樹，魯 恒.無人機影像的同名點匹配[J].光電工程，2012，39(3):19-24.Tang M，Li Y S，Lu H.Correspondence points matching of UAV images[J].Opto-Electronic Engineering，2012，39(3):19-24.

[7]魯 恒，李永樹，何 敬，等.一種基于特征點的無人機影像自動拼接方法[J].地理與地理信息科學(xué)，2010，26(5):16-19.Lu H，Li Y S，He J，et al.A automatic mosaic method in unmanned aerial vehicle images based on feature points[J].Geography and Geo-Information Science，2010，26(5):16-19.

[8]Zhang N，Wu X L.Lossless compression of color mosaic images[J].IEEE Transactions on Image Processing，2006，15(6):1379-1388.

[9]Lucchese L，Leorin S，Cortelazzo G M.Estimation of two-dimensional affine transformations through polar curve matching and its application to image mosaicking and remote-sensing data registration[J].IEEE Transactions on Image Processing，2006，15(10):3008-3019.

[10]林曉帆，林立文，鄧 濤.基于SURF描述子的遙感影像配準[J].計算機工程，2010，36(12):216-218.Lin X F，Lin L W，Deng T.Remote sensing image registration based on SURF descriptor[J].Computer Engineering，2010，36(12):216-218.

[11]孫卜郊，周東華.基于NCC的存在旋轉(zhuǎn)的圖像匹配方法[J].傳感器與微系統(tǒng)，2008，27(5):43-45.Sun B J，Zhou D H.Rotated image registration method based on NCC[J].Transducer and Microsystem Technologies，2008，27(5):43-45.

[12]Li H，Manjunath B S，Mitra S K.A contour-based approach to multisensor image registration[J].IEEE Transactions on Image Processing，1995，4(3):320-334.