CN113256653B - 一种面向高层地物的异源高分遥感影像配准方法 - Google Patents
一种面向高层地物的异源高分遥感影像配准方法 Download PDFInfo
- Publication number
- CN113256653B CN113256653B CN202110570103.1A CN202110570103A CN113256653B CN 113256653 B CN113256653 B CN 113256653B CN 202110570103 A CN202110570103 A CN 202110570103A CN 113256653 B CN113256653 B CN 113256653B
- Authority
- CN
- China
- Prior art keywords
- image
- rise
- shadow
- ground object
- feature
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 58
- 230000009466 transformation Effects 0.000 claims abstract description 59
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 claims abstract description 17
- 238000001514 detection method Methods 0.000 claims abstract description 15
- 238000012216 screening Methods 0.000 claims abstract description 14
- 238000012937 correction Methods 0.000 claims abstract description 10
- 239000013598 vector Substances 0.000 claims description 26
- 238000000605 extraction Methods 0.000 claims description 11
- 230000008569 process Effects 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 8
- 230000011218 segmentation Effects 0.000 claims description 8
- 230000003044 adaptive effect Effects 0.000 claims description 5
- 238000001914 filtration Methods 0.000 claims description 5
- 238000013507 mapping Methods 0.000 claims description 5
- 230000000877 morphologic effect Effects 0.000 claims description 4
- 238000013519 translation Methods 0.000 claims description 4
- 230000001174 ascending effect Effects 0.000 claims description 3
- 230000035945 sensitivity Effects 0.000 claims description 3
- 238000012360 testing method Methods 0.000 claims description 3
- 230000003287 optical effect Effects 0.000 description 11
- 238000013461 design Methods 0.000 description 5
- 238000002474 experimental method Methods 0.000 description 5
- 230000008859 change Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 230000008030 elimination Effects 0.000 description 3
- 238000003379 elimination reaction Methods 0.000 description 3
- 238000003384 imaging method Methods 0.000 description 3
- 230000005855 radiation Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000002457 bidirectional effect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 241000282461 Canis lupus Species 0.000 description 1
- UDHXJZHVNHGCEC-UHFFFAOYSA-N Chlorophacinone Chemical compound C1=CC(Cl)=CC=C1C(C=1C=CC=CC=1)C(=O)C1C(=O)C2=CC=CC=C2C1=O UDHXJZHVNHGCEC-UHFFFAOYSA-N 0.000 description 1
- 241000282320 Panthera leo Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000005260 corrosion Methods 0.000 description 1
- 230000007797 corrosion Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- YOBAEOGBNPPUQV-UHFFFAOYSA-N iron;trihydrate Chemical compound O.O.O.[Fe].[Fe] YOBAEOGBNPPUQV-UHFFFAOYSA-N 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000011524 similarity measure Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/14—Transformations for image registration, e.g. adjusting or mapping for alignment of images
- G06T3/147—Transformations for image registration, e.g. adjusting or mapping for alignment of images using affine transformations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/13—Edge detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
- G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
- G06T7/337—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving reference images or patches
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/90—Determination of colour characteristics
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20036—Morphological image processing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20112—Image segmentation details
- G06T2207/20164—Salient point detection; Corner detection
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种面向高层地物的异源高分遥感影像配准方法，步骤为，分别对参考影像和待配准影像进行阴影检测和高层地物筛选；分别对参考影像和待配准影像提取相位一致性特征点；对提取的相位一致性特征点集进行配准；对待配准影像运用仿射变换式，并采用双线性插值完成粗配准；在粗配准基础上，进一步采用小三角面元微分纠正，实现精配准。
Description
技术领域
本发明属于图像处理领域，特别涉及了一种异源高分遥感影像配准方法。
背景技术
遥感图像配准是指将两幅或多幅来自于不同视角、不同时间或者不同传感器的同一场景的遥感图像进行几何校准的过程，其目的是将参考图像和待配准图像进行对准，获得更为全面的图像描述供观察或进一步处理。近年来，由于传感器的日益丰富，遥感数据逐渐向着多角度、多尺度、多波段、多类型的方向发展，异源遥感影像配准已经成为了配准领域的研究热点。
与普通影像配准不同，异源高分辨率遥感影像往往存在更为严峻的局部变形问题，在配准过程中往往由于成像条件差异，导致图像之间地物的辐射特性不同，地物的描述特征不一定保持稳定。地物在不同成像模式下也可能存在较大的几何和辐射差异。其中，由于拍摄条件不同而易发生相对形变的高层地物所产生的相对视察偏移尤为严重，因其高度之间存在差异，导致其特征间不能服从同一空间变换，很可能造成大量特征点的丢失与同名点对的不合理剔除，而由于其结构、颜色特征保留较好，是配准重点区域，因此如若全部剔除，剩余过少同名点对又将大大影响配准精度。同时，由于太阳光的照射，高层地物周围往往存在大量阴影区域，而异源遥感影像间不同的拍摄角度导致阴影区域大小、形状也存在较大差异，使得对特征点所构造的特征描述符不够准确，为同名点匹配带来很大干扰。
为此，学者们已经开展了针对异源影像配准的大量研究工作。其方法大致上可以分为两类，即基于特征的配准方法与基于区域的配准方法。在基于特征的配准方法中，常用的特征有点、边缘、轮廓、区域等特征，提取的特征应对异源配准影像应对噪声、地物辐射特性和分辨率应保持稳定。Lowe等设计尺度不变特征变换(SIFT)算法，在异源影像上展现一定抗噪声干扰鲁棒性。Bay等设计加速稳健(SURF)算法，在保持尺度不变的基础上加快异源影像配准速度。Ye等依据相位特征提出了带有方向信息的相位一致直方图描述子(Histogram of Orientated Phase Congruency，HOPC)，展现了良好的多源匹配效果，但由于采用模板匹配，耗时较长。Li等基于约束点特征，设计了应用于红外与光学异源影像之间的配准方法。在基于区域的配准方法中，通常使用归一化互相关系数(Normalized CrossCorrelation Coefficient，NCCC)、互信息(Mutual Information，MI)、梯度等统计信息来完成配准。Inglada等采用归一化互相关系数和互信息两种相似性测度，对异源光学影像进行了实验，针对异源影像之间的角度、尺度偏差问题等仍需要做出一定改良。信息测度对于异源遥感影像配准已证实具备较高的鲁棒性，Pual等基于互信息并采用SPSA优化以完成光学与SAR异源影像的配准。Shams等依据梯度信息确定初始参数,接着基于互信息完成精配准，针对角度、尺度偏差具有一定的鲁棒性。Yan等引入方向梯度距离直方图配合数据驱动灰狼优化算法，以此应对多模态影像非线性强度差异的配准问题。
上述诸多方法都在不同方面做出了有效改进，尽管如此，这些方法还是大多忽略了高层地物上严重相对视察偏移所导致的特征点稀少与同名点对误差等问题，同时阴影干扰的排除也未在大多方法中得到有效实施。
发明内容
为了解决上述背景技术提到的技术问题，本发明提出了一种面向高层地物的异源高分遥感影像配准方法，显著提高异源遥感影像的配准精度。
为了实现上述技术目的，本发明的技术方案为：
一种面向高层地物的异源高分遥感影像配准方法，包括以下步骤：
(1)分别对参考影像和待配准影像进行阴影检测和高层地物筛选；
(2)分别对参考影像和待配准影像提取相位一致性特征点；
(3)对步骤(2)提取的相位一致性特征点集进行配准；
(4)对待配准影像运用仿射变换式，并采用双线性插值完成粗配准；在粗配准基础上，进一步采用小三角面元微分纠正，实现精配准。
进一步地，在步骤(1)中，进行阴影检测的方法如下：
将彩色的RGB图像转换为HSV图像，在获得色调H、饱和度S和明度V三个分量的基础上，进一步使用迭代法获得最佳分割阈值，并将每一个像素灰度与最佳分割阈值进行比较，获得阴影像素；接着对图像进行形态学闭运算，将细微阴影区域连通化，在此基础上，对阴影区域进行筛除，筛除阴影面积和阴影长宽比不在设定值范围内的阴影区域，得到正确阴影区域。
进一步地，在步骤(1)中，采用种子点区域生长分割算法，将阴影区域沿阴影方向平移获得带有相对高度信息的种子点，并进行区域生长分割，筛选出相对应的高层地物对象。
进一步地，在步骤(2)中，使用Log-Gabor滤波器在频域对图像进行不同尺度和方向的滤波处理，利用滤波后图像的幅值和相位信息，计算出图像各像素的相位一致性值：
上式中，PC(x,y)表示影像上像点(x,y)相位一致性的强度幅值，Wo(x,y)是基于频率分布的权重系数加权项，Aso(x,y)为像点在Log-Gabor滤波器特定尺度s和方向o上的振幅；ΔΦso(x,y)为相位偏移；T为噪声阈值；ε为避免除数为0的常数；
符号表示值为正时取本身，否则取0；
对于每一个特定Log-Gabor滤波器方向o，将所有尺度下的卷积结果都带入上式，计算出每个特定方向o下的相位一致性测度，并进一步计算得出相位一致最大矩M和相位一致最小矩m，两者分别用于边缘提取与角点提取：
a＝∑(PC(o)cos(o))2
b＝2∑(PC(o)cos(o))(PC(o)sin(o))
c＝∑(PC(o)sin(o))2
其中，PC(o)为特定方向o下的相位一致性测度。
进一步地，分别提取高层地物与低层地物上所有像素相位一致最小矩m，使用最大类间方差自适应方法得到高层地物初始阈值Yg0与低层地物初始阈值Yd0；
自适应计算高层地物更新阈值Yg1，其步骤如下：
(a)将高层地物上所有像素相位一致最小矩值按升序排列，剔除最前端等于或接近0的像素，对保留像素标号1至N；
(b)寻找到第一个大于等于初始阈值Yg0的像素Yk，统计其所在序列分位数
(c)将此序列分位数
乘以0.85，并设置高层地物更新阈值Yg1为序列中此分位数处的相位一致最小矩值，并基于此提取高层地物上的特征点；
对于低层地物，将影像分块，并统计每块上的低层地物特征点数量占总数量的比例，对于占比大于1/4的影像块，以最终各块内特征点数占总数量的比例和最小为目标函数f，以占比排序不变和最终占比不低于原来占比1/2为限制条件，自适应计算低层地物各影像块的更新阈值
并基于此提取低层地物上的特征点：
上式中，ki表示初始影像块内特征点数量占总数量比例，ki'表示最终各块内特征点数量占总数量比例，n表示初始占比大于1/4的影像块数量，S'表示最终低层地物特征点总数，Yi表示各影像块内所有像素相位一致最小矩值；
将低层地物上提取的特征点与高层地物上提取的特征点合并为特征点集。
进一步地，步骤(3)的具体过程如下：
(301)对于每一个提取的特征点，在周围9种图块内，基于各方向相位一致性测度构建特征描述符；统计各图块内阴影面积大小，引入阴影面积加权特征向量距离以排除阴影干扰；采用双向匹配法，完成同名点对匹配，获得匹配点对集；(302)采用RANSAC随机一致性检验，并针对高层地物上的同名点对设计变换误差自适应惩罚因子，以降低高层地物空间变化差异对映射方程的影响；
滤除错误点对并估计仿射变换式参数，确定空间仿射变换式的形式：
其中，(X,Y)和(X′,Y′)分别为待配准影像和参考影像的坐标，λ为尺度缩放因子，θ为影像的相对旋转角，(c,r)为图像在二维平面的相对平移量。
进一步地，在步骤(301)中，所述阴影面积加权特征向量距离如下：
上式中，dist'表示阴影面积加权特征向量距离，disti表示各图块初始特征向量距离，ui表示与相对阴影面积大小呈负相关的系数权重，
分别表示参考影像和待配准影像对应图块的特征向量，rm表示与初始特征向量距离对应的相对阴影面积大小，ri表示第i对图块的相对阴影面积大小。
进一步地，在步骤(302)中，通过变换误差自适应惩罚因子修正高层地物上同名点对变换误差：
E'＝E×F
上式中，E'表示修正后的高层地物上同名点对变换误差，E表示高层地物上同名点对变换误差，F表示变换误差自适应惩罚因子，P为惩罚因子系数，q表示惩罚因子关于相对高度的灵敏度，h1、h2分别表示两图像上同一高层地物相对高度，(xa,ya)、(xb,yb)分别表示两图像同名点对坐标，H表示仿射变换式。
进一步地，在步骤(4)中，采用最近距离算法在粗配准后影像上构建局部小三角面元，每个小三角形都是唯一最简单形状，在每个小三角面元内，逐个进行一次多项式的纠正：
根据三角形的各顶点坐标确定系数a0、a1、a2、b0、b1、b2。
采用上述技术方案带来的有益效果：
1、本发明通过阴影检测并结合区域分割技术，筛选出高层地物，在此基础上对不同类型地物自适应筛选特征点，提高高层地物关键特征点数及特征点整体分布均匀性。
2、本发明在特征点匹配中引入阴影面积加权特征向量，以降低阴影区域对特征点周围纹理结构相似度的影响，提高同名点对匹配准确度。
3、本发明在仿射变换式解算阶段针对高层地物点对，基于地物相对高度设计自适应惩罚因子，在充分利用所有点对间变换关系的条件下，降低高层地物空间变化差异对映射方程的影响权重，提高最终配准精度。
附图说明
图1是本发明的方法流程图；
图2是本发明中获取种子点示意图；
图3是本发明中9种图块类型示意图；
图4是南京区域异源光学影像实验结果对比图；
图5是西安区域异源光学影像实验结果对比图。
具体实施方式
以下将结合附图，对本发明的技术方案进行详细说明。
本发明设计了一种面向高层地物的异源高分遥感影像配准方法，如图1所示，包括4个步骤：(1)阴影检测与高层地物筛选；(2)相位一致性特征点提取；(3)相位一致性特征点集配准；(4)图像变换。
(1)阴影检测与高层地物筛选
利用开源opencv工具包，基于HSV彩色空间变换，在参考影像与待配准影像上进行阴影检测。同时依据所检测阴影区域的面积、长宽比进行筛除，并计算阴影对象沿阴影方向平均跨度，作归一化处理，记录为对应高层地物相对高度。经检验阴影检测效果良好。在此基础上采用种子点区域生长分割算法，将阴影区域沿阴影方向平移获得带有相对高度信息的种子点，并进行区域生长分割，筛选出相对应的高层地物对象。
高层地物周围往往由于太阳光照射，周围存在大量阴影，其可用于提供高层地物位置、大小、高度等信息，因此首先对影像采用HSV空间变换进行阴影检测，进而利用形态学闭运算处理获得阴影检测结果图像。
HSV模型是基于色调(H)、饱和度(S)和明度(V)的六棱锥颜色空间。按照对应关系式，将彩色的RGB图像转换为HSV图像，并根据HSV模型阴影区域的特点进行阴影检测。
在获得三个分量的基础之上，进一步使用迭代法获得最佳分割阈值，并将每一个像素灰度与最佳分割阈值进行比较，获得阴影像素。接着对图像进行形态学闭运算，将细微阴影区域连通化。在此基础之上，对阴影区域进行筛除，筛除阴影面积小于阈值Smin、阴影面积大于阈值Smax、阴影长宽比小于阈值Cmin、阴影长宽比大于阈值Cmax的阴影区域。对于上述四个阈值，分别统计所有阴影区域面积直方图和长宽比直方图，并获得每个直方图第一个长间隔左侧有效数值与第一个长间隔右侧有效数值，基于此将Smin设置为50，Smax设置为1000，Cmin设置为0.08，Cmax设置为12。最终提取得到正确阴影区域。
对所得的阴影检测结果图利用相位一致性最大矩获取阴影边缘，再以阴影连通域为单位,利用Ransac方法对每个连通域进行直线检测，将长度大于10的直线予以保留，并计算统计所有直线角度，获取阴影边缘直线角度直方图，选取其中最大峰值处为阴影方向β，次大方向α。对每幅影像上的阴影区域，在α方向上进行腐蚀运算，以略微扩大搜索范围，其中结构性元se设置为ster1('line’,a,α),a值取为6，然后将阴影区域A1沿阴影方向β平移b像素，b取值略微大于a*sinα，A1平移后得到新阴影A2，选取A2中不属于A1的像素点为种子点A：
A＝A2-A2∩A1
获取种子点的过程如图2所示。
统计初始阴影区域A1上沿阴影方向归一化跨度作为对应高层地物相对高度h，赋予种子点，使其带有相对高度信息。对得到的各高层地物种子点，利用区域增长分割算法筛选出高层地物。在区域增长中，从每一个种子点出发，逐步迭代选择周围与其具有相似属性的邻点并合并。对于种子点A，搜索到其邻点d，若d满足下式则将d加入区域U，以d作为新种子点进行迭代搜索。
U(t)＝{{d∈U}∪{lowthred＜d＜highthred}}
其中lowthred和highthred分别为门限控制的上下阈值。种子点区域生长停止，即获得了高层地物筛选结果。
(2)相位一致性特征点提取
考虑到异源高分辨率遥感影像在灰度强度上往往存在非线性变化，亮度上也存在差异性，因此引入相位一致性对影像进行转换，并计算相位一致最大矩与相位一致最小矩来完成特征点提取。
相位一致性是一种使用傅里叶谐波分量描述信号局部强度的特征，它利用特定方向o不同尺度s的Log Gabor滤波器，将信号分解为频率域下与o和s相对应的傅里叶谐波分量，并加权叠合。相位一致性特征具有强度不变特性，可抵抗异源影像中的亮度差异与灰度非线性差异。
生成相位一致性图时，使用Log-Gabor滤波器在频域对图像进行不同尺度和方向的滤波处理，利用滤波后图像的幅值和相位信息，计算出图像各像素的相位一致性值：
上式中，PC(x,y)表示影像上像点(x,y)相位一致性的强度幅值，Wo(x,y)是基于频率分布的权重系数加权项，Aso(x,y)为像点在Log-Gabor滤波器特定尺度s和方向o上的振幅；ΔΦso(x,y)为相位偏移；T为噪声阈值；ε为避免除数为0的常数；
符号表示值为正时取本身，否则取0。
在此基础之上，对于每一个特定Log-Gabor滤波器方向，将所有尺度下的卷积结果都带上式，即可计算出每个特定方向o下的相位一致性测度，则可进一步计算得出相位一致最大矩M与相位一致最小矩m，两者分别用于边缘提取与角点提取：
a＝∑(PC(o)cos(o))2
b＝2∑(PC(o)cos(o))(PC(o)sin(o))
c＝∑(PC(o)sin(o))2
基于高层地物筛选结果，在高层地物与低层地物中分别采用阈值自适应特征点提取策略，提高高层地物上关键特征点数量，同时减少密度过高的低层地物中特征点数，提高特征点总体分布均匀性。最终在参考影像与待配准影像中分别提取出特征点集。
在高层地物筛选结果基础之上，对整幅图像计算相位一致最大矩M、最小矩m，由于角点相较于边缘点，特征信息更加明显集中，因此主要依据最小矩m提取角点，而最大矩M则作为筛选角点的限制条件，即剔除不满足一定最大矩阈值的角点。
分别提取高层地物与低层地物上所有像素相位一致最小矩值，使用最大类间方差自适应方法得到高层地物初始阈值Yg0与低层地物初始阈值Yd0。为提高层地物上关键特征点数量以更细致得描述其复杂纹理特征，进一步自适应计算高层地物更新自适应阈值Yg1，其步骤如下：步骤1：将高层地物上所有像素相位一致最小矩值按升序排列，剔除最前端等于或接近0的像素，对保留像素标号1至N；步骤2：寻找到第一个大于等于初始阈值Yg0的像素Yk，统计其所在序列分位数
步骤3：将此序列分位数乘以0.85，并设置更新阈值Yg1为序列中此分位数处相位一致最小矩值，并基于此提取高层地物上满足条件的特征点。
对于低层地物，则需减小密度过高处的特征点数量以避免特征点冗余现象。为此将影像分为9块，并统计每块上的低层地物特征点数量占总数量S比例ki，对于占比大于1/4的n块影像块，以各块内特征点数占总数量S'比例ki'和最小为目标函数，以占比排序不变和最终占比不低于原来占比1/2为限制条件，自适应降低筛选阈值，定下上述限制条件的目的是为保证各块纹理复杂度与其中特征点数呈正相关的前提。最终基于各图块更新自适应阈值
在低层地物上提取特征点，并与高层地物上特征点合并为特征点集，具体过程如下：
上式中，ki表示初始影像块内特征点数量占总数量比例，ki'表示最终各块内特征点数量占总数量比例，n表示初始占比大于1/4的影像块数量，S'表示最终低层地物特征点总数，Yi表示各影像块内所有像素相位一致最小矩值。
综上所述，设计阈值自适应特征点提取策略能够有效提高高层地物上特征点数量以及特征点总体分布均匀性，且依据影像本身属性自适应选取，很好地排除了异源影像差异所带来的干扰。
(3)相位一致性特征点集配准
步骤1：点集间对应关系评估。对于每一个提取所得的特征点，在周围9种图块内，基于各方向相位一致性测度PC(o)构建特征描述符。在此基础之上，统计各图块内阴影面积大小，基于此引入阴影面积加权特征向量距离，以排除阴影干扰。采用双向匹配法，完成同名点对匹配，获得匹配点对集。基于相位一致测度的特征描述符，具有光照不变性与灰度不变性，可排除异源遥感影像拍摄条件不同所带来的干扰。
由于对于异源遥感影像，其呈像时间、角度不同，其建筑物周围阴影的大小、形状、边缘轮廓也不同。在获得同名点对时，是根据特征点周围纹理、结构信息相似性来匹配的，阴影会为此带来较大干扰，但若简单采用阴影剔除算法进行消除，则又会丢失大量纹理特征，且阴影边缘干扰依然存在。因此在对高层地物的特征点进行描述和匹配时有必要设计策略对此类干扰进行剔除。
本发明引入了阴影面积加权特征向量距离dist'。如图3所示，筛选特征点周围9种类型图块，每个图块大小都为32×32像素，对其分别构建相位一致特征描述符，同时计算出各图块内阴影面积大小。
为避免高层地物筛选可能存在的些许误差，即本该正确匹配的点对，但两点在异源影像上分别处于不同类地物上，因此依然计算两影像间任意两点之间匹配程度，但同时为降低运算复杂度，分以下两类处理：(1)若所匹配的点对都为低层地物点，则只选择以特征点为中心的图块1计算特征向量距离；(2)若所匹配的点对上有任意一点在所在影像上被确认为高层地物上特征点，则基于各图块内阴影面积大小，设计加权特征向量距离。其步骤如下：
(a)对于两点间9对对应图块，各计算之间初始特征向量距离disti。
(b)将两点间对应图块内阴影面积相加，并除以2倍图块总面积，得到与初始特征向量距离对应的9个相对阴影面积大小rm，其值都在0至1之间。
(c)计算出与相对阴影面积大小呈负相关的系数权重ui，乘以对应图块初始特征距离并相加，得到最终阴影面积加权特征向量距离dist'：
上式中，dist'表示阴影面积加权特征向量距离，disti表示各图块初始特征向量距离，ui表示与相对阴影面积大小呈负相关的系数权重，
分别表示参考影像和待配准影像对应图块的特征向量，rm表示与初始特征向量距离对应的相对阴影面积大小，ri表示第i对图块的相对阴影面积大小。
此特征距离中，阴影面积和最小的图块间距离所占权重最大，此距离所受阴影干扰最小，最为合理真实，而阴影面积较大的则相反。利用此加权特征距离，能够在综合考虑特征点周围各类图块匹配程度的条件下，极大程度得降低阴影干扰，从而有效提高正确同名点对数量，进而提高匹配精度。
基于上述特征向量距离，对影像上特征点之间进行双向匹配，并判定当特征点最小距离与次小距离之比小于阈值t，且双向满足此条件时，则认为两个特征点是一对匹配点对，最终获得匹配点集。
步骤2：空间变换确定。采用RANSAC随机一致性检验，并针对高层地物上的同名点对设计了变换误差自适应惩罚因子，从而降低高层地物空间变化差异对映射方程的影响。最终滤除错误点对并估计仿射变换式参数，确定空间变换式，空间仿射变换式的形式如下：
其中，(X,Y)和(X′,Y′)分别为待配准影像和参考影像的坐标，λ为尺度缩放因子，θ为影像的相对旋转角，(c,r)为图像在二维平面的相对平移量。
异源高分辨率遥感影像中，高层地物上同名点对无法与其余地物上点对服从一致空间变换，若是采取大多现有方法中统一估计仿射变换式参数的策略，将会带来较大误差。若分别在高层地物与低层地物上解算变换式，则又会由于同名点对数量的降低及高层地物上相对形变的存在，变换式都不准确。
本发明在解算仿射变换式时引入了自适应惩罚因子F以获取正确空间变换关系。对于高层地物上的同名点对，所在地物相对高度越高，则对其变换误差的合理容忍度就越高，因为其空间变换本身就与正确变换关系之间存在与相对高度呈正相关的误差。若不添加惩罚因子，则在解算仿射参数时以所有内点变换均方根误差最小为目标，会构造出更迁就于高层地物上空间变换关系的变换式，必定降低配准精度。因此对于高层地物上同名点对变换误差E，添加自适应惩罚因子进行修正，其值在0-1之间，与高层地物所对应的相对高度h呈负相关：
E'＝E×F
上式中，E'表示修正后的高层地物上同名点对变换误差，E表示高层地物上同名点对变换误差，F表示变换误差自适应惩罚因子，P为惩罚因子系数，q表示惩罚因子关于相对高度的灵敏度，h1、h2分别表示两图像上同一高层地物相对高度，(xa,ya)、(xb,yb)分别表示两图像同名点对坐标，H表示仿射变换式。
在迭代确定初始变换模型中选取8对同名点对，并对变换误差进行惩罚因子修正，基于最小二乘原理，以误差均方根最小为目标函数估算初始变换参数。再进一步基于修正误差筛选内点，并代入所有内点，再次进行最小二乘，确认每次迭代过程中的最终变换参数。通过引入变换误差自适应惩罚因子，既充分考虑了所有地物上的同名点对变换关系，又大大降低了高层地物空间变化差异对映射方程的影响权重。
(4)图像变换
对原始待配准影像运用仿射变换式，并采用双线性插值完成粗配准。在粗配准基础之上，进一步采用小三角面元微分纠正，实现精配准。采用最近距离算法在粗配准后影像上构建局部小三角面元，每个小三角形都是唯一最简单形状。在每个小三角面元内，逐个进行一次多项式的纠正：
逐步依据三角形的各顶点坐标确定系数a0、a1、a2、b0、b1、b2，再对待配准影像进行高精度矫正，完成精配准。
分别采用本发明方法、SIFT算法、SURF算法与HOPC算法对两组数据进行完整配准实验，并从量化指标与视觉分析两方面对比最终实验结果。两组数据如下：
1、使用ZY_3(资源3号)，在2018年所拍摄的南京区域异源光学影像数据，其分辨率为2.1m，处于RGB波段。使用GF_2(高风2号)，在2020年所拍摄的南京对应区域光学影像数据，其分辨率为0.81m，处于RGB波段。将资源3号影像作为参考图像，高风2号影像作为待配准图像。
2、用GF_2(高风2号)，在2018年所拍摄的西安区域异源光学影像数据，其分辨率为0.81m，处于RGB波段。使用GF_1(高分1号)，在2020年所拍摄的西安对应区域光学影像数据，其分辨率为2m，处于RGB波段。将高分1号影像作为参考图像，高风2号影像作为待配准图像。
南京区域异源光学影像对比实验结果如表1、图4所示。西安区域异源光学影像对比实验结果如表2、图5所示。分析结果可得，本发明方法在所有图像配准中均展现出最佳性能。相较于南京区域影像，西安区域影像上存在大量汽车、集装箱等斑点噪声信息，因此各项配准指标都有所下降，但同时影像上也存在较多高层地物及阴影区域，相较其他算法，本发明方法对其特殊处理以提取得到最多特征点对数，展现出最佳配准性能。SIFT算法虽然具备尺度不变特性，但其在对异源高分辨率遥感影像的配准中提取了大量冗余特征点，造成匹配率与配准精度下降。SURF算法基于灰度极值信息筛选特征点，在异源影像中提取数量过少，特别是对西安影像数据，直接导致配准失败。HOPC算法模板匹配与本发明一样基于相位一致性提取特征点，能较正确得对异源影像间相似几何结构信息进行配准，但由于其未针对相对形变高层地物采取特殊处理，因此配准精度受到干扰以至低于本发明方法，同时由于其模板遍历校准的过程，其耗时远远大于本发明方法。
表1
表2
实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。
Claims (7)
1.一种面向高层地物的异源高分遥感影像配准方法，其特征在于，包括以下步骤：
(1)分别对参考影像和待配准影像进行阴影检测和高层地物筛选；
(2)分别对参考影像和待配准影像提取相位一致性特征点；
(3)对步骤(2)提取的相位一致性特征点集进行配准；
(4)对待配准影像运用仿射变换式，并采用双线性插值完成粗配准；在粗配准基础上，采用小三角面元微分纠正，实现精配准；
在步骤(2)中，使用Log-Gabor滤波器在频域对图像进行不同尺度和方向的滤波处理，利用滤波后图像的幅值和相位信息，计算出图像各像素的相位一致性值：
上式中，PC(x,y)表示影像上像点(x,y)相位一致性的强度幅值，Wo(x,y)是基于频率分布的权重系数加权项，Aso(x,y)为像点在Log-Gabor滤波器特定尺度s和方向o上的振幅；ΔΦso(x,y)为相位偏移；T为噪声阈值；ε为避免除数为0的常数；
符号表示值为正时取本身，否则取0；
对于每一个特定Log-Gabor滤波器方向o，将所有尺度下的卷积结果都带入上式，计算出每个特定方向o下的相位一致性测度，并计算得出相位一致最大矩M和相位一致最小矩m，两者分别用于边缘提取与角点提取：
a＝∑(PC(o)cos(o))2
b＝2∑(PC(o)cos(o))(PC(o)sin(o))
c＝∑(PC(o)sin(o))2
其中，PC(o)为特定方向o下的相位一致性测度；
步骤(3)的具体过程如下：
(301)对于每一个提取的特征点，在周围9种图块内，基于各方向相位一致性测度构建特征描述符；统计各图块内阴影面积大小，引入阴影面积加权特征向量距离以排除阴影干扰；采用双向匹配法，完成同名点对匹配，获得匹配点对集；
(302)采用RANSAC随机一致性检验，并针对高层地物上的同名点对设计变换误差自适应惩罚因子，通过变换误差自适应惩罚因子修正高层地物上同名点对变换误差，以降低高层地物空间变化差异对映射方程的影响；
滤除错误点对并估计仿射变换式参数，确定空间仿射变换式的形式：
其中，(X,Y)和(X′,Y′)分别为待配准影像和参考影像的坐标，λ为尺度缩放因子，θ为影像的相对旋转角，(c,r)为图像在二维平面的相对平移量。
2.根据权利要求1所述面向高层地物的异源高分遥感影像配准方法，其特征在于，在步骤(1)中，进行阴影检测的方法如下：
将彩色的RGB图像转换为HSV图像，在获得色调H、饱和度S和明度V三个分量的基础上，使用迭代法获得最佳分割阈值，并将每一个像素灰度与最佳分割阈值进行比较，获得阴影像素；接着对图像进行形态学闭运算，将细微阴影区域连通化，在此基础上，对阴影区域进行筛除，筛除阴影面积和阴影长宽比不在设定值范围内的阴影区域，得到正确阴影区域。
3.根据权利要求1所述面向高层地物的异源高分遥感影像配准方法，其特征在于，在步骤(1)中，采用种子点区域生长分割算法，将阴影区域沿阴影方向平移获得带有相对高度信息的种子点，并进行区域生长分割，筛选出相对应的高层地物对象。
4.根据权利要求1所述面向高层地物的异源高分遥感影像配准方法，其特征在于，分别提取高层地物与低层地物上所有像素相位一致最小矩m，使用最大类间方差自适应方法得到高层地物初始阈值Yg0与低层地物初始阈值Yd0；
自适应计算高层地物更新阈值Yg1，其步骤如下：
(a)将高层地物上所有像素相位一致最小矩值按升序排列，剔除最前端等于或接近0的像素，对保留像素标号1至N；
(b)寻找到第一个大于等于初始阈值Yg0的像素Yk，统计其所在序列分位数
(c)将此序列分位数
乘以0.85，并设置高层地物更新阈值Yg1为序列中此分位数处的相位一致最小矩值，并基于此提取高层地物上的特征点；
对于低层地物，将影像分块，并统计每块上的低层地物特征点数量占总数量的比例，对于占比大于1/4的影像块，以最终各块内特征点数占总数量的比例和最小为目标函数f，以占比排序不变和最终占比不低于原来占比1/2为限制条件，自适应计算低层地物各影像块的更新阈值
并基于此提取低层地物上的特征点：
上式中，ki表示初始影像块内特征点数量占总数量比例，ki'表示最终各块内特征点数量占总数量比例，n表示初始占比大于1/4的影像块数量，S'表示最终低层地物特征点总数，Yi表示各影像块内所有像素相位一致最小矩值；
将低层地物上提取的特征点与高层地物上提取的特征点合并为特征点集。
5.根据权利要求1所述面向高层地物的异源高分遥感影像配准方法，其特征在于，在步骤(301)中，所述阴影面积加权特征向量距离如下：
上式中，dist'表示阴影面积加权特征向量距离，disti表示各图块初始特征向量距离，ui表示与相对阴影面积大小呈负相关的系数权重，
分别表示参考影像和待配准影像对应图块的特征向量，rm表示与初始特征向量距离对应的相对阴影面积大小，ri表示第i对图块的相对阴影面积大小。
6.根据权利要求1所述面向高层地物的异源高分遥感影像配准方法，其特征在于，在步骤(302)中，通过变换误差自适应惩罚因子修正高层地物上同名点对变换误差：
E'＝E×F
上式中，E'表示修正后的高层地物上同名点对变换误差，E表示高层地物上同名点对变换误差，F表示变换误差自适应惩罚因子，P为惩罚因子系数，q表示惩罚因子关于相对高度的灵敏度，h1、h2分别表示两图像上同一高层地物相对高度，(xa,ya)、(xb,yb)分别表示两图像同名点对坐标，H表示仿射变换式。
7.根据权利要求1所述面向高层地物的异源高分遥感影像配准方法，其特征在于，在步骤(4)中，采用最近距离算法在粗配准后影像上构建局部小三角面元，每个小三角形都是唯一最简单形状，在每个小三角面元内，逐个进行一次多项式的纠正：
根据三角形的各顶点坐标确定系数a0、a1、a2、b0、b1、b2。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110570103.1A CN113256653B (zh) | 2021-05-25 | 2021-05-25 | 一种面向高层地物的异源高分遥感影像配准方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110570103.1A CN113256653B (zh) | 2021-05-25 | 2021-05-25 | 一种面向高层地物的异源高分遥感影像配准方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113256653A CN113256653A (zh) | 2021-08-13 |
| CN113256653B true CN113256653B (zh) | 2023-05-09 |
Family
ID=77184161
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110570103.1A Active CN113256653B (zh) | 2021-05-25 | 2021-05-25 | 一种面向高层地物的异源高分遥感影像配准方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113256653B (zh) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114565653B (zh) * | 2022-03-02 | 2023-07-21 | 哈尔滨工业大学 | 一种存在旋转变化和尺度差异的异源遥感图像匹配方法 |
| CN115049708B (zh) * | 2022-04-12 | 2023-04-07 | 南京雷电信息技术有限公司 | 基于lsd直线检测与模板匹配的sar图像配准方法 |
| CN115876785B (zh) * | 2023-02-02 | 2023-05-26 | 苏州誉阵自动化科技有限公司 | 一种用于产品缺陷检测的视觉识别系统 |
| CN116385502B (zh) * | 2023-03-09 | 2024-04-19 | 武汉大学 | 一种基于几何约束下区域搜索的图像配准方法 |
| CN117853487B (zh) * | 2024-03-07 | 2024-05-14 | 浙江合丰科技有限公司 | 一种基于图像处理技术的fpc连接器裂纹检测方法及系统 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103077527A (zh) * | 2013-02-05 | 2013-05-01 | 湖北工业大学 | 一种稳健的多源卫星遥感影像配准方法 |
| CN103514606A (zh) * | 2013-10-14 | 2014-01-15 | 武汉大学 | 一种异源遥感影像配准方法 |
| CN104167003A (zh) * | 2014-08-29 | 2014-11-26 | 福州大学 | 一种遥感影像的快速配准方法 |
| CN105261014A (zh) * | 2015-09-30 | 2016-01-20 | 西南交通大学 | 一种多传感器遥感影像匹配方法 |
| CN111476251A (zh) * | 2020-03-26 | 2020-07-31 | 中国人民解放军战略支援部队信息工程大学 | 一种遥感影像匹配方法及装置 |
| CN112634335A (zh) * | 2020-12-25 | 2021-04-09 | 清华大学 | 面向非线性辐射畸变的稳健遥感影像特征点对提取方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107563438B (zh) * | 2017-08-31 | 2019-08-30 | 西南交通大学 | 一种快速鲁棒的多模态遥感影像匹配方法和系统 |
-
2021
- 2021-05-25 CN CN202110570103.1A patent/CN113256653B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103077527A (zh) * | 2013-02-05 | 2013-05-01 | 湖北工业大学 | 一种稳健的多源卫星遥感影像配准方法 |
| CN103514606A (zh) * | 2013-10-14 | 2014-01-15 | 武汉大学 | 一种异源遥感影像配准方法 |
| CN104167003A (zh) * | 2014-08-29 | 2014-11-26 | 福州大学 | 一种遥感影像的快速配准方法 |
| CN105261014A (zh) * | 2015-09-30 | 2016-01-20 | 西南交通大学 | 一种多传感器遥感影像匹配方法 |
| CN111476251A (zh) * | 2020-03-26 | 2020-07-31 | 中国人民解放军战略支援部队信息工程大学 | 一种遥感影像匹配方法及装置 |
| CN112634335A (zh) * | 2020-12-25 | 2021-04-09 | 清华大学 | 面向非线性辐射畸变的稳健遥感影像特征点对提取方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113256653A (zh) | 2021-08-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113256653B (zh) | 一种面向高层地物的异源高分遥感影像配准方法 | |
| CN109409292B (zh) | 基于精细化特征优化提取的异源图像匹配方法 | |
| CN110097093B (zh) | 一种异源图像精确匹配方法 | |
| CN111260731B (zh) | 一种棋盘格亚像素级角点自适应检测的方法 | |
| CN104574347B (zh) | 基于多源遥感数据的在轨卫星图像几何定位精度评价方法 | |
| CN111079556A (zh) | 一种多时相无人机视频图像变化区域检测及分类方法 | |
| CN107564017B (zh) | 一种城市高分遥感影像阴影检测及分割方法 | |
| CN112017223B (zh) | 基于改进SIFT-Delaunay的异源图像配准方法 | |
| CN103679714A (zh) | 一种基于梯度互相关的光学和sar图像自动配准方法 | |
| CN102800099B (zh) | 多特征多级别的可见光与高光谱图像高精度配准方法 | |
| CN113012234B (zh) | 基于平面变换的高精度相机标定方法 | |
| Chen et al. | Robust affine-invariant line matching for high resolution remote sensing images | |
| CN105405138B (zh) | 基于显著性检测的水面目标跟踪方法 | |
| CN110222661B (zh) | 一种用于运动目标识别及跟踪的特征提取方法 | |
| CN116612165A (zh) | 一种针对大视角差异sar图像的配准方法 | |
| CN115601407A (zh) | 一种红外与可见光图像配准方法 | |
| CN111311596A (zh) | 一种基于改进lbp特征的遥感影像变化检测方法 | |
| CN105631860B (zh) | 基于局部排序方向直方图描述子的图像同名点提取方法 | |
| CN106250687B (zh) | 去扁化ipp的沉积砾石圆度计算方法 | |
| CN116758266A (zh) | 一种指针式仪表的读数方法 | |
| CN115035326B (zh) | 一种雷达图像与光学图像精确匹配方法 | |
| CN114565653B (zh) | 一种存在旋转变化和尺度差异的异源遥感图像匹配方法 | |
| CN115423851A (zh) | 一种基于os-sift的可见光-sar图像配准算法 | |
| CN115511928A (zh) | 多光谱图像的匹配方法 | |
| CN115588033A (zh) | 基于结构提取的合成孔径雷达与光学图像配准系统及方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |