本发明属于信号处理技术领域，特别涉及一种地形高程自适应的降维图像配准方法，适用于干涉合成孔径雷达insar图像的配准。
背景技术：
干涉合成孔径雷达insar是一般sar功能的延伸和发展，它利用多个接收天线或单个天线多次观测得到的回波数据进行干涉处理；由于insar能够增强sar获取信息的能力，在军事、国民经济的各个领域有广阔的应用前景，因而获得了广泛重视。在星载insar高程测量中难以采用单航过方式，原因是航迹与场景的距离比较远，所以在星载insar干涉高程测量中目前广泛采用双航过的工作方式，即雷达卫星两次以基本相同的方向，不同的航迹通过同一地区的上空，将两次成像过程中雷达的回波数据和系统参数记录下来进行干涉处理；合成孔径雷达干涉测量地形高程时需要对同一场景获得两幅相干性很高的sar图像，因此需要进行图像配准的步骤，使配准后的图像中的同一位置的像素对应地面同一小块区域，以保证两幅图像的相干性。配准的精度直接影响到两幅图像的相干性。
以一个航过成像(称为主图像)的各像素为基准，将其像素对应的地面样点在另一个航过图像(称为辅图像)中的位置找出来，并将这一位置用偏移量表示，再根据偏移量，对辅图像的像素做位移，使其与主图像中同样位置的像素能对应于地面的同一点。配准处理分为粗配准和精配准；粗配准对辅图像做整体位移，不考虑局部性，精配准要考虑局部性，且要对辅图像进行插值，以达到高配准精度。
图像的配准首先要确定一幅图像中的每一个像素在另外一幅图像中的位置，即配准的二维偏移量；在基线较短的情况下，对于一块比较小的平坦的成像场景，只要将一副像进行简单的二维平移就可以将两幅图像对应像素配准；而对于基线较长的或地面起伏剧烈的情况下，不同方位、距离的像素进行配准需要不同的平移量，这使得图像配准过程的复杂度大大增加，不能用简单的图像平移和插值来实现。
已有的insar图像配准方法主要可分为两大类方法，第一类方法首先利用基于图像数据或图像特征的方法估计控制点处配准偏移量，然后根据二维二元低阶多项式模型计算其他像素处的配准偏移量，该第一类方法没有考虑地形高程对配准偏移量的影响，当基线较长或观测场景内地形复杂度较高时，利用二维二元低阶多项式模型拟合配准偏移量的精度较低，利用控制点处配准偏移量估计多项式系数时，控制点数量和位置分布对系数估计精度的影响较大；第二类方法首先利用几何配准法计算每一像素的初始配准偏移量，然后利用控制点处高精度配准偏移量根据二维二元一阶多项式模型校正初始配准偏移量，从而获取全场景高精度配准偏移量；该第二类方法需要进行逐像素几何配准处理，运算效率较低。
技术实现要素：
针对上述现有技术存在的不足，本发明的目的在于提出一种地形高程自适应的降维图像配准方法，该种地形高程自适应的降维图像配准方法能够追踪主辅图像间配准偏移量随地形高程的变化，是一种提高配准偏移量拟合精度的地形高程自适应的降维图像配准方法。
为达到上述技术目的，本发明采如下技术方案予以实现。
一种地形高程自适应的降维图像配准方法，包括以下步骤：
步骤1，确定主图像和辅图像，然后确定主辅图像间的配准偏移量；
步骤2，根据所述主辅图像间的配准偏移量，建立配准偏移量拟合模型，得到配准偏移量拟合模型中的8个待确定拟合系数；
步骤3，求解配准偏移量拟合模型中的8个待确定拟合系数，进而确定最终配准偏移量拟合模型；
步骤4，根据最终配准偏移量拟合模型，得到辅图像中每个像素点相对于主图像的配准偏移量，然后对辅图像进行重采样，进而得到精确配准的图像。
本发明与现有的技术相比具有以下优点：
第一，本发明中考虑了地形高程对配准偏移量的影响，提高了配准偏移量的拟合精度；在实际运用中，星载双航过的方式中会有复杂的情况出现，不仅有效基线长度较长，而且两次飞行的航迹不完全平行，具有一个很小的夹角；因此，长基线、地面起伏剧烈和航迹不平行都是图像配准中经常出现的非理想因素；本发明中提出的包含地形高程项的二维三元一阶多项式模型可以有效的解决地形高程对主辅图像配准偏移量的影响，从而提高配准偏移量的拟合精度。
第二，本发明中无须进行逐个像素点的几何配准处理，提高了运算效率：传统的几何配准法需要计算每一像素的初始配准偏移量，然后利用控制点处高精度配准偏移量根据二维二元一阶多项式模型校正初始配准偏移量，从而获取全场景高精度配准偏移量；而在本发明中提出的二维三元一阶多项式模型，只需要通过某一个像素点的配准偏移量计算就可以确定模型中的未知系数，其他像素处可以通过确定好的模型直接进行配准偏移量的步骤，减轻了运算的繁琐。
第三，本发明中通过降维处理将需要利用控制点处配准偏移量估计的多项式系数减少为两个常数项系数，提高了系数估计的稳健性；以往在利用基于图像数据或图像特征的方法估计控制点处配准偏移量，然后根据二维二元低阶多项式模型计算其他像素处的配准偏移量的方法中，利用控制点处配准偏移量估计多项式系数时，控制点数量和位置分布对系数估计精度的影响较大；在本发明中提出的二维三元一阶多项式模型，只需要用以往的方法确定两个常数项系数，其他系数都已经通过计算得出，所以控制点选取带来的影响对本发明中提出的模型较小，提高了系数估计的稳健性。
附图说明
下面结合附图和具体实施方式对本发明作进一步详细说明。
图1为本发明的一种地形高程自适应的降维图像配准方法实现流程图；
图2(a)为sar幅度图；
图2(b)为在图2(a)方框所示区域内选择12个特显点以定量评估配准精度示意图；
图2(c)为观测场景对应的srtmdem示意图；
图3为采用本发明方法配准后主辅图像间的干涉相位图；
图4为二阶多项式方法和本发明方法配准误差统计结果图；
图5(a)为二阶多项式方法的方位向配准误差示意图；
图5(b)为二阶多项式方法的距离向配准误差示意图；
图5(c)为本发明方法的方位向配准误差示意图；
图5(d)为本发明方法的距离向配准误差示意图。
具体实施方式
参照图1，为本发明的一种地形高程自适应的降维图像配准方法实现流程图；其中所述地形高程自适应的降维图像配准方法，包括以下步骤：
步骤1，确定星载雷达，所述星载雷达的检测区域内存在目标，且所述星载雷达对其检测区域进行两次观测后，得到两幅观测图像；两次观测下得到的两幅观测图像由于时间间隔的影响在距离向和方位向上会有像素偏差。
在两幅观测图像中任意选取其中一幅观测图像作为主图像，另外一幅观测图像作为辅图像，主图像和辅图像分别包含h个像素点；h为大于0的正整数。
选取主图像中任意一个像素点，记为主图像中像素点，所述主图像中像素点坐标为(am,rm)，由式(1)-(3)计算wgs84坐标系下目标在地表上的三维位置t，t＝(tx,ty,tz)t：
r＝|p-t|(1)
其中，r表示星载雷达到目标的斜距，t表示方位时间，p表示wgs84坐标系下星载雷达的位置矢量，p＝(px,py,pz)t，px表示wgs84坐标系下星载雷达在x轴的位置分量，py表示wgs84坐标系下星载雷达在y轴的位置分量，pz表示wgs84坐标系下星载雷达在z轴的位置分量，t表示wgs84坐标系下目标在地表上的三维位置，t＝(tx,ty,tz)t，tx表示wgs84坐标系下目标在地表上x轴的分量，ty表示wgs84坐标系下目标在地表上y轴的分量，tz表示wgs84坐标系下目标在地表上z轴的分量，fd表示多普勒频率，v表示wgs84坐标系下星载雷达的速度矢量，v＝(vx,vy,vz)t，vx表示wgs84坐标系下星载雷达在x轴的速度分量，vy表示wgs84坐标系下星载雷达在y轴的速度分量，vz表示wgs84坐标系下星载雷达在z轴的速度分量，(·)t表示转置，λ表示星载雷达发射信号的波长；h表示目标高程，可以通过dem数字高程模型得到；re表示地球赤道半径，rp表示地球极地半径；wgs84坐标系是一种国际上采用的地心坐标系。
然后根据wgs84坐标系下目标在地表上的三维位置t，计算得到与主图像中像素点坐标(am,rm)对应的辅图像中像素点坐标(as,rs)，其表达式为：
as＝(t-t0)·fp(4)
其中，t0表示辅图像中方位向上第一个像素点对应的方位时间，fp表示脉冲重复频率，r0表示辅图像中距离向第一个像素点对应的星载雷达斜距，fs′表示脉冲采样频率，c表示光速，r表示星载雷达到目标的斜距，t表示方位时间。
主图像中像素点坐标(am,rm)与所述与主图像中像素点坐标(am,rm)对应的辅图像中像素点坐标(as,rs)的差值，为主辅图像间的配准偏移量。
具体地，步骤1中所述的配准偏移量计算模型，就是需要确定某一像素点在主辅图像上各自对应的坐标，具体实施如下：
(1a)通过正向定位获取主图像像素对应的目标三维位置，正向定位的意思是将成像参数(如多普勒频率等)和轨道数据(卫星位置等)代入到距离-多普勒算法模型，求解地球上场景目标的三维位置，即从图像像素点入手求解目标三维位置的过程称为正向定位；
(1b)通过反向定位获取上述目标在辅图像上的坐标，反向定位的意思是利用目标三维位置、成像参数和轨道数据计算得到该目标点对应在图像上的像素点二维坐标。
步骤2，根据所述主辅图像间的配准偏移量，确定配准偏移量拟合模型：在主辅图像间基线较长或观测场景内地形复杂度较高时，需要考虑地形高程对配准偏移量的影响，式(6)和式(7)用包含地形高程项的二维三元一阶多项式模型拟合主辅图像间的坐标关系，进而得到配准偏移量拟合模型，并得到配准偏移量拟合模型的8个待确定拟合系数，分别为第一待确定拟合系数d0、第二待确定拟合系数d1、第三待确定拟合系数d2、第四待确定拟合系数d3、第五待确定拟合系数g0、第六待确定拟合系数g1、第七待确定拟合系数g2和第八待确定拟合系数g3；所述配准偏移量拟合模型为：
as＝d0+d1am+d2rm+d3h(6)
rs＝g0+g1am+g2rm+g3h(7)
具体地，传统配准方法假设目标高程对配准偏移量的影响可以忽略，用含有方位和距离坐标项的二维二元低阶多项式模型拟合配准偏移量随主图像方位和距离坐标的变化，但是在本方法中，考虑地形高程对配准偏移量的影响，用包含地形高程项h的二维三元一阶多项式模型拟合主辅图像间的坐标关系。
步骤3，估计配准偏移量拟合模型中的8个待确定拟合系数：需要确定配准偏移量拟合模型中的8个待确定拟合系数，利用主辅图像的成像参数和轨道数据计算6个一次项系数，分别为第二待确定拟合系数d1、第三待确定拟合系数d2、第四待确定拟合系数d3、第六待确定拟合系数g1、第七待确定拟合系数g2和第八待确定拟合系数g3。
具体地，首先推导目标三维位置相对于主图像方位时间、雷达斜距和目标高程的一阶偏导数，将式(1)～式(3)中目标三维位置分别对方位时间、雷达斜距和目标高程求一阶偏导数得到：
其中，
再推导目标对应辅图像方位时间和雷达斜距相对于目标三维位置的一阶偏导数，将式(1)和(2)中方位时间和雷达斜距分别对目标三维位置求一阶偏导数得到：
其中，下标s表示辅图像，as表示辅图像对应的星载雷达加速度矢量，ts表示辅图像对应的方位时间，ps表示辅图像对应的星载雷达位置矢量，vs表示辅图像对应的星载雷达速度矢量，rs表示确定辅图像时星载雷达到目标的斜距，t表示wgs84坐标系下目标在地表上的三维位置，t＝(tx,ty,tz)t，fd,s表示确定辅图像时星载雷达的多普勒频率。
结合上述推导结果得到6个一次项系数估计值，分别为第二拟合系数估计值
其中，第二拟合系数估计值
之后在主图像内任意选取l个像素点作为控制点，0<l<h，进而得到l个控制点，然后采用基于图像数据或图像特征的方法，得到l个控制点分别相对于主图像的配准偏移量，主图像像素点坐标(am,rm)已知，使用主图像像素点坐标(am,rm)分别与l个控制点的配准偏移量相加，进而得到l个控制点各自对应的辅图像像素点坐标，将l个控制点各自对应的辅图像像素点坐标和主图像像素点坐标(am,rm)代入式(6)和式(7)，进而得到l组系数估计值，每组系数估计值包括第一待确定拟合系数估计值和第五待确定拟合系数估计值，所述l组系数估计值数值近似；选取l组系数估计值中每l个数据的中间值，分别对应作为第一拟合系数估计值
具体地，步骤3中估计包含地形高程项的二维三元一阶多项式模型中的未知系数，包括如下步骤：
(3a)根据正向定位模型，推导目标三维位置相对于主图像方位时间、雷达斜距和目标高程的一阶偏导数。
(3b)根据反向定位模型，推导目标对应辅图像方位时间和雷达斜距相对于目标三维位置的一阶偏导数。
(3c)结合上述两步结果得到目标在辅图像上的二维坐标相对于其在主图像上的二维坐标和高程的一阶偏导数，即可得到第二拟合系数估计值
(3d)在主图像上选取控制点，采用基于图像数据或图像特征的方法估计控制点处配准偏移量；
(3e)补偿控制点处配准偏移量的一次分量，利用剩余配准偏移量得到第一拟合系数估计值
(3f)计算其他像素处配准偏移量，得到最终配准偏移量拟合模型，从而让主辅图像像素对应于地面同一成像单元。
所述最终配准偏移量拟合模型为：
为了降低控制点对系数估计的影响，提高系数估计的稳健性，本发明给出了利用主辅图像的成像参数和轨道数据计算一次项系数的公式；在实测数据处理过程中，利用观测场景中心处目标对应的主辅图像的成像参数和轨道数据计算d1～d3和g1～g3。
步骤4，根据步骤3得到的最终配准偏移量拟合模型确定主辅图像剩余h-l个像素点的相对偏移量：将主图像中剩余h-l个像素点所在区域均匀划分为m个控制块，将每个控制块的中心像素点坐标代入最终配准偏移量拟合模型中，得到对应像素点在辅图像中的坐标，主辅图像的坐标差值即为对应相对偏移量；进而得到m个相对偏移量；m为大于0的正整数。
然后对m个相对偏移量进行三次样条插值(通常采用10倍插值)，得到辅图像剩余h-l个像素点中每个像素点相对于主图像的配准偏移量；最后根据辅图像剩余h-l个像素点中每个像素点相对于主图像的配准偏移量，以及l个控制点分别相对于主图像的配准偏移量，对辅图像进行重采样，重采样后的结果为精确配准的图像，至此完成了主辅图像间的配准。
本发明能够提高配准偏移量拟合精度，提高了估计的稳健性，完成了主辅图像间的配准，可用于星载insar的图像配准处理。
下面给出一个实际运用中的例子，对本发明进一步说明：
利用terrasar-x分别于2010年3月8日和2010年3月19日对兰州地区进行重复航过观测任务时获取的两幅sar复图像验证本说明方法，主要系统参数如表1所示。取2010年3月8日获取的sar复图像为主图像，2010年3月19日获取的sar复图像为辅图像。
sar幅度图如图2(a)所示，在图2(a)方框所示的区域内选择12个特显点以定量评估配准精度，如图2(b)所示，图2(c)为观测场景对应的srtm(shuttleradartopographymission)dem雷达地形测绘的数字高程，已将其投影到主图像斜距平面，观测场景幅宽为15km×15km，观测场景内目标高程变化695m；利用图1所示的处理流程对主辅图像进行配准处理，通过最后处理得出的干涉相位图来验证本发明方法的有效性。
表1terrasar-x主要系统参数
表1中列出了terrasar-x雷达卫星的主要参数，通过雷达载波频率可以得到雷达的波长0.031m，脉冲重复频率即方位向参数，脉冲采样频率即距离向参数，方位向参数和距离向参数就是实际中的成像参数。
具体主辅图像配准实施步骤：
(1)针对主图像上某一像素(am,rm)，将多普勒频率、方位向参数、距离向参数和terrasar-x雷达卫星的实时位置及速度代入式(1)～式(3)所示的距离-多普勒模型，求解目标三维位置，再计算目标三维位置对应在辅图像上的二维坐标(as,rs)。
(2)列出包含地形高程项的二维三元一阶多项式模型利用观测场景中心处目标对应的主辅图像的成像参数和terrasar-x雷达卫星轨道数据根据式(15)～式(20)计算d1～d3和g1～g3。
(3)在主图像内均匀选择100个控制点，控制点处配准偏移量采用实相关函数法估计，再利用剩余配准偏移量估计d0和g0。
(4)用拟合模型计算两幅图像中的任何一个场景目标位置的相对偏移量，根据相对偏移量对整幅辅图像进行插值重采样(采用10倍插值)，得到精确配准的图像，至此完成了主辅图像间的配准。
配准后主辅图像间的干涉相位图如图3所示；利用本发明方法配准后的主辅图像可以获得清晰的干涉相位图，验证了本发明方法的有效性。
为了验证本发明方法的精确性和稳健性，利用mento-carlo实验对比分析本文方法与二阶多项式方法(利用二维二元二阶多项式模型拟合配准偏移量，并且多项式系数均利用控制点处配准偏移量估计)在不同控制点数量下的配准性能。
具体比较过程如下：
(1)mento-carlo实验次数为100次，在主图像内均匀选择900个控制点，控制点处配准偏移量采用实相关函数法估计；
(2)每次实验从上述900个控制点中随机选择100个控制点估计多项式拟合系数；
(3)配准误差利用图2(b)所示的12个特显点进行估计。对配准后主辅图像上的同名特显点进行200倍插值，将峰值间的距离向和方位向偏移量作为该特显点处的配准误差。图4给出了二阶多项式方法和本发明方法配准误差均方根值随控制点数量的变化。
由图4清晰可知，二阶多项式方法配准误差随控制点数量变化较大，而本说明方法配准误差随控制点数量变化较小；由表1可知，主辅图像间的垂直基线大于水平基线，导致目标高程对距离向配准偏移量的影响大于方位向配准偏移量，所以二阶多项式方法的距离向配准误差大于方位向配准误差；本发明方法可以追踪配准偏移量随目标高程的变化，所以本发明方法的方位向和距离向配准误差近似相同，并且可以满足0.1像素的配准精度要求；上述结果验证了本发明方法的精确性和稳健性。
为了更直观证明本发明方法相对于二阶多项式方法的精确性和稳健性，图5(a)、图5(b)、图5(c)和图5(d)给出了控制点数量为10个时，100次试验的12个特显点处的配准误差分布；图5(a)、图5(b)、图5(c)和图5(d)分别为二阶多项式方法的方位向和距离向配准误差、本发明方法的方位向和距离向配准误差。
由图5(a)、图5(b)、图5(c)和图5(d)可知，二阶多项式方法的方位向和距离向配准误差最大值分别为3.47像素和3.45像素，而本发明方法的方位向和距离向配准误差最大值分别为0.17像素和0.30像素。
总而言之，本发明方法公开了一种地形高程自适应的降维图像配准方法，解决了长基线或复杂场景星载insar图像的配准问题；其实现步骤是：
构建配准偏移量拟合模型，并利用观测场景中心处目标对应的主辅图像的成像参数和轨道数据计算其中6个未知系数；在主图像上选取一定数量的控制点，采用基于图像数据或图像特征的方法估计控制点处配准偏移量；估计剩余的2个常数项拟合系数；计算其他像素处配准偏移量，以此偏移量对所有辅图像重采样，完成主辅图像配准；本发明方法不仅能够精确拟合主辅图像间的配准偏移量，而且通过降维处理将需要利用控制点处配准偏移量估计的多项式系数减少为两个常数项系数，提高了系数估计的稳健性。
综上所述，仿真实验验证了本发明的正确性，有效性和可靠性。
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围；这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。