基于方向相位稠密特征的多传感器遥感影像配准方法和系统

基于方向相位稠密特征的多传感器遥感影像配准方法和系统 A multi-sensor remote sensing registration method and system based on dense feature of orientated phase first-author 叶沅鑫，1985年生，男，副教授，研究方向为遥感图像处理与分析。E-mail： yeyuanxin110@163.com 叶沅鑫，1985年生，男，副教授，研究方向为遥感图像处理与分析。E-mail： yeyuanxin110@163.com yeyuanxin110@163.com 西南交通大学 地球科学与环境工程学院, 成都 611756 西南交通大学 地球科学与环境工程学院, 成都 611756 Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, Chengdu 611756, China Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, Chengdu 611756, China 针对多传感器遥感影像间显著几何畸变和灰度差异造成的配准困难问题，本文提出了一种基于结构相似性的快速精确配准方法。为了构建多传感器影像间的稳健结构特征，引入光照和对比度不变性的相位一致性模型，利用相位一致性特征值和特征方向构建了一种逐像素的三维结构特征描述符—方向相位稠密特征（Dense Feature of Orientated Phase， DFOP），该特征描述符通过捕捉影像的几何结构分布，能够有效抵抗多传感器影像间的灰度差异。然后，基于模板匹配的策略，将DFOP描述符变换到频率域，并进行快速子像素精度匹配。此外，在DFOP的基础上，本文研制了一种快速鲁棒的多传感器遥感影像配准系统。通过利用多种地貌类型的多传感器遥感影像对所提出的方法进行测试，实验结果表明，相较于其它基于灰度或结构特征的匹配方法，DFOP能够获得更高的匹配正确率，且所开发的配准系统（下载链接1）也优于商业软件ENVI以及ERDAS中配准模块的配准精度。 To solve the problem of registration difficulty caused by significant geometric distortion and gray differences between multi-sensor remote sensing images, this paper proposes a fast and accuracy registration method based on structural similarity between images, where the phase congruency model with illumination and contrast invariances is introduced to construct the robust structural features descriptors of images. Firstly, the intensity and orientation of phase congruency are used to build a pixel-wise three-dimensional structural feature representation named DFOP (Dense Feature of Orientated Phase), which can effectively resist the grayscale difference between multi-sensor images by capturing geometric structures of images. Next, the DFOP feature descriptor is transformed into the frequency domain, and the single-step DFT approach is used to achieve the fast matching with a subpixel accuracy by employing a template matching scheme. In addition, we also developed a fast and robust automatic multi-sensor remote sensing image registration system based on proposed DFOP. Finally, the proposed method and registration system has been validated using multiple pairs of multi-sensor remote sensing images (including optical, LIDAR and SAR) covering different scenes, the results show that the proposed DFOP achieves higher correct matching rate, and the developed registration system outperforms the registration module of ENVI and ERDAS in registration accuracy. Our system is available at https://github.com/yeyuanxin110/Remote-Sensing-Image-Registration-system.git https://github.com/yeyuanxin110/Remote-Sensing-Image-Registration-system.git 多传感器遥感影像 影像配准 相位一致性 DFOP 影像配准系统 multi-sensor remote sensing images image registration phase congruency dense feature of orientated phase images registration system 1　引言 1 航空航天和遥感成像技术的快速发展为人们提供了不同传感器的遥感影像，而对这些多传感器遥感影像集成应用的基础是高精度的几何配准。影像配准是将从不同传感器、不同时间或不同视角下获得的同一区域的多张影像统一到同一坐标系下并进行对齐和叠加的过程（ Zitová和Flusser， 2003 Zitová和Flusser， 2003 ）。作为遥感影像应用中重要的预处理步骤，影像的配准精度会直接影响到后续应用，如变化检测（ Bruzzone和Bovolo， 2013 Bruzzone和Bovolo， 2013 ）、影像融合（ Schmitt 等， 2017 Schmitt 等， 2017 ）以及影像拼接（LI等， 2015）等。由于多传感器遥感影像间成像机理的不同所引起的几何畸变和非线性灰度差异，如何实现多传感器遥感影像间的快速精确配准一直是学术界和工业界研究的难点。 遥感影像配准方法主要可分为基于特征的方法和基于区域的方法。基于特征方法的主要步骤包括特征提取、特征描述、特征匹配和几何校正，即首先需要在影像间提取具有高重复率的显著稳健特征，之后利用特征之间的相关性确定最佳匹配并建立变换关系。常用的特征包括点特征（ Huo等， 2012 Huo等， 2012 ）、线特征（ Li和Leung， 2007 Li和Leung， 2007 ）和面特征（ Goncalves等， 2011 Goncalves等， 2011 ）。具有尺度和旋转不变性的尺度不变特征变换（Scale Invariant Feature Transform， SIFT）（ Lowe， 2004 Lowe， 2004 ）被广泛应用于遥感影像的配准，SIFT虽能抵抗影像间的旋转和尺度差异，但是由于SIFT是基于像素梯度信息进行特征表达，在多传感器影像配准中难以获得稳定的同名点。为了克服多传感器影像间的非线性灰度差异， Dellinger等（2015） Dellinger等（2015） 通过采用比值梯度代替灰度差梯度对SIFT进行改进，提出了SAR-SIFT，在一定程度上减小了斑点噪声的影响。 Xiang等（2018） Xiang等（2018） 通过使用多尺度指数加权均值比率（multiscale ratio of exponentially weighted averages， ROEWA）和多尺度Sobel分别在SAR和光学影像上计算影像的梯度信息，并在尺度空间检测Harris特征点，提出了OS-SIFT算法。虽然上述方法对SIFT算法做出了一定程度上的改进，但对具有非线性灰度差异的多传感器遥感影像并不具有较好的稳健性。 基于区域的方法是通过使用固定大小的模板窗口，然后以某种相似性测度为准则确定最优匹配，并以模板窗口的中心点像素作为同名点。常用的相似性测度有灰度差平方和（Sum of Squared Differences， SSD）（ Hisham等， 2015 Hisham等， 2015 ）、归一化互相关（Normalized Cross Correlation， NCC）（ Martinez等， 1999 Martinez等， 1999 ）以及互信息（Mutual Information， MI）（ Suri和Reinartz， 2010 Suri和Reinartz， 2010 ）等。由于这些相似性测度对非线性灰度变化较为敏感，在多传感器遥感影像配准中表现不够理想。如 Lu等（2020） Lu等（2020） 使用结构张量描述影像的几何结构属性，并构建特征描述符进行匹配。 Ye等（2017） Ye等（2017） 使用具有光照和对比度不变性的相位一致性构建了相位一致性方向直方图（Histogram of Orientated Phase Congruency， HOPC），但由于HOPC采用稀疏采样的方式构建特征描述符，难以非常精确的刻画影像的结构细节，其匹配精度会受到一定程度的影响。鉴于此， Ye等（2019） Ye等（2019） 构建了一种逐像素表达的稠密结构特征描述符—方向梯度通道特征（Channel Features of Orientated Gradients， CFOG），CFOG通过逐像素的结构特征表达提高了匹配精度，但由于CFOG采用的是影像的梯度信息，而梯度信息容易受影像间噪声和非线性灰度差异的影响。 随着深度学习方法在计算机视觉任务的巨大进步，国内外学者将深度学习引入到了影像配准。 Nan等（2019） Nan等（2019） 提出一种对Siamese网络进行优化，提取多模态影像共有特征的遥感影像匹配方法。 Zhu等（2019） Zhu等（2019） 提出了一种基于双分支深度学习网络的特征学习方法，将图像匹配问题转化为两类分类问题。 Lan等（2021） Lan等（2021） 通过引入D2-Net特征提取的思想，提出了一种基于深度学习特征的匹配方法CMM-Net。虽然利用深度学习进行影像配准逐渐成为了研究热点，但基于深度学习的方法需要大量数据进行训练，因而限制了它们在工程实践中的广泛应用。更多遥感影像配准方法请参考 Feng等（2021） Feng等（2021） 发表的最新综述性论文。 综上所述，为解决多传感器遥感影像间因几何畸变和非线性灰度差异造成配准困难的问题，本文借助于几何结构特征能够有效刻画多传感器影像间共有属性的性质，首先利用具有光照和对比度不变性的相位一致性模型构建了一种逐像素表达的三维结构特征描述符——方向相位稠密特征（Dense Feature of Orientated Phase， DFOP），然后采用模板匹配的策略设计了一种基于DFOP的子像素自动配准方法，实现了多传感器遥感影像间的高精度配准。此外，本文基于所设计的自动配准方法开发了一套多传感器遥感影像自动配准系统，在利用DFOP进行匹配的同时，引入遥感影像的地理信息或有理多项式参数（rational polynomial coefficient， RPC）进行几何空间约束，旨在解决多传感器影像间的几何畸变和灰度差异造成的配准困难问题。本文主要贡献如下： 1：构建了一种基于方向相位特征的三维结构特征描述符DFOP，该描述符能够有效抵抗多传感器遥感影像间的非线性灰度差异，并在此基础上设计了一种快速的子像素匹配方法； 2：研制了一套多传感器遥感影像自动配准系统，该系统利用遥感影像自身的地理空间信息进行几何粗配准，然后采用DFOP进行精确配准，可有效实现多传感器遥感影像间的快速精确配准。 2　基于DFOP描述符的快速匹配 1 2.1　DFOP描述符 2 DFOP特征描述符是受CFOG的启发提出的。CFOG是采用梯度方向和梯度幅值构建的三维特征描述符（ Ye 等， 2019 Ye 等， 2019 ），而由于CFOG特征描述符是采用逐像素的方式表达影像的几何结构信息，在一定程度上能够抵抗影像间的光照和对比度变化，但CFOG采用的是影像的梯度信息，易受多传感器遥感影像间非线性灰度差异的影响。因此，本文利用对影像间非线性灰度差异抵抗性更强的相位一致性特征值和特征方向构建DFOP特征描述符，具体步骤如下： （1） 计算相位一致性特征值： Kovesi等利用Log-Gabor滤波改进了相位一致性的计算模型（ Kovesi， 1999 Kovesi， 1999 ），即对Log-Gabor进行傅里叶逆变换获得其在空间域中的偶对称滤波器 <math id="M1"><msub><mrow><mi>L</mi></mrow><mrow><mi>e</mi><mi>v</mi><mi>e</mi><mi>n</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824376&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824375&type= 5.33400011 3.72533321 和奇对称滤波器 <math id="M2"><msub><mrow><mi>L</mi></mrow><mrow><mi>o</mi><mi>d</mi><mi>d</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824426&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824425&type= 4.91066647 3.72533321 ，通过将影像 <math id="M3"><mi>I</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824424&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824423&type= 10.15999985 3.64066648 在不同尺度 n 和方向 o 下分别与偶对称滤波器和奇对称滤波器进行卷积运算，获得对应的响应分量 <math id="M4"><msub><mrow><mi>e</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824382&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824381&type= 12.44600010 4.74133301 和 <math id="M5"><msub><mrow><mi>r</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824384&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824383&type= 12.27666664 4.74133301 。 <math display="block" id="M6"><mo stretchy="false">[</mo><msub><mrow><mi>e</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo>,</mo><msub><mrow><mi>r</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mo>=</mo><mo stretchy="false">[</mo><mi>I</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mi mathvariant="normal">*</mi><msubsup><mrow><mi>L</mi></mrow><mrow><mi>e</mi><mi>v</mi><mi>e</mi><mi>n</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msubsup><mo>,</mo><mi>I</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mi mathvariant="normal">*</mi><msubsup><mrow><mi>L</mi></mrow><mrow><mi>o</mi><mi>d</mi><mi>d</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msubsup><mo stretchy="false">]</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824386&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824385&type= 67.64866638 5.07999992 式中，*表示卷积运算。 影像 <math id="M7"><mi>I</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824424&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824423&type= 10.15999985 3.64066648 在尺度 n 和方向 o 处的振幅 <math id="M8"><msub><mrow><mi>A</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824390&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824389&type= 13.03866673 4.74133301 和相位 <math id="M9"><msub><mrow><mi>ϕ</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824416&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824415&type= 13.37733269 4.74133301 可以通过以下公式获得： <math display="block" id="M10"><msub><mrow><mi>A</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo>=</mo><mroot><mrow><msub><mrow><mi>e</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><msup><mrow><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><mo>+</mo><msub><mrow><mi>r</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><msup><mrow><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mrow><mrow/></mroot></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824394&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824393&type= 49.44533539 5.67266655 <math display="block" id="M11"><msub><mrow><mi>ϕ</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo>=</mo><mi mathvariant="normal">a</mi><mi mathvariant="normal">r</mi><mi mathvariant="normal">c</mi><mi mathvariant="normal">t</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">n</mi><mfenced separators="|"><mrow><mfrac><mrow><msub><mrow><mi>r</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mrow><mi>e</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></mrow></mfrac></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824396&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824395&type= 41.74066544 10.32933331 之后通过对影像在多个方向、多尺度下的相位值进行整合，获得影像的相位一致性，其计算表达式为： <math display="block" id="M12"><mi>P</mi><mi>C</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><munder><mo movablelimits="false">∑</mo><mrow><mi>n</mi></mrow></munder><mrow><munder><mo movablelimits="false">∑</mo><mrow><mi>o</mi></mrow></munder><mrow><msub><mrow><mi>W</mi></mrow><mrow><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mfenced open="⌊" close="⌋" separators="|"><mrow><msub><mrow><mi>A</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mi mathvariant="normal">Δ</mi><msub><mrow><mi>ϕ</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo>-</mo><mi>T</mi></mrow></mfenced></mrow></mrow></mrow><mrow><munder><mo movablelimits="false">∑</mo><mrow><mi>n</mi></mrow></munder><mrow><munder><mo movablelimits="false">∑</mo><mrow><mi>o</mi></mrow></munder><mrow><msub><mrow><mi>A</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo>+</mo><mi>ξ</mi></mrow></mrow></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824398&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824397&type= 75.26866913 14.64733315 式中， <math id="M13"><mi>P</mi><mi>C</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824400&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824399&type= 13.29266548 3.64066648 表示点 <math id="M14"><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824537&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824536&type= 8.21266651 3.64066648 处的相位一致性值， <math id="M15"><msub><mrow><mi>W</mi></mrow><mrow><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824404&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824403&type= 12.95400047 4.74133301 为给定频率扩展的权重因子， T 为噪声补偿值， <math id="M16"><mi>ξ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824407&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824405&type= 1.69333339 3.30200005 是一个避免分母为零的常数， <math id="M17"><mfenced open="⌊" close="⌋" separators="|"><mrow/></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824409&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824408&type= 4.14866638 3.72533321 表示取正运算符。 <math id="M18"><msub><mrow><mi>ϕ</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824416&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824415&type= 13.37733269 4.74133301 是一个相位偏差函数，其表达式定义为： <math display="block" id="M19"><mtable columnalign="left"><mtr><mtd><mi mathvariant="normal">Δ</mi><msub><mrow><mi>ϕ</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo>=</mo><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mo stretchy="false">(</mo><msub><mrow><mi>ϕ</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo>-</mo><msub><mrow><mover><mrow><mi>ϕ</mi></mrow><mo>¯</mo></mover></mrow><mrow><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo>-</mo></mtd></mtr><mtr><mtd><mfenced open="|" close="|" separators="|"><mrow><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mo stretchy="false">(</mo><msub><mrow><mi>ϕ</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo>-</mo><mover><mrow><msub><mrow><mi>ϕ</mi></mrow><mrow><mi>o</mi></mrow></msub></mrow><mo>¯</mo></mover><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow></mfenced></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824414&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824413&type= 59.18199921 11.93799973 式中， <math id="M20"><msub><mrow><mi>ϕ</mi></mrow><mrow><mi>n</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824416&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824415&type= 13.37733269 4.74133301 是影像在尺度 n 和方向 o 下的相位值， <math id="M21"><msub><mrow><mover><mrow><mi>ϕ</mi></mrow><mo>¯</mo></mover></mrow><mrow><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824418&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824417&type= 12.78466606 4.91066647 表示对应方向所有尺度上的加权均值相位。 （2） 计算相位一致性特征方向： 尽管利用相位一致性特征值可以获得影像的结构边缘信息，但其只考虑了像素的特征值，而忽视了其对应的特征方向信息，因而对影像的区域几何结构特征描述十分有限。类似于梯度方向，相位一致性方向表示影像特征变化最快的方向，在构建特征描述符中具有重要的作用。 本文根据奇对称滤波器 <math id="M22"><msub><mrow><mi>L</mi></mrow><mrow><mi>o</mi><mi>d</mi><mi>d</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824426&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824425&type= 4.91066647 3.72533321 的性质，即Log-Gabor函数傅里叶逆变换对应的奇对称滤波器近似于一种导数差分模板，影像 <math id="M23"><mi>I</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824424&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824423&type= 10.15999985 3.64066648 与奇对称滤波器的卷积结果表示了影像在对应方向上的能量变化。根据 式（4） 式（4） ，影像的相位一致性使用了多个方向的奇对称滤波器，因此这里将影像 <math id="M24"><mi>I</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824424&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824423&type= 10.15999985 3.64066648 与奇对称滤波器 <math id="M25"><msub><mrow><mi>L</mi></mrow><mrow><mi>o</mi><mi>d</mi><mi>d</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824426&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824425&type= 4.91066647 3.72533321 在多个方向的卷积结果分别投影到水平方向和垂直方向，获得在水平方向对应的能量 a 和垂直方向对应的能量 b 。之后通过计算它们之间的反正切值，即获得相位一致性特征方向，计算公式如下： <math display="block" id="M26"><mtable columnalign="left"><mtr><mtd><mi>a</mi><mo>=</mo><munder><mo movablelimits="false">∑</mo><mrow><mi>θ</mi></mrow></munder><mrow><mo stretchy="false">(</mo><mi>o</mi><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mi>c</mi><mi>o</mi><mi>s</mi><mi>θ</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mi>b</mi><mo>=</mo><munder><mo movablelimits="false">∑</mo><mrow><mi>θ</mi></mrow></munder><mrow><mo stretchy="false">(</mo><mi>o</mi><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mi>s</mi><mi>i</mi><mi>n</mi><mi>θ</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mi>ϕ</mi><mo>=</mo><mi mathvariant="normal">a</mi><mi mathvariant="normal">r</mi><mi mathvariant="normal">c</mi><mi mathvariant="normal">t</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">n</mi><mo stretchy="false">(</mo><mi>b</mi><mo>,</mo><mi>a</mi><mo stretchy="false">)</mo></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824428&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824427&type= 27.51666832 20.40466690 式中， <math id="M27"><mi>o</mi><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824430&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824429&type= 6.60400009 3.47133350 表示在方向 <math id="M28"><mi>θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824432&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824431&type= 1.77800000 2.70933342 上的奇对称滤波器卷积结果，并将特征方向限制在 <math id="M29"><mfenced open="[" close="]" separators="|"><mrow><mn mathvariant="normal">0</mn><mo>,</mo><mi>π</mi></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824440&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824439&type= 9.05933285 3.89466691 ，以解决多传感器遥感影像间灰度反向的问题（Zhang等， 2008）。 （3） 构建DFOP描述符：Ye等的研究表明，逐像素的特征描述符能够获得影像的精细结构特征（ Ye等， 2019 Ye等， 2019 ）。因此，本文利用相位一致性特征方向和特征值构建一种逐像素的三维结构特征描述符，即首先将相位一致性特征方向划分为多个方向特征层，然后将能够描述影像几何结构属性的相位一致性特征值按照一定的权重划分到不同的方向特征层。如 图1 图1 所示，以为 <math id="M30"><msup><mrow><mn mathvariant="normal">0</mn></mrow><mrow><mo>°</mo></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824436&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824435&type= 2.79399991 2.96333337 起点，以 <math id="M31"><msup><mrow><mn mathvariant="normal">30</mn></mrow><mrow><mo>°</mo></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824438&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824437&type= 4.65666676 2.96333337 为间隔，将步骤（2）中计算的相位一致性特征方向 <math id="M32"><mfenced open="[" close="]" separators="|"><mrow><mn mathvariant="normal">0</mn><mo>,</mo><mi>π</mi></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824440&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824439&type= 9.05933285 3.89466691 划分为6个等分方向区间，并以等分方向区间的中心构建一个方向层（第一层的方向为 <math id="M33"><msup><mrow><mn mathvariant="normal">15</mn></mrow><mrow><mo>°</mo></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824442&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824441&type= 4.48733330 2.96333337 ），则这6个等分方向区间构成了6个不同的方向层。随后，我们可以直接将相位一致性特征值根据其特征方向投影到对应的方向层，但直接根据特征方向划分没有考虑到特征值对相邻方向层的贡献，会导致对影像的结构信息描述不够准确。如当某一像元的特征方向为 <math id="M34"><mn mathvariant="normal">30</mn><mo>.</mo><msup><mrow><mn mathvariant="normal">1</mn></mrow><mrow><mo>°</mo></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824444&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824443&type= 7.11199999 2.96333337 ，特征值为100时，若采用直接根据特征方向对特征值投影的方式，该像元的特征值需要投影到方向层2中，但其特征方向与第1层和第2层的距离分别为 <math id="M35"><mn mathvariant="normal">15</mn><mo>.</mo><msup><mrow><mn mathvariant="normal">1</mn></mrow><mrow><mo>°</mo></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824446&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824445&type= 7.11199999 2.96333337 和 <math id="M36"><mn mathvariant="normal">14</mn><mo>.</mo><msup><mrow><mn mathvariant="normal">9</mn></mrow><mrow><mo>°</mo></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824448&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824447&type= 7.28133297 2.96333337 ，即该像元的特征方向在第1层与第2层中的方向层中的贡献是接近的。因此，我们首先采用距离加权的方式，将每个像元的特征方向根据其与相邻方向层距离分别分配对应权重到方向层中。在这里，由于 <math id="M37"><mn mathvariant="normal">0</mn><mo>~</mo><msup><mrow><mn mathvariant="normal">15</mn></mrow><mrow><mo>°</mo></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824450&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824449&type= 8.38199997 2.96333337 只与第1层相邻，其对应的特征值在第1层的权重为1，同理， <math id="M38"><msup><mrow><mn mathvariant="normal">165</mn></mrow><mrow><mo>°</mo></mrow></msup><mo>~</mo><msup><mrow><mn mathvariant="normal">180</mn></mrow><mrow><mo>°</mo></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824452&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824451&type= 14.98600006 2.96333337 对应的特征值在第6层的权重也为1。 随后，根据图像像元在每个方向层中的权重，将像元的相位一致性特征值计算到对应方向层中的像元位置，形成逐像素的具有几何结构特征的三维稠密特征描述符。为了减少影像噪声的影响，采用二维高斯滤波对三维描述符的每一层进行高斯滤波，之后利用一维高斯滤波沿方向层（以下简称为Z轴）进行平滑处理。最后，为了进一步减小多传感器遥感影像间的非线性灰度差异问题，通过对滤波后的特征描述符沿Z轴进行归一化处理，形成稠密的三维特征描述符。由于提出的特征描述符是基于影像的相位一致性特征方向和特征值，且是逐像素进行特征表达，可精确地描述了影像的形状和几何结构属性，整个构建流程如 图2 图2 所示。 2.2　子像素快速匹配策略 2 由于所构建的DFOP描述符是逐像素的三维特征描述符，其数据量较大。若直接采用传统的相似性度量进行匹配，则存在匹配速度较慢的问题，因此本文采用单步离散傅里叶变换方法（Guizar等， 2008）对方向相位稠密特征进行子像素快速匹配，整个匹配流程如 图3 图3 所示，其流程可分为以下两步： （1） 计算DFOP描述符间的整像素平移：假定 <math id="M39"><mi>f</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824477&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824476&type= 13.37733269 3.64066648 和 <math id="M40"><mi>g</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824479&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824478&type= 13.63133430 3.64066648 分别为参考影像和输入影像方向相位通道特征描述符，特征描述符间仅存在平移关系，即 <math id="M41"><mi>f</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo stretchy="false">)</mo><mo>=</mo><mi>g</mi><mo stretchy="false">(</mo><mi>x</mi><mo>-</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>,</mo><mi>y</mi><mo>-</mo><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>,</mo><mi>z</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824467&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824466&type= 46.14333344 4.74133301 。根据傅里叶变换的性质，两个函数在空间域中的平移等于其在频率域中的线性相位差（ Foroosh等， 2002 Foroosh等， 2002 ），在频率域下两者的关系可以表示为： <math id="M42"><mi>F</mi><mo stretchy="false">(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo stretchy="false">)</mo><mo>=</mo><mi>G</mi><mo stretchy="false">(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo stretchy="false">)</mo><mi mathvariant="normal">e</mi><mi mathvariant="normal">x</mi><mi mathvariant="normal">p</mi><mo stretchy="false">(</mo><mo>-</mo><mn mathvariant="normal">2</mn><mi>π</mi><mi>i</mi><mo stretchy="false">(</mo><mi>u</mi><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>+</mo><mi>v</mi><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo stretchy="false">)</mo><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824469&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824468&type= 66.88666534 4.74133301 ，它们的归一化互功率谱可以写为： <math display="block" id="M43"><mtable columnalign="left"><mtr><mtd><mi>Q</mi><mo stretchy="false">(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mi>F</mi><mo stretchy="false">(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo stretchy="false">)</mo><mi>G</mi><msup><mrow><mo stretchy="false">(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo stretchy="false">)</mo></mrow><mrow><mi mathvariant="normal">*</mi></mrow></msup></mrow><mrow><mfenced open="|" close="|" separators="|"><mrow><mi>F</mi><mo stretchy="false">(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo stretchy="false">)</mo><mi>G</mi><msup><mrow><mo stretchy="false">(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo stretchy="false">)</mo></mrow><mrow><mi mathvariant="normal">*</mi></mrow></msup></mrow></mfenced></mrow></mfrac><mo>=</mo></mtd></mtr><mtr><mtd><mi mathvariant="normal">e</mi><mi mathvariant="normal">x</mi><mi mathvariant="normal">p</mi><mo stretchy="false">(</mo><mo>-</mo><mn mathvariant="normal">2</mn><mi>π</mi><mi>i</mi><mo stretchy="false">(</mo><mi>u</mi><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>+</mo><mi>v</mi><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mover accent="true"><mi mathvariant="bold-italic">k</mi><mo>⃗</mo></mover></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824471&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824470&type= 50.29199982 16.84866714 式中， <math id="M44"><mi>F</mi><mo stretchy="false">(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824473&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824472&type= 14.81666565 3.47133350 和 <math id="M45"><mi>G</mi><mo stretchy="false">(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824475&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824474&type= 14.73200035 3.47133350 分别为 <math id="M46"><mi>f</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824477&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824476&type= 13.37733269 3.64066648 和 <math id="M47"><mi>g</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824479&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824478&type= 13.63133430 3.64066648 的傅里叶变换形式，*表示复共轭， <math id="M48"><mover accent="true"><mi mathvariant="bold-italic">k</mi><mo>⃗</mo></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824481&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824480&type= 1.77800000 3.89466691 是一个3维单位矢量。通过对归一化互功率谱 <math id="M49"><mi>Q</mi><mo stretchy="false">(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824491&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824490&type= 14.81666565 3.64066648 进行傅里叶逆变换可以得到一个脉冲函数 <math id="M50"><mi>δ</mi><mo stretchy="false">(</mo><mi>x</mi><mo>-</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>,</mo><mi>y</mi><mo>-</mo><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>,</mo><mi>z</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824485&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824484&type= 28.19400024 4.74133301 ，特征描述符间的整像素平移量 <math id="M51"><mo stretchy="false">(</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>,</mo><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824487&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824486&type= 10.58333302 4.74133301 可由脉冲函数峰值的位置计算得到。 （2） 计算DFOP描述间的子像素平移： 单步离散傅里叶变换方法通过在初始同名点邻域位置（邻域半径 R =1.5像素）进行上采样互相关获取同名点的子像素精度。根据设定的子像素精度因子 <math id="M52"><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824505&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824504&type= 2.37066650 3.72533321 ，以及归一化互功率谱 <math id="M53"><mi>Q</mi><mo stretchy="false">(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824491&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824490&type= 14.81666565 3.64066648 ，构建两个大小分别为 <math id="M54"><mi mathvariant="bold">k</mi><mi mathvariant="bold">e</mi><msub><mrow><mi mathvariant="bold">r</mi></mrow><mrow><mi>c</mi></mrow></msub><mo stretchy="false">(</mo><mn mathvariant="normal">1.5</mn><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub><mo>,</mo><mi>N</mi><mo stretchy="false">)</mo><mo>,</mo><mi mathvariant="bold">k</mi><mi mathvariant="bold">e</mi><msub><mrow><mi mathvariant="bold">r</mi></mrow><mrow><mi>r</mi></mrow></msub><mo stretchy="false">(</mo><mi>M</mi><mo>,</mo><mn mathvariant="normal">1.5</mn><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824493&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824492&type= 45.29666519 4.57200003 的频率域核矩阵。之后将核矩阵与归一化互功率谱进行矩阵相乘获得子像素邻域影像 <math id="M55"><mi>M</mi><mi>a</mi><mi>t</mi><mo stretchy="false">(</mo><mn mathvariant="normal">1.5</mn><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub><mo>,</mo><mn mathvariant="normal">1.5</mn><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824503&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824502&type= 24.63800049 4.57200003 ，表示如下： <math display="block" id="M56"><mtable columnalign="left"><mtr><mtd><mi>M</mi><mi>a</mi><mi>t</mi><mo stretchy="false">(</mo><mn mathvariant="normal">1.5</mn><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub><mo>,</mo><mn mathvariant="normal">1.5</mn><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub><mo stretchy="false">)</mo><mo>=</mo></mtd></mtr><mtr><mtd><mi mathvariant="bold">k</mi><mi mathvariant="bold">e</mi><msub><mrow><mi mathvariant="bold">r</mi></mrow><mrow><mi>c</mi></mrow></msub><mo stretchy="false">(</mo><mn mathvariant="normal">1.5</mn><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub><mo>,</mo><mi>N</mi><mo stretchy="false">)</mo><mo>×</mo><mi mathvariant="bold-italic">Q</mi><mo stretchy="false">(</mo><mi>N</mi><mo>,</mo><mi>M</mi><mo stretchy="false">)</mo><mo>×</mo><mi mathvariant="bold">k</mi><mi mathvariant="bold">e</mi><msub><mrow><mi mathvariant="bold">r</mi></mrow><mrow><mi>r</mi></mrow></msub><mo stretchy="false">(</mo><mi>M</mi><mo>,</mo><mn mathvariant="normal">1.5</mn><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub><mo stretchy="false">)</mo></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824497&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824496&type= 62.82266998 9.99066734 式中， <math id="M57"><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824499&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824498&type= 2.37066650 3.72533321 表示设定的子像素精度因子，1.5为邻域半径， M 和 N 分别为模板影像的行和列， <math id="M58"><mi>Q</mi><mo stretchy="false">(</mo><mi>N</mi><mo>,</mo><mi>M</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824501&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824500&type= 13.20800018 3.64066648 表示归一化互功率谱的第一层。之后通过将子像素邻域影像 <math id="M59"><mi>M</mi><mi>a</mi><mi>t</mi><mo stretchy="false">(</mo><mn mathvariant="normal">1.5</mn><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub><mo>,</mo><mn mathvariant="normal">1.5</mn><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824503&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824502&type= 24.63800049 4.57200003 的峰值位置与子像素精度因子 <math id="M60"><msub><mrow><mi>k</mi></mrow><mrow><mi>f</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824505&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824504&type= 2.37066650 3.72533321 相除得到同名点处的子像素平移量 <math id="M61"><mo stretchy="false">(</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>,</mo><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824507&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824506&type= 10.58333302 4.74133301 ，最后将整像素平移量 <math id="M62"><mo stretchy="false">(</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>,</mo><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824509&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824508&type= 10.58333302 4.74133301 与子像素平移量 <math id="M63"><mo stretchy="false">(</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>,</mo><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824511&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824510&type= 10.58333302 4.74133301 相加得到特征描述符间的子像素级平移量 <math id="M64"><mo stretchy="false">(</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo>,</mo><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824513&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824512&type= 10.58333302 4.74133301 。 3　多传感器遥感影像配准系统流程 1 上述描述的匹配方法可以有效地抵抗图像间的非线性灰度差异，为了处理图像间的几何畸变（如旋转和尺度差异），本文利用卫星图像自带的地理信息或有理多项式参数进行空间约束和几何粗纠正，并在此基础上研制了一套多传感器遥感影像配准系统。配准流程图如 图4 图4 所示，主要包括粗配准和精配准两个阶段： 粗配准阶段： （1） 选择粗配准模式：所研发的系统提供了两种配准模式，分别为基于有理函数模型（rational function model， RFM）（ Pan等， 2019 Pan等， 2019 ）的配准模式和基于地理坐标信息的配准模式。读取参考影像和输入影像，根据输入影像是否经过地理校正然后采用对应的粗配准模式。 （2） 特征点提取：由于FAST算子（ Rosten和Drummond， 2006 Rosten和Drummond， 2006 ）已被证明具有快速提取特征点的优势，本文首先使用FAST算子在参考影像上采用分块的策略提取特征点。分块的目的是避免提取的特征点在整幅影像上分布不均匀，产生聚簇现象，如 图5 图5 （a）所示。根据提取的总特征点数量，将参考影像划分为 <math id="M65"><mi>M</mi><mo>×</mo><mi>M</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824518&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824517&type= 10.32933331 2.62466669 的互不重叠的区域，然后分别在每块区域采用FAST算子进行特征点提取，之后根据响应值的大小从大到小选择前 <math id="M66"><mi>k</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824520&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824519&type= 1.69333339 2.62466669 个点作为特征点。 （3） 区域粗配准：由于参考影像与输入影像间往往存在尺度和旋转差异，因此对匹配模板提取DFOP描述符前，首先需要消除模板影像间的尺度和旋转差异。根据步骤（1）中的粗配准模式，即若输入影像经过地理校正，则直接根据特征点的地理坐标 <math id="M67"><mi>P</mi><mo stretchy="false">(</mo><mi>G</mi><mi>e</mi><msub><mrow><mi>o</mi></mrow><mrow><mi>x</mi></mrow></msub><mo>,</mo><mi>G</mi><mi>e</mi><msub><mrow><mi>o</mi></mrow><mrow><mi>y</mi></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824525&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824524&type= 20.23533440 4.57200003 直接预测其在输入影像上的初始位置 <math id="M68"><mi>P</mi><mo stretchy="false">(</mo><mi>X</mi><mo>,</mo><mi>Y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824531&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824530&type= 12.36133385 3.47133350 ，然后以点 <math id="M69"><mi>P</mi><mo stretchy="false">(</mo><mi>X</mi><mo>,</mo><mi>Y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824531&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824530&type= 12.36133385 3.47133350 为中心，根据参考影像和输入影像的分辨率比值确定一定尺寸的模板影像，然后建立几何变换关系消除模板影像间的旋转和尺度差异形成粗配准。若输入影像为携带RPC文件的卫星影像，则依据RFM模型预测特征点在输入影像上的图像坐标 <math id="M70"><mi>P</mi><mo stretchy="false">(</mo><mi>X</mi><mo>,</mo><mi>Y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824531&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824530&type= 12.36133385 3.47133350 ，然后根据影像间的分辨率比值关系确定模板影像，最后对模板影像进行局部几何校正实现影像间的区域粗配准，具体操作见文献（ 王蒙蒙等， 2022 王蒙蒙等， 2022 ）。 精配准阶段： （4） 构建DFOP描述符：根据1.1节所提出的DFOP描述符构建方法，对形成粗配准的模板影像构建DFOP特征描述符。 （5） 子像素快速匹配：采用1.2节所提出的快速匹配方法，将特征描述符变换到频率域空间，采用单步离散傅里叶方法实现同名点的快速匹配。 （6） 对参考影像上的特征点进行遍历匹配，即重复（3）~（5）步骤，实现所有特征点的快速匹配。 （7） 粗差剔除：由于多传感器影像间因视角、大气条件等差异造成的影响，使得匹配获取的同名点对不可避免地存在错误同名点对，因此需要对匹配结果进一步优化，保证匹配结果的精度和可靠性。根据本文方法配准模式的不同，本文方法在粗差剔除和几何校正部分也分为了两种模式。即若输入影像包含地理空间信息，则直接采用随机抽样一致性算法（Random Sample Consensus， RANSAC）（ Fischler和Bolles， 1987 Fischler和Bolles， 1987 ）建立透视变换模型进行粗差剔除对初始匹配结果进行优化，然后根据优化后的同名点用最小二乘法求解透视变换模型系数，最后根据模型系数对输入影像进行几何校正实现影像间的配准，如下式所示： <math display="block" id="M71"><mtable columnalign="left"><mtr><mtd><msup><mrow><mi>x</mi></mrow><mrow><mi mathvariant="normal">'</mi></mrow></msup><mo>=</mo><mfrac><mrow><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mi>x</mi><mo>+</mo><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mi>y</mi><mo>+</mo><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub></mrow><mrow><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">7</mn></mrow></msub><mi>x</mi><mo>+</mo><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">8</mn></mrow></msub><mi>y</mi><mo>+</mo><mn mathvariant="normal">1</mn></mrow></mfrac></mtd></mtr><mtr><mtd><msup><mrow><mi>y</mi></mrow><mrow><mi mathvariant="normal">'</mi></mrow></msup><mo>=</mo><mfrac><mrow><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">4</mn></mrow></msub><mi>x</mi><mo>+</mo><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">5</mn></mrow></msub><mi>y</mi><mo>+</mo><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">6</mn></mrow></msub></mrow><mrow><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">7</mn></mrow></msub><mi>x</mi><mo>+</mo><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">8</mn></mrow></msub><mi>y</mi><mo>+</mo><mn mathvariant="normal">1</mn></mrow></mfrac></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824533&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824532&type= 27.51666832 19.30400085 <math display="block" id="M72"><mfenced open="|" close="|" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd><mtd><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mtd><mtd><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">4</mn></mrow></msub></mtd><mtd><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">5</mn></mrow></msub></mtd><mtd><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">6</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">7</mn></mrow></msub></mtd><mtd><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">8</mn></mrow></msub></mtd><mtd><mn mathvariant="normal">1</mn></mtd></mtr></mtable></mrow></mfenced><mo>≠</mo><mn mathvariant="normal">0</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824535&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824534&type= 22.60600090 13.03866673 式中， <math id="M73"><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824537&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824536&type= 8.21266651 3.64066648 和 <math id="M74"><mo stretchy="false">(</mo><msup><mrow><mi>x</mi></mrow><mrow><mi mathvariant="normal">'</mi></mrow></msup><mo>,</mo><msup><mrow><mi>y</mi></mrow><mrow><mi mathvariant="normal">'</mi></mrow></msup><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824539&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824538&type= 9.82133293 3.97933316 分别表示同名点在输入影像和参考影像上的像素坐标， <math id="M75"><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>.</mo><mo>.</mo><mo>.</mo><msub><mrow><mi>l</mi></mrow><mrow><mn mathvariant="normal">8</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824541&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824540&type= 7.28133297 3.72533321 为透视变换模型系数。若输入影像不包含地理空间信息，但携带RPC参数的卫星影像，则采用RANSAC根据RFM模型对初始匹配结果进行粗差剔除。但由于RFM模型通常包含系统误差，因此本文在粗差剔除时，在RFM模型像方中添加仿射变换模型对其进行改正（ Grodecki和Dial， 2003 Grodecki和Dial， 2003 ），改正后的模型可以写为： <math display="block" id="M76"><mtable columnalign="left"><mtr><mtd><mi>X</mi><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mi>X</mi><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mi>Y</mi><mo>=</mo><mi mathvariant="normal">R</mi><mi mathvariant="normal">F</mi><msub><mrow><mi mathvariant="normal">M</mi></mrow><mrow><mi>x</mi></mrow></msub><mo stretchy="false">(</mo><mi>P</mi><mo>,</mo><mi>L</mi><mo>,</mo><mi>H</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mi>Y</mi><mo>+</mo><msub><mrow><mi>b</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>+</mo><msub><mrow><mi>b</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mi>X</mi><mo>+</mo><msub><mrow><mi>b</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mi>Y</mi><mo>=</mo><mi mathvariant="normal">R</mi><mi mathvariant="normal">F</mi><msub><mrow><mi mathvariant="normal">M</mi></mrow><mrow><mi>y</mi></mrow></msub><mo stretchy="false">(</mo><mi>P</mi><mo>,</mo><mi>L</mi><mo>,</mo><mi>H</mi><mo stretchy="false">)</mo></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824543&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824542&type= 55.03333664 9.82133293 式中， <math id="M77"><mo stretchy="false">(</mo><mi>X</mi><mo>,</mo><mi>Y</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824545&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824544&type= 9.56733322 3.47133350 为输入影像的图像坐标， <math id="M78"><mo stretchy="false">(</mo><mi>P</mi><mo>,</mo><mi>L</mi><mo>,</mo><mi>H</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824547&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824546&type= 13.46199989 3.47133350 为参考影像上同名点的地理空间坐标， <math id="M79"><mi mathvariant="normal">R</mi><mi mathvariant="normal">F</mi><msub><mrow><mi mathvariant="normal">M</mi></mrow><mrow><mi>x</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824549&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824548&type= 8.46666718 3.80999994 和 <math id="M80"><mi mathvariant="normal">R</mi><mi mathvariant="normal">F</mi><msub><mrow><mi mathvariant="normal">M</mi></mrow><mrow><mi>y</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824551&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824550&type= 8.46666718 3.80999994 分别表示地理空间坐标到影像像素坐标 X 和 Y 的变换， <math id="M81"><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>.</mo><mo>.</mo><mo>.</mo><msub><mrow><mi>b</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824553&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824552&type= 8.72066593 3.72533321 为添加的仿射变换模型系数。在粗差剔除后，使用RFM模型和添加的仿射变化模型对输入的卫星影像进行几何校正和输出保存，实现输入影像和参考影像间的高精度快速配准。 （8） 影像输出与保存：对配准后的输入影像进行输出保存。 4　实验与分析 1 为了验证本文方法和系统的有效性，首先将所提出的匹配算法与其它经典的模板匹配算法进行对比分析。所对比的匹配方法包括了基于灰度相似性的NCC和MI算法，以及基于结构特征的HOPC和CFOG算法。然后将自主研制的配准系统与目前著名的商业遥感影像处理软件ENVI和ERDAS中的配准模块进行配准实验对比。实验数据平台处理器是Inter（R） Core（TM） i7-10700K CPU 3.79GHz，RAM是32GB，Windows10 X64操作系统。 4.1　实验数据 2 本文利用10组多传感器遥感影像 （包括了光学、SAR和LiDAR等， 图6 图6 ） 对本文方法和研发系统的有效性进行测试和评价。其中，第1~8组影像用于测试DFOP描述符的匹配性能，由于采用的是模板匹配的方法，因此在测试前已经对每组影像进行了粗配准，使影像间仅存在平移关系。第9~10组影像用于测试所研发系统的配准性能，每组实验影像间不仅存在旋转、尺度和时相差异，且存在非线性灰度差异。所使用的10组影像涵盖了多种地貌类型，如结构信息丰富的城镇、畸变较大的山地等，实验数据的详细信息如 表1 表1 所示。 对比方案 对比类型 实验编号 传感器 影像大小（像素） 地面采样间隔（米） 时间（年-月） 16 方向性能对比 8 光学-SAR 2 1 GF-2 1500×1500 5 2018-04 GF-3 1500×1500 5 2018-03 2 2 ZY-3 1500×1500 5 2017-01 GF-3 1500×1500 5 2018-03 2 3 Sentinel-1 1500×1500 10 2018-10 Sentinel-2 1500×1500 10 2018-8 2 4 GF-1 1500×1500 8 2019-12 GF-3 1500×1500 8 2017-05 6 光学-光学 2 5 GF-1 1500×1500 16 2020-10 GF-6 1500×1500 16 2020-11 2 6 Google earth 1500×1500 1 2019-11 Ikonos 1500×1500 1 2000-02 2 7 Google earth 1500×1500 0.5 2020-01 SuperView-1 1500×1500 0.5 2020-05 2 LIDAR-光学 2 8 LIDAR 524×524 2.5 2012-06 Airborne visible 524×524 2.5 2012-06 4 软件模块对比 2 光学-光学 2 9 Google earth 17808×18973 1 2012-03 Quickbird 6819×7193 2.5 2002-01 2 光学-SAR 2 10 GF-3 7751×7393 5 2018-03 GF-2 7300×6908 4 2018-04 4.2　匹配性能对比分析 2 在对比实验中，这里采用匹配正确率（匹配正确率=正确匹配点对/总匹配点对）和匹配时间作为评价指标对五种算法进行定量对比分析。对比方案如下：首先使用FAST算子在参考影像上采用分块的策略提取200个均匀分布的特征点，并根据匹配模板的大小在输入影像上确定［-20， 20］的搜索区域，之后使用不同尺寸的匹配模板进行模板窗口滑动匹配，并采用投影变换作为约束条件把误差在2个像素的匹配点视为正确匹配点。前8组实验的匹配正确率随模板大小变化的结果如 图7 图7 所示， 图8 图8 显示了本文方法在模板大小为 <math id="M82"><mn mathvariant="normal">105</mn><mo>×</mo><mn mathvariant="normal">105</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824570&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824569&type= 15.66333294 2.70933342 像素时的正确匹配点。 size of <math id="M84"><mn mathvariant="normal">105</mn><mo>×</mo><mn mathvariant="normal">105</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824570&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824569&type= 15.66333294 2.70933342 pixels for all the test cases. 从 图7 图7 的实验结果可以看出，在第1~4组的光学与SAR影像匹配实验数据中，直接依赖于灰度相关性的NCC方法在4组实验中均未能取得令人满意的结果，这是因为NCC对影像间的灰度只具有线性不变性，而难以处理具有非线性灰度差异的遥感影像。MI方法的匹配结果虽优于NCC，但较于HOPC、CFOG以及本文提出的DFOP，其匹配正确率依然较低。本文方法在实验结果中的表现优于CFOG，其原因是DFOP是采用相位一致性构建的逐像素的结构特征描述符，而CFOG采用的影像梯度。相较于梯度，相位一致性对影像间的对比度和非线性灰度变化具有更强的稳健性。另外，同样采用相位一致性的HOPC的匹配正确率也低于DFOP，这是因为HOPC虽采用相位一致性构建特征描述符，但其构建特征描述符时采用稀疏抽样的方式（ Ye等，2017 Ye等，2017 ），而DFOP采用的是逐像素的特征表达方式，可以更精确捕获影像的结构细节。在实验4结果中（ 图7 图7 （d）），由于实验数据的地形为山区，HOPC、CFOG和DFOP的匹配正确率均有所下降，其原因是山区影像内部存在显著的局部畸变，导致两幅影像间提取的结构特征描述符存在一定程度的形变。而在实验5中，每种方法的匹配正确率均处于同一水平，从实验数据分析来看，实验5的数据分别是GF-1和GF-6，二者均为中分辨率光学卫星影像，实验数据影像内部几何结构特征丰富，数据间几乎不存在时相差异和非线性灰度差异，因此，每种方法的匹配正确率均达到了90%。而同样是光学影像，但在影像间传感器不同的实验6和实验7中，每种方法的匹配正确率却未能达到实验5中的匹配正确率结果。从影像数据信息来看，而在实验6中，影像间虽不存在时相差异，但影像内部的几何结构特征较少，造成了匹配正确率下降的结果。实验7中的影像地形为山区，影像内部存在局部畸变的影响，致使匹配正确率下降，尽管如此，DFOP的匹配正确率依然高于其它方法。实验8中LIDAR与光学影像的匹配结果同样证明了DFOP的优越性。通过以上多组多传感器数据实验，证明了所提出的DFOP对遥感影像间的非线性灰度差异具有更强的稳健性。 各方法的匹配耗时如 图9 图9 所示，由于MI需要计算每个匹配窗口的联合直方图，它的计算复杂度较高，在所有方法中最为耗时。相比于MI，NCC的计算效率具有一定程度的提高，但是随着模板窗口的增大，其计算时间增长剧烈。对于采用模板几何结构特征的HOPC、CFOG以及DFOP三种方法，由于影像的相位一致性计算过程比影像的梯度信息复杂，因此在构建特征描述符时，由相位一致性信息构建的HOPC和DFOP比采用影像梯度信息的CFOG所花费的时间较多，所以它们的匹配时间较长。对于HOPC和DFOP而言，由于HOPC是在空间域中利用特征描述符的归一化相关系数进行匹配，而DFOP是将特征描述符转换到频率空间，利用快速傅里叶变换进行匹配加速，其计算效率要明显优于HOPC，尤其是当模板窗口逐渐增大时，DFOP的计算优势更加明显。 4.3　配准系统对比 2 通过上述实验的比较分析可知，DFOP描述符对多传感器遥感影像间的非线性灰度差异具有更强的鲁棒性，所提出的匹配方法也能够实现影像间的快速稳健匹配。于是，本文基于所提出的DFOP研制了一套多传感器遥感影像配准系统，并将其与目前著名的商业遥感影像处理软件ENVI和ERDAS中的配准模块进行配准实验对比，其中，ENVI采用MI作为相似性测度进行匹配，ERDAS采用NCC作为相似性测度，并利用了影像金字塔进行加速匹配。使用的数据为第9~10组实验数据，其中，第9组中的Google earth和第10组中的GF-3均为几何校正后的影像，并将它们作为每组实验中的参考影像。而第9组中的输入影像为带有地理坐标信息的QuickBird影像，与参考影像间存在尺度、时相等差异。第10组实验中的输入影像为未经过几何校正但带有RPC参数文件的GF-2卫星影像。实验设置如下：在每组实验中设定匹配200对同名点，误差剔除阈值设置为2个像素，匹配模板大小均设置为 <math id="M85"><mn mathvariant="normal">85</mn><mo>×</mo><mn mathvariant="normal">85</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824575&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24824574&type= 11.93799973 2.70933342 。采用人工刺点的方式在参考影像和配准后的输入影像间 30对同名点作为检查点，并分别从匹配正确率和均方根误差（root mean squared error， RMSE）两个方面对所开发的系统与其它软件中的配准模块进行对比评价。使用的软件版本分别为ENVI5.3、ERDAS2014，本文所研制的配准系统基于Qt5.0采用C++语言和GDAL开源库实现。 实验9和实验10在各软件配准模块的正确匹配点分布对比如 图10 图10 所示，其中，×表示对应软件中不包含的配准模式。 表2 表2 显示了实验9和实验10的匹配正确率以及RMSE。从实验的结果可以看出，ENVI和ERDAS软件中的配准模块在实验9和实验10中的匹配正确率和RMSE均较低。由于实验9影像间的时相差异近10年，匹配难度较大，ENVI的匹配正确率仅为19.5%，而由于ERDAS采用金字塔的匹配策略，因而其匹配正确率高于ENVI，匹配正确率达到了36%，但所研制系统的匹配正确率为66%，远高于ENVI和ERDAS中的配准模块。在实验10中，由于ERDAS不具基于卫星RPC参数的配准模式，这里以×表示配准失败。所开发的配准系统采用DFOP描述符对影像的几何结构进行特征描述和表达，因而对实验10中光学和SAR影像间的非线性灰度差异具有更强的鲁棒性，匹配正确率为64.33%，远远高于ENVI的10.67%。实验9和实验10配准后的棋盘格及其放大子图分别如图 11 11 、 12 12 所示。从棋盘格子图及其RMSE可以看出，由于本文系统匹配的正确同名点更多，所建立的几何变换关系更加稳健，因而配准后的精度更高。 注： ×表示对应软件中不包含的配准模式（即基于RPC的配准） 2 实验编号 2 ENVI 2 ERDAS 2 本文系统 匹配正确率 RMSE 匹配正确率 RMSE 匹配正确率 RMSE 9 19.5% 1.62 36% 1.37 66% 1.09 10 10.67% 8.35 × × 64.33% 0.97 4　结论 1 针对多传感器遥感影像间几何畸变和灰度差异造成的配准困难问题，本文利用相位一致性构建了一种能够表达影像精细结构特征的特征描述符DFOP，并采用模板匹配的策略在频率域空间设计了一种快速稳健的子像素匹配方法。为了验证DFOP描述符的性能和表现，将其与MI、NCC、HOPC以及CFOG等进行对比与分析，结果表明：DFOP在多组多传感器遥感影像中获得了最高的匹配正确率。此外，本文基于所提出的匹配方法研制了一套多传感器遥感影像配准系统，并通过两组实验证明了所研制的系统配准结果优于ENVI和ERDAS软件中的配准模块。 本文所构建的DFOP描述符虽能够精细表达影像的结构特征，但其不具有旋转和尺度不变性，此外，所研制的配准系统是基于输入影像需具备地理坐标信息或遥感卫星影像携带RPC参数信息文件。因此，在今后的研究中，将针对描述符的旋转和尺度问题展开研究。 参考文献（ References ） Bruzzone L and Bovolo F 2013 . A Novel Framework for the Design of Change-Detection Systems for Very-High-Resolution Remote Sensing Images . Proceedings of the IEEE . 101 ( 3 ): 609 – 630 [ DOI: 10.1109/JPROC.2012.219-7169 http://dx.doi.org/10.1109/JPROC.2012.219-7169 ] A Novel Framework for the Design of Change-Detection Systems for Very-High-Resolution Remote Sensing Images Dellinger F , Delon J , Gousseau Y , Michel J and Tupin F . 2015 . SAR-SIFT: A SIFT-Like Algorithm for SAR Images . IEEE Transactions on Geoscience and Remote Sensing , 53 ( 1 ): 453 – 466 [ DOI: 10.1109/TGRS.2014.2323552 http://dx.doi.org/10.1109/TGRS.2014.2323552 ] SAR-SIFT: A SIFT-Like Algorithm for SAR Images Feng R T , Shen H F , Bai J J and Li X H . 2021 . Advances and Opportunities in Remote Sensing Image Geometric Registration: A Systematic Review of State-of-the-Art Approaches and Future Research Directions [J]. IEEE Geoscience and Remote Sensing Magazine , 9 ( 4 ): 120 - 142 [ DOI: 10.1109/MGRS.2021.3081763 http://dx.doi.org/10.1109/MGRS.2021.3081763 ]. Advances and Opportunities in Remote Sensing Image Geometric Registration: A Systematic Review of State-of-the-Art Approaches and Future Research Directions Fischler M A and Bolles R C . 1987 . Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography . Readings in Computer Vision , 726 - 740 [ DOI: 10.1016/B978-0-08-051581-6.50070-2 http://dx.doi.org/10.1016/B978-0-08-051581-6.50070-2 ] Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography Foroosh H , Zerubia J B and Berthod M . 2002 . Extension of phase correlation to subpixel registration . IEEE Transactions on Image Processing , 11 ( 3 ): 188 - 200 [ DOI: 10.1109/83.9889-53 http://dx.doi.org/10.1109/83.9889-53 ] Extension of phase correlation to subpixel registration Goncalves H , Corte-Real L and Goncalves J A . 2011 . Automatic Image Registration Through Image Segmentation and SIFT [J]. IEEE Transactions on Geoscience and Remote Sensing , 49 ( 7 ): 2589 - 2600 [ DOI: 10.1109/TGRS.2011.21-09389 http://dx.doi.org/10.1109/TGRS.2011.21-09389 ] Automatic Image Registration Through Image Segmentation and SIFT Grodecki J and Dial G . 2003 . Block adjustment of high-resolution satellite images described by rational polynomials . Photogrammetric Engineering and Remote Sensing , 69 ( 1 ): 59 - 68 [ DOI: 10.14318/PERS.69.1.59 http://dx.doi.org/10.14318/PERS.69.1.59 ] Block adjustment of high-resolution satellite images described by rational polynomials Guizar-Sicairos M , Thurman S T and Fienup J R . 2008 . Efficient subpixel image registration algorithm . Optics Letters , 33 ( 2 ): 156 - 158 [ DOI: 10.1364/OL.33.000156 http://dx.doi.org/10.1364/OL.33.000156 ] Efficient subpixel image registration algorithm Hisham M B , Yaakob S N , Raof R A A , Nazren A A and Wafi N . 2015 . Template Matching Using Sum of Squared Difference and Normalized Cross Correlation // 2015 IEEE Student Conference on Research and Development (SCOReD) . Kuala Lumpur, Malaysia . 100 – 104 [ DOI: 10.-1109/ SCORED.2015.7449303 http://dx.doi.org/10.-1109/SCORED.2015.7449303 ] Template Matching Using Sum of Squared Difference and Normalized Cross Correlation Huo C L , Pan C H , Huo L G and Zhou Z X . 2012 . Multilevel SIFT Matching for Large-Size VHR Image Registration [J]. IEEE Geoscience and Remote Sensing Letters , 9 ( 2 ): 171 - 175 [ DOI: 10.1109/LGRS.2011.216349-1 http://dx.doi.org/10.1109/LGRS.2011.216349-1 ] Multilevel SIFT Matching for Large-Size VHR Image Registration Kovesi P . 1999 . Image Features from Phase Congruency . Journal of Computer Vision Research , 1 ( 3 ): 1 - 26 [ DOI: 10.1.1.54. 5658 http://dx.doi.org/10.1.1.54.5658 ] Image Features from Phase Congruency Lan C Z , Lu W J , Yu J M and Xu Q . 2021 . Deep learning algorithm for feature matching of cross modality remote sensing images . Acta Geodaetica et Cartographica Sinica , 2021, 50 ( 2 ): 189 - 202 Deep learning algorithm for feature matching of cross modality remote sensing images 蓝朝桢 , 卢万杰 , 于君明 , 徐青 . 2021 . 异源遥感影像特征匹配的深度学习算法 . 测绘学报 , 50 ( 2 ): 189 - 202 [ DOI: 10.11947/j.AGCS.2021.202-00048 http://dx.doi.org/10.11947/j.AGCS.2021.202-00048 ] 异源遥感影像特征匹配的深度学习算法 Li X H , Hui N , Shen H F , Fu Y J and Zhang L P . 2015 . A Robust Mosaicking procedure for High Spatial Resolution Remote Sensing Images [J]. ISPRS Journal of Photogrammetry and Remote Sensing , 109 : 108 - 125 [ DOI: 10.1016/j.isprsjprs.2015.09.009 http://dx.doi.org/10.1016/j.isprsjprs.2015.09.009 ] A Robust Mosaicking procedure for High Spatial Resolution Remote Sensing Images Li Z H and Leung H . 2007 . Contour-based Multisensor Image Registration with Rigid Transformation // 2007 10th International Conference on Information Fusion . Quebec, QC, Canada [ DOI: 10.1109/ICIF:2007.4408109 http://dx.doi.org/10.1109/ICIF:2007.4408109 ] Contour-based Multisensor Image Registration with Rigid Transformation Lowe D G . 2004 . Distinctive Image Features from Scale-Invariant Keypoints . International Journal of Computer Vision , 60 ( 2 ): 91 – 110 [ DOI: 10.1023/B:VISI.0000029664. 99615.94 http://dx.doi.org/10.1023/B:VISI.0000029664.99615.94 ] Distinctive Image Features from Scale-Invariant Keypoints Lu J Z , Hu M Q , Dong J , Han S L and Su A . 2020 . A novel dense descriptor based on structure tensor voting for multi-modal image matching . Chinese Journal of Aeronautics . 33 ( 9 ): 2408 - 2419 [ DOI: 10.1016/j.cja.2020.02.002 http://dx.doi.org/10.1016/j.cja.2020.02.002 ] A novel dense descriptor based on structure tensor voting for multi-modal image matching Martinez A , Garcia J and Abad F . 1999 . A Correlation-Symbolic Approach to Automatic Remotely Sensed Image Rectification // IEEE 1999 International Geoscience and Remote Sensing Symposium . Hamburg, Germany. IEEE: 336 – 338 [ DOI: 10.1109/IGARSS.1999.773490 http://dx.doi.org/10.1109/IGARSS.1999.773490 ] A Correlation-Symbolic Approach to Automatic Remotely Sensed Image Rectification Nan K , Qi H and Ye Y X . 2019 . A template matching method of multimodal remote sensing images based on deep convolutional feature representation . Acta Geodaetica et Cartographica Sinica . 48 ( 6 ): 727 - 736 . A template matching method of multimodal remote sensing images based on deep convolutional feature representation 南轲 , 齐华 , 叶沅鑫 . 2019 . 深度卷积特征表达的多模态遥感影像模板匹配方法 . 测绘学报 , 48 ( 6 ): 727 - 736 [ DOI: 10.10.1194-7/j.AGCS.2019.20180432 http://dx.doi.org/10.10.1194-7/j.AGCS.2019.20180432 ] 深度卷积特征表达的多模态遥感影像模板匹配方法 Pan X C , Jiang T , Yu A Z , Wang X and Zhang Y . 2019 . Geo-positioning of remote sensing images with reference image . Journal of Remote Sensing , 23 ( 4 ): 673 - 684 Geo-positioning of remote sensing images with reference image 潘雪琛 , 姜挺 , 余岸竹 , 王鑫 , 张一 . 2019 . 基准影像数据辅助遥感影像几何定位 . 遥感学报 , 23 ( 4 ): 673 - 684 [ DOI: 10.11834/jrs.20197490 http://dx.doi.org/10.11834/jrs.20197490 ] 基准影像数据辅助遥感影像几何定位 Schmitt M , Tupin F and Zhu X X . 2017 . Fusion of SAR and Optical Remote Sensing Data — Challenges and Recent Trends // IEEE International Geoscience and Remote Sensing Symposium (IGARSS) . Fort Worth, TX : I EEE: 5458 – 5461 [ DOI: 10.1109/IGARSS.2017.8128239 http://dx.doi.org/10.1109/IGARSS.2017.8128239 ] Fusion of SAR and Optical Remote Sensing Data — Challenges and Recent Trends Rosten E and Drummond T . 2006 . Machine Learning for High-Speed Corner Detection // 9 th European Conference on Computer Vision , Graz, Austria: 430 - 443 [ DOI: 10.1007/-11744023_34 http://dx.doi.org/10.1007/-11744023_34 ] Machine Learning for High-Speed Corner Detection Suri S and Reinartz P . 2010 . Mutual-Information-Based Registration of TerraSAR-X and Ikonos Imagery in Urban Areas . IEEE Transactions on Geoscience and Remote Sensing . 48 ( 2 ): 939 – 949 [ DOI: 10.1109/TGRS.2-009.2034-842 http://dx.doi.org/10.1109/TGRS.2-009.2034-842 ] Mutual-Information-Based Registration of TerraSAR-X and Ikonos Imagery in Urban Areas Wang M M , Ye Y X , Zhu B and Zhang G . 2022 . An Automatic Registration Method for Optical and SAR Images Based on Spatial Constraint and Structure Features . Geomatics and Information Science of Wuhan University , 47 ( 1 ): 141 - 148 An Automatic Registration Method for Optical and SAR Images Based on Spatial Constraint and Structure Features 王蒙蒙 , 叶沅鑫 , 朱柏 , 张过 . 2022 . 基于空间几何约束和结构特征的光学影像与SAR影像配准方法 . 武汉大学学报信息科学版 , 47 ( 1 ): 141 - 148 [ DOI: 10.13202/j.whugis201903-54 http://dx.doi.org/10.13202/j.whugis201903-54 ] 基于空间几何约束和结构特征的光学影像与SAR影像配准方法 Xiang Y M , Wang F and You H J . 2018 . OS-SIFT: A Robust SIFT-Like Algorithm for High-Resolution Optical-to-SAR Image Registration in Suburban Areas . IEEE Transactions on Geoscience and Remote Sensing , 56 ( 6 ): 3078 – 3090 [ DOI: 10.1109/TGRS.2018.2790483 http://dx.doi.org/10.1109/TGRS.2018.2790483 ] OS-SIFT: A Robust SIFT-Like Algorithm for High-Resolution Optical-to-SAR Image Registration in Suburban Areas Ye Y X , Bruzzone L , Shan J , Bovolo F and Zhu Q . 2019 . Fast and Robust Matching for Multimodal Remote Sensing Image Registration . IEEE Transactions on Geoscience and Remote Sensing . 57 ( 11 ): 9059 – 9070 [ DOI: 10.1109/TGR-S.2019.2924684 http://dx.doi.org/10.1109/TGR-S.2019.2924684 ] Fast and Robust Matching for Multimodal Remote Sensing Image Registration Ye Y X , Shan J , Bruzzone L and Shen L . 2017 . Robust Registration of Multimodal Remote Sensing Images Based on Structural Similarity . IEEE Transactions on Geoscience and Remote Sensing , 55 ( 5 ): 2941 – 2958 [ DOI: 10.1109/TGRS.2017.2656380 http://dx.doi.org/10.1109/TGRS.2017.2656380 ] Robust Registration of Multimodal Remote Sensing Images Based on Structural Similarity Z Y , C Z G and X Y . 2008 . Multi-spectral remote image registration based on SIFT . Electronics Letters , 44 ( 2 ): 107 - 108 [ DOI: 10.1049/el:20082477 http://dx.doi.org/10.1049/el:20082477 ] Multi-spectral remote image registration based on SIFT Zhu H , Jiao L C , Ma W P , Liu F and Zhao W . 2019 . A Novel Neural Network for Remote Sensing Image Matching . IEEE Transactions on Neural Networks and Learning Systems , 30 ( 9 ): 2853 – 2865 [ DOI: 10.1109/TNNLS.2018. 288757 http://dx.doi.org/10.1109/TNNLS.2018.288757 ] A Novel Neural Network for Remote Sensing Image Matching Zitová B and Flusser J . 2003 . Image Registration Methods: A Survey . Image and Vision Computing . 21 ( 11 ): 977 – 1000 [ DOI: 10.1016/S0262-8856(03)00137-9 http://dx.doi.org/10.1016/S0262-8856(03)00137-9 ] Image Registration Methods: A Survey 10.11834/jrs.20221765 P237 A 叶沅鑫，王蒙蒙，杨超，喻智睿，葛旭明.XXXX.基于方向相位稠密特征的多传感器遥感影像配准方法和系统.遥感学报，XX（XX）： 1-13 叶沅鑫，王蒙蒙，杨超，喻智睿，葛旭明.XXXX.基于方向相位稠密特征的多传感器遥感影像配准方法和系统.遥感学报，XX（XX）： 1-13 YE Yuanxin，WANG Mengmeng，YANG Chao，YU Zhirui，GE Xuming. XXXX. A multi-sensor remote sensing registration method and system based on dense feature of orientated phase. National Remote Sensing Bulletin， XX（XX）：1-13 YE Yuanxin，WANG Mengmeng，YANG Chao，YU Zhirui，GE Xuming. XXXX. A multi-sensor remote sensing registration method and system based on dense feature of orientated phase. National Remote Sensing Bulletin， XX（XX）：1-13 国家自然科学基金（编号：41971281） Supported by Support by National Natural Science Foundation of China (No. 41971281) 2022-03-04 遥感学报