技术领域
[0001] 本
发明属于地表三维形变监测技术领域，具体涉及一种融合SAR与光学偏移量技术的地表三维形变计算方法及系统，可应用于
地震、滑坡、火山等地质灾害监测领域的研究中。
背景技术
[0002] 目前，遥感技术能够有效应用于针对诸如地震、滑坡、火山、地裂缝等地质灾害的调查、监测及预警工作中。其中基于SAR数据和光学数据的两类技术作为遥感领域重要的两个分支，在地质灾害研究中均充当了重要的
角色。
[0003] 针对光学的应用主要通过目视解译或多个波段数据的组合变换进而提取表征地质灾害体特征的因素，此类方法皆是对灾害体展开定性的研究，利用光学偏移量技术可以实现地质灾害所引起地表形变的定量分析，然而传统的方法往往只能获取地表二维形变。
[0004] 基于SAR数据的地表形变监测技术主要包括
合成孔径雷达干涉测量(InSAR)技术及SAR偏移量技术等，目前利用SAR数据的合成孔径雷达干涉测量(InSAR)及偏移量技术能够有效的应用于地表三维形变场的恢复，但为了获取研究区域地表三维形变，传统的方法至少需要三个来自于不同雷达成像几何结构的SAR数据。当可用SAR卫星数据源较少，造成该方法在实际应用中不足以满足三维形变场计算条件时，这些技术便无法开展。
发明内容
[0005] 针对
现有技术的
缺陷和不足，本发明给出了一种融合SAR与光学偏移量技术的地表三维形变计算方法及系统，解决现有的地表三维形变解算在SAR数据源不足时使用受限的问题。
[0006] 为了实现上述目的，本发明采用如下技术方案予以实现：
[0007] 融合SAR与光学偏移量技术的地表三维形变计算方法，包括以下步骤：
[0008] 步骤1，获取待监测区域形变前后的两景SAR数据，将两景SAR
数据处理为两景单视复数影像，根据两景单视复数影像计算待监测区域地表的方位向形变DAZI和视线向形变DLOS；
[0009] 步骤2，获取待监测区域形变前后的两景光学数据，两景光学数据与两景SAR数据
覆盖区域相同，根据两景光学数据计算地表的南北向形变DNS和东西向形变DEW；
[0010] 步骤3，对步骤1获得的方位向形变DAZI、视线向形变DLOS和步骤2获得的南北向形变DNS和东西向形变DEW进行重
采样，使其具有相同的空间
分辨率；
[0011] 步骤4，通过公式(1)，获得地表的真实垂直向形变Dup、真实南北向形变Dns和真实东西向形变Dew；
[0012]
[0013] 式中，θ表示卫星入射角，表示卫星飞行方位角。
[0014] 具体的，所述的步骤1中利用SAR偏移量
跟踪方法计算待监测区域地表的方位向形变DAZI和视线向形变DLOS。
[0015] 具体的，所述的步骤2包括：
[0016] 步骤2.1，对于获取的两景光学数据，判断两景光学数据的覆盖范围大小是否相同，若相同，进行步骤2.2；若不同，裁剪出两景数据相同的公共区域，对两景公共区域进行步骤2.2；
[0017] 步骤2.2，使用基于ENVI
软件的COSI-Corr光学偏移量计算工具获得地表的南北向形变量DNS和东西向形变量DEW。
[0018] 优选的，所述的步骤2还包括对南北向形变DNS和东西向形变DEW的去除误差过程：根据得到的南北向形变DNS和东西向形变DEW掩模形变区域，在远离形变区域的
位置选取若干
像素点，根据这些像素点的形变值及地理位置利用一次拟合多项式拟合整个监测区域的轨道误差趋势面，然后将原始形变图与此误差趋势面相减，得到去除误差后的南北向形变DNS和东西向形变DEW。
[0019] 优选的，所述的步骤4中，在使用公式(1)计算Dup、Dns和Dew时，采用赫尔默特方差分量估计的方法
迭代来最终确定SAR与光学两类形变观测值的权，进而求得地表的真实垂直向形变量Dup、真实南北向形变量Dns和真实东西向形变量Dew。
[0020] 本发明还公开了融合SAR与光学偏移量技术的地表三维形变计算系统，该系统包括：
[0021] SAR偏移量计算模
块，用于获取待监测区域形变前后的两景SAR数据，将两景SAR数据处理为两景单视复数影像，根据两景单视复数影像，利用SAR偏移量跟踪方法计算待监测区域地表的方位向形变DAZI和视线向形变DLOS；
[0022] 光学偏移量计算模块，用于获取待监测区域形变前后的两景光学数据，两景光学数据与SAR偏移量计算模块中两景SAR数据覆盖区域相同，根据两景光学数据计算地表的南北向形变DNS和东西向形变DEW；
[0023] 重采样模块，用于对SAR偏移量计算模块获得的方位向形变DAZI、视线向形变DLOS和光学偏移量计算模块获得的南北向形变DNS和东西向形变DEW进行重采样，使其具有相同的空间分辨率；
[0024] 地表三维形变计算模块，用于通过公式(1)，获得地表的真实垂直向形变Dup、真实南北向形变Dns和真实东西向形变Dew；
[0025]
[0026] 式中，θ表示卫星入射角，表示卫星飞行方位角。
[0027] 具体的，所述的SAR偏移量计算模块中利用SAR偏移量跟踪方法计算待监测区域地表的方位向形变DAZI和视线向形变DLOS。
[0028] 具体的，所述的光学偏移量计算模块包括：
[0029] 对于获取的两景光学数据，判断两景光学数据的覆盖范围大小是否相同，若相同，使用基于ENVI软件的COSI-Corr光学偏移量计算工具获得地表的南北向形变量DNS和东西向形变量DEW；若不同，裁剪出两景数据相同的公共区域，再使用基于ENVI软件的COSI-Corr光学偏移量计算工具生成地表的南北向形变量DNS和东西向形变量DEW。
[0030] 进一步的，所述的光学偏移量计算模块还包括误差去除模块，用于去除南北向形变DNS和东西向形变DEW中的误差，具体为：根据得到的南北向形变DNS和东西向形变DEW掩模形变区域，在非形变区域的位置选取若干像素点，根据这些像素点的形变值及地理位置利用一次拟合多项式模拟整个待监测区域的轨道误差趋势面，然后将原始形变图与此误差趋势面相减，得到的去除误差后的南北向形变DNS和东西向形变DEW。
[0031] 优选的，所述的地表三维形变计算模块中，在使用公式(1)计算Dup、Dns和Dew时，采用赫尔默特方差分量估计的方法迭代来最终确定SAR与光学两类形变观测值的权，进而求得地表的真实垂直向形变量Dup、真实南北向形变量Dns和真实东西向形变量Dew。
[0032] 与现有技术相比，本发明的有益效果是：
[0033] 本发明通过融合SAR与光学偏移量技术得到一种用于计算地表三维形变的方法，该方法的优点在于利用SAR与光学数据，通过偏移量技术获得地表形变，然后通过建立地表真实三维形变求解模型，且通过赫尔默特方差分量估计的方法最终确定SAR与光学两类形变观测值的权重，最终能够有效地获得地表三维真实形变，解决了现有仅依赖SAR数据的地表三维形变解算在SAR数据源不足时无法恢复地表三维形变的问题，为地表三维形变计算提供了一种新的方法，此工作对于研究地表形变具有重要意义。
附图说明
[0034] 图1为本发明
实施例记载的地表三维形变计算方法
流程图。
[0035] 图2为实施例记载的凯库拉地震沿方位向形变监测结果。
[0036] 图3为实施例记载的凯库拉地震沿视线向形变监测结果。
[0037] 图4为实施例记载的凯库拉地震在南北向上的形变监测结果，图中沿北方向形变为正，沿南方向形变为负。
[0038] 图5为实施例记载的凯库拉地震在东西向上的形变监测结果，图中沿东方向形变为正，沿西方向形变为负。
[0039] 图6为利用本发明所得凯库拉地震地表三维形变场；图中黑色箭头表示地表
水平向形变，箭头方向为水平形变方向，长度表示形变大小，
颜色变化表示垂直向形变。
[0040] 以下结合
说明书附图和具体实施方式对本发明做具体说明。
具体实施方式
[0041] 以下给出本发明的具体实施方式，需要说明的是本发明并不局限于以下具体实施例中，凡在本
申请技术方案
基础上做的等同变换均落入本发明的保护范围。
[0042] 如图1所示，本发明实施例提供的一种融合SAR与光学偏移量技术的地表三维形变计算方法包括以下步骤：
[0043] 步骤1，获取待监测区域的两景SAR数据，如地震前后同一区域的两景数据；将两景SAR数据处理为两景单视复数影像，利用SAR偏移量跟踪方法根据两景单视复数影像计算待监测区域地表的方位向形变DAZI和视线向形变DLOS。
[0044] SAR偏移量跟踪方法的具体思路是：
选定其中一景作为主影像，将另一景从影像配准重采样到该主影像，使主影像和从影像具有相同的空间几何结构；利用SAR幅度信息对两景SAR影像进行粗配准，计算主、从影像之间的初始偏移量，然后在粗配准的基础上，选取一定大小的搜索窗口并计算主从影像之间的归一化互相关系数ρ(x,y)，如式(2)所示，[0045]
[0046] 式中，x和y为主影像坐标，xs和ys为两个方向上的偏移量，m和n为窗口大小，fm和fs分别为主影像幅度和从影像幅度。
[0047] 当互相关系数ρ(x,y)达到最大时便可获得精确的子像素级的偏移量，最终获得地表的方位向形变DAZI和视线向形变DLOS。
[0048] 步骤2，获取待监测区域形变前后的两景光学数据，两景光学数据与步骤1中两景SAR数据获取时间点大致相同、覆盖区域相同。根据两景光学数据计算地表的光学偏移量，即南北向形变DNS和东西向形变DEW；
[0049] 需要注意的是，在利用两景光学数据计算地表的光学偏移量时，首先判断两景光学数据的覆盖范围大小是否相同，若相同，则使用基于ENVI软件的COSI-Corr光学偏移量计算工具生成地表的南北向形变量DNS和东西向形变量DEW；若不同，需要先裁剪出两景数据相同的公共区域，对两景公共区域再使用基于ENVI软件的COSI-Corr光学偏移量计算工具生成地表的南北向形变量DNS和东西向形变量DEW。
[0050] 对于得到的南北向形变量DNS和东西向形变量DEW，由于一些光学影像数据没有经过严格的正射校正和几何校正，所以在计算所得形变场中存在明显的轨道误差，因此需要去除该误差，具体为：根据得到的南北向形变DNS和东西向形变DEW掩模形变区域，在非形变区域的位置选取若干像素点，根据这些像素点的形变值及地理位置利用一次拟合多项式拟合整个待监测区域的轨道误差趋势面，然后将原始形变图与此误差趋势面相减，得到去除误差后的南北向形变DNS和东西向形变DEW。
[0051] 步骤3，对步骤1获得的方位向形变DAZI、视线向形变DLOS和步骤2获得的南北向形变DNS和东西向形变DEW进行重采样，使其具有相同的空间分辨率，并裁剪为相同的范围。
[0052] 另外，若光学偏移量结果和SAR偏移量结果的
坐标系不一致时，需要将两者转化为相同的坐标系，如光学偏移量结果通常为UTM坐标系，而经过地理编码后的SAR偏移量结果往往为WGS 84地理坐标系，因此首先需要将光学偏移量结果由UTM坐标系转为与SAR结果相同的地理坐标系。
[0053] 步骤4，根据SAR卫星成像几何关系建立方位向形变和视线向形变与地表真实三维形变之间的关系式(3)和(4)，
[0054]
[0055]
[0056] 式中，Dup表示地表的真实垂直向形变，Dns表示真实南北向形变，Dew表示真实东西向形变；
[0057] 再结合光学所得南北向和东西向形变，可建立SAR与光学偏移量联合求解得到地表三维形变的观测方程，如式(1)所示。
[0058]
[0059] 式中，θ表示卫星入射角，表示卫星飞行方位角。
[0060] 最后通过公式(1)，获得地表的真实垂直向形变Dup、真实南北向形变Dns和真实东西向形变Dew。
[0061] 在使用公式(1)计算Dup、Dns和Dew时，采用赫尔默特方差分量估计的方法通过迭代定权，最终确定SAR与光学形变观测值在计算过程中的权重，并最终求得地表三维形变场Dup、Dns和Dew。具体步骤如下：
[0062] 步骤1)，首先确定初始方差，通过给定一个遍历窗口，分别计算方位向、视线向、东西向及南北向形变量文件中每个中心像元在该窗口内的方差作为初始方差，利用式(5)权计算公式计算初始权，式中c为任意给定常数，σi2为第i类观测方差，pi即为权值。令计算所得SAR与光学两类形变观测值的权分别为p1和p2；
[0063]
[0064] 步骤2)，令式(4)中L1＝[DLOS DAZI]T，L2＝[DNS DEW]T， A为系数矩阵，A＝[A1 A2]T，其中，
[0065]
[0066]
[0067] 将式(4)改写为误差方程,如式(8)所示，
[0068]
[0069] 由式(8)可得法方程为
[0070]
[0071] 式中，
[0072] 步骤3)，最终建立残差平方和 和 之间的关系式如式(10)所示，其中，S矩阵如式(11)所示，
[0073]
[0074]
[0075] 根据公式(10)求得两类观测值迭代后单位权方差，计算公式如式(12)所示，[0076]
[0077] 步骤4)，根据所求得两类观测值迭代后单位权方差，利用公式(13)再次定权，式中，k为迭代次数，
[0078]
[0079] 然后重复步骤2)和步骤3)，进行迭代，直至两类观测值迭代后的单位权方差相等迭代停止，最终求得 即求得地表三维形变。
[0080] 本发明的实施例还提供了一种融合SAR与光学偏移量技术的地表三维形变计算系统，该系统包括：
[0081] SAR偏移量计算模块，用于获取待监测区域形变发生前后的两景SAR数据，将两景SAR数据处理为两景单视复数影像，利用SAR偏移量跟踪方法根据两景单视复数影像中计算待监测区域地表的方位向形变DAZI和视线向形变DLOS。SAR偏移量跟踪方法的具体思路同上文记载。
[0082] 光学偏移量计算模块，用于获取待监测区域形变发生前后的两景光学数据，两景光学数据与SAR偏移量计算模块中两景SAR数据覆盖区域相同，根据两景光学数据，利用光学偏移量技术计算地表的南北向形变DNS和东西向形变DEW；
[0083] 需要注意的是，在利用两景光学数据计算地表的光学偏移量时，首先判断两景光学数据的覆盖范围大小是否相同，若相同，则使用基于ENVI软件的COSI-Corr光学偏移量计算工具生成地表的南北向形变量DNS和东西向形变量DEW；若不同，需要先裁剪出两景数据的公共区域，对两景公共区域再使用基于ENVI软件的COSI-Corr光学偏移量计算工具生成地表的南北向形变量DNS和东西向形变量DEW。
[0084] 对于得到的南北向形变量DNS和东西向形变量DEW，由于一些光学影像数据没有经过严格的正射校正和几何校正，所以在计算所得形变场中存在明显的轨道误差，因此还设置了误差去除模块，误差去除模块用于通过掩模得到的南北向形变DNS和东西向形变DEW形变区域，在非形变区域的位置选取若干像素点，根据这些像素点的形变值及地理位置利用一次拟合多项式拟合整个待监测区域的轨道误差趋势面，然后将原始形变图与此误差趋势面相减，得到去除误差后的南北向形变DNS和东西向形变DEW。
[0085] 重采样模块，用于对SAR偏移量计算模块获得的方位向形变DAZI、视线向形变DLOS和光学偏移量计算模块获得的南北向形变DNS和东西向形变DEW进行重采样，使其具有相同的空间分辨率；
[0086] 另外，若光学偏移量结果和SAR偏移量结果的坐标系不一致时，需要设置一个坐标系转化模块，该模块用于将两个不同的坐标系转化为相同的坐标系，如光学偏移量结果通常为UTM坐标系，而经过地理编码后的SAR偏移量结果往往为WGS 84地理坐标系，因此首先需要将光学偏移量结果由UTM坐标系转为与SAR结果相同的地理坐标系。
[0087] 地表三维形变计算模块，用于通过公式(1)，获得地表的真实垂直向形变Dup、真实南北向形变Dns和真实东西向形变Dew；
[0088]
[0089] 式中，θ表示卫星入射角，表示卫星飞行方位角。
[0090] 在地表三维形变计算模块中，使用公式(1)计算Dup、Dns和Dew时，采用赫尔默特方差分量估计的方法通过迭代定权，最终确定SAR与光学形变观测值在计算过程中的权重，并最终求得地表三维形变场Dup、Dns和Dew。具体步骤同上文中的步骤1)至步骤4)。
[0091] 实施例
[0092] 本实施例选择位于新西兰东北部的凯库拉县地区作为实验区域，该地区于2016年11月14日发生了MW7.8级大地震，地震发生之后，地表破裂明显，且形变量级较大。实验选用的SAR数据为2016年2月20日和2016年12月10日覆盖该地区的两景地震发生前后Sentinel-
1数据，光学数据为2016年2月19日和2017年2月23日覆盖该地区的两景地震发生前后Sentinel-2数据，由于本实施例仅研究该地震同震三维形变的计算，且同震形变量级远远大于震间形变，因此震后的SAR与光学数据获取时间
节点相差时间内地表发生的震间形变可忽略不计，即可视为数据获取时间大致相同。
[0093] 实验过程
[0094] 步骤1，首先对于获取的两景Sentinel-1数据进行预处理生成单视复数影像(SLC)，然后选定2016年2月20日影像作为主影像，将另一景从影像配准重采样到该主影像，使主从影像具有相同的空间几何结构。然后利用SAR偏移量技术，设置搜索窗口为64*64，搜索步长为8*8，计算方位向和视线向二维形变场，分别如图2和图3所示。图中Azimuth为卫星飞行方向，即为方位向，LOS为沿雷达视线方向，即为视线向，远离卫星方向形变量为负，靠近卫星方向形变量为正。
[0095] 步骤2，获取与SAR数据覆盖相同地区的2016年2月19日和2017年2月23日两景Sentinel-2数据，首先进行数据预处理，包括
辐射定标和大气校正，最终提取第8波段进行光学偏移量的计算。由于数据覆盖范围大小不同，因此需裁剪出两景数据的公共区域。然后利用基于ENVI软件的COSI-Corr光学偏移量计算工具计算光学偏移量，设置搜索窗口大小为32*32，移动步长为8*8，掩膜
阈值为0.9，迭代次数为2次，计算获得地表在南北向和东西向上的形变。由于一些光学影像数据没有经过严格的正射校正和几何校正，所以在计算所得形变场中存在明显的轨道误差，因此在非形变区域的位置选取若干像素点，根据这些像素点的形变值及地理位置利用一次拟合多项式模拟整个研究区域的轨道误差趋势面，然后将原始形变图与此误差趋势面相减，达到去除轨道误差的目的，最终获得地表在南北向和东西向上的形变如图4和图5所示。
[0096] 步骤3，对SAR与光学偏移量结果进行重采样，使二者具有相同的空间分辨率和坐标系统。由于Sentinel-2光学偏移量的结果为UTM坐标系，因此首先需要将光学偏移量结果由UTM坐标系转为WGS 84下的地理坐标系，然后将SAR与光学结果重采样为相同的空间分辨率，并且裁剪为相同的范围，最终导出为TIFF格式文件。
[0097] 步骤4，根据式(1)，通过步骤1)至步骤4)获得地表的真实垂直向形变Dup、真实南北向形变Dns和真实东西向形变Dew，如图6所示为地表的三维形变场。由图6中矢量箭头的方向及垂直方向的形变可以看出地震造成各板块之间相互作用，引起各板块的地表发生不同的形变，从而能够有效地解释地震所造成板块运动的机理。
