高分三號(hào)順軌影像快速匹配并行計(jì)算方法

耿 勇 勝，洪 中 華*，童 小 華，劉 世 杰，馮 永 玖，張 云，周 汝 雁，韓 彥 嶺，潘 海 燕，王 靜，楊 樹 瑚，徐 利 軍

(1.上海海洋大學(xué)信息學(xué)院，上海 201306；2.同濟(jì)大學(xué)測(cè)繪與地理信息學(xué)院，上海 200092)

0 引言

星載SAR影像廣泛應(yīng)用于圖像融合與拼接、三維重建、變化檢測(cè)等方面[1-4]，而SAR影像匹配是這類應(yīng)用的重要基礎(chǔ)。星載SAR影像匹配主要有基于強(qiáng)度的方法和基于特征的方法[5,6]，前者需要計(jì)算大量的相關(guān)性矩陣，時(shí)間復(fù)雜度高，后者包含特征提取、特征描述子構(gòu)建和相似性準(zhǔn)則匹配等過程，其速度更快、精度更高。SIFT(Scale Invariant Feature Transform)方法是具有代表性的特征匹配方法之一[7]，具有良好的尺度不變性和抗噪性[8]，已廣泛應(yīng)用于各種遙感影像匹配中[9,10]。然而，SAR影像不同于光學(xué)影像，存在相干斑噪聲，導(dǎo)致SIFT算法難以提取高質(zhì)量的特征點(diǎn)。針對(duì)這一難題，Dellinger等[11]提出了SAR-SIFT算法，通過圖像強(qiáng)度指數(shù)加權(quán)算子(ROWEA)和梯度影像對(duì)數(shù)計(jì)算以降低相干斑噪聲的影響，但由于SAR-SIFT采用CPU串行計(jì)算處理SAR衛(wèi)星影像效率較低，無法滿足地震等自然災(zāi)害發(fā)生后快速響應(yīng)場(chǎng)景之需。

針對(duì)遙感影像匹配效率不高的難題，當(dāng)前處理方法主要有4種：分塊匹配方法、特征尺度降維法、CPU并行計(jì)算、GPU并行計(jì)算。分塊匹配策略的核心是減少影像匹配中非重疊區(qū)域的冗余計(jì)算，如楊佳賓等[12]通過無人機(jī)影像的POS信息提取重疊區(qū)域以減少Dense SIFT快速拼接的冗余計(jì)算，但該類方法僅采用分塊和CPU串行的方式效率仍不高；魏天華[13]采用特征尺度降維和多核DSP方法提高了SAR-SIFT匹配效率，但該方法是通過降低金字塔維數(shù)損失匹配精度以提高效率；利用多核CPU并行計(jì)算可加倍加快處理效率，在影像匹配中應(yīng)用更為廣泛，如程棟等[14]提出利用OpenMP基于多核CPU并行計(jì)算對(duì)SAR影像匹配中的互相關(guān)模板匹配進(jìn)行加速(加速比為8)，宋偉等[15]利用OpenMP基于多核CPU并行計(jì)算對(duì)無人機(jī)影像匹配方法進(jìn)行加速(加速比優(yōu)于2)，但CPU核數(shù)有限導(dǎo)致該類方法效率提升并不明顯；GPU核數(shù)遠(yuǎn)大于CPU核數(shù)，更適用于影像匹配效率提升，如許川佩等[16]利用OpenCL基于GPU并行加速計(jì)算提高SIFT算法圖像匹配效率(加速比優(yōu)于10)，張杰等[17]提出利用CUDA基于GPU對(duì)特征提取中的圖像金字塔和梯度計(jì)算進(jìn)行加速(加速比為40)，肖漢等[18]利用基于CUDA的GPU并行計(jì)算對(duì)SIFT算法加速(加速比約為30)。

綜合上述研究可以發(fā)現(xiàn)，目前針對(duì)SAR影像SIFT匹配加速的研究較多，且取得了明顯進(jìn)展，但專門針對(duì)SAR-SIFT的并行計(jì)算方法未見報(bào)道。因此，本文提出一種基于CPU-GPU協(xié)同并行計(jì)算的SAR-SIFT快速匹配方法，以期提高SAR影像匹配速度。

1 SAR影像快速匹配并行計(jì)算方法

本文提出的SAR-SIFT快速匹配方法主要步驟為：首先對(duì)SAR影像進(jìn)行分塊預(yù)處理并提取主副影像間的重疊區(qū)域，然后采用Fork/Join模式的CPU-GPU協(xié)同并行計(jì)算以提升SAR-SIFT匹配效率。

1.1 SAR影像分塊預(yù)處理與重疊區(qū)域提取

針對(duì)星載SAR影像文件大導(dǎo)致匹配效率低問題，本文先對(duì)影像進(jìn)行預(yù)處理，求出重疊區(qū)域，再進(jìn)行后續(xù)匹配。1)分別利用主副影像有理函數(shù)模型參數(shù)(RPCs)計(jì)算影像4個(gè)角點(diǎn)的經(jīng)緯度坐標(biāo)，然后通過求交判斷該影像是否與其他影像相交；2)將每幅影像依次設(shè)為主影像并進(jìn)行分塊，求出與相交副影像的重疊區(qū)域，以減少影像分塊時(shí)不重疊區(qū)域的冗余計(jì)算；3)將主影像分成N×N塊，并求得每一子塊4個(gè)角點(diǎn)的經(jīng)緯度坐標(biāo)，進(jìn)而求得主副影像塊之間的重疊區(qū)域；4)采用RPCs反算出主副影像塊中重疊區(qū)域的行列坐標(biāo)，而后對(duì)主副影像塊中的重疊區(qū)域進(jìn)行匹配。

1.2 基于SAR-SIFT算法的連接點(diǎn)提取

SAR-SIFT算法在特征尺度空間和特征描述方面對(duì)SIFT進(jìn)行了改進(jìn)，提高了SAR影像的匹配精度，但隨著數(shù)據(jù)規(guī)模增大，計(jì)算效率逐漸降低。SAR-SIFT生成尺度空間響應(yīng)函數(shù)R(x,y,σ)計(jì)算如下[11]：

(1)

R(x,y,σ)=Det(CSH(x,y,σ))-t·Trace(CSH(x,y,σ))

(2)

首先，通過尺度空間響應(yīng)函數(shù)R(x,y,σ)在每個(gè)尺度空間上尋找3×3范圍內(nèi)的極大值進(jìn)行特征點(diǎn)提取，其中響應(yīng)函數(shù)值大于0.8。在確定特征點(diǎn)坐標(biāo)后，通過特征點(diǎn)鄰域的點(diǎn)經(jīng)過高斯加權(quán)求得特征點(diǎn)主方向，然后采用GLOH[20]方法構(gòu)建特征描述子。利用特征描述子初匹配后，采用RANSAC篩選出正確的匹配點(diǎn)對(duì)，即連接點(diǎn)。在連接點(diǎn)提取過程中，通過Fork/Join模式對(duì)金字塔影像進(jìn)行分層并行處理，該并行模式在CPU端和GPU端分別采用OpenMP和CUDA技術(shù)對(duì)SAR-SIFT匹配方法進(jìn)行加速。

1.3 基于OpenMP和CUDA的SAR-SIFT協(xié)同并行計(jì)算方法

本文提出的SAR-SIFT并行計(jì)算方法主要技術(shù)流程如圖1所示。

圖1 SAR-SIFT協(xié)同并行計(jì)算流程Fig.1 Framework of collaborative parallel computing for SAR-SIFT

(1)在獲取輸入影像后，采用GPU構(gòu)建特征金字塔，以充分發(fā)揮GPU并行處理和浮點(diǎn)運(yùn)算力的優(yōu)勢(shì)；同時(shí)利用CPU對(duì)SAR-SIFT協(xié)同并行計(jì)算，其遵循Fork/Join并發(fā)模式。在構(gòu)建8層特征金字塔時(shí)，由OpenMP編譯指令分配8個(gè)獨(dú)立線程處理對(duì)應(yīng)子任務(wù)中的響應(yīng)函數(shù)。在每個(gè)子任務(wù)中，影像從內(nèi)存?zhèn)魅腼@存中，由GPU完成SAR-HRRIS響應(yīng)函數(shù)的計(jì)算。該計(jì)算首先對(duì)影像做高斯模糊操作并計(jì)算其梯度幅值和梯度方向，然后計(jì)算像素點(diǎn)協(xié)方差矩陣，最后計(jì)算當(dāng)前層的SAR-HARRIS響應(yīng)函數(shù)。主線程等待8個(gè)工作線程處理完畢后一起將該層的SAR-HRRIS響應(yīng)函數(shù)返回CPU內(nèi)存中，該過程中僅需將影像傳入顯存，由GPU統(tǒng)一計(jì)算特征尺度，以減少顯存與內(nèi)存之間數(shù)據(jù)傳輸耗時(shí)問題。

(2)完成特征金字塔構(gòu)建后，利用OpenMP在CPU端并行提取特征點(diǎn)和構(gòu)建特征描述子。在提取特征點(diǎn)時(shí)，由OpenMP編譯指令分配8個(gè)工作線程用于提取每一層的特征點(diǎn)；在構(gòu)建特征描述子時(shí)，由OpenMP分配N個(gè)工作線程用于特征描述子構(gòu)建；匹配階段應(yīng)用knnmatch方法(基于歐氏距離匹配)[21]，采用RANSAC[22]三倍中誤差篩除錯(cuò)誤的匹配點(diǎn)。這種CPU-GPU協(xié)同并行計(jì)算方法充分利用了計(jì)算機(jī)硬件性能，可顯著提高SAR影像匹配效率。

2 實(shí)驗(yàn)與分析

2.1 實(shí)驗(yàn)數(shù)據(jù)與環(huán)境

本文實(shí)驗(yàn)數(shù)據(jù)為2021年5月26日采集的覆蓋云南省大理市漾濞縣的5景高分三號(hào)L1A級(jí)全極化條帶影像，該影像工作頻段為C頻段，方位向空間分辨率為8 m，距離向空間分辨率為6～9 m，幅寬為30 km。漾濞縣地勢(shì)北高南低，海拔最高為4 122 m，最低為1 174 m，總面積為1 954 km2，其中山區(qū)面積占98.4%。

實(shí)驗(yàn)硬件環(huán)境：CPU I-9-10900X，顯卡 RTX2080，內(nèi)存 64 GB，顯存 11 GB；軟件環(huán)境：Windows 10 64位操作系統(tǒng)、Visual Studio 2017、Opencv 3.4.11、GDAL 2.3.1、Eigen 3、CUDA Toolkit 10.1。所有實(shí)驗(yàn)在同一環(huán)境下運(yùn)行，采用C++編寫代碼，在x64debug模式下記錄程序運(yùn)行時(shí)間。

算法參數(shù)設(shè)置如下：在構(gòu)建金字塔時(shí)，金字塔層數(shù)設(shè)為8；在提取特征點(diǎn)時(shí)，SAR-HARRIS響應(yīng)函數(shù)的閾值設(shè)為0.8；在構(gòu)建特征描述子時(shí)，在提取的所有特征點(diǎn)中保留響應(yīng)值前4 000個(gè)特征點(diǎn)用于構(gòu)建特征描述子匹配，從而提高匹配速度；在匹配階段，歐氏距離匹配閾值設(shè)為0.8；最后保存每塊影像上最優(yōu)的10個(gè)連接點(diǎn)。

2.2 匹配結(jié)果分析

采用本文提出的SAR-SIFT快速匹配方法在云南漾濞震區(qū)5景高分三號(hào)順軌影像的連接點(diǎn)提取結(jié)果如圖2所示，其中重疊區(qū)域中的紫色十字符號(hào)代表相鄰兩景影像的連接點(diǎn)坐標(biāo)，在左右兩邊對(duì)應(yīng)的高分三號(hào)匹配影像中，紫線代表連接點(diǎn)連線。由圖2可知，影像重疊區(qū)域約占每景影像的20%，提取的連接點(diǎn)分布較均勻且全部落入重疊區(qū)域內(nèi)。本文采用RMSE[23]作為匹配精度指標(biāo)，得到A、B、C、D 4個(gè)區(qū)域的RMSE分別為0.3076、0.2673、0.2478、0.2569，均為亞像元精度，滿足圖像拼接需求。

圖2 高分三號(hào)遙感影像匹配結(jié)果Fig.2 Matching results of GF-3 remote sensing images

2.3 匹配效率分析

為驗(yàn)證本文提出的基于重疊區(qū)域分塊匹配策略和CPU-GPU協(xié)同并行計(jì)算在SAR-SIFT匹配方法中的有效性，設(shè)計(jì)了兩組對(duì)比實(shí)驗(yàn)：1)在分塊匹配策略中，設(shè)計(jì)了未分塊以及從4×4分塊到10×10分塊8個(gè)對(duì)比方案；在CPU-GPU協(xié)同并行計(jì)算方案中，設(shè)計(jì)了CPU單線程、CPU多線程、GPU并行、CPU-GPU協(xié)同并行4種對(duì)比方案，其中GPU并行是指在CPU串行執(zhí)行每一層金字塔構(gòu)建時(shí)采用GPU并行計(jì)算。本文重點(diǎn)關(guān)注星載SAR影像匹配效率，采用加速比作為量化參數(shù)[24]。實(shí)驗(yàn)結(jié)果如表1所示，其中耗時(shí)包括影像預(yù)處理、通訊消耗和特征點(diǎn)提取、特征描述子構(gòu)建、特征匹配所需時(shí)間，加速比1為未分塊影像匹配耗時(shí)除以不同分塊影像匹配最低耗時(shí)，加速比2為單線程影像匹配耗時(shí)除以多線程影像匹配耗時(shí)。2)將本文方法與OpenMP加速[15]、分塊匹配[25]、減層降維[13]不同分塊策略進(jìn)行對(duì)比，結(jié)果如表2所示。

由表1可知：1)影像匹配耗時(shí)隨著CPU線程數(shù)增加而遞減，當(dāng)線程數(shù)大于5時(shí)耗時(shí)減少并不明顯，主要原因是5個(gè)線程只能同時(shí)處理8層金字塔影像中的5層，導(dǎo)致后3層金字塔影像處理時(shí)出現(xiàn)2個(gè)線程等待，當(dāng)線程數(shù)增至8時(shí)，加速比為5.74；GPU并行方案可大幅提升處理效率，加速比達(dá)12.67；CPU八線程+GPU協(xié)同并行計(jì)算在所有計(jì)算方案中耗時(shí)最少，加速比達(dá)14.44。2)同一線程數(shù)的大部分情況下，6×6分塊方案的加速比最高。3)在所有計(jì)算方案中，分塊方案比未分塊方案快2倍以上。4)在分塊方案中，匹配次數(shù)隨著影像分塊數(shù)增加而增加，進(jìn)而增加耗時(shí)，而未分塊影像只需匹配一次，耗時(shí)較少；此外，GPU并行計(jì)算時(shí)分塊影像需多次傳入GPU顯存，導(dǎo)致通訊時(shí)間變長(zhǎng)、加速比降低，但分塊可減少冗余計(jì)算，使加速比下降不明顯。5)“6×6分塊+CPU八線程+GPU并行計(jì)算”為最優(yōu)方案，其相對(duì)于未分塊CPU串行計(jì)算方案的加速比達(dá)45.15。由表2中的加速比可知，本文方法在5景高分三號(hào)順軌影像匹配中效率最高，加速比為8.90，驗(yàn)證了本文方法的有效性。

表1 影像不同分塊和不同并行計(jì)算方案的實(shí)驗(yàn)結(jié)果Table 1 Experimental results of different image blocks by using various parallel computing schemes

表2 不同分塊、并行策略耗時(shí)對(duì)比Table 2 Time-consuming comparison of different image blocks and parallel strategies

3 結(jié)語

本文提出一種基于重疊區(qū)域分塊匹配策略和CPU-GPU協(xié)同計(jì)算的SAR-SIFT特征匹配方法，采用分塊匹配策略以減少非重疊區(qū)域的冗余計(jì)算，并基于Fork/Join模式CPU-GPU協(xié)同并行計(jì)算以提高SAR-SIFT特征匹配效率。為驗(yàn)證本文方法的有效性，在云南漾濞縣5景高分三號(hào)影像上開展影像匹配實(shí)驗(yàn)，并與當(dāng)前主流的分塊匹配策略、OpenMP多線程、減層降維等方法進(jìn)行對(duì)比，結(jié)果表明本文方法效率最高，加速比為8.90，為解決SAR影像匹配效率低問題提供了一種解決方案，為進(jìn)一步提升大區(qū)域SAR影像匹配性能提供了有效策略和方法。