2. 四川省应急测绘与防灾减灾工程技术研究中心, 四川 成都 610041;
3. 教育部轨道交通安全协同创新中心, 四川 成都 610031;
4. 高速铁路运营安全空间信息技术国家地方联合工程实验室, 四川 成都 610031
2. Sichuan Engineering Research Center for Emergency Mapping & Disaster Reduction, Chengdu 610041, China;
3. Collaborative Innovation Center for Rail Transport Safety, Chengdu 610031, China;
4. State-Province Joint Engineering Laboratory of Spatial Information Technology for High-speed Railway Safety, Chengdu 610031, China
随着空间技术的发展，SAR(synthetic aperture radar)成像系统得到了越来越广泛的应用。将SAR影像与光学影像进行融合处理能够更好地获取地物信息，有助于目标识别、变化检测等，其中SAR影像与光学影像匹配是十分重要的一步。因此SAR影像与光学影像匹配技术成为目前的研究热点[1, 2, 3, 4, 5]。
SAR影像存在严重的相干斑噪声干扰，并且与光学影像之间存在非线性灰度差异，导致许多匹配方法难以直接用于SAR影像与光学影像匹配。为了克服这个问题，研究人员提出了许多改进算法，其中应用最广泛的是基于局部特征的匹配方法。研究发现SIFT(scale invariant feature transform)方法[6]检测结果中的噪声点绝大多数来自高斯差分尺度空间的第1阶影像，因此可以跳过第1阶影像，从第2阶开始进行特征检测。同时，当参考影像与待匹配影像的相对旋转关系已知时，无须计算特征方向，直接为每个特征点分配一个特征方向，从而改善匹配效果[7, 8, 9, 10]。一些方法通过设置影像强度阈值来获取空间分布均匀的特征点，并通过增大特征区域来提高匹配的正确率[11]。这类方法通过增大特征区域提高了特征描述符的稳健性，但是也降低了特征描述符的可区分性。还有一些方法采用各向异性尺度空间代替高斯尺度空间抑制噪声影响，或通过多视处理平滑噪声，或采用先滤波再匹配的方法[12, 13, 14]。这些方法在抑制噪声影响的同时降低了影像的分辨率，造成影像信息丢失，不利于后续影像处理。
以上方法主要针对SAR影像相干斑噪声进行改进，没有考虑SAR影像与光学影像之间的非线性灰度差异，因而难以获得可靠的匹配结果。针对这个问题，本文从特征检测、描述和匹配3个方面出发，提出一种可靠的SAR影像与光学影像匹配方法。首先，对影像进行粗纠正，消除影像旋转和分辨率差异；其次，通过对数变换改进相位一致性方法，提取对SAR影像相干斑噪声稳健的特征点；然后，基于高斯伽玛型(Gaussian-Gamma-shaped，GGS)双边窗口算子[15]计算影像边缘强度图，在SAR影像和光学影像上获取一致性较强的目标轮廓特征，并构造梯度方向直方图模式的特征描述符；最后，联合特征描述符和几何约束条件，实现SAR影像与光学影像匹配。本文方法整体流程如图 1所示。由于本文特征检测方法通过相位一致性最大最小矩进行特征判别，对于任意像素，如果最大矩大于阈值，则判断该像素为边缘点，在这些边缘点中，如果最小矩也大于阈值，则判断该边缘点为角点，其本质是在边缘图上进行特征点检测，检测到的特征点位于影像边缘上[16]。因此，本文方法设计的“先提取特征点，然后在边缘强度图上进行特征点描述和匹配”的方案不会产生特征点位于非边缘区不利于特征描述的问题。
|
| 图 1 本文方法整体流程图 Fig. 1 Flow chart of the proposed method |
为了便于匹配，本文方法在影像精确匹配之前，采用文献[1]的方法，利用遥感影像空间参考中的地理坐标信息，建立影像之间的几何关系，自动计算参考影像4个角点在待匹配影像上的同名点，生成4对控制点，通过式(1)的透视变换模型实现影像粗纠正
式中，l1、l2、…、l8为透视变换模型参数，且
相位一致性特征检测方法加入了噪声抑制过程，对影像噪声具有较强的稳健性[16]，但原始的相位一致性特征检测方法假设不同影像区域的噪声响应值相等，因此，其噪声抑制只适用于加性噪声。对于SAR影像上的乘性噪声，相位一致性方法的噪声抑制效果不足。因为对SAR影像进行对数变换以后影像噪声近似服从高斯分布[17]，本文利用相位一致性方法进行特征检测时，首先将输入影像进行对数变换，然后计算特征响应值，并通过逆变换在原始影像上提取特征点，具体过程如下。
(1) 对输入影像进行对数变换，将影像乘性相干斑噪声转换为高斯加性噪声。
(2) 在对数变换后的影像上计算每个像素点的相位一致性度量值PCln[16, 18]
式中，An(x)为信号上点x处第n个傅里叶分量的幅值；ϕ(x)为相位角；ϕ(x)为整体平均相位角；W(x)为频率传播的权重；ε是常数，其目的是避免分母为0；T为影像噪声的一个估计量；符号“$\left\lfloor {} \right\rfloor $”表示当其内部计算结果大于0时，输出结果为内部计算值本身，否则输出结果为0。
(3) 计算Iln(x,y)上每个像素点的相位一致性最小矩mln(x,y)
式中 式中，PCln(x,y,θ)表示方向θ上的相位一致性度量值，求和运算是计算各个方向(6个方向)的相位一致性累加值。
(4) 对mln(x,y)进行对数变换的逆变换，得到逆变换结果m(x,y)。
(5) 对m(x,y)进行归一化处理，将相位一致性最小矩归一化到[0, 1]区间内
(6) 对每个像素的归一化相位一致性最小矩$\tilde m$(x,y)进行非极大值抑制，并设置最小矩阈值进行特征点判别：如果归一化相位一致性最小矩为局部极大值，且大于阈值，则认为该像素点为特征点。
1.3 基于Gaussian-Gamma-shaped边缘强度图的特征描述由于1.1节中影像粗纠正整体上消除了参考影像与待匹配影像之间的尺度和旋转变化，因此在特征描述时可以对影像上所有特征设置统一的特征区域大小和特征方向。假设影像分辨率为λ，设置特征区域半径为r个像素
特征描述时，由于SAR影像噪声干扰以及SAR影像与光学影像之间灰度差异较大，如图 2所示，使得SIFT这类基于梯度的方法往往不能获得稳定的匹配结果。
|
| 图 2 SAR影像与光学影像示例 Fig. 2 Example of SAR and optical images |
从图 2中可以看到，SAR影像与光学影像地物轮廓相似性较高。因此，本文先提取影像边缘信息，再基于边缘信息进行特征点描述和匹配。SAR影像边缘信息提取方法主要有两大类：一类是先对影像进行滤波处理，再采用传统光学影像边缘检测方法进行边缘提取，常用的滤波方法有Lee滤波算法[19]、Gamma-Map滤波算法[20](GM)、双边滤波算法[21](BF)、非局部均值滤波算法[22](NLm)、PPB[23]算法等；另一类是直接在SAR影像进行顾及噪声的边缘提取，常用的方法有矩形双边窗口(Rect)比值算子[24, 25, 26]和GGS双边窗口比值算子[15]。本文采用如图 3(a)所示的模拟SAR影像对上述两类方法进行测试。对于第1类方法，分别采用Lee、GM、BF、NLm和PPB滤波器对SAR影像进行滤波，滤波后的影像均采用Canny算子进行边缘检测。对于第2类方法，分别采用基于Rect窗口和GGS窗口的比值算子进行边缘检测。测试结果如图 3(b)—图 3(h)所示。
|
| 图 3 边缘检测结果 Fig. 3 Edge detection results |
从图 3所示结果可以看出，基于Lee滤波器、GM滤波器以及BF滤波器滤波后影像的边缘检测结果中包含大量的噪声点，因此，这3种方法不适合用于SAR影像边缘信息检测。后面4种方法的检测结果优于前面3种方法。其中，相对于GGS比值算子而言，NLm滤波器、PPB滤波器和Rect比值算子检测结果中存在更多错误边缘点。另外，NLm滤波和PPB滤波的时间效率非常低。综上所述，GGS比值算子在SAR影像边缘检测中表现出最优的性能。因此，本文采用GGS比值算子计算参考影像和待匹配影像的边缘强度图，提取边缘信息，具体步骤如下[15]。
(1) 计算每个像素双边GGS窗口中两个区域的局部均值
式中，θ=0、π${1 \over P}$、π${2 \over P}$、…、π${{P - 1} \over P}$是为了能够检测出多个方向的边缘特征将双边GGS窗口进行旋转的方向值；P为划分的方向个数；WUθ、WLθ是方向为θ的窗口函数
水平方向的双边窗口函数WU和WL定义如下
式中，参数σx>1控制窗口长度，参数α>1控制窗口的宽度，参数β>0控制双边窗口之间的间隔距离。
(2) 计算每个像素的最小均值比值R(x,y)
式中，m1(x,y|θp)和m2(x,y|θp)为θp方向上GGS窗口中两个区域的局部均值。
(3) 计算每个像素的边缘强度响应值ESM(x,y)=1－R(x,y)。
图 2所示的SAR影像与光学影像的GGS边缘强度图如图 4所示。从图 4可以看出，虽然SAR影像受到严重的噪声干扰并且与光学影像的灰度差异非常大，但是它们的GGS边缘强度图在整体上具有很强的相似度。得到边缘强度图以后，本文在边缘强度图上利用梯度方向直方图(histogram of oriented gradient，HOG)结构计算特征描述符，流程如图 5所示。
|
| 图 4 GGS边缘强度图 Fig. 4 GGS edge strength map |
|
| 图 5 特征描述符构造流程图 Fig. 5 Flow chart of feature descriptor construction |
由于影像边缘强度图纹理不丰富，基于边缘强度图的特征描述符中最近邻的点可能并非同名点。为了提高特征匹配率，本文采用n-NN的方法，对参考影像上的每个特征点选择n对候选同名特征。在候选同名特征集合中，基于特征描述符相似性和特征几何约束条件选择种子点，并以种子点为控制基础，计算最终的匹配结果，具体过程如下。
(1) 对参考影像上的H个特征点，计算候选匹配特征集合中所有n×H对初始匹配的特征描述符距离，并选择其中描述符距离最近且满足影像x和y方向上坐标偏差小于阈值(dx,dy)的K对特征作为种子匹配对。
(2) 对每一对种子匹配Mi(Fssar,Fsopt),(i=1,2,…,K)，设置对应的匹配集合setAlli为空集，并将该对种子匹配归入匹配集合setAlli中，setAlli中的特征对标记为(Fasar,Faopt)i。
(3) 对于候选匹配特征集合中任意一对特征(Fjsar,Fjopt)，首先对特征点之间的影像坐标偏差进行阈值判断。如果大于阈值(dx,dy)，则认为该对特征不是同名特征，进入下一对特征的判别，否则，计算该对特征与setAlli中每对特征的距离差值和方向差值，如果(Fjsar,Fjopt)与setAlli中超过95%的特征对之间的距离和方向关系满足式(11)，则将该特征对归入setAlli中
式中，|·|表示欧氏距离函数；ts和tθ分别为距离比值阈值和方向阈值。
(4) 完成候选特征集合中所有特征对的遍历后，保存setAlli。
(5) 选择其他种子匹配迭代进行步骤(2)至步骤(4)，直至完成K次迭代运算。迭代结束以后，setAlli,(i=1,2，…,K)中特征对数最多的一组结果即为几何约束后的匹配结果。
经过以上操作以后，初始候选特征匹配集合中大量的错误匹配特征将被剔除，再基于式(1)的透视变换模型采用RANSAC方法进行错误匹配剔除，获得最终的匹配结果。
2 试验及结果分析将本文方法分别与SIFT[6]、Schwind[7]、Suri[8]、Fan[10]方法进行对比分析以验证本文方法的有效性。在本文试验中，对数变换相位一致性特征检测中最小矩阈值设置为0.3；GGS边缘强度图计算中双边窗口长度σx=3，宽度α=2，窗口间隔β=1，每个像素点窗口所取的方向个数P=8；基于几何约束的匹配方法中对于参考影像上的每个特征点所选择的候选特征个数H=20，种子匹配特征对数K=10，坐标距离阈值dx=10，dy=10，距离比值阈值ts=0.2，方向阈值tθ=5°。
2.1 试验数据本文试验数据如图 6和图 7所示，其中SAR影像分别为COSMO-Skymed(CS)影像(X波段、HH极化)、RADARSAT-2(R2)影像(C波段、HH极化)、TerraSAR-X(TX)影像(X波段、HH极化)，光学影像为资源三号全色影像(ZY-3)。
|
| 图 6 试验数据1 Fig. 6 Experimental dataset 1 |
|
| 图 7 试验数据2 Fig. 7 Experimental dataset 2 |
试验以正确匹配特征个数作为各种方法性能的评价指标，统计结果见表 1。
| 试验影像像对 | SIFT | Suri | Schwind | Fan | 本文方法 | |
| SAR&SAR | 图 6(a)&(b) | 3 | 5 | 4 | 21 | 44 |
| 图 6(b)&(c) | 146 | 202 | 418 | 1457 | 426 | |
| 图 7(b)&(c) | 1 | 3 | 7 | 8 | 25 | |
| SAR&光学 | 图 6(a)&(d) | 1 | 1 | 3 | 2 | 93 |
| 图 6(b)&(d) | 1 | 1 | 3 | 12 | 99 | |
| 图 6(c)&(d) | 1 | 2 | 4 | 8 | 143 | |
| 图 7(a)&(d) | 0 | 0 | 2 | 7 | 18 | |
| 图 7(b)&(d) | 0 | 0 | 2 | 17 | 26 | |
| 图 7(c)&(d) | 2 | 1 | 4 | 1 | 25 | |
从表 1试验结果可以看出，对于SAR影像与SAR影像匹配，除了在第2对影像上Fan的方法获得的正确匹配特征个数最多以外，另外两对影像上都是本文方法获得了最多的正确匹配特征。所有的方法在第2对影像上都能获得较多的正确匹配特征，这主要是因为第2对影像的影像分辨率更高，在参考影像和待匹配影像上检测出了更多的特征点，即在更多的特征中匹配出了更多的同名特征。另外，虽然其他方法表现的性能不如本文方法，但是仍然都能获得一定数量的同名特征，主要是因为在参考影像与待匹配影像上均存在大量的相干斑噪声，参考影像上的噪声点有可能在待匹配影像上对应的位置也存在噪声点，即在匹配结果中可能存在同名噪声点。同时，由于参考影像与待匹配影像均为SAR影像，影像上相同地物目标之间的灰度值差异较小，因此其他几种匹配方法都能够获得一定数量的正确匹配特征。
对于SAR影像与光学影像匹配，在所有试验中本文方法都获得了最多的正确匹配特征。其他匹配方法在多对影像上都匹配失败，主要原因是：
(1) 在SAR影像上检测出的大量噪声点在光学影像上不存在同名噪声点，这些噪声点对匹配过程造成极大的干扰，导致出现大量的错误匹配，正确匹配特征个数大幅减少。
(2) 参考影像与待匹配影像灰度值差异变化非常大，甚至灰度值相对大小发生逆转，这种非单调的灰度值变化直接导致SIFT描述符中梯度计算结果发生变化，导致本组试验中其他几种匹配方法失败。在本文方法中，SAR影像与光学影像的GGS边缘强度图之间具有较强的相似度，因此在GGS边缘强度图的基础上计算特征描述符能够获得较好的匹配结果。
本文方法对SAR影像与光学影像的匹配结果如图 8所示。
|
| 图 8 本文方法匹配结果 Fig. 8 Matching results of the proposed method |
本文在影像匹配结果的基础上进行影像配准，并人工选择控制点计算均方根误差(RMSE)来衡量配准精度。根据表 1所示匹配结果，在图 6(b)和图 6(c)所示影像上所有匹配方法都能够获得4对以上同名点，因此本文选择该对影像作为试验数据，配准精度见表 2。从表 2所示结果可以看出，基于本文匹配方法进行影像配准的精度最高，坐标值偏差d小于1个像素。配准结果如图 9所示。
| 像素 | |||||
| 试验影像像对 | SIFT | Suri | Schwind | Fan | 本文方法 |
| 图 6(b)&(c) | 1.58 | 1.47 | 1.37 | 1.61 | 0.98 |
|
| 图 9 基于本文匹配方法进行影像配准的结果 Fig. 9 Image registration based on the matching result of the proposed method |
本文针对SAR影像相干斑噪声及其与光学影像的灰度差异导致现有匹配方法难以获得可靠匹配结果的问题，提出了一种基于GGS双边窗口边缘强度图的匹配方法。本文方法的创新之处在于利用GGS双边窗口算子对SAR影像和光学影像能够获得相似影像轮廓的特性，构建了相似性较强的特征描述符，同时，基于遥感影像地理空间参考信息对影像进行粗纠正，并设计几何约束条件与特征描述符相似性相结合进行影像匹配。以上改进使得本文方法能够较好地处理SAR影像噪声干扰和SAR影像与光学影像的灰度差异问题，对于SAR影像与SAR影像之间的匹配以及SAR影像与光学影像之间的匹配都能够获得较好的匹配结果。但是，本文方法基于影像先验信息进行粗纠正，在地形起伏严重的区域匹配效果不理想。后续研究工作将尝试解决地形起伏区域的SAR影像与光学影像匹配问题。
| [1] | 岳春宇, 江万寿. 几何约束和改进SIFT的SAR影像和光学影像自动配准方法[J]. 测绘学报, 2012, 41(4):570-576. YUE Chunyu, JIANG Wanshou. An Automatic Registration Algorithm for SAR and Optical Images Based on Geometry Constraint and Improved SIFT[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(4):570-576. |
| [2] | 李雨谦, 皮亦鸣, 王金峰. 基于水平集的SAR图像与光学图像的配准[J]. 测绘学报, 2010, 39(3):276-282. LI Yuqian, PI Yiming, WANG Jinfeng. The Registration between SAR and Optical Image Based on Level Set[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(3):276-282. |
| [3] | 史磊, 李平湘, 杨杰. 利用SIFT与粗差探测进行SAR影像配准[J]. 武汉大学学报(信息科学版), 2010, 35(11):1296-1299. SHI Lei, LI Pingxiang, YANG Jie. SAR Imagery Registration Based on SIFT and Data Snooping[J]. Geomatics and Information Science of Wuhan University, 2010, 35(11):1296-1299. |
| [4] | 王峰, 尤红建, 傅兴玉. 应用于SAR图像配准的自适应SIFT特征均匀分布算法[J]. 武汉大学学报(信息科学版), 2015, 40(2):159-163. WANG Feng, YOU Hongjian, FU Xingyu.Auto-adaptive Well-distributed Scale-invariant Feature for SAR Images Registration[J]. Geomatics and Information Science of Wuhan University, 2015, 40(2):159-163. |
| [5] | 李锡军, 关泽群, 盛文彬. 一种基于共性特征的SAR图像与卫星光学图像匹配方法[J]. 湘潭大学自然科学学报, 2010, 32(1):112-116. LI Xijun, GUAN Zequn, SHENG Wenbin. SAR Image and Satellite Optical Image Matching Method Based on Common Characteristics[J]. Natural Science Journal of Xiangtan University, 2010, 32(1):112-116. |
| [6] | LOWE D G. Distinctive Image Features from Scale-invariant Keypoints[J]. International Journal of Computer Vision, 2004, 60(2):91-110. |
| [7] | SCHWIND P, SURI S, REINARTZ P,et al. Applicability of the SIFT Operator to Geometric SAR Image Registration[J]. International Journal of Remote Sensing, 2010, 31(8):1959-1980. |
| [8] | SURI S, SCHWIND P, UHL J,et al. Modifications in the SIFT Operator for Effective SAR Image Matching[J]. International Journal of Image and Data Fusion, 2010, 1(3):243-256. |
| [9] | FAN Bin, WU Fuchao, Hu Zhanyi. Aggregating Gradient Distributions into Intensity Orders:A Novel Local Image Descriptor[C]//Proceedings of IEEE Conference on Computer Vision and Pattern Recognition.Providence:IEEE,2011:2377-2384. |
| [10] | FAN Bin,HUO Chunlei,PAN Chunhong,et al.Registration of Optical and SAR Satellite Images by Exploring the Spatial Relationship of the Improved SIFT[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 10(4):657-661. |
| [11] | LIU Lining,WANG Yunhong,WANG Yiding.SIFT Based Automatic Tie-point Extraction for Multitemporal SAR Images[C]//International Workshop on Education Technology and Training and International Workshop on Geoscience and Remote Sensing.Shanghai:IEEE, 2008:499-503. |
| [12] | WESSEL B, HUBER M, ROTH A. Registration of Near Real-time SAR Images by Image-to-image Matching[C]//Proceedings of Photogrammetric Image Analysis:PIA07.Munich:[s.n.], 2007:179-184. |
| [13] | YU Xiaoping,LIU Tong,LI Pingxiang,et al.The Application of Improved SIFT Algorithm in High Resolution SAR Image Matching in Mountain Areas[C]//Proceedings of International Symposium on Image and Data Fusion.Tengchong:IEEE,2011:1-4. |
| [14] | WANG Shanhu, YOU Hongjian, FU Kun. BFSIFT:A Novel Method to Find Feature Matches for SAR Image Registration[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(4):649-653. |
| [15] | SHUI Penglang, CHENG Dong. Edge Detector of SAR Images Using Gaussian-Gamma-shaped Bi-windows[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(5):846-850. |
| [16] | KOVESI P. Phase Congruency Detects Corners and Edges[C]//Proceedings of the 7th Digital Image Computing:Techniques and Applications. Sydney:[s.n.], 2003:309-318. |
| [17] | XIE Hua, PIERCE L E, ULABY F T. Statistical Properties of Logarithmically Transformed Speckle[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(3):721-727. |
| [18] | 陈敏, 朱庆, 朱军, 等. 多光谱遥感影像亮度空间相位一致性特征点检测[J]. 测绘学报, 2016, 45(2):178-185. DOI:10.11947/j.AGCS.2016.20150030. CHEN Min, ZHU Qing, ZHU Jun, et al. Interest Point Detection for Multispectral Remote Sensing Image Using Phase Congruency in Illumination Space[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(2):178-185. DOI:10.11947/j.AGCS.2016.20150030. |
| [19] | LEE J S. Speckle Suppression and Analysis for Synthetic Aperture Radar Images[J]. Optical Engineering, 1986, 25(5):636-643. |
| [20] | ARSENAULT H H, APRIL G. Properties of Speckle Integrated with a Finite Aperture and Logarithmically Transformed[J]. Journal of the Optical Society of America, 1976, 66(11):1160-1163. |
| [21] | PARIS S, DURAND F. A Fast Approximation of the Bilateral Filter Using a Signal Processing Approach[C]//Proceedings of the 9th European Conference on Computer Vision. Graz:Springer,2006:568-580. |
| [22] | BUADES A, COLL B, MOREL J M. A Review of Image Denoising Algorithms, with a New One[J]. Multiscale Modeling & Simulation, 2006, 4(2):490-530. |
| [23] | DELEDALLE C A,DENIS L,TUPIN F.Iterative Weighted Maximum Likelihood Denoising with Probabilistic Patch-based Weights[J]. IEEE Transactions on Image Processing, 2009, 18(12):2661-2672. |
| [24] | OLIVER C J, BLACKNELL D, WHITE R G. Optimum Edge Detection in SAR[J]. IEEE Proceedings:Radar, Sonar and Navigation, 1996, 143(1):31-40. |
| [25] | GERMAIN O, REFREGIER P. Edge Location in SAR Images:Performance of the Likelihood Ratio Filter and Accuracy Improvement with an Active Contour Approach[J]. IEEE Transactions on Image Processing, 2001, 10(1):72-78. |
| [26] | SCHOU J, SKRIVER H, NIELSEN A A, et al. CFAR Edge Detector for Polarimetric SAR Images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(1):20-32. |