CN115115678A - InSAR干涉测量复数图像高精度配准方法 - Google Patents
InSAR干涉测量复数图像高精度配准方法 Download PDFInfo
- Publication number
- CN115115678A CN115115678A CN202210297307.7A CN202210297307A CN115115678A CN 115115678 A CN115115678 A CN 115115678A CN 202210297307 A CN202210297307 A CN 202210297307A CN 115115678 A CN115115678 A CN 115115678A
- Authority
- CN
- China
- Prior art keywords
- image
- points
- registration
- insar
- point
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 176
- 230000009466 transformation Effects 0.000 claims abstract description 43
- 238000011156 evaluation Methods 0.000 claims abstract description 37
- 238000004364 calculation method Methods 0.000 claims abstract description 35
- 230000001427 coherent effect Effects 0.000 claims abstract description 16
- 238000012952 Resampling Methods 0.000 claims abstract description 14
- 238000011426 transformation method Methods 0.000 claims abstract description 11
- 230000008569 process Effects 0.000 claims description 26
- 238000005070 sampling Methods 0.000 claims description 22
- 230000000694 effects Effects 0.000 claims description 18
- 239000013598 vector Substances 0.000 claims description 18
- 238000012545 processing Methods 0.000 claims description 15
- 238000007689 inspection Methods 0.000 claims description 14
- 238000005305 interferometry Methods 0.000 claims description 14
- 230000008859 change Effects 0.000 claims description 10
- 238000001914 filtration Methods 0.000 claims description 9
- 230000006835 compression Effects 0.000 claims description 6
- 238000007906 compression Methods 0.000 claims description 6
- 230000007423 decrease Effects 0.000 claims description 3
- 238000009795 derivation Methods 0.000 claims description 3
- 238000009499 grossing Methods 0.000 claims description 3
- 238000005286 illumination Methods 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 230000004044 response Effects 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 238000001514 detection method Methods 0.000 abstract description 4
- 238000011439 discrete element method Methods 0.000 abstract description 4
- 230000000750 progressive effect Effects 0.000 abstract description 2
- 238000000638 solvent extraction Methods 0.000 abstract description 2
- 238000005516 engineering process Methods 0.000 description 18
- 238000010586 diagram Methods 0.000 description 13
- 238000011161 development Methods 0.000 description 12
- 230000018109 developmental process Effects 0.000 description 12
- 238000012544 monitoring process Methods 0.000 description 11
- 238000002474 experimental method Methods 0.000 description 5
- 238000013507 mapping Methods 0.000 description 3
- 230000003287 optical effect Effects 0.000 description 3
- 238000010587 phase diagram Methods 0.000 description 3
- 206010040844 Skin exfoliation Diseases 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000035618 desquamation Effects 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 230000004927 fusion Effects 0.000 description 2
- 230000010365 information processing Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000008034 disappearance Effects 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000011524 similarity measure Methods 0.000 description 1
- 238000003809 water extraction Methods 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/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
-
- 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
- G06T2207/10044—Radar image
-
- 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/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
-
- 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/20212—Image combination
- G06T2207/20224—Image subtraction
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Image Analysis (AREA)
Abstract
本申请结合InSAR图像自身特点分步骤递进提出两种提高配准精度的方法，第一阶段采用基于改进特征协同变换方法进行InSAR图像配准，采用多项式拟合求出主从图像分别在距离和方位向的偏移多项式,然后进行重采样和干涉计算得到干涉图，对主从图像进行关键点的检验和配准，分块将图像分成若干大小适当区域进行特征点匹配；第二阶段在改进特征协同变换方法的基础上结合基于附加SPS量能评价的非线性InSAR配准法来提高图像局部配准的精度，消除配准之后相干计算得出的干涉图中的SPS点,基于大量奇异旋转点对图像进行非线性配准；两种方法具有逐步递进效果，生成质量较高的干涉图和DEM，图像配准精度高，错误率低，具有很大的实际运用价值。
Description
技术领域
本申请涉及一种InSAR复数图像配准方法，特别涉及一种InSAR干涉测量复数图像高精度配准方法，属于InSAR图像配准技术领域。
背景技术
合成孔径雷达干涉测量InSAR是一种快速发展的大地测量技术，能全天候获得高精度、连续覆盖的地面高程和地表形变信息。近年来，InSAR己在地形测绘(地面高程模型构建)、环境监测、灾害监测评估、能源资源勘查(油气田开采、矿藏资源开采、地下水抽采等)等相关领域得到了广泛应用，并取得了一系列重要成果。高分辨率、多化、多站、多组网SAR系统的长足进展，新型机载SAR和地基SAR系统逐渐成熟,显著增强了地形测绘和地表形变监测能力。其中多站InSAR、多基线InSAR技术推动了高精度地形测绘的发展，而时序InSAR技术的发展显著提高了地表形变测量的精度和可靠性,，在地球科学研发应用中发挥的作用越来越大。
但由于受到数据源质量、数量和相关参数信息的限制，目前大多数InSAR成像干涉数据尚未大规模应用，急需更大的发展，缺少通过对有限的SAR复数图像数据的研究，实现InSAR数据处理中的图像配准技术和后续的解缠工作，最后大致恢复成出DEM的方法。
图像配准是数字图像处理技术的核心、难点和热点，伴随着计算机视觉与信息处理技术的发展，图像匹配方法的复杂程度也越来越高，从开始简单的模板匹配，发展至角度旋转的图像匹配，再发展至基于特征点提取的匹配方法，这一系列发展使图像匹配性能得到了很大提升。
配准方法的分类可以依据不同的准则，依据变换模型的复杂程度对配准方法进行分类,并归纳了配准技术的实现步骤：特征三维、相似度测量、搜索三维和策略。
对于InSAR图像的配准方法，目前用的比较多的是先确定控制点，再用控制点拟合多项式，同时对待配准的数据进行坐标系转换以及重采样,目前比较常用的配准方法是模板或模式匹配中的相关系数法,虽然有配准的精度不能满足InSAR对数据的要求,但其结果可以提供有效的初始值给图像后续进一步的四配。图像配准是图像处理领域的核心，如果配准只使用点特征，就很难满足目前各种方法的要求，现有技术的方法已被广泛应用于遥感图像处理，计算机视觉，模式识别等。基于特征点方法是一个热点和难点问题,而这类主要方法包括Harris配准方法以及最大稳定极值区域(MSER)配准方法等，这些方法被成功的应用于光学图像配准。由于SAR图像和光学图像之前不同的机制,但如果将基于特征的配准方法直接使用在SAR图像,在图像配准中的一些缺陷将不能避免，如低精度，高错误率。
InSAR图像配准中存在较大的形变,同一特征灰度图像之间存在较大差异，基于灰度的图像配准方法在SAR图像配准中是无效的。到目前为止，有许多图像配准方法，例如，基于图像边缘特征点的同名点匹配法，但这种方法非常复杂。纹理驱动的SAR图像自动配准方法可用于具有更多的亮点或者特征线的SAR图像配准，但这些特征在正常情况下都不能被发现。基于不变矩的图像自动配准方法最初需要在封闭区域中提取图像特征，然后配准相同的名称点，但由于SAR图像的信噪比低，封闭区域的特征点不能从SAR图像处理过程中获得。结合小波迭代细化和信息融合的同源的SAR图像配准方法，在SAR图像失真小的时候是有效的，但很难配准具有不同的分辨率的图像。现有技术的普通图像配准方法在主图像与副图像之间存在较大变形、旋转和缩放等问题时，不再适用于实时SAR图像配准。
为了合成SAR干涉图,需要通过观察两个相同的位置，即“主”和“从”地图得到的两个复数图像，InSAR能够在相差微小的时间差甚至同时获得同一被监测目标地点的两幅甚至更多的图像，在得到干涉图这一部消除了干扰噪声，使得干涉图上的信息只与所监测的目标地点的位置有关，从而充分利用获得的SAR图像中的数据信息，与其它测量地形形变的方法相比，其区别在于能准确、经济、高效、快速的解决大范围区域地形与形变中棘手的问题，与其它光学图像对比，InSAR突出的优点在于：一是InSAR能够适应各种天气气候及工作时段，监测范围广泛，能够较快的获得信息数据；二是InSAR不仅可以运用多种方法来获得监测难度大的目标地点的地形信息,还可以用来做高质量的地表形变监测；三是InSAR可提供的SAR图像数据时间跨度大，且其监测地面高程变化的精度可达到毫米级甚至更高，可用比其它方法更少的花费来获得连续的、高精度的观测信息，可以用来对地面沉降的持续观测。
然而，有许多影响InSAR数据质量的关键因素，其中一个重要因素是InSAR图像配准，决定InSAR技术的鲁棒性和精度，从而决定它的应用和未来发展，提高InSAR图像配准的准确性将促进InSAR技术的发展，扩大其实际应用范围。
综合来看，现有技术的InSAR复数图像配准存在明显不足，其主要缺陷和设计难点包括：
第一，对于InSAR图像的配准方法，现有技术一般是先确定控制点，再用控制点拟合多项式，同时对待配准的数据进行坐标系转换以及重采样,配准采用相关系数法,精度不能满足InSAR对数据的要求,图像配准只使用点特征，很难满足要求，现有技术基于特征点的方法一般只能应用于光学图像配准，由于SAR图像和光学图像之前不同的机制,如果将基于特征的配准方法直接使用在SAR图像,在图像配准中的一些缺陷将不能避免，如低精度，高错误率；而且由于受到数据源质量、数量和相关参数信息的限制，现有技术缺少通过对有限的SAR复数图像数据的配准，无法实现InSAR图像配准技术和后续的解缠工作，无法生成高质量DEM；
第二，InSAR图像配准中存在较大的形变,同一特征灰度图像之间存在较大差异，现有技术基于灰度的图像配准方法在SAR图像配准中是无效的；基于图像边缘特征点的同名点匹配法非常复杂；纹理驱动的SAR图像自动配准方法可用于具有更多的亮点或特征线的SAR图像配准，但这些特征在正常情况下都不能被发现；由于SAR图像的信噪比低，基于不变矩的图像自动配准方法封闭区域的特征点不能从SAR图像处理过程中获得；结合小波迭代细化和信息融合的同源的SAR图像配准方法，但很难配准具有不同的分辨率的图像；现有技术图像配准方法在主图像与副图像之间存在较大变形、旋转和缩放等问题时，不再适用于实时SAR图像配准；图像配准直接决定InSAR技术的鲁棒性和精度，决定InSAR的应用和未来发展，当前InSAR图像配准的低准确性制约了InSAR技术的发展，制约了其实际应用范围；
第三，现有技术缺少多项式拟合求主从图像分别在距离和方位向的偏移多项式的方法,无法重采样和干涉计算得到干涉图，无法提取出InSAR复数图像的幅度，无法分块对所得的主从图像进行关键点的检验和配准，现有技术缺少尺度三维图像极值检验、InSAR图像特征点定位、InSAR方向角度确定、图像特征点InSAR描述符、SAR图像配准的可靠方法，不能充分采用图像的相位和幅值信息，无法在高效率与配准结果的精度之间平衡，配准不准确导致严重的失相干，获得的干涉相位图中条纹边缘存在大量的噪声，条纹模糊，需要轨道参数等额外信息，同时配准所需的时间长，而且精度低；
第四，现有技术往往受到干涉图上大量的奇异点SPS所影响，而这些SPS点就是来自于精度不高的配准结果，目前没有特别有效消除这类奇异点的方法，现有技术处理图像配准过程中SPS点的方法效果不理想，现有技术缺少附加SPS量能评价的非线性InSAR配准方法，缺少SPS的量能作为配准评价标准，无法消除配准之后相干计算得出的干涉图中的SPS点,无法基于大量的奇异旋转点对图像进行非线性配准，导致生成的干涉图质量较低，得到的DEM信噪比低。
发明内容
针对现有技术的不足，本申请结合InSAR图像自身特点分步骤递进提出两种有效提高配准精度的方法，第一阶段是采用基于改进特征协同变换方法进行InSAR图像配准，采用多项式拟合求出主从图像分别在距离和方位向的偏移多项式,然后进行重采样和干涉计算得到干涉图，提取出InSAR复数图像的幅度，对主从图像进行关键点的检验和配准，分块将图像分成若干大小适当区域进行特征点匹配；第二阶段是在改进特征协同变换方法的基础上结合基于附加SPS量能评价的非线性InSAR配准法来提高图像局部配准的精度，通过消除配准之后相干计算得出的干涉图中的SPS点,基于大量的奇异旋转点对图像进行非线性配准；两种方法具有逐步递进提高InSAR复数图像配准精度的效果，逐步递进生成质量较高的干涉图和DEM，图像配准精度高，错误率低，具有很大的实际运用价值。
为实现以上技术效果，本申请所采用的技术方案如下：
InSAR干涉测量复数图像高精度配准方法，结合InSAR图像自身的特点分步骤递进提出两种有效提高配准精度的方法，第一阶段是采用基于改进特征协同变换方法进行InSAR图像配准，采用多项式拟合求出主从图像分别在距离和方位向的偏移多项式,然后进行重采样和干涉计算得到干涉图，提取出InSAR复数图像的幅度，对所得的主从图像进行关键点的检验和配准，采取分块处理方式，将图像分成若干大小适当区域进行特征点匹配；第二阶段是在改进特征协同变换方法的基础上结合基于附加SPS量能评价的非线性InSAR配准法来提高图像局部配准的精度，采用SPS的量能作为配准评价标准，通过消除配准之后相干计算得出的干涉图中的SPS点,基于大量的奇异旋转点对图像进行非线性配准，与改进特征协同变换配准方法结合提高InSAR图像配准之后所生成的干涉图的质量；本申请提出的两种方法具有逐步递进提高InSAR复数图像配准精度的效果，逐步递进生成质量较高的干涉图和DEM；
第一阶段，基于改进特征协同变换的InSAR图像配准方法：一是尺度三维图像极值检验，二是InSAR图像特征点的定位，三是InSAR方向角度的确定，四是图像特征点的InSAR描述符，五是改进特征协同变换的SAR图像配准；
第二阶段，基于附加SPS量能评价的非线性InSAR配准法：一是SPS点的设置与求解，二是基于SPS量能评价的InSAR配准操作，三是附加SPS的非线性InSAR图像精配准。
InSAR干涉测量复数图像高精度配准方法，进一步的，尺度三维图像极值检验：构建尺度三维高斯金字塔模型，首先对InSAR图像进行不同尺度的高斯模糊和降采样，得出一组大小不同的图像，从大到小，由下至上构成一个金字塔形状结构，一共有O组，每一组又分为S层，初始图像放在金字塔的底层，即第一层，每一次降采样后得到的图像作为金字塔的下一层，金字塔一共有n层，其总层数是由初始InSAR图像的大小及最高层图像的大小决定，其计算式为：
n＝log2{min(M，N)}-t，t∈[0，log2{min(M，N)} 式1
m,n为初始图像的大小，t为最高层图像的最小维数对数值，把目标图像的尺度三维L(x,y,σ)定义为：
L(x，y，σ)＝G(x，y，σ)*I(x，y) 式2
G(x,y,σ)是尺度可变高斯函数，I(x,y)为初始图像：
式中(x,y)是坐标位置，参数σ表示图像的平滑程度，不同σ值代表图像的不同尺度，表现图像的不同特征，式2的I(x,y)的分辨率是无限的，σ＝0，但尺度为零的图像不存在，图像的初始尺度设定为0.5，由小尺度σ1生成大尺度σ2的过程为：
采用高斯三维差分在不同尺度下取图像关键点，对相邻高斯三维尺度的图像相减，得到式5：
D(x,y,σ)为响应值图像，k为相邻尺度三维倍数，对高斯金字塔中上下相邻的图像求差值，得到高斯差分图，所求关键点由高斯差分三维尺度的局部极值点组成，求解关键点的第一步是对单组内各相邻的图像之间进行对比，为求出高斯三维差分函数的极值点,在高斯三维差分金字塔中除去部分像素值过小的点，这些像素之间的对比度低，是非稳定点，所有像素点都要和其周围的点进行比较,判断其图像域与尺度域是否与相邻点不同，正中心的检验点和同尺度及相邻尺度一共26个像素点比较，保证在全局三维空间都能对极值点进行成功检验，从每组中提取出x个尺度的极值点，每组高斯差分金字塔有x+2层，而每组高斯金字塔则有x+3层图像，3≤x≤5。
InSAR干涉测量复数图像高精度配准方法，进一步的，InSAR图像特征点的定位：对离散三维所提取的极值点进行拟合，将求出的离散三维极值点通过插值求出真实极值点，对D(x)进行尺度三维展开,得到：
式中的X＝(x,y,σ)T，通过对式6的求导使方程等于零，算得X的偏移量为：
其对应极值点的方程的值为：
InSAR干涉测量复数图像高精度配准方法，进一步的，InSAR方向角度的确定：采用每个特征点所在层的高斯尺度图像内的局部结构来得出其方向基准，该高斯尺度图像的尺度σ已知，并且该尺度是相对于高斯金字塔所在组的基准层的尺度，局部结构是在高斯尺度图像中以特征点为中心，以r为半径的区域内计算所有像素梯度的幅角和幅值，半径r为：
r＝3×1.5σ 式9
其中σ是相对于所在组的基准层的高斯尺度图像的尺度，像素梯度的幅值和幅角的计算公式为：
式中L的尺度是每一个关键点分别所在的尺度，完成图像关键点检验，所有关键点含有三类信息：位置、尺度和方向，用这些信息得出对应的改进特征协同变换特征区域；
在计算得出特征点梯度值后,对结果进行统计,并用直方图的形式表现特征点领域内像素每一个梯度方向所对应的幅值大小，用高斯函数对直方统计图进行平滑处理消除噪声，直方图的最大值处所指的梯度方向就是关键点的梯度方向。
InSAR干涉测量复数图像高精度配准方法，进一步的，图像特征点的InSAR描述符：将关键点邻域进行分区，根据每一区域内的梯度信息得到一个单独的向量，该向量是对其对应领域内图像信息的唯一描述，确定特征点描述符的具体过程为：首先定位一块求解描述子区域，在特征点所在尺度图像上进行对应特征点描述，将特征点周围划为d×d个小块，每一小块都是一个边长为3σ的正方形，即在每一小块的每一条边上都有3σ个像素点，其中σ是特征点对应高斯金字塔组的最低端层的尺度大小，特征点周围区域的边长是3σ(d+1)，在特征点周围的那块区域总共有3σ(d+1)×3σ(d+1)个像素点；
把获得的特征点周围邻域以特征为中点旋转到对应特征点的主方向,为使得经过旋转之后的正方形全部被包含进区域内，使这个区域的半径为正方形中点到对角点的距离，即对角线的1/2，即式12：
旋转后的关键点的坐标变为：
式中[x’,y’]T是进行坐标轴旋转之后的坐标；
把区域里全部的采样点都分散到各自的子邻域,把子邻域里的梯度按值给其分散到8个角度上,求出它的值,采样点在区域旋转之后被分配到半径为r的圆内的一个d×d的子块内，求出会改变子块内的采样点的梯度及方向，分至8个不同方向，进行旋转操作后的采样点在子方格块内的坐标用式14表示：
将式14所得的采样点线性插值，计算其对每个种子点的贡献值，求得特征向量后，对其做归一化操作消除光照改变对结果的影响，根据特征点尺度完成对描述向量排序后，生成改进特征协同变换特征描述向量。
InSAR干涉测量复数图像高精度配准方法，进一步的，改进特征协同变换的SAR图像配准：采用多项式拟合求出主从图像分别在距离和方位向的偏移多项式,然后进行重采样和干涉计算得到干涉图，提取出SAR复数图像的幅度，生成参考图像和输入图像的幅度图，对所得的主从图像进行关键点的检验和配准，待配准的主从图像间的像素偏移量,在距离向和方位向上没有太大的偏移和尺度变化，即主图像的特征点在辅图像中对应位置的附近，将图像分成若干大小适当区域进行特征点匹配，得到每一个特征点在另一幅图像中对应的同名点后，按照式15设置合适的临界值除去由于误差导致的错误点：
U＝{Δx，Δy|abs(Δx-median{Δx})＜R，abs(Δy-median{Δy)＜R} 式15
式中U为改进特征协同变换匹配到的点集合，(△x,△y)为同名的的偏移量，R为临界值，如果匹配点满足式15则认为是正确匹配点，除去误配点之后获得成功匹配点；
根据最小二乘法对这些检测到的匹配点按照式16进行四参数二次多项式拟合,得到主从图像的偏移多项式参数，用生成的多项式对从图像进行重采样，生成其分别对应的相干图，在获得一定数目的正确配准点对之后进行进一步的配准点对精化,对所有的配准点对进行均方根误差排序，选取最小的16个点作为精细配准点。
InSAR干涉测量复数图像高精度配准方法，进一步的，SPS点的设置与求解：通过减少SPS点的数量来消除局部失真最后提高DEM的质量；
首先将主从图像进行配准，然后对主从图像产生的干涉图I采用式17进行16点均值滤波：
式17中的各个参数解释：分别将主从图像上的8×2个像素看成一块全局M(m,n)和S(m,n),每一对相对应的块生成16点集干涉图中的一个像素点，而且在图中M(p,q；m,n)表示的是在M(m,n)中从上至下第p行第q列的像素，16点均值滤波通过求解干涉图上16个像素的平均值降低噪声，即产生的相位图在方位角方向有八倍的压缩，在距离向上有两倍的压缩﹐然后求出它的SPS点，再采用局部非线性配准插值方法处理之后的主从图像。
InSAR干涉测量复数图像高精度配准方法，进一步的，基于SPS量能评价的InSAR配准操作：得出SPS点数量进行基于SPS点量能的局部非线性配准,经过八次迭代插值计算的主从图像为1/8像素坐标系统，而初始图像为1像素坐标系统，I/8(m,n)表示1/8像素坐标系统产生干涉图的像素，I(m,n)是1像素坐标系统产生的干涉图的像素，(m,n)定义的是在干涉图全局图像中的定位，(p,q)定义的是在1/8像素坐标系统局部的位置，将主从图像分块成每块拥有(8×8)×(2×8)＝1024个子像素,然后进行运算得到I/8(m,n)，将主从图像上这些拥有1024个像素块称为M/8(m,n)和S/8(m,n)；
首先做一个1/8像素坐标系统的局部干涉图，像素值为I/8(m,n)其值和16点均值滤波之后的干涉图中像素值I(m,n)关联，其关系式为：
式18中的M/8(p,q；m,n)表示的是M/8(m,n)块中第p行第q列的像素，S/8(p,q；m,n)同理，如果在块I/8(m,n)，I/8(m+1,n)，I/8(m,n+1)，I/8(m+1,n+1)中存在一个SPS点，闭环旋转插值的值是非零，即：
|(θ/8(m+1，n)-θ/8(m,n))+(θ/8(m，n)-θ/8(m，n+1))+(θ/8(m+1，n+1)-θ/8(m+1，n))+(θ/8(m，n+1)-θ/8(m+1，n+1))|≥2π 式19
其中θ/8(m,n)＝arg(I/8(m,n)).即在I/8(m,n)-I/8(m+1,n+1)这四个区域内所对应的主从图像内有局部的失真，移动从图像S/8(m,n)-S/8(m+1,n+1)中的一块或多块以减少干涉图中的失真，1/8像素长度的距离移动S/8(m,n)中的点S/8(p,q；m,n)产生S/8 new(p,q；m,n)，然后检查生成干涉图中SPS点的数量是否有随着这一操作减少，如果没有减少，继续将剩下的三个块分别进行移动操作。
InSAR干涉测量复数图像高精度配准方法，进一步的，在采用迭代干涉运算去除SPS点之后,仍然有部分SPS点无法去除，同时对4|(2×2)个块进行移动，将其视为一个大的块，对S/8(m,n)，S/8(m,n+1)，S/8(m+1,n)，S/8(m+1,n+1)进行相同的移动操作，将这四块视为一个大的块,进行上述的相干运算的到干涉图,进一步抹去SPS点,如果还是无法达到要求，采用9|(3×3)更大的块来进行运算，查验最后的效果；
如果上述移动都不能使的SPS点减少，将从图像向其它七个方向移动，然后对于余下的SPS点，尝试着移动2/8像素长度的距离，以前面的方式操作，如果移动距离高达8/8＝1个像素距离之后删除所有的SPS点，开始删除下一个干涉图的SPS点，提高InSAR的配准精度。
与现有技术相比，本申请的创新点和优势在于：
第一，本申请结合InSAR图像自身特点分步骤递进提出两种有效提高配准精度的方法，第一阶段是采用基于改进特征协同变换方法进行InSAR图像配准，采用多项式拟合求出主从图像分别在距离和方位向的偏移多项式,然后进行重采样和干涉计算得到干涉图，提取出InSAR复数图像的幅度，对所得主从图像进行关键点的检验和配准，分块将图像分成若干大小适当区域进行特征点匹配；第二阶段是在改进特征协同变换方法的基础上结合基于附加SPS量能评价的非线性InSAR配准法来提高图像局部配准的精度，通过消除配准之后相干计算得出的干涉图中的SPS点,基于大量的奇异旋转点对图像进行非线性配准，与改进特征协同变换配准方法结合提高InSAR图像配准之后所生成的干涉图的质量；两种方法具有逐步递进提高InSAR复数图像配准精度的效果，逐步递进生成质量较高的干涉图和DEM，图像配准精度高，错误率低，具有很大的实际运用价值；
第二，本申请提出基于改进特征协同变换的InSAR图像配准方法，采用多项式拟合求出主从图像分别在距离和方位向的偏移多项式,然后重采样和干涉计算得到干涉图，提取InSAR复数图像的幅度，对主从图像进行关键点检验和配准，分块将图像分成若干大小适当区域进行特征点匹配，通过尺度三维图像极值检验、InSAR图像特征点定位、InSAR方向角度确定、图像特征点InSAR描述符、SAR图像配准，充分利用图像的相位和幅值信息，在高效率与配准结果的精度之间平衡，配准后获得干涉图条纹更加清晰，条纹边缘的斑点噪声相对较少，改进特征协同变换进行配准时不需要轨道参数等额外信息，同时改进特征协同变换配准所需的时间比较短，而且拥有更高的精度；
第三，本申请提出用一种基于附加SPS量能评价的非线性InSAR配准法，结合改进特征协同变换配准方法提高InSAR图像配准精度，通过SPS点的设置与求解、基于SPS量能评价的InSAR配准操作、附加SPS的非线性InSAR图像精配准，采用SPS点的量能作为主从图像以局部和非线性方式配准过程中的评价标准,可以生成质量更高的干涉图，这是InSAR图像配准的一个附加技巧，产生更高质量的干涉图,提高信噪比，该方法在提高InSAR复数图像配准精度方面具备有效性；
第四，本申请InSAR能够在相差微小的时间差甚至同时获得同一被监测目标地点的两幅甚至更多图像，得到干涉图消除了干扰噪声，使得干涉图上的信息只与所监测的目标地点的位置有关，从而充分利用获得的SAR图像中的数据信息，能准确、经济、高效、快速的解决大范围区域地形与形变中棘手问题，能够适应各种天气气候及工作时段，监测范围广，能较快获得信息数据；提供的SAR图像数据时间跨度大，可用更少的花费来获得连续的、高精度的观测信息，针对影响InSAR数据质量关键因素的图像配准，本申请进行了一系列优化改进，InSAR技术的鲁棒性和精度大幅提高，促进InSAR技术的发展，扩大其实际应用范围。
附图说明
图1是改进特征协同变换的InSAR图像配准过程示意图。
图2是生成参考图像和输入图像的主从图像示意图。
图3是基于改进特征协同变换的InSAR图像配准效果图。
图4是相干系数法InSAR图像配准效果图。
图5是改进特征协同变换精配准点对结果示意图。
图6是两种配准法的干涉相位图进行局部放大示意图。
图7是本申请去平地相位后的干涉结果图。
图8是截取部分的条纹及相位解缠结果示意图。
图9是通过最小二乘法解缠之后的整幅图恢复出DEM模型示意图。
图10是干涉图和SPS点配准时主从图像移动变换示意图。
图11是基于附加SPS量能评价的非线性InSAR配准法的流程图。
图12是本申请的块与干涉图像像素之间的关系图。
图13是基于SPS量能评价的图像符号及状态定义图。
图14是InSAR配准操作移动从影像过程示意图。
图15是SPS量能评价的InSAR配准干涉图及SPS点分布图。
图16是基于附加SPS量能评价的非线性InSAR配准后得到的DEM。
图17是SPS量能评价的InSAR配准前后平均信噪比和峰值信噪比对比图。
具体实施方法
下面结合附图，对本申请提出的InSAR干涉测量复数图像高精度配准方法的技术方案进行进一步的描述，使本领域的技术人员能更好的理解本申请并能够予以实施。
合成孔径雷达干涉测量(InSAR)是在遥感应用领域广泛获得地表高程及形变信息的重要技术，它不但能适应各种天气及工作时段、还能产生质量很高的数字高程模型(DEM)以满足人们日益增长的需求。采用InSAR可获得大范围内的地表高程和形变信息,在对地表高程值监测、地面活动监测、海洋洋流监测方面取得了大量成果,显示了巨大的潜力。
InSAR数据处理过程包括InSAR复数图像数据配准、干涉图的生成和去平地效应等相关处理、干涉相位的解缠绕，最后高程计算地理编码处理等。其中InSAR的图像配准过程是将主从复数图像中的公共区域的像素点一一对应起来，然后获得整个公共区域的干涉相位,该步骤结果的准确性直接影响生成的干涉图的质量以及InSAR后续的数据处理工作，所以InSAR复数图像配准是最关键、最基本的一个步骤，一般都要达到像元级别的配准精度，后续生成的干涉图和DEM才会有较好的质量。但伴随着计算机视觉与信息处理技术的发展，现有技术的图像配准方法已经不能达到生成高质量干涉图和DEM的标准。
本申请结合InSAR图像自身的特点分步骤递进提出两种有效提高配准精度的方法，第一阶段是采用基于改进特征协同变换方法进行InSAR图像配准，第二阶段是采用改进特征协同变换方法结合基于附加SPS量能评价的非线性InSAR配准法来提高图像局部配准的精度，采用SPS的量能作为配准评价标准，作为图像配准过程中的一个附加技巧，与改进特征协同变换配准方法结合可以有效提高InSAR图像配准之后所生成的干涉图的质量。为了证明两阶段方法在提高InSAR复数图像配准方面的有效性，分别采用传统相干系数法、基于改进特征协同变换的InSAR图像配准方法、基于附加SPS量能评价的非线性InSAR配准法进行实验，通过对比生成的干涉图的质量和信噪比得出有效的结论，最终的实验结果表明本申请提出的两种方法都能提高InSAR复数图像配准的精度，且具有递进增强的效果，并且生成质量较高的干涉图和DEM。
一、基于改进特征协同变换的InSAR图像配准方法
InSAR现有技术的一些图像配准方法不能充分采用图像的相位和幅值信息，无法在高效率与配准结果的精度之间平衡，本申请提出采用改进特征协同变换方法进行InSAR复数图像配准，并通过对比实验证明改进特征协同变换InSAR复数图像配准的可行性与有效性。
(一)尺度三维图像极值检验
构建尺度三维高斯金字塔模型，首先对InSAR图像进行不同尺度的高斯模糊和降采样，得出一组大小不同的图像，从大到小，由下至上构成一个金字塔形状结构，一共有O组，每一组又分为S层，初始图像放在金字塔的底层，即第一层，每一次降采样后得到的图像作为金字塔的下一层，金字塔一共有n层，其总层数是由初始InSAR图像的大小及最高层图像的大小决定，其计算式为：
n＝log2{min(M，N)}-t、t∈[0，log2{min(M，N)} 式1
m,n为初始图像的大小，t为最高层图像的最小维数对数值，把目标图像的尺度三维L(x,y,σ)定义为：
L(x，y，σ)＝G(x，y，σ)*I(x，y) 式2
G(x,y,σ)是尺度可变高斯函数，I(x,y)为初始图像：
式中(x,y)是坐标位置，参数σ表示图像的平滑程度，不同σ值代表图像的不同尺度，表现图像的不同特征，式2的I(x,y)的分辨率是无限的，σ＝0，但尺度为零的图像不存在，图像的初始尺度设定为0.5，由小尺度σ1生成大尺度σ2的过程为：
采用高斯三维差分在不同尺度下取图像关键点，对相邻高斯三维尺度的图像相减，得到式5：
D(x,y,σ)为响应值图像，k为相邻尺度三维倍数，对高斯金字塔中上下相邻的图像求差值，得到高斯差分图，所求关键点由高斯差分三维尺度的局部极值点组成，求解关键点的第一步是对单组内各相邻的图像之间进行对比，为求出高斯三维差分函数的极值点,在高斯三维差分金字塔中除去部分像素值过小的点，这些像素之间的对比度低，是非稳定点，所有像素点都要和其周围的点进行比较,判断其图像域与尺度域是否与相邻点不同，正中心的检验点和同尺度及相邻尺度一共26个像素点比较，保证在全局三维空间都能对极值点进行成功检验，从每组中提取出x个尺度的极值点，每组高斯差分金字塔有x+2层，而每组高斯金字塔则有x+3层图像，3≤x≤5。
(二)InSAR图像特征点的定位
在改进特征协同变换过程中，需要得到特征点准确的尺度与位置信息，即需要达到亚像素级的准确度，所以接下来必须对离散三维所提取的极值点进行拟合，将求出的离散三维极值点通过插值求出真实极值点，对D(x)进行尺度三维展开,得到：
式中的X＝(x,y,σ)T，通过对式6的求导使方程等于零，算得X的偏移量为：
其对应极值点的方程的值为：
(三)InSAR方向角度的确定
采用每个特征点所在层的高斯尺度图像内的局部结构来得出其方向基准，该高斯尺度图像的尺度σ已知，并且该尺度是相对于高斯金字塔所在组的基准层的尺度，局部结构是在高斯尺度图像中以特征点为中心，以r为半径的区域内计算所有像素梯度的幅角和幅值，半径r为：
r＝3×1.5σ 式9
其中σ是相对于所在组的基准层的高斯尺度图像的尺度，像素梯度的幅值和幅角的计算公式为：
式中L的尺度是每一个关键点分别所在的尺度，完成图像关键点检验，所有关键点含有三类信息：位置、尺度和方向，用这些信息得出对应的改进特征协同变换特征区域。
在计算得出特征点梯度值后,对结果进行统计,并用直方图的形式表现特征点领域内像素每一个梯度方向所对应的幅值大小，为消除噪声对梯度方向角度的影响，用高斯函数对直方统计图进行平滑处理，直方图的最大值处所指的梯度方向就是关键点的梯度方向。
(四)图像特征点的InSAR描述符
将关键点邻域进行分区，根据每一区域内的梯度信息得到一个单独的向量，该向量是对其对应领域内图像信息的唯一描述，确定特征点描述符的具体过程为：首先定位一块求解描述子区域，在特征点所在尺度图像上进行对应特征点描述，将特征点周围划为d×d个小块，每一小块都是一个边长为3σ的正方形，即在每一小块的每一条边上都有3σ个像素点，其中σ是特征点对应高斯金字塔组的最低端层的尺度大小，特征点周围区域的边长是3σ(d+1)，在特征点周围的那块区域总共有3σ(d+1)×3σ(d+1)个像素点。
为使特征点能保持旋转不变特性，把获得的特征点周围邻域以特征为中点旋转到对应特征点的主方向,为使得经过旋转之后的正方形全部被包含进区域内，使这个区域的半径为正方形中点到对角点的距离，即对角线的1/2，即式12：
旋转后的关键点的坐标变为：
式中[x’,y’]T是进行坐标轴旋转之后的坐标；
把区域里全部的采样点都分散到各自的子邻域,把子邻域里的梯度按值给其分散到8个角度上,求出它的值,采样点在区域旋转之后被分配到半径为r的圆内的一个d×d的子块内，求出会改变子块内的采样点的梯度及方向，分至8个不同方向，进行旋转操作后的采样点在子方格块内的坐标用式14表示：
将式14所得的采样点线性插值，计算其对每个种子点的贡献值，求得特征向量后，对其做归一化操作消除光照改变对结果的影响，根据特征点尺度完成对描述向量排序后，生成改进特征协同变换特征描述向量。
(五)改进特征协同变换的SAR图像配准
InSAR配准的两幅图像没有太大的位移量和尺度大小的改变，主从图像上的关键点位置差距不太大，用分块匹配模式，提高运算效率和准确率，改进特征协同变换方法的匹配过程如图1所示。
采用改进特征协同变换配准完成后，采用多项式拟合求出主从图像分别在距离和方位向的偏移多项式,然后进行重采样和干涉计算得到干涉图。
提取出SAR复数图像的幅度，生成参考图像和输入图像的幅度图，如图2所示。
对所得的主从图像进行关键点的检验和配准，待配准的主从图像间的像素偏移量,在距离向和方位向上没有太大的偏移和尺度变化，即主图像的特征点在辅图像中对应位置的附近，因此采取分块处理方式，将图像分成若干大小适当区域进行特征点匹配，这样既可以大幅提高配准效率，又可以降低误配率，得到每一个特征点在另一幅图像中对应的同名点后，按照式15设置合适的临界值除去由于误差导致的错误点：
U＝{Δx，Δy|abs(Δx-median{Δx})＜R，abs(Δy-median{Δy})＜R} 式15
式中U为改进特征协同变换匹配到的点集合，(△x,△y)为同名的的偏移量，R为临界值，如果匹配点满足式15则认为是正确匹配点，除去误配点之后获得成功匹配点。
根据最小二乘法对这些检测到的匹配点按照式16进行四参数二次多项式拟合,得到主从图像的偏移多项式参数。
分别用本申请方法和相干系数法生成的多项式对从图像进行重采样，生成其分别对应的相干图，相干系数越高，干涉数据的质量越好，将配准后的主、从图像应用式生成相干图。
如图3和图4所示。相干图中越亮的位置表示该处的相干系数越高，即干涉相位的质量越好，而较暗的位置可能是由相位噪声引起的干涉相位不连续或地形起伏剧烈引起的干涉相位不连续，从两幅相干图可以看出，改进特征协同变换匹配得到的相干图其高亮区域比相干系数的多，同样也可以统计对比两种方法得到的相干系数值，做出相关系数分布直方图来对比来评价两种匹配方法的精度。
从上面两张图可以看出，用相干系数法配准得到的相干系数值大量的分布在0.5到0.7之前，而由改进特征协同变换配所多得到的的相干系数值则大量的分布在0.7到0.9直接，相干系数越高，配准效果越好，后续得到的干涉图质量也会越高,所以通过对比之后可以得出改进特征协同变换配准效果要比传统方法好，相干系数法在配准过程中的不准确会导致严重的失相干。
在获得一定数目的正确配准点对之后进行进一步的配准点对精化,对所有的配准点对进行均方根误差排序，选取最小的16个点作为精细配准点,检测出的配准点如图5。
(六)改进特征协同变换配准方法和相干系数法实验对比
1.生成干涉图与去平地效应
通过将相干系数和改进特征协同变换配准后的主从图像的复数数据进行共辄相乘，并提取出干涉相位，便得到了干涉相位图，相位值通过灰度值显示，其中一个白色的区域代表相位值为π，而黑色的区域代表相位值为-π。
分别将两种配准法的干涉相位图进行局部放大,如图6所示。比较两幅干涉相位图的条纹发现：在相干系数法获得的干涉相位图中，条纹边缘存在大量的噪声，使条纹较为模糊；而改进特征协同变换方法配准后获得干涉图，其中的干涉条纹更加清晰，条纹边缘的斑点噪声相对较少，配准的结果较好。
这些干涉图都是由平坦的地面与地面高程信息叠加在一起产生,因此得到干涉相位图之后，进行去平地效应才能得到真实的相位，平地相位通过估算条纹的频率，然后根据式16去除：
式16中，xn为图像强度，
为条纹频率估算值，列方向为距离向；m,n分别代表图像行列方向的位置，m＝0,1,…,M-1,n＝0,1,…,N-1。
去平地相位后的干涉图如图7所示，对比图6可以发现，初始的干涉相位图的条纹几乎都是由平地相位产生,它平行于方位向，延伸着距离向一级一级地分布,与平地相位长生的原理相符合，一旦平地相位去除，便可以清楚的从干涉图看出地势和高程变化的趋势。
由于配准过程中，失配准的点在干涉相位图上会构成相位噪声，这些相位噪声会对后续相位解缠绕造成误差，并且传播到整幅相位图中，严重影响相位解缠绕结果的可靠性。因此尽可能滤除相位噪声点，提高图像的信噪比。
2.相位解缠对比
干涉图的相位值是经过对复数图像的相位反正切函数运算得到的,只是干涉的相位主值，需要对相位恢复才能得到真实的相位值，初始干涉图的分辨率较高，对整幅干涉图进行解缠非常耗时，且不方便对解缠结果进行对比分析,于是从初始干涉图中截取部分干涉条纹图采用最小二乘法进行相位解缠。由于InSAR对图像配准精度要求很高，配准的精度直接影响干涉图的质量，从而对解缠结果影响很大，分别使用枝切法解缠，对最终结果进行比较。截取部分的条纹及相位解缠结果如图8所示：
3.修复DEM模型
在得到真实的干涉相位值后，去平地相位后的干涉相位值与地表高程成正比，解缠后的相位图反应地表各点间的相对高度关系,即地形变化，根据通过最小二乘法解缠之后的整幅图恢复出DEM模型，如图9。
改进特征协同变换方法配准后InSAR图像相对于传统的相干系数法拥有质量更高的干涉条纹，而且，改进特征协同变换进行配准时不需要轨道参数等额外信息，同时改进特征协同变换配准所需的时间比较短，而且拥有更高的精度。
二、基于附加SPS量能评价的非线性InSAR配准法
要得到一个准确的数字高程模型DEM就必须拥有质量较高的干涉图，但其结果往往受到干涉图上大量的奇异点SPS所影响，而这些SPS点就是来自于精度不高的配准结果，然而目前没有特别有效消除这类奇异点的方法，虽然现有技术有针对如何处理图像配准过程中SPS点的方法，但效果不理想。本申请提出用一种基于附加SPS量能评价的非线性InSAR配准法，结合改进特征协同变换配准方法提高InSAR图像配准精度，这种方法可以生成质量更高的干涉图，以产生精度更高的DEM,并通过实验表明该方法可以得到信噪较高的DEM。
实际干涉图包含大量SPS非零旋转点，绝大多数SPS会导致不能恢复出一个高精度的DEM，SPS点的产生有两个起源：一是违反抽样定理，由于地形信息过于粗糙而导致相位混叠现象而产生；二是由于存在主和从观察轨道之间的差异，在用两个稍微不同的视角观测同一目标地点时其反射、散射、和透视值不同,导致干涉图有局部相位畸变。通常配准是通过大尺度的图像振幅互相关最大化而小尺度的将复振幅相关性最大化，在求解相关性时在一定范围内求平均，以减少所包含的噪声，但大区域的求平均操作会丢失一些在图像配准中用来消除失真的点，这种代价会使得图像配准的精度受到限制，反射、散射和透视差异导致的局部变形，在互相关过程中产生SPS点。本申请采用SPS点的量能作为主从图像以局部和非线性方式配准过程中的评价标准,产生更高质量的干涉图,这是InSAR图像配准的一个附加技巧，提高信噪比。
(一)SPS点的设置与求解
通过截取相位干涉图中的一小块进行实验分析，发现对干涉图上所截取的一小块部分对应主从图像进行处理可以有效减少SPS点数量；图10显示了干涉图和SPS点是如何随着图像配准时从图像在与主图像最佳的配准位置左右以1/8像素的整数倍距离移动而产生变换的，实验表明，许多SPS移动，出现或者消失，但干涉图的条纹只会产生小的变化，在干涉图中白色的SPS点表示顺时针的非零旋转点，而黑色的代表逆时针的非零旋转点。
采用以上方法的干涉图上SPS点的出现和消失比起全局相关计算产生的结果更加局部，在传统配准方法中增加局部的非线性配准过程所得到的干涉图全局图形前后变化并不会太大，但这种局部精配准可以提高局部的干涉图质量，在一个理想化没有畸变的干涉图中不存在SPS点。本申请通过减少SPS点的数量来消除局部失真最后提高DEM的质量。
图11是基于附加SPS量能评价的非线性InSAR配准法的流程图，首先将主从图像进行配准，然后对主从图像产生的干涉图I采用式17进行16点均值滤波：
式17中的各个参数通过图12进一步解释：分别将主从图像上的8×2个像素看成一块全局M(m,n)和S(m,n),每一对相对应的块生成在图12(a)中那样的16点集干涉图中的一个像素点，而且在图中M(p,q；m,n)表示的是在M(m,n)中从上至下第p行第q列的像素，16点均值滤波通过求解干涉图上16个像素的平均值降低噪声，即产生的相位图在方位角方向有八倍的压缩，在距离向上有两倍的压缩﹐然后求出它的SPS点，再采用局部非线性配准插值方法处理之后的主从图像。
(二)基于SPS量能评价的InSAR配准操作
得出SPS点数量进行基于SPS点量能的局部非线性配准,结合图12进行说明，经过八次迭代插值计算的主从图像为1/8像素坐标系统，而初始图像为1像素坐标系统。图13列出了各名词的定义,I/8(m,n)表示1/8像素坐标系统产生干涉图的像素，I(m,n)是1像素坐标系统产生的干涉图的像素，(m,n)定义的是在干涉图全局图像中的定位，(p,q)定义的是在1/8像素坐标系统局部的位置，将主从图像分块成每块拥有(8×8)×(2×8)＝1024个子像素,然后进行运算得到I/8(m,n)，将主从图像上这些拥有1024个像素块称为M/8(m,n)和S/8(m,n)。
首先做一个1/8像素坐标系统的局部干涉图，像素值为I/8(m,n)其值和16点均值滤波之后的干涉图中像素值I(m,n)关联，其关系式为：
式18中的M/8(p,q；m,n)表示的是M/8(m,n)块中第p行第q列的像素，S/8(p,q；m,n)同理，如果在块I/8(m,n)，I/8(m+1,n)，I/8(m,n+1)，I/8(m+1,n+1)中存在一个SPS点，如图12(c)所示，闭环旋转插值的值是非零，即：
|(θ/8(m+1，n)-θ/8(m,n))+(θ/8(m，n)-θ/8(m，n+1))+(θ/8(m+1，n+1)-θ/8(m+1，n))+(θ/8(m，n+1)-θ/8(m+1，n+1))|≥2π 式19
其中θ/8(m,n)＝arg(I/8(m,n)).即在I/8(m,n)-I/8(m+1,n+1)这四个区域内所对应的主从图像内有局部的失真，移动从图像S/8(m,n)-S/8(m+1,n+1)中的一块或多块以减少干涉图中的失真。如图12中所示，1/8像素长度的距离移动S/8(m,n)中的点S/8(p,q；m,n)产生S/8 new(p,q；m,n)，然后检查生成干涉图中SPS点的数量是否有随着这一操作减少，如果没有减少，继续按照图14所示将剩下的三个块分别进行移动操作。
在采用上面所述的迭代干涉运算去除SPS点之后,会发现仍然有部分SPS点无法去除，同时对4|(2×2)个块进行移动，将其视为一个大的块，对S/8(m,n)，S/8(m,n+1)，S/8(m+1,n)，S/8(m+1,n+1)进行相同的移动操作，将这四块视为一个大的块,进行上述的相干运算的到干涉图,进一步抹去SPS点,如果还是无法达到要求，采用9|(3×3)更大的块来进行运算，查验最后的效果。
如果上述移动都不能使的SPS点减少，将从图像按照图14所示，向其它七个方向移动，然后对于余下的SPS点，尝试着移动2/8像素长度的距离，以前面的方式操作，如果移动距离高达8/8＝1个像素距离之后删除所有的SPS点，开始删除下一个干涉图的SPS点，提高InSAR的配准精度。
(三)附加SPS的非线性InSAR图像精配准
对图7中截取画面按照上面的步骤进行实验，在采用基于附加SPS量能评价的非线性InSAR配准法前的局部干涉图以及提取出来的SPS点如图15(a)所示，右边的图是其SPS点，总数一共是1016个，在进行上述实验步骤第一次的移动计算之后,抹去了59％的SPS点，剩下了401个SPS点，如图15(b)，进行第二轮计算之后剩下了325个SPS点，接下来的单块移动计算中，没有任何的SPS点得到消除抹去，所以移动接下来的三个块，即上面所说的4个块一起移动进行后续计算，得出来的干涉图仅剩下185个SPS点，如图15(d),甚至还可以9个块一起移动，得出来的干涉图其SPS点更少。
总之通过实验结果可以看出，采用基于附加SPS量能评价的非线性InSAR配准法可以有效的减少干涉图上面的SPS点，实例中达到了81.91％的擦除率，极大的提高了干涉图的质量。采用基于附加SPS量能评价的非线性InSAR配准法之后的全局干涉图，与图7的初始干涉图进行直观对比，发现采用了基于附加SPS量能评价的非线性InSAR配准法之后的干涉图条纹明显更加清晰。
选取与图7同一局部进行相位解缠,最后通过与最小二乘法进行相位解缠之后，最后恢复出高程图16。
得到DEM之后通过对比采用基于附加SPS量能评价的非线性InSAR配准法前后的平均信噪比和峰值信噪比来直观的体现其作用,如图17。图17可以直观的显示基于附加SPS量能评价的非线性InSAR配准法可以同时提高DEM的平均信噪比和峰值信噪比，说明了基于附加SPS量能评价的非线性InSAR配准法在提高InSAR复数图像配准精度方面的有效性。
Claims (9)
1.InSAR干涉测量复数图像高精度配准方法，其特征在于，结合InSAR图像自身的特点分步骤递进提出两种有效提高配准精度的方法，第一阶段是采用基于改进特征协同变换方法进行InSAR图像配准，采用多项式拟合求出主从图像分别在距离和方位向的偏移多项式,然后进行重采样和干涉计算得到干涉图，提取出InSAR复数图像的幅度，对所得的主从图像进行关键点的检验和配准，采取分块处理方式，将图像分成若干大小适当区域进行特征点匹配；第二阶段是在改进特征协同变换方法的基础上结合基于附加SPS量能评价的非线性InSAR配准法来提高图像局部配准的精度，采用SPS的量能作为配准评价标准，通过消除配准之后相干计算得出的干涉图中的SPS点,基于大量的奇异旋转点对图像进行非线性配准，与改进特征协同变换配准方法结合提高InSAR图像配准之后所生成的干涉图的质量；本申请提出的两种方法具有逐步递进提高InSAR复数图像配准精度的效果，逐步递进生成质量较高的干涉图和DEM；
第一阶段，基于改进特征协同变换的InSAR图像配准方法：一是尺度三维图像极值检验，二是InSAR图像特征点的定位，三是InSAR方向角度的确定，四是图像特征点的InSAR描述符，五是改进特征协同变换的SAR图像配准；
第二阶段，基于附加SPS量能评价的非线性InSAR配准法：一是SPS点的设置与求解，二是基于SPS量能评价的InSAR配准操作，三是附加SPS的非线性InSAR图像精配准。
2.根据权利要求1所述InSAR干涉测量复数图像高精度配准方法，其特征在于，尺度三维图像极值检验：构建尺度三维高斯金字塔模型，首先对InSAR图像进行不同尺度的高斯模糊和降采样，得出一组大小不同的图像，从大到小，由下至上构成一个金字塔形状结构，一共有O组，每一组又分为S层，初始图像放在金字塔的底层，即第一层，每一次降采样后得到的图像作为金字塔的下一层，金字塔一共有n层，其总层数是由初始InSAR图像的大小及最高层图像的大小决定，其计算式为：
n＝log2{min(M，N)}-t，t∈[0，log2{min(M，N)} 式1
m,n为初始图像的大小，t为最高层图像的最小维数对数值，把目标图像的尺度三维L(x,y,σ)定义为：
L(x，y，σ)＝G(x，y，σ)*I(x，y) 式2
G(x,y,σ)是尺度可变高斯函数，I(x,y)为初始图像：
式中(x,y)是坐标位置，参数σ表示图像的平滑程度，不同σ值代表图像的不同尺度，表现图像的不同特征，式2的I(x,y)的分辨率是无限的，σ＝0，但尺度为零的图像不存在，图像的初始尺度设定为0.5，由小尺度σ1生成大尺度σ2的过程为：
采用高斯三维差分在不同尺度下取图像关键点，对相邻高斯三维尺度的图像相减，得到式5：
D(x,y,σ)为响应值图像，k为相邻尺度三维倍数，对高斯金字塔中上下相邻的图像求差值，得到高斯差分图，所求关键点由高斯差分三维尺度的局部极值点组成，求解关键点的第一步是对单组内各相邻的图像之间进行对比，为求出高斯三维差分函数的极值点,在高斯三维差分金字塔中除去部分像素值过小的点，这些像素之间的对比度低，是非稳定点，所有像素点都要和其周围的点进行比较,判断其图像域与尺度域是否与相邻点不同，正中心的检验点和同尺度及相邻尺度一共26个像素点比较，保证在全局三维空间都能对极值点进行成功检验，从每组中提取出x个尺度的极值点，每组高斯差分金字塔有x+2层，而每组高斯金字塔则有x+3层图像，3≤x≤5。
3.根据权利要求1所述InSAR干涉测量复数图像高精度配准方法，其特征在于，InSAR图像特征点的定位：对离散三维所提取的极值点进行拟合，将求出的离散三维极值点通过插值求出真实极值点，对D(x)进行尺度三维展开,得到：
式中的X＝(x,y,σ)T，通过对式6的求导使方程等于零，算得X的偏移量为：
其对应极值点的方程的值为：
4.根据权利要求1所述InSAR干涉测量复数图像高精度配准方法，其特征在于，InSAR方向角度的确定：采用每个特征点所在层的高斯尺度图像内的局部结构来得出其方向基准，该高斯尺度图像的尺度σ已知，并且该尺度是相对于高斯金字塔所在组的基准层的尺度，局部结构是在高斯尺度图像中以特征点为中心，以r为半径的区域内计算所有像素梯度的幅角和幅值，半径r为：
r＝3×1.5σ 式9
其中σ是相对于所在组的基准层的高斯尺度图像的尺度，像素梯度的幅值和幅角的计算公式为：
式中L的尺度是每一个关键点分别所在的尺度，完成图像关键点检验，所有关键点含有三类信息：位置、尺度和方向，用这些信息得出对应的改进特征协同变换特征区域；
在计算得出特征点梯度值后,对结果进行统计,并用直方图的形式表现特征点领域内像素每一个梯度方向所对应的幅值大小，用高斯函数对直方统计图进行平滑处理消除噪声，直方图的最大值处所指的梯度方向就是关键点的梯度方向。
5.根据权利要求1所述InSAR干涉测量复数图像高精度配准方法，其特征在于，图像特征点的InSAR描述符：将关键点邻域进行分区，根据每一区域内的梯度信息得到一个单独的向量，该向量是对其对应领域内图像信息的唯一描述，确定特征点描述符的具体过程为：首先定位一块求解描述子区域，在特征点所在尺度图像上进行对应特征点描述，将特征点周围划为d×d个小块，每一小块都是一个边长为3σ的正方形，即在每一小块的每一条边上都有3σ个像素点，其中σ是特征点对应高斯金字塔组的最低端层的尺度大小，特征点周围区域的边长是3σ(d+1)，在特征点周围的那块区域总共有3σ(d+1)×3σ(d+1)个像素点；
把获得的特征点周围邻域以特征为中点旋转到对应特征点的主方向,为使得经过旋转之后的正方形全部被包含进区域内，使这个区域的半径为正方形中点到对角点的距离，即对角线的1/2，即式12：
旋转后的关键点的坐标变为：
式中[x’,y’]T是进行坐标轴旋转之后的坐标；
把区域里全部的采样点都分散到各自的子邻域,把子邻域里的梯度按值给其分散到8个角度上,求出它的值,采样点在区域旋转之后被分配到半径为r的圆内的一个d×d的子块内，求出会改变子块内的采样点的梯度及方向，分至8个不同方向，进行旋转操作后的采样点在子方格块内的坐标用式14表示：
将式14所得的采样点线性插值，计算其对每个种子点的贡献值，求得特征向量后，对其做归一化操作消除光照改变对结果的影响，根据特征点尺度完成对描述向量排序后，生成改进特征协同变换特征描述向量。
6.根据权利要求1所述InSAR干涉测量复数图像高精度配准方法，其特征在于，改进特征协同变换的SAR图像配准：采用多项式拟合求出主从图像分别在距离和方位向的偏移多项式,然后进行重采样和干涉计算得到干涉图，提取出SAR复数图像的幅度，生成参考图像和输入图像的幅度图，对所得的主从图像进行关键点的检验和配准，待配准的主从图像间的像素偏移量,在距离向和方位向上没有太大的偏移和尺度变化，即主图像的特征点在辅图像中对应位置的附近，将图像分成若干大小适当区域进行特征点匹配，得到每一个特征点在另一幅图像中对应的同名点后，按照式15设置合适的临界值除去由于误差导致的错误点：
U＝{Δx，Δy|abs(Δx-median{Δx})＜R，abs(Δy-median{Δy})＜R} 式15
式中U为改进特征协同变换匹配到的点集合，(△x,△y)为同名的的偏移量，R为临界值，如果匹配点满足式15则认为是正确匹配点，除去误配点之后获得成功匹配点；
根据最小二乘法对这些检测到的匹配点按照式16进行四参数二次多项式拟合,得到主从图像的偏移多项式参数，用生成的多项式对从图像进行重采样，生成其分别对应的相干图，在获得一定数目的正确配准点对之后进行进一步的配准点对精化,对所有的配准点对进行均方根误差排序，选取最小的16个点作为精细配准点。
7.根据权利要求1所述InSAR干涉测量复数图像高精度配准方法，其特征在于，SPS点的设置与求解：通过减少SPS点的数量来消除局部失真最后提高DEM的质量；
首先将主从图像进行配准，然后对主从图像产生的干涉图I采用式17进行16点均值滤波：
式17中的各个参数解释：分别将主从图像上的8×2个像素看成一块全局M(m,n)和S(m,n),每一对相对应的块生成16点集干涉图中的一个像素点，而且在图中M(p,q；m,n)表示的是在M(m,n)中从上至下第p行第q列的像素，16点均值滤波通过求解干涉图上16个像素的平均值降低噪声，即产生的相位图在方位角方向有八倍的压缩，在距离向上有两倍的压缩﹐然后求出它的SPS点，再采用局部非线性配准插值方法处理之后的主从图像。
8.根据权利要求1所述InSAR干涉测量复数图像高精度配准方法，其特征在于，基于SPS量能评价的InSAR配准操作：得出SPS点数量进行基于SPS点量能的局部非线性配准,经过八次迭代插值计算的主从图像为1/8像素坐标系统，而初始图像为1像素坐标系统，I/8(m,n)表示1/8像素坐标系统产生干涉图的像素，I(m,n)是1像素坐标系统产生的干涉图的像素，(m,n)定义的是在干涉图全局图像中的定位，(p,q)定义的是在1/8像素坐标系统局部的位置，将主从图像分块成每块拥有(8×8)×(2×8)＝1024个子像素,然后进行运算得到I/8(m,n)，将主从图像上这些拥有1024个像素块称为M/8(m,n)和S/8(m,n)；
首先做一个1/8像素坐标系统的局部干涉图，像素值为I/8(m,n)其值和16点均值滤波之后的干涉图中像素值I(m,n)关联，其关系式为：
式18中的M/8(p,q；m,n)表示的是M/8(m,n)块中第p行第q列的像素，S/8(p,q；m,n)同理，如果在块I/8(m,n)，I/8(m+1,n)，I/8(m,n+1)，I/8(m+1,n+1)中存在一个SPS点，闭环旋转插值的值是非零，即：
|(θ/8(m+1，n)-θ/8(m，n))+(θ/8(m，n)-θ/8(m，n+1))+(θ/8(m+1，n+1)-θ/8(m+1，n))+(θ/8(m，n+1)-θ/8(m+1，n+1))|≥2π 式19
其中θ/8(m,n)＝arg(I/8(m,n)).即在I/8(m,n)-I/8(m+1,n+1)这四个区域内所对应的主从图像内有局部的失真，移动从图像S/8(m,n)-S/8(m+1,n+1)中的一块或多块以减少干涉图中的失真，1/8像素长度的距离移动S/8(m,n)中的点S/8(p,q；m,n)产生S/8 new(p,q；m,n)，然后检查生成干涉图中SPS点的数量是否有随着这一操作减少，如果没有减少，继续将剩下的三个块分别进行移动操作。
9.根据权利要求8所述InSAR干涉测量复数图像高精度配准方法，其特征在于，在采用迭代干涉运算去除SPS点之后,仍然有部分SPS点无法去除，同时对4|(2×2)个块进行移动，将其视为一个大的块，对S/8(m,n)，S/8(m,n+1)，S/8(m+1,n)，S/8(m+1,n+1)进行相同的移动操作，将这四块视为一个大的块,进行上述的相干运算的到干涉图,进一步抹去SPS点,如果还是无法达到要求，采用9|(3×3)更大的块来进行运算，查验最后的效果；
如果上述移动都不能使的SPS点减少，将从图像向其它七个方向移动，然后对于余下的SPS点，尝试着移动2/8像素长度的距离，以前面的方式操作，如果移动距离高达8/8＝1个像素距离之后删除所有的SPS点，开始删除下一个干涉图的SPS点，提高InSAR的配准精度。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210297307.7A CN115115678A (zh) | 2022-03-24 | 2022-03-24 | InSAR干涉测量复数图像高精度配准方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210297307.7A CN115115678A (zh) | 2022-03-24 | 2022-03-24 | InSAR干涉测量复数图像高精度配准方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115115678A true CN115115678A (zh) | 2022-09-27 |
Family
ID=83325323
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210297307.7A Pending CN115115678A (zh) | 2022-03-24 | 2022-03-24 | InSAR干涉测量复数图像高精度配准方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115115678A (zh) |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107102333A (zh) * | 2017-06-27 | 2017-08-29 | 北京航空航天大学 | 一种星载InSAR长短基线融合解缠方法 |
| CN111896955A (zh) * | 2020-08-06 | 2020-11-06 | 武汉大学 | 一种船载sar交轨干涉处理方法 |
| CN111998766A (zh) * | 2020-08-31 | 2020-11-27 | 同济大学 | 一种基于时序InSAR技术的地表形变反演方法 |
| CN112017223A (zh) * | 2020-09-11 | 2020-12-01 | 西安电子科技大学 | 基于改进SIFT-Delaunay的异源图像配准方法 |
| US20200394780A1 (en) * | 2017-06-15 | 2020-12-17 | The University Of Nottingham | Land deformation measurement |
| CN112711021A (zh) * | 2020-12-08 | 2021-04-27 | 中国自然资源航空物探遥感中心 | 一种多分辨率InSAR交互干涉时序分析方法 |
-
2022
- 2022-03-24 CN CN202210297307.7A patent/CN115115678A/zh active Pending
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20200394780A1 (en) * | 2017-06-15 | 2020-12-17 | The University Of Nottingham | Land deformation measurement |
| CN107102333A (zh) * | 2017-06-27 | 2017-08-29 | 北京航空航天大学 | 一种星载InSAR长短基线融合解缠方法 |
| CN111896955A (zh) * | 2020-08-06 | 2020-11-06 | 武汉大学 | 一种船载sar交轨干涉处理方法 |
| CN111998766A (zh) * | 2020-08-31 | 2020-11-27 | 同济大学 | 一种基于时序InSAR技术的地表形变反演方法 |
| CN112017223A (zh) * | 2020-09-11 | 2020-12-01 | 西安电子科技大学 | 基于改进SIFT-Delaunay的异源图像配准方法 |
| CN112711021A (zh) * | 2020-12-08 | 2021-04-27 | 中国自然资源航空物探遥感中心 | 一种多分辨率InSAR交互干涉时序分析方法 |
Non-Patent Citations (3)
| Title |
|---|
| 王青松;瞿继双;黄海风;余安喜;董臻;: "联合实、复相关函数的干涉SAR图像配准方法", 测绘学报, no. 04, 15 August 2012 (2012-08-15) * |
| 石晓进;张云华;: "基于Fourier-Mellin变换和相干系数法的重复轨道干涉SAR图像配准新方法", 电子与信息学报, no. 04, 15 April 2009 (2009-04-15) * |
| 程海琴;陈强;刘国祥;杨莹辉;: "基于相干曲面移动拟合的SAR影像高精度配准方法", 测绘科学, no. 05, 27 December 2012 (2012-12-27) * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109633648B (zh) | 一种基于似然估计的多基线相位估计装置及方法 | |
| CN111208512B (zh) | 一种基于视频合成孔径雷达的干涉测量方法 | |
| CN106249236B (zh) | 一种星载InSAR长短基线图像联合配准方法 | |
| CN107341781B (zh) | 基于改进相位一致性特征矢量底图匹配的sar影像校正方法 | |
| CN113311433B (zh) | 一种质量图和最小费用流结合的InSAR干涉相位两步解缠方法 | |
| CN105425216A (zh) | 基于图像分割的重复航过极化InSAR图像配准方法 | |
| Zhang et al. | Satellite SAR geocoding with refined RPC model | |
| CN109100719A (zh) | 基于星载sar影像与光学影像的地形图联合测图方法 | |
| CN115457022B (zh) | 基于实景三维模型正视影像的三维形变检测方法 | |
| CN106204539A (zh) | 一种基于形态学梯度的反演城区建筑物沉降的方法 | |
| CN111562575B (zh) | 一种用于地面沉降的监测方法 | |
| CN112051572A (zh) | 一种融合多源sar数据三维地表形变监测方法 | |
| CN108022259A (zh) | 干涉sar复图像配准方法和系统 | |
| CN115540908A (zh) | 基于小波变换的InSAR干涉条纹匹配方法 | |
| CN110988876B (zh) | 闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质 | |
| CN108562900A (zh) | 一种基于高程校正的sar图像几何配准方法 | |
| CN113238228B (zh) | 基于水准约束的三维地表形变获取方法、系统及装置 | |
| CN109166084B (zh) | 一种基于邻近点梯度关系的sar几何畸变定量模拟方法 | |
| CN107544069B (zh) | 基于平面近似模型的多基线相位解缠绕方法 | |
| CN115546264A (zh) | 一种星载InSAR图像精细配准与立体测量方法 | |
| CN115115678A (zh) | InSAR干涉测量复数图像高精度配准方法 | |
| CN109035312A (zh) | 一种dem辅助的sar图像高精度配准方法 | |
| CN115979207A (zh) | 一种围垦工程沉降监测方法、系统、设备和存储介质 | |
| CN115113202A (zh) | 一种基于二维高斯模型的干涉相位迭代解缠方法 | |
| CN105116410A (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 |