面向城區寬基線立體像對視角變化的結構自適應特征點匹配

陳 敏，朱 慶，何海清，嚴少華，趙怡濤

1. 西南交通大學地球科學與環境工程學院，四川 成都 611756； 2. 東華理工大學測繪工程學院，江西 南昌 330013

影像匹配是遙感影像處理與應用的關鍵步驟之一。經過數十年的發展，研究人員提出了許多適用于不同類型影像的匹配方法[1-5]。現有影像匹配方法大體上可分為兩類：基于灰度的匹配方法和基于特征的匹配方法[6-7]。基于灰度的匹配方法能夠獲得亞像素級精度，但對影像灰度變化和幾何變形比較敏感[8]。相對而言，基于特征的匹配方法能夠較好地克服基于灰度的匹配方法對影像灰度變化和幾何變形穩健性不足的問題。隨著SIFT(scale invariant feature transform)算法[9]的成功，基于特征的匹配方法受到了越來越多的關注[10-11]。

針對影像視角變化，研究人員通過模擬影像仿射或投影空間，并在模擬空間進行特征匹配，獲得了較好的匹配結果[12]。但這類方法時間效率較低，在實際應用中受到限制。在攝影測量領域，高精度POS(position and orientation system)數據通常被用來輔助大視角變化影像(如傾斜影像)的匹配，即利用POS信息對影像進行粗糾正，整體上降低影像幾何變形的影響，再采用傳統方法進行特征點匹配[13-17]。這類方法能夠在一定程度上改善影像匹配的效果，但全局變換難以準確描述影像之間的局部幾何變形。針對這個問題，將整幅影像分成多個子區域，分別對子區域進行特征點檢測和匹配，可以克服全局影像幾何糾正的不足，增加匹配點的數量[18-19]。對于無高精度POS數據的情況，可通過初匹配獲取一定數量的匹配點來估計立體像對之間的幾何變換模型，進而對影像進行粗糾正[19-21]。基于影像粗糾正的方法雖然能夠改善影像匹配效果，但仍然存在以下問題：①影像糾正只能在一定程度上緩解平面場景的幾何變形，城區影像由于存在顯著的遮擋問題，對于位于視差不連續處的特征點(如建筑物角點或邊緣附近的特征點)，無論經過影像全局幾何糾正還是分區域處理，都難以在同名點之間獲得影像內容一致的特征區域，進而產生錯誤匹配；②通過影像初匹配進行幾何糾正的方法依賴于初始匹配結果，對于存在大視角變化的城區寬基線影像，現有方法難以獲得可靠的初匹配結果，導致最終匹配結果不可靠。

為此，本文面向無高精度POS信息的城區寬基線影像，針對傳統方法為同名點計算的特征區域影像內容不一致，導致同名點特征描述符相似度低、匹配失敗的問題，充分挖掘特征點鄰域幾何結構信息，提出結構自適應的特征區域計算方法，在不同視角情況下獲得影像內容一致性的同名特征區域和相似特征描述符，并設計可靠的特征點匹配算法，提高正確匹配點對的數量和匹配正確率。

1 本文方法

城區影像場景多為人工建筑物，能夠提取大量的點特征和直線特征。本文方法通過挖掘點特征與直線特征之間的幾何關系來實現城區寬基線影像特征點匹配：首先，對立體像對提取點特征和直線特征，利用特征點鄰域內的直線特征表達特征點的幾何結構方向信息，為特征點計算結構自適應的不變特征區域和特征描述符，通過雙向匹配策略獲取可靠的初匹配點，并估計立體像對核線幾何關系；其次，針對部分特征點難以僅利用幾何結構方向信息構建視角不變特征區域的問題，聯合特征點幾何結構方向信息和核線幾何約束，構建視角不變特征區域進行特征匹配；最后，設計特征匹配擴展算法增加匹配點數量，提高特征匹配率，并通過RANSAC(random sample consensus)算法[22]剔除錯誤匹配。本文方法的整體流程如圖1所示，其中斜體標記步驟為本文方法的關鍵步驟。

圖1 本文方法整體流程Fig.1 Flow chart of the proposed method

1.1 結構自適應特征點雙向匹配

結構自適應的特征點雙向匹配是本文方法的初匹配步驟，目的是獲取一定數量的匹配點并估計立體像對核線幾何關系，其算法流程如圖2所示。

圖2 結構自適應特征點雙向匹配流程Fig.2 Flow chart of the structure adaptive bidirectional feature point matching method

具體匹配方法如下：

(1) 利用直線特征表達特征點的幾何結構方向信息。如圖3所示，以特征點pi為中心，確定大小為m×m的局部鄰域Ri，提取與Ri相交的非平行直線特征。以特征點為原點，以與直線特征平行的方向向量表示該特征點的幾何結構方向。在每個方向上，以鄰域內與該方向一致且最長的直線特征的長度，作為該幾何結構方向的向量長度。

(2) 挖掘幾何結構方向信息構建特征點支撐區域。根據幾何結構方向數量差異，將特征點分以下3種情況進行處理：

1) 如果特征點具有3個及以上幾何結構方向，對其中任意兩個滿足夾角α∈[θ,π-θ]的幾何結構方向，分別在對應的幾何結構方向上尋找顯著點。其中，在每個方向上尋找顯著點的方法如圖4所示[23]。圖4中，pi為特征點，Oi表示pi的一個幾何結構方向，為Oi確定一個長度為|Oi|+S，寬度為2S的向量支撐區域(其中，|Oi|表示向量Oi的模，參數S用于控制向量支撐區域的尺寸)。如果向量支撐區域內存在直線特征且與該幾何結構方向的交點也位于向量支撐區域內，則認為該交點是該幾何結構方向上的一個顯著點。由兩個方向上的顯著點與特征點構成平行四邊形區域，即為特征點支撐區域(圖5(a))。如果在一個幾何結構方向上存在多個顯著點，則每個顯著點分別用于構建特征點支撐區域，得到多個支撐區域。

2) 如果特征點具有兩個幾何結構方向，首先沿各個方向尋找顯著點，然后以特征點為中心，在幾何結構方向的反方向確定對稱點作為虛擬顯著點，最后分別由顯著點、虛擬顯著點和特征點確定支撐區域(圖5(b))。

3) 如果特征點少于兩個幾何結構方向，則無法利用上述方法構建視角不變支撐區域。后續1.3節的特征點匹配擴展步驟將對這類特征點進行處理。

圖4 幾何結構方向顯著點的確定Fig.4 Stable point determination on astructure orientation

圖5所示為構建特征點支撐區域示意圖。圖5(a)所示為特征點具有3個幾何結構方向的情況，由于O1的向量支撐區域內沒有直線特征，即該方向上沒有顯著點，因此無法與其他幾何結構方向一起構建特征點支撐區域。O2和O3的向量支撐區域內分別存在直線特征且交點也位于向量支撐區域內，因此可以在O2和O3所夾范圍內確定一個特征點支撐區域(藍色虛線標記區域)。圖5(b)所示為特征點具有兩個幾何結構方向的情況，獲得顯著點p1和p2以后，分別在幾何結構方向的反方向確定對稱點q1和q2作為虛擬顯著點，形成兩個特征點支撐區域(橙色和藍色虛線標記區域)。分別由顯著點和虛擬顯著點確定特征點支撐區域可以保證當特征點位于視差不連續區域時至少有一個支撐區域具備視角不變性。在圖5(b)所示情況中，顯著點確定的支撐區域(橙色虛線標記區域)的影像內容不具備視角不變性，而虛擬顯著點確定的支撐區域(藍色虛線標記區域)的影像內容具備視角不變性。

圖3 特征點幾何結構方向的表達Fig.3 Interest point structure orientation

圖5 特征點支撐區域Fig.5 Interest point support region

(3) 基于特征點支撐區域計算特征區域和特征描述符。如果一個特征點存在多個支撐區域，將該特征點視作多個不同的特征點分別分配一個支撐區域。將平行四邊形支撐區域歸一化得到正方形特征區域。在支撐區域歸一化時，為所有特征點設置相同的特征區域尺寸Tr×Tr，并分別將支撐區域中特征點及其對角線頂點映射到正方形特征區域的左下角和右上角頂點。由支撐區域與特征區域4個頂點的對應關系計算單應性矩陣進行特征區域歸一化。將特征區域劃分為16個子區域，在每個子區域內計算8個方向的梯度方向直方圖，得到128維特征描述符。最后，對特征描述符進行歸一化處理，提高特征描述符對影像光照變化的穩健性[9]。

通過本文提出的特征區域歸一化處理，一方面可以避免建筑物同一頂(側)面上不同位置的角點因為對應于同一個支撐區域而產生錯誤匹配，另一方面能夠消除影像旋轉變化的影響。如圖6所示，特征點p1和p2雖然支撐區域相同，但是經過本文方法得到的特征區域和特征描述符具有較強的可區分性，能夠避免錯誤匹配。雖然特征點p1和p3所在的影像存在旋轉變化，但是本文提出的特征區域歸一化方法能夠消除影像旋轉變化，得到相似的特征區域和特征描述符，提高特征匹配率。

圖6 特征區域和特征描述符計算過程Fig.6 Feature region and descriptor computation

(4) 利用雙向的NNDR(nearest neighbor distance ratio)匹配策略[24]進行特征點匹配，并結合RANSAC算法剔除錯誤匹配，得到初匹配集合MSet1和基礎矩陣F。

雖然初匹配能夠正確匹配部分特征點，但是該方法只對滿足以下條件的特征點有效：特征點具有至少兩個滿足夾角約束條件并且能夠獲得顯著點(圖4所示方法)的幾何結構方向。該條件導致初匹配步驟獲得的匹配點數量有限。為了提高匹配點數量，本文方法在初匹配之后分別設計雙重核線約束的結構自適應特征點匹配(1.2節)和特征點匹配擴展(1.3節)算法。前者用于處理具有兩個及兩個以上幾何結構方向的未匹配特征點；后者用于處理經過前面兩個匹配步驟仍然未匹配成功的特征點。

1.2 雙重核線約束的結構自適應特征點匹配

(1)

圖7 核線約束的特征點支撐區域Fig.7 Feature support region computation based on epipolar geometric constraint

獲得參考影像特征點及其候選同名特征點的支撐區域以后，采用1.1節提出的方法計算特征區域和特征描述符，并基于NNDR匹配策略[25]得到特征點匹配集合MSet2。

1.3 特征點匹配擴展

根據特征點是否位于某個已匹配的特征點支撐區域內，將未匹配的特征點分成兩類(第1類：是；第2類：否)，并分別進行匹配：

1.3.1 結合幾何約束和特征描述符相似性匹配第1類特征點

在幾何約束方面，如圖8所示，假設(pi,qi)為一對已匹配的特征點，圖中實線平行四邊形區域為其支撐區域。X為參考影像上位于點pi支撐區域內的一個未匹配特征點，由搜索影像上所有位于特征點qi支撐區域內的未匹配特征點構成X的候選匹配點集合CX。由于特征點支撐區域是基于直線特征確定的局部區域，可以近似為平面區域。由仿射幾何可知，如果特征點X與特征點Y∈CX為一對同名點，則|XA|/|XB|=|YE|/|YF|，且|XC|/|XD|=|YG|/|YH|。其中，點A、B、C、D分別為過特征點X且與點pi的支撐區域的邊平行的直線與點pi的支撐區域的交點；點E、F、G、H分別為過特征點Y且與點qi的支撐區域的邊平行的直線與點qi的支撐區域的交點。

圖8 第1類未匹配特征點匹配擴展Fig.8 The matching expansion for the unmatched point in the first class

計算特征描述符時，首先，根據X在pi的支撐區域中的位置確定X的支撐區域：如果|XA|≤|XB|，則點B被視作一個顯著點，否則點A被視作顯著點；如果|XC|≤|XD|，則點D被視作第2個顯著點，否則點C被視作第2個顯著點。由點X與兩個顯著點共同確定特征點X的支撐區域；然后，根據顯著點的對應關系(A→E、B→F、C→G和D→H)，確定候選匹配點對應的顯著點及支撐區域；最后，采用1.1節所述方法計算特征點X及其候選匹配點的特征區域和特征描述符。

按式(2)計算參考影像特征點與所有候選同名特征點的相似性度量值sim(X,Y)，并尋找最相似的候選匹配點，如果其相似性度量值大于閾值Tsim，則認為該特征點與參考特征點為一對匹配點。完成第1類特征點匹配以后，得到匹配集合MSet3。

sim(X,Y)=

(2)

式中，τ為仿射不變量閾值；DescX和DescY分別為特征點X及其候選匹配點Y的特征描述符。

1.3.2 基于單應變換匹配第2類特征點

首先，基于前面所有匹配結果{MSeti,i=1,2,3}估計影像之間的單應矩陣H；然后，為參考影像上所有第2類未匹配特征點確定以特征點為中心的正方形特征區域，基于單應變換為搜索影像上所有第2類未匹配特征點確定以特征點為中心的四邊形特征區域，并將四邊形特征區域歸一化為正方形區域；最后，計算所有未匹配特征點的特征描述符，并在核線約束下通過NNDR匹配策略得到匹配集合MSet4。利用RANSAC算法對所有匹配結果{MSeti,i=1,…,4}剔除錯誤匹配，得到最終匹配結果。

2 試驗及結果分析

2.1 試驗數據

為了驗證本文算法的有效性，分別采用6對典型的局部影像塊、一組原始大小的三視航空傾斜影像和一組原始大小的三視無人機影像進行特征點匹配試驗。如圖9所示：(a)和(b)所示分別為平房和廣場區域傾斜像對，其中參考影像為下視影像，搜索影像為斜視影像；(c)所示為無人機影像對；(d)—(f)所示均為傾斜影像，其中(d)和(e)中參考影像為下視影像，搜索影像為斜視影像，(f)中參考影像和搜索影像均為斜視影像；(g)所示為三視航空傾斜影像；(h)所示為三視無人機影像。試驗數據詳細信息見表1。

表1 試驗數據詳細信息

2.2 參數設置

本文方法相關參數設置為：構建特征點鄰域幾何結構方向向量時，特征點鄰域大小m×m=11×11像素；鄰域幾何結構方向向量支撐區域參數s=20像素；鄰域幾何結構方向向量夾角閾值θ=10°；歸一化特征區域尺寸Tr×Tr=65×65像素；點到核線的距離閾值Te=20像素；特征匹配擴展算法中，仿射不變量閾值τ=0.3，特征描述符相似性閾值Tsim=0.65。在本文試驗中，Harris算子[27]和LSD算子[28]分別被用于提取點特征和直線特征。本文試驗中所有對比方法的參數均按原文獻作者推薦的參數值進行設置。本文所有試驗均在相同平臺環境(Windows 10，Intel Core i7 3.6 GHz，RAM 32 GB)下完成。

2.3 試驗結果分析

2.3.1 局部影像像對匹配結果分析

本文試驗首先基于局部影像對(圖9(a)—(f))進行兩兩匹配，將本文方法與多種匹配算法進行對比分析以驗證本文方法對典型影像區域的有效性。對比方法包括：分別將Harris-Affine算子[25]、Hessian-Affine算子[25]、MSER算子[26]和DoG算子[9]與SIFT特征描述符和NNDR匹配策略組合而成的4種特征匹配算法(HarAff、HesAff、MSER和SIFT算法)、ASIFT算法[12]，以及基于影像粗糾正的ISIFT算法[19-20]。

在這部分試驗中，以正確匹配特征對數和匹配正確率(正確匹配特征對數/總匹配對數)為評價指標，統計結果如圖10所示。其中，通過人工檢查的方式來統計正確匹配特征對數。

從圖10(a)所示的正確匹配對數可以看出，HarAff、HesAff、MSER和SIFT 4種方法在所有像對上都只能獲得少量的匹配對，表明這4種方法難以適用于視角變化較大且存在遮擋問題的城區寬基線影像。相對于上述4種方法，ASIFT和ISIFT算法通過改進匹配策略來提高算法的穩健性。其中，ASIFT算法通過模擬影像仿射空間來消除影像之間的幾何變形，獲得了更好的匹配效果。影像對2中場景深度變化不顯著，ASIFT算法模擬的仿射空間能夠較好地擬合影像局部區域，因此ASIFT算法在影像對2上獲得了最多的匹配對。但是ASIFT方法中模擬的仿射空間不連續，在影像視角變化和場景深度變化較大的影像對4、5和6上，許多局部區域沒有被模擬的仿射空間所覆蓋，即難以通過模擬仿射空間來消除這些局部區域的幾何變形，因此ASIFT方法在這3對影像上獲得的特征對數較少。此外，ASIFT算法為特征點分配規則特征區域的方法難以適用于地物遮擋和視差不連續的情況。ISIFT算法也是基于影像模擬和粗糾正的思想，但是ISIFT算法只對整幅影像進行一次模擬，當參考影像與待匹配影像視角變化大且影像場景深度變化較大時，整幅影像之間不服從同一個全局變換模型，ISIFT方法獲得的幾何變換模型只能糾正影像中的部分區域。因此ISIFT方法雖然能夠改善SIFT方法的匹配結果，但是其改善程度有限。此外，ISIFT算法依賴于初匹配的結果。如圖中所示影像對5的匹配結果，ISIFT算法因為初匹配結果難以準確估計影像之間的幾何變換模型，最終匹配失敗。

相對于以上方法，本文方法在除影像對2以外的所有5對影像上都得到了最多的正確匹配。這主要得益于本文方法能夠根據特征點鄰域結構自適應地獲取影像內容一致的特征區域。無論特征點位于平面區域或是視差不連續區域，本文方法得到的同名特征區域之間都具有較高的相似度，更容易在特征匹配過程中被正確識別出來。此外，本文方法中的特征匹配擴展算法有助于獲得更多的匹配對。

從圖10(b)所示的匹配正確率可以看出，本文方法和ASIFT算法的匹配正確率優于其他方法。當影像初匹配能夠獲得一定數量的正確匹配用于估計影像幾何變換模型時，ISIFT算法也能獲得較高的匹配正確率，但是當影像視角變化導致無法通過初匹配來估計影像幾何變換模型時，ISIFT算法將匹配失敗。

2.3.2 完整三視影像匹配結果分析

除了采用局部影像像對進行算法驗證以外，本文試驗還利用兩組三視影像測試本文方法的匹配性能。鑒于SIFT算法的廣泛應用以及ASIFT算法對影像視角變化的穩健性，這部分試驗將本文方法與SIFT算法和ASIFT算法進行對比分析。在具體實施時，考慮到原始SIFT算法和ASIFT算法在處理較大尺寸影像時計算內存開銷非常大，且算法時間效率極低，本文采用更加高效的GPU版本的SIFT算法(SIFTGPU[29])以及下采樣模式的ASIFT算法(先將原始影像下采樣為800×600像素大小的影像進行匹配，再將匹配結果反算回原始影像[30])。此部分試驗以三度重疊匹配數量和匹配效率為評價指標。3種方法在兩組三視影像上匹配的統計結果如下表2所示，三度重疊匹配如圖11和圖12所示。由于ASIFT方法的三度重疊匹配數量為0，因此圖11中只列出SIFTGPU方法和本文方法的結果。

表2 三視影像匹配統計結果

圖11 三視影像1的三度重疊匹配結果Fig.11 Three-time overlapped matches on three-view image dataset 1

圖12 三視影像2的三度重疊匹配結果Fig.12 Three-time overlapped matches on three-view image dataset 2

從表2統計結果可以看出，在三度重疊匹配數量方面，本文方法在兩組數據上獲得的三度重疊匹配數量都遠超SIFTGPU和ASIFT算法。尤其在三視影像1上，影像之間視角變化大，且影像場景為密集建筑區域，大量特征點位于視差不連續的邊緣附近，SIFTGPU和ASIFT算法幾乎匹配失敗，而本文方法仍然能夠獲得463個三度重疊匹配；在算法時間效率方面，SIFTGPU算法的時間效率最高。本文方法在分步匹配中利用初匹配估計同名核線來約束后續匹配過程，時間效率優于ASIFT算法，但相對于SIFTGPU而言，運算效率仍然較低。

此外，表2統計結果顯示ASIFT方法獲得的三度重疊匹配少于SIFTGPU方法，尤其在三視影像1上的三度重疊匹配數量為0。對ASIFT方法在該數據集上的匹配結果進行仔細檢查發現：ASIFT方法在三視影像集1中三對立體像對之間獲得的兩兩匹配數量分別為108對、33對和222對，如圖13所示。其中，影像1和影像3由于視角差異太大(圖13(b))，ASIFT方法獲得的匹配點數量非常少，直接影響了三度重疊匹配的數量。

圖13 ASIFT方法(下采樣模式)在三視影像1上的兩兩匹配結果Fig.13 Matches of ASIFT (down-sampling mode) on image pairs in three-view image dataset 1

為了進一步驗證ASIFT算法在三視影像1上的匹配效果，使用普通模式(直接在原始影像上進行匹配)的ASIFT方法對三視影像1進行匹配試驗。具體實施時，采用OpenCV中的ASIFT算子，分別在三視影像1中的3幅影像上提取了2 523 709、3 746 247和3 402 273個特征點。數百萬個特征點進行盲匹配和窮舉搜索帶來了巨大的時間開銷(約67個小時)，然而三度重疊匹配數量仍然是0。3幅影像兩兩匹配的結果如圖14所示。從圖14(b)可以看出，對于視角變化非常大的影像1和影像3，ASIFT算法獲得的匹配點非常少，與圖13(b)的結果一致。此外，從圖14(a)和圖14(c)所示匹配結果可以看出，雖然ASIFT方法在兩兩影像之間能夠獲得一些匹配點，但是匹配點在多視影像上的重復率非常低。

圖14 ASIFT方法(普通模式)在三視影像1上的兩兩匹配結果Fig.14 Matches of ASIFT (general mode) on image pairs in three-view image dataset 1

綜上所述，本文方法在匹配效果方面優于SIFTGPU方法和ASIFT方法，在匹配時間效率方面優于ASIFT方法，但低于SIFTGPU方法。筆者將在后續研究中通過算法和程序優化提高本文方法的時間效率。

3 結 論

本文針對城區寬基線影像視角變化導致傳統方法難以為同名點計算影像內容一致的特征區域和相似特征描述符，進而導致匹配失敗的問題，提出了一種結構自適應的特征點匹配方法。本文方法的創新之處在于利用城區影像點特征與直線特征的幾何關系定義了特征點幾何結構方向信息，構建了結構自適應的特征區域和特征描述符，在此基礎上設計特征匹配和擴展算法，實現了城區寬基線影像的可靠匹配。以上改進使得本文方法能夠較好地處理因影像視角變化導致的幾何變形和遮擋問題，對于大視角變化的城區寬基線影像能夠獲得較好的匹配結果。但是，由于本文方法在特征匹配時利用了粗略的核線約束，因此匹配結果中的誤匹配都滿足該約束條件。在后續剔除誤匹配時，部分誤匹配難以通過(基于基礎矩陣的)RANSAC算法來剔除。此外，本文方法的運算效率仍需進一步提高。后續工作將研究如何有效剔除錯誤匹配，并通過算法和程序優化提高本文方法的運算效率。