时间序列干涉SAR的形变量测量方法及SAR系统与流程

时间序列干涉sar的形变量测量方法及sar系统
技术领域
1.本发明涉及一种合成孔径雷达(sar)检测技术，更具体地涉及一种时间序列干涉sar的形变量测量方法及采用其的sar系统。


背景技术：

2.sar(synthetic aperture radar)，即合成孔径雷达，是一种主动式的对地观测系统，可安装在飞机、卫星、宇宙飞船等飞行平台上，全天时、全天候对地实施观测、并具有一定的地表穿透能力。因此，sar系统在灾害监测、环境监测、海洋监测、资源勘查、农作物估产、测绘和军事等方面的应用上具有独特的优势，可发挥其他遥感手段难以发挥的作用，因此越来越受到世界各国的重视。
3.而卫星雷达时序干涉技术(ts
‑
insar)则是一种利用卫星雷达多时相大数据干涉相干原理，实现800公里太空探测地面毫米形变的遥感技术，是目前唯一的一种高精度、全天时、全天候、全覆盖、全自动的面监测技术手段，一次探测范围平均可达1500平方公里以上，可为国土、交通、测绘、铁道、电力、水利、市政、减灾、住建、能源等行业提供精准形变测量、状态监控、异常预警、灾后定损等服务，高效解决地理信息监测市场痛点。
4.该形变监控测量技术具备高精度、全天时、全天候、全覆盖、全自动等特性，但由于卫星本身的运动及太空俯拍视角，如何使测量的形变更加精确就是现在迫切需要解决的技术问题。


技术实现要素：

5.有鉴于此，本发明的主要目的在于提出一种时间序列干涉sar的形变量测量方法及采用其的sar系统，以期至少部分地解决上述技术问题。
6.为了实现上述目的，本发明第一个方面提供了一种时间序列干涉sar 的形变量测量方法，包括以下步骤：
7.步骤1：对sar影像集中的所有影像进行数据预处理，从而选取出主影像和辅影像，并将所有主、辅影像前置滤波和计算干涉相位，生成干涉图；
8.步骤2：从步骤1得到的干涉图的干涉相位中去除平地和地形相位，生成差分干涉相位，逐像元计算生成差分干涉图；
9.步骤3：对步骤2得到的差分干涉图的差分干涉相位进行时间和空间域的线性变形相位估计，得到每个点目标的时间序列变形相位；
10.步骤4：依据雷达波长参数，进行相位转形变的计算，由此得到所述sar 影像的形变量测量值。
11.本发明第二个方面还提供了一种合成孔径雷达系统，包括发射机和接收机，以及处理器和存储器，所述存储器用于存储可执行程序，其中当所述可执行程序被所述处理器执行时，所述处理器执行如上所述的形变量测量方法。
12.本发明第三个方面还提供了一种存储介质，存储有可执行程序，其中当所述可执
行程序被执行时，实现如上所述的形变量测量方法。
13.基于上述技术方案可知，本发明的时间序列干涉sar的形变量测量方法及采用其的sar系统相对于现有技术至少具有如下有益效果之一：
14.本发明的方法形变测量精度高，计算步骤简洁高效，且可以克服微小的地表形变测量由于时间跨度的增加和两幅图像多普勒质心频率差异变大所导致的时间去相关的难题，计算结果稳定；
15.本发明可以排除噪声和阴影等因素的影响，匹配结果中剔除了误匹配现象，保障了影像配准的准确性；
16.本发明通过影像偏移量粗估计和精估计，能够极大地提高影像配准的精确度和形变计算的准确率；
17.本发明不仅对干涉图的差分干涉相位进行时间和空间域的线性变形相位估计，还在需要时进行非线性变形相位估计，可以去除大气、噪声等残余相位，得到每个点目标的时间序列变形相位；
18.本发明可以利用已有的水准、gps等高精度控制点数据(同期观测的变形量)来修正基准，保证计算的准确性。
附图说明
19.图1是本发明的时间序列干涉sar的形变量测量方法的方框流程图；
20.图2是本发明实施例1的ts
‑
insar技术的数据处理的工艺流程图；
21.图3是本发明多时相影像精配准模块轨道影像初始偏移估计示意图；
22.图4是本发明的影像偏移粗估计步骤中的锚点分布示意图；
23.图5是本发明的影像偏移粗估计的主辅影像选取示意图；
24.图6是本发明的影像偏移粗估计的算法流程图；
25.图7是本发明的采用均值方差法剔除粗差点的流程示意图。
具体实施方式
26.在对于具体实施例的介绍过程中，对结构、性能、效果或者其他特征的细节描述是为了使本领域的技术人员对实施例能够充分理解。但是，并不排除本领域技术人员可以在特定情况下，以不含有上述结构、性能、效果或者其他特征的技术方案来实施本发明。
27.附图中的流程图仅是一种示例性的流程演示，不代表本发明的方案中必须包括流程图中的所有的内容、操作和步骤，也不代表必须按照图中所显示的的顺序执行。例如，流程图中有的操作/步骤可以分解，有的操作/ 步骤可以合并或部分合并，等等，在不脱离本发明的发明主旨的情况下，流程图中显示的执行顺序可以根据实际情况改变。
28.附图中的框图一般表示的是功能实体，并不一定必然与物理上独立的实体相对应。即，可以采用软件形式来实现这些功能实体，或在一个或多个硬件模块或集成电路中实现这些功能实体，或在不同网络和/或处理单元装置和/或微控制器装置中实现这些功能实体。
29.ts
‑
insar技术是一种对长时间序列sar影像集中的永久散射体进行时间和空间域变形量计算，以提取高精度时序变形信息的干涉测量方法。其数据处理工艺流程包括差分
干涉计算、时间及空间形变量估算等处理流程，能够比较精确地还原卫星测量形变。
30.具体地，本发明公开了一种时间序列干涉sar的形变量测量方法，包括以下步骤：
31.步骤1：对sar影像集中的所有影像进行数据预处理，从而选取出主影像和辅影像，并将所有主、辅影像前置滤波和计算干涉相位，生成干涉图；
32.步骤2：从步骤1得到的干涉图的干涉相位中去除平地和地形相位，生成差分干涉相位，逐像元计算生成差分干涉图；
33.步骤3：对步骤2得到的差分干涉图的差分干涉相位进行时间和空间域的线性变形相位估计，得到每个点目标的时间序列变形相位；
34.步骤4：依据雷达波长参数，进行相位转形变的计算，由此得到所述sar 影像的形变量测量值。
35.其中，步骤1的选取出主影像和辅影像的步骤中，所述主影像的数量为1。
36.其中，所述选择主影像的步骤具体包括：
37.a.计算所有影像像对间的时间和空间基线，生成时间和空间基线分布图；
38.b.选择时间和空间基线居中的一幅影像作为主影像。
39.其中，步骤1的对sar影像集中的所有影像进行数据预处理的步骤例如包括：
40.将所有sar影像相对主影像进行配准、裁剪，并组合生成时间序列干涉图集。
41.其中，将所有sar影像相对主影像进行配准、裁剪的具体步骤例如如下：
42.a.将所有影像相对主影像进行配准，要求已组合好的像对，根据主影像进行配准，并将所有影像裁剪成范围一致的区域；
43.b.将所有数据裁剪成范围一致的区域；
44.c.对所有已配准的干涉像对，按照时间序列分别与主影像进行像对组合，逐像元计算干涉相位，生成时间序列干涉图集。
45.其中，步骤a例如具体要求符合如下规定：
46.选择配准算法，设置配准参数，对所有影像中的每个像对进行配准计算；
47.主辅影像配准时要求方位向和距离向误差均小于0.25个像元，且计算配准多项式的同名点应在整景影像上均匀分布。
48.其中，步骤b具体要求例如符合如下规定：
49.所有配准影像裁剪后的公共区域应大于或等于设计的监测工作范围，如有缺失应及时补充数据；
50.选择配准影像中的公共区域作为insar处理范围，将所有影像裁剪成相同范围的区域。
51.其中，步骤c中主影像例如通过如下步骤选取：
52.对同一地区获取的n+1幅sar图像，依据sar图像的时间基线、空间基线、多普勒中心频率构成的三维空间分布图来选择一幅合适的sar图像作为主影像，将其它图像作为辅影像，组合成n幅干涉图；其中，在充分考虑n+1幅sar图像之间的时间基线、空间基线和多普勒质心频率差三个因素最优化的基础上，建立如下的综合相关系数模型：
[0053][0054]
其中，
[0055][0056]
上述公式中，γ
m
为综合相关系数，k为一常数，n为配准后的sar干涉图数量，和b
⊥
，c
分别表示主影像、从影像对应的空间垂直基线和时间基线，t
k，m
和t
c
分别表示主影像、从影像对应的相对参考时刻的时间累积量，f
k，m
和f
c
分别表示主影像、从影像对应的多普勒质心频率，x为方位向位置，c为光速；当综合相关系数达到最大时，所对应的图像即为选取的公共区域的主影像，简称公共主影像。
[0057]
其中，步骤a、b中将所有sar影像相对主影像进行配准、裁剪的步骤例如具体包括：
[0058]
基于slc参数文件中的轨道信息计算辅影像相对于主影像的初始偏移量；
[0059]
对影像偏移进行粗估计，具体包括：使用主辅slc数据的互相关优化粗略估算距离和方位配准偏移量，并根据偏移量文件粗略估计配准偏移多项式；
[0060]
对影像偏移进行精估计，具体包括：使用主辅slc数据的互相关优化精细估算距离和方位配准偏移量，并根据偏移量文件精细估计配准偏移多项式；
[0061]
基于偏移参数文件中的距离和方位向的偏移多项式，使用二维sinc插值进行slc重采样；
[0062]
将dem与主影像进行配准，并将dem范围裁剪成与主影像一致。
[0063]
其中，所述主辅slc数据的互相关优化粗略估算优选是像素级的。
[0064]
其中，所述影像偏移进行粗估计的步骤例如是通过采用锚点匹配方法来实现。
[0065]
其中，所述影像偏移进行粗估计的步骤例如具体包括：
[0066]
采用基于窗口的搜索方法，即在主影像上截取n
×
n大小的单视复影像数据作为基准窗口，在辅窗口上截取m
×
m大小的单视复影像数据作为搜索窗口，然后在辅影像上移动搜索窗口，以预定的准则来确定主影像上按照一定间隔均匀设置的锚点对应的同名点；其中，影像偏移粗估计单元将影像同名点匹配到1个像素。
[0067]
其中，所述预定的准则例如为采用相关系数法和频谱极大值法进行偏移量粗估计。
[0068]
其中，所述频谱极大值法具体包括：以寻求两幅复数图像的频谱极大值作为配准测度，当两幅复图像精确配准时，形成的干涉条纹质量最高。该步骤的具体衡量标准例如要求满足如下关系式：
[0069]
f＝fft(r
·
s
*
)；
[0070]
其中，f为相对最大频谱，fft()为快速傅立叶变换(复共轭)，r和 s分别为主、辅图像对应像素的复数值。
[0071]
通常，频谱极大值法的基本思路是定义复数干涉图谱的最大模值与其它频率成分的模值之和的比值，即使相对最大频谱为配准过程所使用的评价函数。其具体操作一般是在主图像中取一待匹配点，在其周围以较大的窗口(搜索窗口，m
×
m)截取子图像，记作s1，在辅图像中同一位置取对应点，在其周围取较小的窗口(匹配窗口，n
×
n)截取子图像，记作 s2，其中m>2
×
n，以保证计算的准确性。将s2在s1内按二维方向逐点滑动，并且每滑动一次，计算一次该位置的干涉图频谱，即f＝fft(s1’,s2*)。其中，s1’表示s1与s2同等大小的子图像。由此，最大频谱法实际上是利用了insar影像的相位信息进行配准，即当两幅影像相位
越接近时，对应的频谱值越大，这也被称之为相位相关。
[0072]
其中，所述影像偏移精估计时所述精细估算步骤优选为亚像素级。
[0073]
其中，所述影像偏移进行精估计的步骤例如具体包括：
[0074]
影像偏移精估计在粗估计求取的辅影像相对于主影像偏移量的基础上，配准到子像素级，要求到1/10个像素，本单元采用过采样法和相关系数fft插值方法；
[0075]
将主影像的基准窗口和辅影像搜索窗口的数据内插k倍，然后以此为基准采用相关系数法或频谱极大值法或最小平均波动函数法求解偏移量，以此匹配到亚像素级别。
[0076]
其中，对于误匹配点，例如采用均值方差法或多项式拟合法来剔除。
[0077]
其中，对于影像插值重采样步骤，其是通过锚点的偏移量来拟合全图各个像素的偏移量，再利用插值辅图像来得到新的辅图像，根据偏移量对辅影像重采样。
[0078]
其中，影像插值重采样步骤例如可以采用双线性插值法或四点三次卷积插值法来实现。
[0079]
其中，将dem与主影像配准和裁剪的步骤中，具体步骤例如如下：
[0080]
a.对dem采样成与主影像一致的分辨率；
[0081]
b.将dem与主影像进行配准，配准精度优于0.5个像元；
[0082]
c.依据配准关系式，计算生成dem坐标系到sar影像坐标系的转换查找表；
[0083]
d.依据转换查找表，利用多项式拟合算法，将dem转换到sar影像坐标系，生成影像坐标系下的dem。
[0084]
其中，步骤c中将所有主、辅影像前置滤波，计算干涉相位并生成干涉图的具体步骤例如如下：
[0085]
a.在频率域，截取主、辅影像的公共频带进行前置滤波，生成滤波后的主、辅影像；
[0086]
b.对已经过前置滤波的主辅影像像元对进行复共轭相乘，生成干涉相位值，逐像元计算生成干涉图。
[0087]
其中，步骤c中将所有主、辅影像前置滤波，计算干涉相位并生成干涉图的步骤之后还可以包括：
[0088]
对时间序列干涉图集的像元进行cs点目标筛选，具体步骤如下：
[0089]
a.采用幅度离差指数法或信噪比法进行sar数据的cs点目标的识别；
[0090]
b.依据ks理论中的smirnov检测方法识别ds点；
[0091]
c.将满足上述条件要求的点目标从干涉图集中提取出来，生成 cs/ds点目标的干涉相位序列。
[0092]
其中，步骤2中去除平地和地形相位的步骤例如具体包括：
[0093]
依据空间基线参数和地球椭球体参数，计算平地相位；
[0094]
利用配准后dem，计算地形相位；
[0095]
从干涉相位中去除平地和地形相位，生成差分干涉相位，逐像元计算生成差分干涉图。
[0096]
其中，步骤2中去除平地和地形相位的步骤之后还可以包括如下步骤：
[0097]
目视检查每景差分干涉图，若含有残余干涉条纹超过半个波长，计算空间基线残余相位，并去除。
[0098]
其中，所述计算空间基线残余相位并去除的具体步骤例如如下：
[0099]
a.利用二次曲面模型对差分干涉图进行空间基线粗估计，得到空间基线的粗估计相位；再利用差分干涉图中差分相位减去粗估计相位，得到残余相位；
[0100]
b.利用快速傅立叶变换对残余相位进行估计，得到残余基线相位；
[0101]
c.将空间基线粗估计相位加上残余基线相位，得到改进的空间基线相位；
[0102]
d.利用改正的空间基线相位，对平地相位去除残余平地相位，计算得到改正后的平地相位和干涉图集。
[0103]
其中，步骤2中去除平地和地形相位的步骤之后还可以包括如下步骤：
[0104]
对干涉图的差分干涉相位进行时间和空间域的线性变形相位估计，如有要求还进行非线性变形相位估计，去除大气或噪声残余相位，得到每个点目标的时间序列变形相位。
[0105]
其中，所述线性变形相位估计和非线性变形相位估计的具体计算步骤例如如下：
[0106]
①
将cs点目标相连接构成delauney不规则三角网(dtin，或称冗余网)，依据点间连接关系求解相邻点差分相位差；
[0107]
②
依据空间基线、时间基线关系，建立cs点目标的二维周期图，以此为目标函数使模型相关系数最大化，估算相邻点间的线性变形速率和高程差值；若监测工作设计书仅要求线性变形成果，则直接输出成果垂直向变形量计算，生成地质灾害体速率图；
[0108]
③
从差分干涉相位中去除步骤
②
中两项相位量，得到残余相位；对所述残余相位进行空间域均值滤波，计算得到主影像大气相位；对去除主影像大气相位的残余相位进行空间域低通滤波和时间域高通滤波，得到从影像大气相位，并进一步分解出非线性变形相位；
[0109]
④
将步骤
②
中线性变形相位和步骤
③
中非线性变形相位相加，结合时间基线参数，得到每个cs点目标的时间序列变形相位。
[0110]
其中，步骤4中依据雷达波长参数，经过相位转形变计算之后，还根据需要结合外部辅助数据，将解缠相位换算为视线方向的形变量，然后依据视线向与垂直向夹角，将视线向形变量转换为垂直向。
[0111]
其中，步骤4中经过相位转形变计算之后还包括地理编码的步骤，所述地理编码步骤是利用dem(digital elevation model，数字高程模型)产品进行地理编码，具体步骤如下：
[0112]
a.利用dem坐标系到sar影像坐标系的转换查找表，完成监测成果由sar影像坐标系到大地坐标系的反变换，即对监测成果垂直向形变量进行地理编码；
[0113]
b.集合所有地理编码后的点目标，将变形量的时间单位换算成年，生成年度变形速率，逐像元计算生成地质灾害体速率图。
[0114]
其中，步骤4中经过相位转形变计算之后还包括基准修正的步骤，地理编码后点目标的灾害体速率是利用已有的水准、gps高精度控制点数据来修正基准。
[0115]
其中，所述基准修正的具体步骤为：
[0116]
a.以同步水准测量结果作为基准参考，在临近点上计算点目标变形量与实测值之间差值的平均值，即与实测变形量之间存在的整体偏差值；
[0117]
b.将上一步得到的整体偏差值加入每个点目标的变形值，修正因参考点不统一产生的insar结果变形量的整体偏差，完成基准修正。
[0118]
本发明还公开了一种电子设备，例如一种合成孔径雷达系统，包括发射机和接收
机，以及处理器和存储器，所述存储器用于存储可执行程序，当所述可执行程序被所述处理器执行时，所述处理器执行如上所述的时间序列干涉sar的形变量测量方法。
[0119]
该电子设备(即合成孔径雷达系统)的处理器和存储器部分例如可以以逻辑控制部件或通用计算设备的形式表现，例如可编程逻辑控制器 (plc)、现场可编程门阵列(fpga)、单片机、单板机、台式电脑、笔记本电脑、服务器、计算机网络等形式。其中，处理器可以是一个，也可以是多个且协同工作。本发明也不排除进行分布式处理，即处理器可以分散在不同的实体设备中。本发明的电子设备并不限于单一实体，也可以是多个实体设备的总和。
[0120]
其中，存储器存储有计算机可执行程序，通常是机器可读的代码，该计算机可执行程序可以被所述处理器执行，以使得电子设备能够执行本发明的方法，或者方法中的至少部分步骤。
[0121]
该存储器包括易失性存储器，例如随机存取存储单元(ram)和/或高速缓存存储单元，还可以是非易失性存储器，如只读存储单元(rom)。
[0122]
可选的，该实施方式中，电子设备还包括有i/o接口，其用于电子设备与外部的设备进行数据交换。i/o接口可以为表示几类总线结构中的一种或多种，包括存储单元总线或者存储单元控制器、外围总线、图形加速端口、处理单元或者使用多种总线结构中的任意总线结构的局域总线。
[0123]
本发明的电子设备中还可以包括上述示例中未示出的元件或组件。例如，有些电子设备中还包括有显示屏等显示单元，有些电子设备还包括人机交互元件，例如按钮、键盘等。只要该电子设备能够执行存储器中的计算机可读程序以实现本发明方法或方法的至少部分步骤，均可认为是本发明所涵盖的电子设备。
[0124]
本发明还公开了一种存储介质，存储有可执行程序，当所述可执行程序被执行时，实现如上所述的时间序列干涉sar的形变量测量方法。
[0125]
该存储介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了可读程序代码。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。可读存储介质还可以是可读存储介质以外的任何可读介质，该可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。可读存储介质上包含的程序代码可以用任何适当的介质传输，包括但不限于无线、有线、光缆、rf等等，或者上述的任意合适的组合。
[0126]
可以以一种或多种程序设计语言的任意组合来编写用于执行本发明操作的程序代码，所述程序设计语言包括面向对象的程序设计语言—诸如 python、java、c++等，还包括常规的过程式程序设计语言—诸如c语言、汇编语言或类似的程序设计语言。程序代码可以完全地在用户计算设备上执行、部分地在用户设备上执行、作为一个独立的软件包执行、部分在用户计算设备上部分在远程计算设备上执行、或者完全在远程计算设备或服务器上执行。在涉及远程计算设备的情形中，远程计算设备可以通过任意种类的网络，包括局域网(lan)或广域网(wan)，连接到用户计算设备，或者，可以连接到外部计算设备(例如利用因特网服务提供商来通过因特网连接)。
[0127]
为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。
[0128]
如图2所示，本实施例的时间序列干涉sar的形变量测量方法，主要包含数据预处理、差分干涉计算、时间及空间形变量估算、形变量计算等处理流程。下文将展开分述之。
[0129]
(一)数据预处理
[0130]
⑴
主影像选择
[0131]
本实施例的ts
‑
insar方法选择单一主影像。在满足空间基线和时间基线要求的前提下，sar主影像的选择和像对组合工作步骤如下：
[0132]
a.计算所有影像像对间的时间和空间基线，生成时间和空间基线分布图；
[0133]
b.选择时间和空间基线居中的一景作为主影像。
[0134]
⑵
影像配准、裁剪和组合
[0135]
所有sar影像对主影像进行配准、裁剪，并组合生成时间序列干涉图集。具体步骤如下：
[0136]
a.所有影像对主影像进行配准。已组合好的像对，根据主影像进行配准，并将所有影像裁剪成范围一致的区域。具体应符合如下规定：
[0137]
选择配准算法，设置配准参数，对每个像对进行配准计算。
[0138]
主辅影像配准时要求方位向和距离向误差均小于0.25个像元，且计算配准多项式的同名点应在整景影像上均匀分布。
[0139]
b.将所有数据裁剪成范围一致的区域。裁剪要求应满足：
[0140]
所有配准影像裁剪后的公共区域应大于或等于设计的监测工作范围，如有缺失应及时补充数据。
[0141]
选择配准影像中的公共区域作为insar处理范围，将所有影像裁剪成相同范围的区域。
[0142]
c.对所有已配准的干涉像对，按照时间序列分别与主影像进行像对组合。逐像元计算干涉相位，生成时间序列干涉图集。配准理论如下：
[0143]
主影像选取
[0144]
对同一地区获取的n+1幅sar图像，依据sar图像的时间基线、空间基线、多普勒中心频率构成的三维空间分布图来选择一幅合适的sar图像作为主图像，将其他图像作为辅图像，组合成n幅干涉图。为了得到质量较高的干涉相位，在干涉处理之前，选取一幅时空分布最优的图像作为干涉配准的公共主图像。公共主图像的选取涉及到后续所形成干涉对的相关性，其时间基线、有效空间基线和多普勒质心频率差是影响干涉相关的重要因素。对于干涉形变测量而言，有效空间基线长度越小越好；而对于微小的地表形变测量来说，图像又必须具有较长的时间跨度，而时间跨度的增加极易导致时间去相关；此外两幅图像多普勒质心频率差异越大，去相关程度也越大。因此，在充分考虑n+1幅sar图像之间的时间基线、空间基线和多普勒质心频率差三个因素最优化的基础上，可建立如下的综合相关系数模型：
[0145][0146]
其中，
[0147][0148]
上述公式中，γ
m
为综合相关系数，k为一常数，n为配准后的sar干涉图数量，和b
⊥
，c
分别表示主影像、从影像对应的空间垂直基线和时间基线，t
k，m
和t
c
分别表示主影像、从影像对应的相对参考时刻的时间累积量，f
k，m
和f
c
分别表示主影像、从影像对应的多普勒质心频率，x为方位向位置，c为光速；当综合相关系数达到最大时，所对应的图像即为选取的公共区域的主影像，简称公共主影像。
[0149]
影像偏移初始化
[0150]
图3是本发明的多时相影像精配准模块轨道影像初始偏移估计的示意图，如图3所示，影像偏移初始化单元基于slc参数文件中的轨道信息来计算辅影像相对于主影像的初始偏移量。
[0151]
影像偏移粗估计
[0152]
功能描述：
[0153]
①
使用主辅slc数据的互相关优化粗略(像素级)估算距离和方位配准偏移量。
[0154]
②
根据偏移量文件粗略估计配准偏移多项式。
[0155]
详细设计：
[0156]
影像偏移粗估计单元的目的是将影像同名点匹配到1个像素，本模块拟采用锚点匹配方法，即在主影像上按照一定间隔均匀设置锚点，然后在辅影像上寻找同名点，如图4所示。
[0157]
本实施例中采用基于窗口的搜索方法，即在主影像上截取n
×
n大小的单视复影像数据(基准窗口)，在辅窗口上截取m
×
m大小的单视复影像数据(搜索窗口)，然后在辅影像上移动主窗口，以某一准则例如相关系数最大来确定复影像对应的同名点，如图5所示。
[0158]
本实施例中拟采用相关系数法和频谱极大值法进行影像偏移粗估计。
[0159]
频谱极大值测度以寻求两幅复数图像的频谱极大值作为配准测度，当两幅复图像精确配准时，形成的干涉条纹质量最高。具体计算公式如下：
[0160]
f＝fft(r
·
s
*
)；
[0161]
其中，f为相对最大频谱，fft()为快速傅立叶变换(复共轭)，r和 s分别为主、辅图像对应像素的复数值。
[0162]
通过相关系数最大的原则来进行估计的流程如图6所示。
[0163]
影像偏移精估计
[0164]
功能描述：
[0165]
①
使用主辅slc数据的互相关优化精细(亚像素级)估算距离和方位配准偏移量；
[0166]
②
根据偏移量文件精细估计配准偏移多项式。
[0167]
详细设计：
[0168]
影像偏移精估计在粗估计求取的辅影像相对于主影像偏移量的基础上，配准到子像素级，一般要求到1/10个像素，本单元采用过采样法和相关系数fft插值方法。
[0169]
将主影像的基准窗口和辅影像搜索窗口的数据内插k倍，然后以此为基准采用相关系数法或频谱极大值法或最小平均波动函数法求解偏移量，以此匹配到亚像素级别，当
采用频谱极大值法时，以保障基准窗口大小为 2的幂次方。由于粗配准已经到了一个像素，因此可考虑搜索窗口大小比基准窗口大小多4，此时搜索范围为[
‑
2,2]。
[0170]
假设粗配准时辅影像相对于主影像的偏移量，主影像锚点位置，通过过采样方法计算出的像点偏移量，则通过精配准后主影像锚点在辅影像上的坐标为：
[0171][0172]
过采样法配准的流程与粗配准流程基本类似，只是对原窗口数据进行内插后再利用三种粗配准的方法求解偏移量。
[0173]
受噪声和阴影等因素影响，匹配结果中某些点会存在误匹配现象，必须将这些误匹配点探测出来，保障影像配准的准确性。粗差剔除的方法有均值方差法和多项式拟合法。
[0174]
根据insar公式，
[0175]
则像素偏移量为
[0176]
从理论分析可知，主影像上同一列锚点(距离向像素坐标相同)，其在辅影像上的偏移量方差应该优于1个像素。因此可将偏离均值的程度作为剔除粗差点的依据。同时，整列锚点控制以方差δ小于某一个值l作为依据。
[0177][0178][0179]
均值方差法剔除粗差点的流程图如图7所示。
[0180]
影像插值重采样
[0181]
功能描述：
[0182]
影像插值模块是基于偏移参数文件中的距离和方位向的偏移多项式，使用二维sinc插值进行slc重采样。
[0183]
详细设计：
[0184]
影像插值重采样单元的基本思路就是通过锚点的偏移量来拟合全图各个像素的偏移量，再利用插值辅图像来得到新的辅图像。由于主图像中所有像元对应辅图像的偏移量不可能完全相同，且不可能对其中的每一个像元通过配准找到其偏移量，所以需要从主影像中有限的已配准得到偏移量的像元中，用二阶或多阶多项式拟合主影像锚点像元图像坐标与偏移量之间的关系，得到多项式系数。利用多项式系数，可以计算主影像中各个像元的偏移量，然后根据偏移量对辅影像重采样。
[0185]
在得到n个锚点的精确偏移量之后，采用二阶(或者高阶，此处以二阶为例)多项式拟合参数：
[0186][0187]
利用上式得到的拟合参数，对主影像中的每个像素的坐标计算出一个偏移量，根据此偏移量在辅图像中找到其匹配的位置。
[0188]
当两幅复影像之间的精确偏移量确定后，就需要对辅影像进行重采样 (插值)，使辅影像中像元与主影像中像元一一对齐。由于最近邻插值法不能满足insar干涉处理中子像元配准对插值精度的要求，所以本实施例采用双线性插值法和四点三次卷积插值法。
[0189]
⑶
dem(digital elevation model，数字高程模型)与主影像配准和裁剪
[0190]
将dem与主影像进行配准，并将dem范围裁剪成与主影像一致。具体步骤如下：
[0191]
a.应对dem采样成与主影像一致的分辨率；
[0192]
b.将dem与主影像进行配准，配准精度应优于0.5个像元；
[0193]
c.依据配准关系式，计算生成dem坐标系到sar影像坐标系的转换查找表；
[0194]
d.依据转换查找表，利用多项式拟合算法，将dem转换到sar影像坐标系，生成影像坐标系下的dem。
[0195]
⑷
干涉图相位计算
[0196]
将所有主、辅影像前置滤波，计算干涉相位，生成干涉图。具体步骤如下：
[0197]
a.前置滤波。在频率域，截取主、辅影像的公共频带进行前置滤波，生成滤波后的主、辅影像；
[0198]
b.干涉相位计算。对已经过前置滤波的主辅影像像元对进行复共轭相乘，生成干涉相位值，逐像元计算生成干涉图。
[0199]
⑸
cs点及ds点目标选取
[0200]
对时间序列干涉图集的像元进行cs点目标筛选。具体步骤如下：
[0201]
a.cs点目标识别。sar数据cs点目标的识别宜采用幅度离差指数法、信噪比法等方法。结合监测区地物类型，宜选择一种或多种方法，以提高 cs点目标识别的准确性；
[0202]
b.ds点目标识别。依据ks理论中的smirnov检测方法识别ds点；
[0203]
c.cs/ds点目标干涉相位序列生成。将满足上述条件要求的点目标从干涉图集中提取出来，生成cs/ds点目标的干涉相位序列。
[0204]
(二)差分干涉计算
[0205]
⑴
平地和空间相位去除
[0206]
依据空间基线参数和地球椭球体参数，计算平地相位；利用配准后 dem，计算地形相位。从干涉相位中去除平地和地形相位，生成差分干涉相位，逐像元计算生成差分干涉图。
[0207]
⑵
空间基线精化
[0208]
目视检查每景差分干涉图，若含有残余干涉条纹超过半个波长，计算空间基线残余相位，并去除。具体步骤如下：
[0209]
a.利用二次曲面模型对差分干涉图进行空间基线粗估计，得到空间基线的粗估计相位；再利用差分干涉图中差分相位减去粗估计相位，得到残余相位；
[0210]
b.利用快速傅立叶变换对残余相位进行估计，得到残余基线相位；
[0211]
c.将空间基线粗估计相位加上残余基线相位，得到改进的空间基线相位；
[0212]
d.利用改正的空间基线相位，对平地相位去除残余平地相位，计算得到改正后的平地相位和干涉图集。
[0213]
⑶
时间/空间域变形估算
[0214]
对干涉图的差分干涉相位应进行时间和空间域的线性变形相位估计，如有要求还应进行非线性变形相位估计，去除大气、噪声等残余相位，得到每个点目标的时间序列变形相位。ts
‑
insar的计算步骤如下：
[0215]
①
相邻点间参数估计
[0216]
将cs点目标相连接构成delauney不规则三角网(dtin，或称冗余网)，依据点间连接关系求解相邻点差分相位差。
[0217]
②
线性变形相位和残余高程相位计算
[0218]
依据空间基线、时间基线关系，建立cs点目标的二维周期图，以此为目标函数使模型相关系数最大化，估算相邻点间的线性变形速率和高程差值。若监测工作设计书仅要求线性变形成果，则可直接输出成果垂直向变形量计算，生成地质灾害体速率图。
[0219]
③
非线性变形相位和大气相位计算
[0220]
从差分干涉相位中去除
②
中两项相位量，得到残余相位。对该残余相位进行空间域均值滤波，计算得到主影像大气相位。对去除主影像大气相位的残余相位进行空间域低通滤波和时间域高通滤波，得到从影像大气相位，并进一步分解出非线性变形相位。
[0221]
④
时间序列变形相位计算
[0222]
将
②
步骤中线性变形相位和
③
中非线性变形相位相加，结合时间基线参数，得到每个cs点目标的时间序列变形相位。
[0223]
⑷
形变量计算
[0224]
①
视线向形变量计算与垂直向转换
[0225]
依据雷达波长参数，经过相位转形变计算之后，根据需要结合外部辅助数据，将解缠相位换算为视线方向的形变量，然后依据视线向与垂直向夹角，将视线向形变量转换为垂直向。
[0226]
②
地理编码
[0227]
地理编码的方法可利用dem产品进行地理编码。具体步骤如下：
[0228]
a.利用dem坐标系到sar影像坐标系的转换查找表，完成监测成果由sar影像坐标系到大地坐标系的反变换，即对监测成果垂直向形变量进行地理编码。
[0229]
b.集合所有地理编码后的点目标，将变形量的时间单位换算成年，生成年度变形速率，逐像元计算生成地质灾害体速率图。
[0230]
③
基准修正
[0231]
地理编码后点目标的灾害体速率可以利用已有的水准、gps等高精度控制点数据(同期观测的变形量)修正基准，具体步骤为：
[0232]
a.以同步水准测量结果作为基准参考，在临近点上计算点目标变形量与实测值之间差值的平均值，即与实测变形量之间存在的整体偏差值；
[0233]
b.上一步得到的整体偏差值加入每个点目标的变形值，修正因参考点不统一产生的insar结果变形量的整体偏差，完成基准修正。
[0234]
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。