航空滤光片阵列多光谱图像曲面拟合双阈值配准

航空滤光片阵列多光谱图像曲面拟合双阈值配准 Dual-threshold registration based on surface fitting of aerial filter array multispectral image 1 first-author 李铜哨，1998年生，男，硕士研究生， 研究方向为数字图像处理。E-mail： 2423688595@qq.com 李铜哨，1998年生，男，硕士研究生， 研究方向为数字图像处理。E-mail： 2423688595@qq.com 2423688595@qq.com 1 corresp 孙文邦，1976年生，男，博士， 研究方向为数字图像处理、信息保障理论。E-mail： chwenbang@163.com 孙文邦，1976年生，男，博士， 研究方向为数字图像处理、信息保障理论。E-mail： chwenbang@163.com chwenbang@163.com 1 2 1 空军航空大学 航空作战勤务学院, 长春 130022 空军航空大学 航空作战勤务学院, 长春 130022 College of Aviation Combat Service, Air Force Aviation University, Changchun 130022, China College of Aviation Combat Service, Air Force Aviation University, Changchun 130022, China 93057部队, 双辽 136400 93057部队, 双辽 136400 93057 troops, Shuangliao 136400,China 93057 troops, Shuangliao 136400,China 图像配准过程中，匹配点位置精度是决定图像配准精度的关键。本文针对航空滤光片阵列多光谱图像因各谱段间存在像点位移而使误匹配点剔除比较困难的问题，提出了一种基于匹配点位置差曲面拟合双阈值剔除方法。首先，选取多光谱中间波段图像作为基准图像，利用SIFT算法分别提取基准图像和待配准图像的匹配点；其次，在基准图像匹配点处逐点计算两波段图像匹配点的位置差，构建匹配点Delaunay三角网，利用IDW（反距离加权）算法拟合整幅图像位置差曲面；然后，对位置差曲面进行平滑处理，再分别向上向下平移一定容差范围，构成位置差三维阈值空间；最后，利用位置差三维阈值空间筛选出精确匹配点，并完成图像配准工作。理论分析与实验结果表明：该方法可以有效筛选出航空滤光片阵列多光谱图像高精度匹配点，进而有效提高多光谱图像配准的精度。 The high-precision registration between aerial filter array multispectral images is an important guarantee for the subsequent image processing and application. In the process of image registration, the position accuracy of matching points is the key to determine the accuracy of image registration. However, the objects of different strips in the same band image are acquired at different moments, the image displacement between single-band images is large, and the difference of geometric errors of matching points between topographic undulating areas and flat areas in the image is more obvious, the false matching points cannot be accurately eliminated by the global matrix. Aiming at the problem that it is difficult to eliminate the mismatched points in the multi-spectral image of aerial filter array because of the image point displacement between each spectral segment, a new method of double threshold elimination based on the matching point position difference surface fitting is proposed in this paper.Firstly, the intermediate band image of the filter array multispectral images is selected as the reference image, and the matching points in the reference image and the image to be registered are respectively extracted by the sub-pixel-level SIFT algorithm. secondly, the position difference of the matching points of the two bands was calculated point-by-point at the matching points of the benchmark image. and the Delaunay triangulation network of matching points in the reference image was constructed. Then, the position difference surface is smoothed, and the position differences between the matching points of the reference image and the corresponding matching points of the image to be registered are calculated point by point, and then a certain tolerance range is shifted up and down respectively to form a three-dimensional position difference threshold space. Finally, the accurate matching points are selected by using the three-dimensional threshold space of position difference, and the image registration is completed.The three-band composite image of the algorithm-registered image in this paper has clear features and well-defined details，which can meet the requirements of subsequent data processing and application. In order to further illustrate the effectiveness of the proposed algorithm, two datasets of filter array multispectral image are registered and verified from both qualitative and quantitative aspects. False colour composite images: the composite image processed by the proposed algorithm has no obvious pseudo-edges and the features are clear, while there are still obvious pseudo-edges in the comparison algorithm; difference image gray-scale histogram: among the experiments of two datasets，the difference image histogram curve of the proposed algorithm has the largest shift to the left, which indicates that the image registered by the proposed algorithm has the smallest difference with the reference image and the best registration effect.Theoretical analysis and experimental results show that the dual-threshold pointing algorithm based on matching point difference fitting curved surface can effectively screen high precision matching points of aerial filter array multispectral image, and then effectively improve the accuracy of the image registration. The surface fitting to the position difference of matching points reflects the trend of image point displacement in each region of the image, and effectively eliminates the false matching points around the correct matching points according to the feature that the image displacement of same region is approximately same. 曲面拟合 双阈值 滤光片阵列 多光谱图像 配准 Curved surface fitting Double threshold Filter array Multispectral image Registration 1　引 言 1 近年来，多光谱相机已经成为重要的光电遥感载荷，其中，滤光片阵列多光谱相机因其较低的平台需求和简单的数据预处理，被广泛应用到航空无人机多光谱遥感中。但是，航空滤光片阵列多光谱相机的各谱段图像间容易产生像点位移，特别是地面起伏地区，破坏了波段图像间的空间一致性。为满足滤光片阵列式多光谱数据的应用要求，必须对获得的多光谱影像进行高精度配准。 目前，图像配准的主要算法是基于地物特征匹配点进行，而匹配点的位置精度是决定图像配准精度的关键（ Zhang等，2019 Zhang等，2019 年）。因此，在多光谱图像匹配点提取过程中，均需要对初始提取的匹配点进行严格筛选，以剔除影响配准结果的误匹配点和低精度匹配点（ Li等，2021 Li等，2021 年）。 目前，文献中匹配点筛选的算法有很多，例如随机样本一致性（RANSAC）算法（ Fischler等，1981 Fischler等，1981 年； Maria等，2014 Maria等，2014 年； Wang等，2017 Wang等，2017 年）、欧式距离约束法（ Gong等，2013 Gong等，2013 年； Wang等，2017 Wang等，2017 年； Tang等，2020 Tang等，2020 年）、对极几何估计法（ Zhao等，2018 Zhao等，2018 年）等。还有一些学者提出综合运用多种筛选方法，例如孟进军（ Meng等，2020 Meng等，2020 年）等通过RANSAC算法、最小二乘法进行匹配点平差，得到均方根误差平均值小于1的精匹配点；李芳芳（ Li等，2009 Li等，2009 年）等根据匹配点对几何位置去除周围相差较大的匹配点、双向匹配进一步筛选匹配点，最后使用RANSAC算法得到精匹配点等。 但是，目前匹配点筛选算法主要采用全域阈值算法，对于各波段图像的成像中心为同一点时，各波段图像匹配点几何误差的差异比较小，基于全域阈值剔除图像匹配点中的误匹配点或低精度匹配点的效果比较好。 而对于航空滤光片阵列多光谱图像，图像中匹配点的几何误差大致可分为两类：地面起伏区域存在着较大的像点位移，匹配点的几何误差较大；而地面相对平坦区域匹配点的几何误差较小，基于全局阈值匹配点筛选算法无法有效区分正确匹配点和误匹配点。 因此，本文提出了一种基于图像匹配点位置差的曲面拟合自适应双阈值筛选方法，用于滤光片阵列多光谱图像的配准，尤其是复杂地形地区的多光谱图像配准。 2　滤光片阵列多光谱成像特点 1 2.1　成像原理 2 滤光片阵列多光谱相机是以小画幅方式工作，其在相机光路系统中加入了一个滤光板，滤光板放置在CCD探测器前面，在滤光板上沿着垂直于飞行方向镀制了几个窄带的带通滤光膜，每个带通滤光膜只能通过一个谱段的图像，如 图1 图1 所示（图中显示为8条滤光片）。 多光谱面阵探测器的若干行像元对应一个光谱带，整个面阵探测器划分为若干光谱带，曝光由快门控制，每次曝光获得某一区域目标不同波段的条带图像。然后，通过平台运动即可获得一定区域的多光谱图像数据。 多光谱相机的多波段图像均是在同一阵面上成像，需要通过图像拼接的方法才能够形成一幅单波段图像，条带图像裁切拼接如 图2 图2 所示。 2.2　图像特征 2 大部分多光谱相机各波段是同时成像，如线性渐变滤光片（ Wang等，2016 Wang等，2016 年）、单镜头棱镜分光式多光谱相机（ Huang等，2008 Huang等，2008 年），其多光谱图像间的像点位移几乎一致。 但是，滤光片阵列多光谱相机同一波段内不同条带的地物不是同一时刻获得的，不同波段的同一地物也不是同一时刻获得的。由于，在同一时刻相机所成的影像中，各个波段条带图像的入射角各不相同，而拼接后的单波段图像的入射角基本相同。因此，不同波段条带图像的入射角不同，导致同一地物在不同波段条带图像的位置存在一定差异，即存在像点位移，尤其是在地形起伏地区很容易形成较大的像点位移，如 图3 图3 所示。 地平面上一点B，像点在光谱条带图像1和光谱条带图像2上分别为b1、b2；若B处存在高程点A，A点的像点分别为a1、a2。由 图3 图3 可知，b1、a1与a2、b2不重合，存在一定的像点位移。拼接后得到的单波段图像之间的同名像点原本应是b1、b2，由于像点位移的存在，同名像点b1、b2变为a1、a2，导致多光谱各波段图像上的相同位置对应着不同的地物。并且，地面起伏越明显，图像的同名像点位移就越大，反之地面相对平坦，像点位移越小，几乎重合。 设相机感光平面为4008H×4008V，成像高度为18000米，地面高程差为200米，则成像区域边缘像素最大的像点位移为 图4 图4 所示的 δ h ，计算式如 式（1） 式（1） 所示。 <math display="block" id="M1"><msub><mrow><mi>δ</mi></mrow><mrow><mi>h</mi></mrow></msub><mo>=</mo><mfrac><mrow><mi>h</mi></mrow><mrow><mi>H</mi></mrow></mfrac><msub><mrow><mi>r</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911729&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911724&type= 13.80066681 7.70466709 式中， r 0 为 图4 图4 中 oa 0 ，计算式为： <math display="block" id="M2"><msub><mrow><mi>r</mi></mrow><mrow><mn mathvariant="normal">0</mn><mi mathvariant="normal">m</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi></mrow></msub><mo>=</mo><mi>o</mi><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>=</mo><mn mathvariant="normal">2004</mn><mo>×</mo><mi>d</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911741&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911735&type= 33.35866928 3.80999994 代入 式（1） 式（1） 中有， <math display="block" id="M3"><msub><mrow><mi>δ</mi></mrow><mrow><mi>h</mi><mi mathvariant="normal">m</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi></mrow></msub><mo>=</mo><mfrac><mrow><mn mathvariant="normal">200</mn></mrow><mrow><mn mathvariant="normal">18000</mn></mrow></mfrac><mo>×</mo><mn mathvariant="normal">2004</mn><mo>×</mo><mi>d</mi><mo>≈</mo><mn mathvariant="normal">22.27</mn><mi>d</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911750&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911745&type= 53.67866135 7.70466709 则，空间排列上的第一波段图像和第八波段图像的匹配点最大位置差为44.53个像素，如 式（4） 式（4） 所示。 <math display="block" id="M4"><mi mathvariant="normal">Δ</mi><mi>d</mi><mo>=</mo><mn mathvariant="normal">2</mn><msub><mrow><mi>δ</mi></mrow><mrow><mi>h</mi><mi mathvariant="normal">m</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">44.53</mn><mi>d</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911761&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911753&type= 31.83466530 3.89466691 例如，选取某型滤光片阵列多光谱图像中的两幅图像，利用SIFT算法提取匹配点，再将两幅图像匹配点位置分别采用红色与蓝色标记叠加在其中一幅图像上，如 图5 图5 所示。 图5 图5 左上角小图为地面起伏地区的匹配点对比情况放大图，明显可以看出匹配点位置几何关系存在较大的误差；右下角小图为平坦地区的匹配点对比情况放大图，匹配点几乎重合，几何误差较小。 正是由于匹配点位置存在较大差异，且不同的区域差异不同，给匹配点筛选带来极大的困难。 3　滤光片阵列多光谱图像配准过程 1 滤光片阵列多光谱图像的配准技术过程与其他多光谱图像配准流程基本相同，主要包括特征点提取、匹配点筛选、灰度重采样等几个过程，具体的技术过程如 图6 图6 所示。 3.1　特征点提取 2 多光谱图像配准时，通常选用中间波段图像作为基准图像，例如某型多光谱相机设置有490nm、530nm、580nm、610nm、650nm、730nm、840nm、900nm这样8个波段，则选择610nm或650nm波段图像作为基准图像，然后将其他波段图像分别向基准图像进行配准。 本文选择610nm波段图像作为基准图像，然后采用SIFT算法依次提取基准图像和另一幅波段图像的特征点，并进行特征匹配以得到匹配点对。最后，去除匹配点对中的重复匹配点。 3.2　匹配点筛选 2 3.2.1　常用匹配点筛选算法分析 3 目前，常用的匹配点筛选方法采取由粗到精的策略。粗匹配阶段主要基于RANSAC算法剔除特征点对中的误匹配点对，得到粗匹配点对；精匹配阶段主要利用最小二乘法进行匹配点平差，得到精匹配点对（ Yang等，2019 Yang等，2019 年）。其中，最小二乘法进行匹配点平差步骤如下： （1）计算单应变换矩阵一般通过最小二乘法估计单应矩阵 <math id="M5"><mi>H</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911800&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911797&type= 2.53999996 2.62466669 ，使得到基准图像和待配准图像的任意一对匹配点对 p 1 （ x 1 ，y 1 ） 、p 2 （ x 2 ，y 2 ）满足以下条件： <math display="block" id="M6"><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr><mtr><mtd><mn mathvariant="normal">1</mn></mtd></mtr></mtable></mrow></mfenced><mo>=</mo><mi>H</mi><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mtd></mtr><mtr><mtd><mn mathvariant="normal">1</mn></mtd></mtr></mtable></mrow></mfenced><mo>=</mo><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">11</mn></mrow></msub></mtd><mtd><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">12</mn></mrow></msub></mtd><mtd><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">13</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">21</mn></mrow></msub></mtd><mtd><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">22</mn></mrow></msub></mtd><mtd><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">23</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">31</mn></mrow></msub></mtd><mtd><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">32</mn></mrow></msub></mtd><mtd><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">33</mn></mrow></msub></mtd></mtr></mtable></mrow></mfenced><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mtd></mtr><mtr><mtd><mn mathvariant="normal">1</mn></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911803&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911802&type= 52.49333191 13.80066681 （2）匹配点平差利用计算得到的单应矩阵 H 预测待配准影像中任意匹配点在参考影像中的坐标值。根据预测值与实际匹配点分别计算 x 和 y 方向的差值和均方根误差（Root Mean Square，RMS）。剔除均方根大于阈值的点对，得到正确的精匹配点对。 这种方法对于各波段图像的成像中心为同一点（即各波段图像同时成像）的多光谱相机而言，其图像匹配点的几何位置大致相同，同名像点之间基本满足单个矩阵变换关系，剔除误匹配点或低精度匹配点的效果较好。但对于滤光片阵列多光谱图像而言，由于其波段图像间存在较大的像点位移，匹配点位置存在差异，在地形起伏地区的欧式距离较大，而平坦区域匹配点坐标大致相同。因此，存在像点位移的滤光片阵列多光谱图像无法使用单个矩阵准确描述匹配点之间的位置坐标转换关系。 例如，某一滤光片阵列多光谱图像，采用全局匹配点计算单应矩阵 H ，并通过 式（5） 式（5） 计算待配准图像匹配点在基准图像上预测点的位置，如 图7 图7 所示。 图7 图7 中红色点为待配准图像匹配点的预测点位置，蓝色点为基准图像匹配点位置。从 图7 图7 可以明显看出，预测点位置与基准图像匹配点偏差较大，将会在匹配点平差过程中筛选不出精准匹配点，导致图像配准失败。因此，利用全局匹配点计算的单应矩阵 H 区分正确匹配点和误匹配点的方法并不能用于滤光片阵列多光谱图像。 另外，有学者提出：直接将待配准多光谱图像中匹配点叠加到基准图像中，并利用 式（6） 式（6） 计算待配准图像和基准图像匹配点位置距离差，然后采取全域阈值消除位置距离误差较大的匹配点（ Liang等，2015 Liang等，2015 年）。 <math display="block" id="M7"><mi mathvariant="normal">Δ</mi><msub><mrow><mi>d</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mroot><mrow><msup><mrow><mfenced separators="|"><mrow><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>-</mo><msub><mrow><mi>x</mi></mrow><mrow><mi mathvariant="normal">j</mi></mrow></msub></mrow></mfenced></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><mo>+</mo><msup><mrow><mfenced separators="|"><mrow><msub><mrow><mi>y</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>-</mo><msub><mrow><mi>y</mi></mrow><mrow><mi mathvariant="normal">j</mi></mrow></msub></mrow></mfenced></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mrow><mrow/></mroot></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911818&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911813&type= 45.38133240 7.11199999 式（6） 式（6） 中（ x i ， y i ）为基准图像中匹配点坐标，（ x j ， y j ）为待配准图像中对应的匹配点坐标。 这种方法在一定程度上解决了利用单应变换矩阵和匹配点平差方法筛选匹配点的固有缺陷，有效保留了平坦地区的匹配点，但是容易将地形起伏较大区域的匹配点整体消除，使得地形起伏较大区域无法进行图像配准。不失一般性，采用位置距离差一维数据来说明，如 图8 图8 所示。 从 图8 图8 可以看出，平坦区域的B点与周围匹配点的位置差存在明显差异，应该是误匹配点，采用全域阈值1可以有效地剔除，但是对于地形起伏地区的匹配点，由于成像机理造成匹配点位置差大于平坦地区，采用全域阈值1会将地形起伏地区的匹配点全部剔除。 若要保留地形起伏地区的匹配点只能提高阈值，如采用全域阈值2，则可以很好地保留地形起伏地区的匹配点，但地形平坦地区的误匹配点B则无法剔除。 因此，学者提出通过匹配点的几何位置去除周围相差较大的匹配点（ Li等，2009 Li等，2009 年）。这种方法很好地解决了利用全域阈值筛选匹配点的缺点，能够有效保留图像上地形较高区域和平坦区域上的匹配点，但对于地形坡面处的匹配点筛选问题仍然没有得到很好的解决。 例如 图8 图8 中，A点和D点处在地形起伏区域、B点处在平坦区域，A、B和D三点与周围匹配点的位置差都存在明显差异，应该属于误匹配点，按照匹配点几何位置去除周围相差较大的匹配点的方法可有效地消除A、B和D三个误匹配点。 图8 图8 中C、E处在地形坡面处，从匹配点位置差的变化趋势来看，C和E也应该是误匹配点，但是采取匹配点的几何位置去除周围相差较大的匹配点的方法，却很难消除。 3.2.2　曲面拟合双阈值匹配点筛选算法 3 根据3.2.1节的分析，传统的剔点方式无法有效剔除滤光片阵列多光谱图像中的误匹配点，因此，本文提出曲面拟合双阈值匹配点筛选算法。图像中的不同区域可根据匹配点的距离差拟合出不同的曲面，利用拟合后的曲面对各匹配点进行剔点，这个过程是自适应的，即各个匹配点根据拟合曲面上所对应点的容差范围确定各自的阈值空间。基本思路如 图9 图9 所示。 利用待配准多光谱图像匹配点叠加在基准图像上计算得到的匹配点位置距离差，拟合插值得到距离差曲面，并对拟合曲面进行平滑处理，如 图9 图9 中黑色双点划线所示。然后对平滑后曲面向上平移一定容差单位，得到上限阈值曲面，如 图9 图9 中蓝色实线所示。所有高于上限阈值曲面的匹配点可以认为是误匹配点。利用上限阈值曲面可以很好消除地形平坦地区的误匹配点B、地形起伏地区的误匹配点A，以及地形坡面上的误匹配点C。但是对于地形起伏地区的误匹配点D、坡面误匹配点E则无法消除。 因此，将平滑后的曲面向下平移一定容差单位，得到下限阈值曲面，如 图9 图9 中红色虚线所示。同理，低于下限阈值曲面的匹配点可以认为是误匹配点，故使用下限阈值曲面可消除匹配点D和点E。 上限阈值曲面和下限阈值曲面构成了一个三维阈值空间。简单地说，所有落入三维阈值空间内部的匹配点位置差可以判定为合理范围内的正确匹配点，而落在三维阈值空间外部的匹配点则判定为误匹配点。 曲面拟合双阈值匹配涉及到关键技术主要包括曲面插值拟合、曲面平滑、阈值空间确定、阈值空间判断4个方面。 （1）曲面插值拟合曲面插值拟合算法较多，例如二元三次插值、二元三次样条插值、反距离加权插值等。如 图10 图10 所示的是滤光片阵列多光谱局部图像采取二元三次 4 插值法生成的距离差曲面，色度条颜色代表着拟合曲面上的像点位移大小。 由 图10 图10 可知，距离差插值拟合结果在局部区域出现相邻点的像点位移变化较大，“上下剧烈抖动”的现象，如 图10 图10 中的蓝色虚线圆圈区域，甚至出现了负值，拟合曲面明显不符合像点位移变化趋势，不利于确定三维阈值空间。 因此，本文参考反距离加权插值方法（ Chen等，2019 Chen等，2019 年）生成距离差曲面，使插值结果均在最大值与最小值之间。首先，利用基准图像匹配点位置构建Delaunay三角网，然后在每个三角形内采取反距离加权插值方法填充，如 图11 图11 所示。 图11 图11 中，蓝色点A、B、C为基准图像匹配点中某个Delaunay三角网顶点，而红色点A΄、B΄、C΄为待配准图像匹配点。通过分别计算A与A΄、B与B΄、C与C΄之间的位置距离差，得到 m 1 、 m 2 、 m 3 。若三角形中某一点距顶点A、B和C的距离分别为 d 1 、 d 2 、 d 3 ，则该点的距离差曲面值为： <math display="block" id="M8"><mi mathvariant="normal">Δ</mi><mi>D</mi><mo>=</mo><mfrac><mrow><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mrow></mfrac></mrow><mrow><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub></mrow></mfrac></mrow></mfrac><msub><mrow><mi>m</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>+</mo><mfrac><mrow><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mrow></mfrac></mrow><mrow><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub></mrow></mfrac></mrow></mfrac><msub><mrow><mi>m</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo>+</mo><mfrac><mrow><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub></mrow></mfrac></mrow><mrow><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>d</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub></mrow></mfrac></mrow></mfrac><msub><mrow><mi>m</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911859&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911857&type= 94.40332794 19.21933365 基于Delaunay三角网利用反距离加权插值生成的距离差曲面如 图12 图12 所示，色度条中不同颜色代表着拟合曲面上的像点位移。 difference surface 2）曲面平滑一般情况下，多光谱图像中的误匹配点分布比较稀疏，数量较少。因此，为了防止平滑过程中，误匹配点影响到其他匹配点，本文选择3×3均值滤波器进行平滑处理。平滑效果图如 图13 图13 所示。 从 图13 图13 中可以明显看出，有许多类似于 图8 图8 中的A、B、C误匹配点（图中红色圆圈标记显示）。若转换曲面下方视角，如 图14 图14 所示，则从 图14 图14 中也可以明显看出类似 图9 图9 中的D、E误匹配点（图中红色圆圈标记显示）。 （3）阈值空间确定如 图9 图9 所示，以平滑后距离差曲面为基准，分别上下浮动一定的容差范围，形成三维阈值空间。其中，阈值范围的选取一定程度上决定了图像中误匹配点剔除的效果优劣，若阈值的范围选择过大，则无法有效剔除图像匹配点中的误匹配点。反之，保留下的正确匹配点无法满足多光谱图像配准精度的需求。 因为本文在提取多光谱图像的特征点时采用了亚像素级的SIFT算法，故特征点的提取精度为±0.5像素。 因此，本文在平滑后的距离差曲面加上±0.5构成三维阈值空间。注：若在平坦区域，加上－0.5可能使下限阈值曲面低于0，但不影响后续的阈值的判断，如 图15 图15 所示。 （4）阈值空间判断SIFT算法是亚像素级的特征点提取算法，其特征匹配点坐标一般不为整数。例如， 图16 图16 所示的匹配点A、B和C。图中红色虚边圆点为匹配点距离差值 d i ，蓝色圆点为上限阈值 U i ，青色圆点为下限阈值 D i 。 根据上限阈值曲面分别插值计算各匹配点的上限阈值 U i ，根据下限阈值曲面分别插值计算各匹配点的下限阈值 D i ，然后利用 式（8） 式（8） 进行判断距离差值 Y i 是否为正确匹配点。 <math display="block" id="M9"><msub><mrow><mi>d</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><mi mathvariant="normal">正</mi><mi mathvariant="normal">确</mi><mi mathvariant="normal">匹</mi><mi mathvariant="normal">配</mi><mi mathvariant="normal">点</mi><mtext>   </mtext><msub><mrow><mi>D</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>≤</mo><msub><mrow><mi>d</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>≤</mo><msub><mrow><mi>U</mi></mrow><mrow><mi>i</mi></mrow></msub></mtd></mtr><mtr><mtd><mi mathvariant="normal">错</mi><mi mathvariant="normal">误</mi><mi mathvariant="normal">匹</mi><mi mathvariant="normal">配</mi><mi mathvariant="normal">点</mi><mtext>   </mtext><msub><mrow><mi>D</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>≥</mo><msub><mrow><mi>d</mi></mrow><mrow><mi>i</mi></mrow></msub><mi mathvariant="normal">或</mi><msub><mrow><mi>d</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>≥</mo><msub><mrow><mi>U</mi></mrow><mrow><mi>i</mi></mrow></msub></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911890&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911889&type= 54.69466782 10.15999985 3.3　灰度重采样 2 本文采用Delaunay三角小面元反距离加权的方法进行重采样，完成图像配准。 首先，分别建立基准图像与待配准图像的Delaunay三角网，然后对每一个三角形依次进行灰度重采样，如 图17 图17 所示。 例如，对其中某一对匹配三角形A 1 B 1 C 1 和A 2 B 2 C 2 ，首先，根据基准图像中三角形A 1 B 1 C 1 的顶点坐标（ x 11 ， y 11 ）、（ x 12 ， y 12 ）、（ x 13 ， y 13 ），待配准图像中三角形A 2 B 2 C 2 的顶点坐标（ x 21 ， y 21 ）、（ x 22 ， y 22 ）、（ x 23 ， y 23 ），利用 式（9） 式（9） 计算出基准图像与待配准图像之间的几何位置转换矩阵 H 。 <math display="block" id="M10"><mi>H</mi><mo>=</mo><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">21</mn></mrow></msub></mtd><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">22</mn></mrow></msub></mtd><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">23</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">21</mn></mrow></msub></mtd><mtd><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">22</mn></mrow></msub></mtd><mtd><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">23</mn></mrow></msub></mtd></mtr><mtr><mtd><mn mathvariant="normal">1</mn></mtd><mtd><mn mathvariant="normal">1</mn></mtd><mtd><mn mathvariant="normal">1</mn></mtd></mtr></mtable></mrow></mfenced><msup><mrow><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">11</mn></mrow></msub></mtd><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">12</mn></mrow></msub></mtd><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">13</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">11</mn></mrow></msub></mtd><mtd><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">12</mn></mrow></msub></mtd><mtd><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">13</mn></mrow></msub></mtd></mtr><mtr><mtd><mn mathvariant="normal">1</mn></mtd><mtd><mn mathvariant="normal">1</mn></mtd><mtd><mn mathvariant="normal">1</mn></mtd></mtr></mtable></mrow></mfenced></mrow><mrow><mo>-</mo><mn mathvariant="normal">1</mn></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911895&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911894&type= 54.10200119 14.39333248 然后，对于基准图像三角形内某一点 P 1 （ x 1 ， y 1 ），利用 式（10） 式（10） ，求得其在待配准图像中对应位置 P 2 （ x 2 ， y 2 ）。 <math display="block" id="M11"><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mtd></mtr><mtr><mtd><mn mathvariant="normal">1</mn></mtd></mtr></mtable></mrow></mfenced><mo>=</mo><mi>H</mi><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr><mtr><mtd><mn mathvariant="normal">1</mn></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911897&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=22911896&type= 19.30400085 13.80066681 最后，在待配准图像三角形A 2 B 2 C 2 中利用反距离加权插值法得到 P 2 点的灰度（ Feng等，2018 Feng等，2018 年）， 并将其赋给 P 1 点。如果每一个三角形中每一点依次插值，就完成了灰度插值工作。 4　实验与分析 1 为了验证本文提出算法的可行性，本文采用了某机载滤光片阵列多光谱相机拍摄某地区的山区多光谱图像。其中，滤光片阵列多光谱相机的探测器采用的是加拿大DALSA公司的全帧转移（Full Frame Transfer）型CCD探测器。 选择其中580nm、610nm和650nm三个波段两组图像进行验证实验，实验图像大小为1977像素×2034像素，如 图18 图18 所示。 1,2,3;4,5,6; 由于多光谱相机的成像区域存在地形起伏，580nm、610nm和650nm波段的多光谱图像之间存在着较大的像点位移，如果直接对三个波段的多光谱图像进行假彩色合成，则合成后的假彩色影像中的图像目标重影（伪边缘）十分明显，不利于多光谱图像后期的应用处理。例如数据一的合成结果，如 图19 图19 所示。 选择610nm波段图像作为基准图像，首先，利用SIFT算法提取多光谱图像的特征点（ Sourabh等，2021 Sourabh等，2021 年），对图像特征点进行特征匹配，并去除匹配点中的多余重复点，得到初始匹配点，如 图20 图20 所示。 1,2;3,4; 图20 图20 中，去除重复点后的初始匹配点以红色圆点表示。其中，图 20（a） 20（a） 、 20 20 分别为数据一和数据二的610nm和580nm波段图像初始匹配点在基准图像上的空间分布情况；图 20（b） 20（b） 、 20 20 分别为数据一和数据二的610nm和650nm波段图像初始匹配点在基准图像上的空间分布情况。 基于上述四组初始匹配点分别拟合匹配点距离差插值拟合曲面，如 图21 图21 所示。其中，图 21（a） 21（a） 、 21 21 分别为数据一的拟合曲面，图 21（c） 21（c） 、 21 21 为数据二的拟合曲面，色度条颜色代表拟合曲面上各区域对应的像点位移。 1;2;3;4; 由 图21 图21 可知，同一区域的多光谱图像构建的距离差拟合曲面形状大体一致，但不同波段的图像间仍存在差异；不同区域的多光谱图像所构建的距离差拟合曲面差异较大。 以数据一拟合的曲面为例， 图21（a） 图21（a） 和 图21（b） 图21（b） 中均可以明显看出匹配点中存在许多的误匹配点，如 图21（a） 图21（a） 和 图21（b） 图21（b） 中红色圆圈标记的匹配点的距离差值，红色圆圈内的匹配点与周围匹配点的距离差存在着明显差异。 为了更好地说明本文算法的匹配点筛选性能，选择了目前常用的两种匹配点筛选算法（算法一：RANSAC算法；算法二：（ Liang等，2015 Liang等，2015 年）的控制点选取算法），与本文提出的匹配点筛选算法进行对比实验。算法一、算法二与本文算法处理后的精确匹配点筛选结果如 图22 图22 所示。 1,2,3;4,5,6;7,8,9;10,11,12; 由 图22 图22 可知，算法一（RANSAC算法）处理滤光片阵列多光谱影像时，在剔除误匹配点的同时将大量正确匹配点也剔除掉，导致影像部分区域无匹配点。 算法二（ （Liang等，2015 （Liang等，2015 年）中算法），剔除了大量地形起伏地区的图像匹配点，剔除后的点无法满足配准精度要求。 3种算法的匹配点数量对比情况，如 表1 表1 所示。 2 数据 2 算法 3 610nm和650nm 谱段配准 3 610nm和580nm 谱段配准 初始 点对 精匹配点对 匹配率 初始 点对 精匹配点对 匹配率 3 数据一 算法一 10592 5652 53.4% 12551 10093 80.4% 算法二 10592 3303 31.2% 12551 4012 32.0% 本文算法 10592 10580 99.9% 12551 12540 99.9% 3 数据二 算法一 11666 3266 28.0% 13313 5703 42.8% 算法二 11666 4767 40.9% 13313 4575 34.4% 本文算法 11666 11629 99.7% 13313 13270 99.7% 通过数据比对发现，算法二剔除点的数量较多，部分区域存在无匹配点现象，且较为明显；算法一剔除点数量介于两者之间，但也存在部分区域无匹配点的问题。本文算法保留匹配点数量最多且匹配点分布在影像各处，能够保证影像配准的精度需要。 因此，算法二以匹配点对 x 方向和 y 方向的误差均值作为阈值，对地形起伏地区而言，在剔除误匹配点的同时，像点位移较大的正确点对也会被剔除。这种方法既减少了匹配点对数量，又降低了影像配准精度。 算法一（RANSAC）中，基于全局匹配点选取的样本模型，会把匹配点中的大量正确点剔除（ Shao等，2019 Shao等，2019 年），并不适合滤光片阵列式多光谱图像。 本文算法，筛选匹配点时既考虑了周围点对，又剔除了不符合点对周围区域像点位移变化的点，有效解决了基于局部区域精筛匹配点的难题。 以数据一的610nm和580nm谱段配准为例，剔除的点对如 图24 图24 所示。其中，右上角图为地形起伏地区，右下角图为平坦区域。 图23 图23 中白色点和绿色点是本文算法剔除的点。 1;2;3;4; 从目视角度上看，右上角图中剔除的点对距离差明显大于周围点对，而周围匹配点对距离差近似一致，说明该点是错误点。 对于右下角图中的两对匹配点，左侧几乎重合的点位于平坦区域，其周围匹配点距离差较小，但该点距离差远小于周围点对，而右侧点的距离差远大于周围匹配点距离差，并不符合周围点对距离差的变化趋势，需要剔除掉。 另外，基于以上得到的匹配点对，对四组影像进行三角小面元配准（ Zheng等，2019 Zheng等，2019 年），并分别从基准图像与配准后图像的差异图像灰度直方图和假彩色合成图像两个方面对三种算法进行说明，以比较三种算法的优劣。 差异图像灰度直方图：对四组实验中基准图像与配准后的图像分别进行作差并取绝对值，得到基准图像与配准后图像的差异图像。 求出差异图像局部区域灰度直方图并以曲线形式显示，如 图24 图24 所示。图 24（a） 24（a） 、 24 24 分别为数据一、数据二的610nm波段和配准后580nm波段差异图像局部区域灰度直方图，图 24（b） 24（b） 、 24 24 分别为数据一、数据二的610nm波段和配准后650nm波段差异图像局部区域直方图。 一般情况下，差异图像灰度直方图向左下偏移程度越大，说明差异图像的灰度值越小，图像的配准效果越准确；直方图越向右上偏移，配准效果越差。因此，本文采用差异图像直方图的偏移情况来判断图像配准结果的好坏，进而说明匹配点筛选算法的优劣。 由 图24 图24 的（a）、（b）、（c）和（d）可以看出，两组图像的配准实验中原始影像之间的位置偏移较大，差异图像直方图均是最靠右上侧。而算法一和算法二的差异图像灰度直方图相对原始影像数据明显向左下偏移，说明算法一和算法二能够起到很好的配准效果。 本文算法相对于算法一和算法二的差异图像直方图，又明显向左下偏移，特别是在610nm-650nm图像配准上效果更为明显。说明本文算法配准效果明显优于算法一和算法二，进而说明本文的匹配点筛选算法效果优于算法一和算法二。 假彩色合成图像：对两组实验中配准后的局部影像分别合成假彩色图像，如 图25 图25 所示。 1,2,3,4;5,6,7,8; 从目视效果上看，未进行配准的原始影像合成后的伪彩色图像中地形起伏地区伪边缘较为严重，算法一和算法二处理后的合成图像未解决伪边缘问题，但伪边缘较原始数据有所减轻。 而本文算法处理后的合成影像中，地物清晰，细节更明确，没有明显的伪边缘现象，说明本文算法能够较好的处理滤光片阵列多光谱图像配准问题。 5　结 论 1 针对机载滤光片阵列式多光谱影像准确选取匹配点的难题，提出了一种基于匹配点坐标的自适应曲面拟合双阈值配准方法，即根据匹配点像点位移拟合的曲面剔除在阈值范围外的控制点，降低了滤光片阵列式多光谱影像中误匹配点、低精度匹配点对影像配准结果的影响。通过实验对比，表明本文所提出的算法较目前常用算法有更好的去除误匹配点效果，较好解决了复杂地形的滤光片阵列式多光谱图像精确配准问题，配准后的影像可以满足影像的后期处理及应用需求。 通过对目前用于遥感图像筛选匹配点的算法进行分析与验证，发现使用全局矩阵无法有效剔除像点位移较大的图像中的误匹配点和低精度匹配点，或剔除误匹配点的同时将大量正确匹配点剔除，无法满足图像配准精度的要求。本文算法构建的位置差拟合曲面可以反映整幅图像像点位移的变化趋势，根据同一区域的匹配点像点位移近似一致的特点，可以有效剔除正确匹配点周围的误匹配点。 阈值的选取一定程度上决定了图像中误匹配点剔除的效果优劣，而本文选择0.5阈值取得了较好的实验效果，却并不能保证完全适用所有图像，针对不同类型的图像数据，怎样快速选择较为合适的阈值是下步研究的重点。此外，拟合曲面部分耗时较长，在图像数据较大时，运行速度更为明显，通过GPU加速、优化程序也是下一步的研究工作。 参考文献（ References ） Chen C F and Li Y Y . 2019 . A Fast Global Interpolation Method for Digital Terrain Model Generation from Large LiDAR-Derived Data. Remote Sensing, 11(11) (DOI 10.3390/rs11111324 http://dx.doi.org/10.3390/rs11111324 ) A Fast Global Interpolation Method for Digital Terrain Model Generation from Large LiDAR-Derived Data. Remote Sensing, 11(11) (DOI Fischler M A and Bolles R C . 1981 . Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography . Commun. Assoc. Comp. March , 24 ( 6 ): 381 - 395 . (DOI 10.1145/358669.358692 http://dx.doi.org/10.1145/358669.358692 ) Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography Feng R , Li X , Zou W , and H Shen . 2018 . Registration of multitemporal GF-1 remote sensing images with weighting perspective transformation model . IEEE , 2018. [DOI 10.1109/ICIP.2017.8296685 http://dx.doi.org/10.1109/ICIP.2017.8296685 ] Registration of multitemporal GF-1 remote sensing images with weighting perspective transformation model Gong Z C , Pei J Ha and Xie W X . 2013 . Research of Accurate Registration Method for Multispectral Remote Sensing Images . Journal of Signal Processing , 29 ( 10 ): 1398 - 1406 Research of Accurate Registration Method for Multispectral Remote Sensing Images 龚志成 , 裴继红 , 谢维信 . 2013 . 多光谱遥感卫星图像的精确配准方法研究信号处理 , 29 ( 10 ): 1398 - 1406 [ DOI: 10.3969/j.issn.1003- 0530.2013.10.018 http://dx.doi.org/10.3969/j.issn.1003-0530.2013.10.018 ] 多光谱遥感卫星图像的精确配准方法研究信号处理 Huang Y , Lan Y and Hoffmann W C . 2008 . Use of Airborne Multi-Spectral Imagery in Pest Management Systems . Agricultural Engineering International , 2008 . Use of Airborne Multi-Spectral Imagery in Pest Management Systems Liang J X , Wan X X and Liu Q . 2015 . Research on Visible Broadband Spectral Image Registration Method Based on SIFT Algorithm . Journal of Hunan University of Technology , 29 ( 02 ): 57 - 63+93 Research on Visible Broadband Spectral Image Registration Method Based on SIFT Algorithm 梁金星 , 万晓霞 , 刘强 . 2015 . 基于SIFT算法的可见光宽带光谱图像配准方法研究 . 湖南工业大学学报 , 29 ( 02 ): 57 - 63+93 [ DOI: 10.3969/j.issn.1673-9833.2015.02.011 http://dx.doi.org/10.3969/j.issn.1673-9833.2015.02.011 ] 基于SIFT算法的可见光宽带光谱图像配准方法研究 Li F F , Xiao B L , Jia Y H and Mao X L . 2009 . Improved SIFT Algorithm and Its Application in Automatic Registration of Remotely sensed Imagery . Geomatics and Information Science of Wuhan University , 34 ( 10 ): 1245 - 1249 Improved SIFT Algorithm and Its Application in Automatic Registration of Remotely sensed Imagery 李芳芳 , 肖本林 , 贾永红 , 毛星亮 . 2009 . SIFT算法优化及其用于遥感影像自动配准 . 武汉大学学报(信息科学版) , 34 ( 10 ): 1245 - 1249 [ DOI: CNKI:SUN:WHCH.0.2009-10-024 http://dx.doi.org/CNKI:SUN:WHCH.0.2009-10-024 ] SIFT算法优化及其用于遥感影像自动配准 Li Q , Han G , Liu P , Yang H and Wu J . 2021 . An Infrared-Visible Image Registration Method Based on the Constrained Point Feature . Sensors , 21 ( 4 ): 1188 . (DOI 10.3390/s21041188 http://dx.doi.org/10.3390/s21041188 ) An Infrared-Visible Image Registration Method Based on the Constrained Point Feature Maria V and Konstantinos K . 2014 . Automatic Descriptor-Based Co-Registration of Frame Hyperspectral Data . Remote Sensing , 6 ( 4 ): 3409 - 3426 . (DOI 10.1145/358669.358692 http://dx.doi.org/10.1145/358669.358692 ) Automatic Descriptor-Based Co-Registration of Frame Hyperspectral Data Meng J J , Wu J Q , Lu L L , Li Q T , Zhang Q , Feng S Y and Yan J . 2020 . A Full-Spectrum Registration Method for Zhuhai-1 Satellite Hyperspectral Imagery . Sensors , 20 ( 21 ): 6298 . (DOI 10.3390/s20216298 http://dx.doi.org/10.3390/s20216298 ) A Full-Spectrum Registration Method for Zhuhai-1 Satellite Hyperspectral Imagery Shao C , Zhang C , Fang Z and Yang G . 2019 . A Deep Learning-based Semantic Filter for RANSAC-based Fundamental Matrix Calculation and the ORB-SLAM System . IEEE Access, PP ( 99 ): 1 - 1 . A Deep Learning-based Semantic Filter for RANSAC-based Fundamental Matrix Calculation and the ORB-SLAM System DOI 10 . 1109/ ACCESS .2019. 2962268 1109/ ACCESS Sourabh Paul and Umesh C . Pati . 2021 . A comprehensive review on remote sensing image registration. International Journal of Remote Sensing , 42 ( 14 ). [DOI 10.1109/ICIP.2017.8296685 http://dx.doi.org/10.1109/ICIP.2017.8296685 ] Pati Tang Z J , Ma G Q , Lu J Q , Wang Z , Fu B and Wang Y J . 2020 . Sonar image mosaic based on a new feature matching method . IET Image Processing , 2020, 14(10). [DOI 10.1049/IET-IPR. 2019. 0695 ] Sonar image mosaic based on a new feature matching method Wang G , Li J N , Su Q T , Zhang X F , Lv G H and Wang H G . 2017 . Algorithm based on morphological component analysis and scale-invariant feature transform for image registration . Journal of Shanghai Jiaotong University (Science) , 22 ( 1 ): 99 – 106 . (DOI 10.1007/s12204-017-1807-7 http://dx.doi.org/10.1007/s12204-017-1807-7 ) Algorithm based on morphological component analysis and scale-invariant feature transform for image registration Wang J J , Shi Y and Liu H H . 2017 . Registration of frame hyperspectral images based on UAV . Chinese Journal of Agricultural Resources and Regional Planning , 38 ( 09 ): 63 - 71 Registration of frame hyperspectral images based on UAV 王晶晶 , 史云 , 刘含海 . 2017 . 无人机载框幅式高光谱影像的波段配准研究 . 中国农业资源与区划 , 38 ( 09 ): 63 - 71 [ DOI: 10.7621/cjarrp.1005-9121.20170910 http://dx.doi.org/10.7621/cjarrp.1005-9121.20170910 ] 无人机载框幅式高光谱影像的波段配准研究 Wang Y and Gong Y . 2016 . Design of Multispectral Imaging Spectrometer Using Linear Variable Filter . Laser & Optoelectronics Progress , 53 ( 01 ): 247 - 253 Design of Multispectral Imaging Spectrometer Using Linear Variable Filter 王颖 , 巩岩 . 2016 . 线性渐变滤光片型多光谱成像光谱仪设计 . 激光与光电子学进展 , 53 ( 01 ): 247 - 253 [ DOI: CNKI:SUN:JGDJ. 0.2016-01-035 http://dx.doi.org/CNKI:SUN:JGDJ.0.2016-01-035 ] 线性渐变滤光片型多光谱成像光谱仪设计 Yang H , X Li , Zhao L and Chen . 2019 . A Novel Coarse-to-Fine Scheme for Remote Sensing Image Registration Based on SIFT and Phase Correlation . Remote Sensing , 11 ( 15 ): 1833 . (DOI 10.3390/rs11151833 http://dx.doi.org/10.3390/rs11151833 ) A Novel Coarse-to-Fine Scheme for Remote Sensing Image Registration Based on SIFT and Phase Correlation Zhang Z , Han D Q , Dezert Jean and Yang Y . 2019 . A New Image Registration Algorithm Based on Evidential Reasoning . Sensors , 19 ( 5 ). [DOI 10.3390/s19051091 http://dx.doi.org/10.3390/s19051091 ] A New Image Registration Algorithm Based on Evidential Reasoning Zhao H R and Lu S H . 2018 . Dense High-definition Image Matching Strategy Based on Scale Distribution of Feature and Geometric Constraint . Acta Geodaetica et Cartographica Sinica , 47 ( 06 ): 790 - 798 Dense High-definition Image Matching Strategy Based on Scale Distribution of Feature and Geometric Constraint 赵红蕊 , 陆胜寒 . 2018 . 基于特征尺度分布与对极几何约束的高清影像快速密集匹配方法 . 测绘学报 , 47 ( 06 ): 790 - 798 [ DOI: CNKI:SUN:CHXB.0.2018-06-012 http://dx.doi.org/CNKI:SUN:CHXB.0.2018-06-012 ] 基于特征尺度分布与对极几何约束的高清影像快速密集匹配方法 Zheng J , Xu Q Z , Zhai B and Y W . 2019 . Accurate hyperspectral and infrared satellite image registration method using structured topological constraints. Infrared Physics & Technology, 104 : 103122. [ DOI 10.1016/ j.infrared.2019. 103122 ] Accurate hyperspectral and infrared satellite image registration method using structured topological constraints. Infrared Physics & Technology, 致谢:衷心感谢学术编辑和审稿专家的宝贵意见和建设性建议，以及在论文撰写过程中给予我帮助的人！ 10.11834/jrs.20211432 A 李铜哨，孙文邦，岳广，张星铭，顾子侣.XXXX.航空滤光片阵列多光谱图像曲面拟合双阈值配准.遥感学报，XX（XX）： 1-15 李铜哨，孙文邦，岳广，张星铭，顾子侣.XXXX.航空滤光片阵列多光谱图像曲面拟合双阈值配准.遥感学报，XX（XX）： 1-15 LI Tongshao，SUN Wenbang，YUE Guang，Zhang Xingming，Gu Zilv. XXXX. Dual-threshold registration based on surface fitting of aerial filter array multispectral image. National Remote Sensing Bulletin， XX（XX）：1-15 LI Tongshao，SUN Wenbang，YUE Guang，Zhang Xingming，Gu Zilv. XXXX. Dual-threshold registration based on surface fitting of aerial filter array multispectral image. National Remote Sensing Bulletin， XX（XX）：1-15 全军军事类研究生课题（编号： DSSQ910252018010） Supported by Supported by Military Postgraduate Projects of the Army of China (No.DSSQ910252018010) 2020-06-23 2021-12-16 遥感学报