本发明提供了一种基于Offset Tracking技术的边坡形变时序监测方法，包括：步骤1：基于SAR数据，采用Offset Tracking技术计算目标区域的二维形变；步骤2：结合地形数据和卫星轨道数据，反演目标区域的准三维形变；步骤3：基于地形数据，校验每个测量点的准三维形变测量结果，剔除测量误差点；步骤4：综合多时间段的形变测量结果，计算目标点的时序形变。该发明能够提升Offset Tracking技术的适用性，为技术的工程化应用提供技术支持。
一种基于Offset Tracking技术的边坡形变时序监测方法
技术领域
[0001] 本发明涉及合成孔径雷达技术领域，特别涉及一种基于Offset Tracking技术的边坡形变时序监测方法。
背景技术
[0002] 近年来，国内每年都会发生多起滑坡事故。一般来说，在边坡宏观失稳之前，其表面都会先发生微小位移。因此，为了更深入地研究边坡滑坡的触发机理，并在第一时间准确
地实现临滑预报，我们需要对边坡进行全覆盖，长时间、连续、实时和高精度的形变监测。作
为一种遥感技术手段，星载SAR技术具有监测覆盖范围大、全天时全天候观测等优势。基于
星载SAR数据的InSAR技术和Offset Tracking技术都具备形变监测能力。
[0003] Offset Tracking技术的测量精度比InSAR技术低一个数量级，当图像分辨率为3m时，形变测量精度可达15cm。但是，与InSAR数据处理需要积累大量SAR图像的条件不同，
Offset Tracking技术只需要两幅SAR图像就能开展形变监测，且不存在周期模糊现象，即
使在干涉失相干的条件下，也能在方位向和距离向上获取几米至数十米的二维剧烈形变信
息。
[0004] 但是，在Offset Tracking技术的工程化应用过程中，还存在两个问题：(1)Offset Tracking技术获取的是距离向和方位向二维形变，而不是空间三维形变，难以较直观地展
示空间形变特征；(2)Offset Tracking技术是采用计算窗口相关性的方式来测量形变，当
窗口的选择不适当时，可能导致形变反演错误。
[0005] 针对这两个问题，本发明首先基于SAR数据，采用Offset Tracking技术计算目标区域的二维形变（距离向和方位向）；然后，结合地形数据和卫星轨道数据，反演目标的准三
维形变；接着，基于地形数据，校验每个测量点的准三维形变测量结果，剔除测量误差点；最
后，综合多时间段的形变测量结果，计算目标点的时序形变。本发明能够提升Offset
Tracking技术的适用性，为技术的工程化应用提供技术支持。
发明内容
[0006] 本发明提供一种基于Offset Tracking技术的边坡形变时序监测方法，能够提升Offset Tracking技术的适用性，为技术的工程化应用提供技术支持。
[0007] 本发明提供一种基于Offset Tracking技术的边坡形变时序监测方法，包括：
[0008] 步骤1：基于SAR数据，采用Offset Tracking技术计算目标区域的二维形变；
[0009] 步骤2：结合地形数据和卫星轨道数据，反演目标区域的准三维形变；
[0010] 步骤3：基于地形数据，校验每个测量点的准三维形变测量结果，剔除测量误差点；
[0011] 步骤4：综合多时间段的形变测量结果，计算目标点的时序形变。
[0012] 在一种可能实现的方式中，步骤1中，基于SAR数据，采用Offset Tracking技术计算目标区域的二维形变，包括：
[0013] 基于SAR数据，并利用影像互相关算法，对形变前后的SAR图像进行粗配准和空间重采样；
[0014] 选用滑动窗口，计算子像元间的相关性；
[0015] 提取对应方位向和距离向的偏移量，同时拟合并去除由两幅图像成像时间和空间不同所带来的系统误差偏移值；
[0016] 获得对应方位向和距离向的地表形变信息；
[0017] 其中，二维形变指的是方位向和距离向的地表形变信息。
[0018] 在一种可能实现的方式中，步骤2中，结合地形数据和卫星轨道数据，反演目标区域的准三维形变，包括：
[0019] 以目标区域中心点为坐标原点，建立空间直角坐标系，其中，x轴正方向指向正东，y轴正方向指向正北，z轴正方向指向天空；
[0020] 基于Offset Tracking技术的形变反演，获取目标区域中的目标点在距离向和方位向的二维形变，其中，距离向的形变向量dRa和方位向的形变向量dAz，在空间直角坐标系
中分别表达为：
[0021] dRa=[xRa,yRa,zRa]=|dRa|[uxRa,uyRa,uzRa]=|dRa|uRa
[0022] dAz=[xAz,yAz,zAz]=|dAz|[uxAz,uyAz,uzAz]=|dAz|uAz
[0023] 其中，uRa为距离向形变的单位向量，uAz为方位向形变的单位向量，|dRa|为距离向测得的形变值，|dAz|为方位向测得的形变值，uxRa为距离向形变单位向量uRa在x轴的分量，
uyRa为距离向形变单位向量uRa在y轴的分量，uzRa为距离向形变单位向量uRa在z轴的分量；xRa
为距离向形变向量dRa在x轴的分量，yRa为距离向形变向量dRa在y轴的分量，zRa为距离向形变
向量dRa在z轴的分量；uxAa为方位向形变单位向量uAz在x轴的分量，uyAa为方位向形变单位向
量uAz在y轴的分量，uzAa为方位向形变单位向量uAz在z轴的分量；xAa为方位向形变向量dAz在x
轴的分量，yAa为方位向形变向量dAz在y轴的分量，zAa为方位向形变向量dAz在z轴的分量；
[0024] 基于地形数据，计算每个目标点位置处的坡度角θ和坡向角φ，并在每个目标点的临近区域建立一个滑动平面，其中，所述滑动平面的方程表示为：
[0025] f(θ,φ)x+g(θ,φ)y+z=0
[0026] 式中，f(θ,φ)和g(θ,φ)为滑动平面的方程系数，它们由坡度角θ和坡向角φ决定，根据空间几何知识，f(θ,φ)和g(θ,φ)通过以下代数方程求解：
[0027]
[0028] 基于上述理论，定义目标点的三维形变向量d=［x,y,z］，基于计算的二维形变以及目标点形变位于滑动平面的基本理论假设，建立以下方程组：
[0029]
[0030] 通过解算上述方程组，可获取目标点的准三维形变向量d，进而反演得到目标区域的准三维形变。
[0031] 在一种可能实现的方式中，步骤3：基于地形数据，校验每个测量点的准三维形变测量结果，剔除测量误差点，包括：
[0032] 基于三维形变反演结果进行第一校验；
[0033] 基于SAR几何照射关系进行第二校验；
[0034] 根据第一校验结果、第二校验结果，剔除测量误差点。
[0035] 在一种可能实现的方式中，基于三维形变反演结果进行第一校验，包括：
[0036] 获取目标点的坡度角θ和坡向角φ；
[0037] 计算所述目标点沿坡向下的梯度矢量S，其数学表达式为：
[0038] S=[tanθcosφ,tanθsinφ,‑1]
[0039] 所述目标点的准三维形变量d与梯度矢量S的夹角φ由下式计算：
[0040]
[0041] 将所述夹角φ的绝对值与夹角阈值φTh进行比较，根据比较结果对测量误差点进行第一剔除。
[0042] 在一种可能实现的方式中，基于SAR几何照射关系进行第二校验，包括：
[0043] 基于地形数据和卫星轨道数据，获取每个目标点的对应的雷达方向向量aLOS和本地法向量aS；
[0044] 确定雷达LOS方向向量aLOS和本地法向量aS的夹角关系，判断所述目标点是否处于阴影和叠掩状态；
[0045] 若是，对测量误差点进行第二剔除。
[0046] 在一种可能实现的方式中，确定雷达LOS方向向量aLOS和本地法向量aS的夹角关系，判断所述目标点是否处于阴影和叠掩状态，包括：
[0047] 计算雷达LOS方向向量aLOS和本地法向量aS的夹角φ;
[0048] 当夹角φ大于90°时，判定所述目标点在SAR图像中处于阴影状态；
[0049] 当夹角φ的绝对值小于φTh时，判定所述目标点在SAR图像中处于叠掩状态。
[0050] 在一种可能实现的方式中，步骤4中，综合多时间段的形变测量结果，计算目标点的时序形变，包括：
[0051] 在目标区域内建立规则网格，并对网格上的点进行插值处理；
[0052] 根据插值处理结果，获取形变序列结果；
[0053] 形变序列结果，获取目标区域不同目标点的形变演化特征，并计算获取所述目标点的时序形变。
[0054] 本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明
书、权利要求书、以及附图中所特别指出的结构来实现和获得。
[0055] 下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
[0056] 附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：
[0057] 图1基于OffsetTracking技术的边坡形变时序监测方法流程图；
[0058] 图2由坡度和坡向角构成的目标点滑动平面示意图；
[0059] 图3目标点阴影状态的判断方法；
[0060] 图4目标点叠掩状态的判断方法；
[0061] 图5基于OffsetTracking技术直接反演的三维形变结果；
[0062] 图6误差点剔除后的三维形变结果；
[0063] 图7某滑坡距离向空间累积形变量；
[0064] 图8某滑坡方位向空间累积形变量；
[0065] 图9某滑坡2018.11.10 2019.10.07准三维形变示意图；~
[0066] 图10为对应位置在距离向的空间累积形变历史曲线；
[0067] 图11为对应位置在方位向的空间累积形变历史曲线。
具体实施方式
[0068] 以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
[0069] 本发明提供一种基于Offset Tracking技术的边坡形变时序监测方法，如图1所示，包括：
[0070] 步骤1：基于SAR数据，采用Offset Tracking技术计算目标区域的二维形变；
[0071] 步骤2：结合地形数据和卫星轨道数据，反演目标区域的准三维形变；
[0072] 步骤3：基于地形数据，校验每个测量点的准三维形变测量结果，剔除测量误差点；
[0073] 步骤4：综合多时间段的形变测量结果，计算目标点的时序形变。
[0074] 该实施例中，在进行步骤1‑4的过程中，是基于DEM数据（地形高程数据）实现的。
[0075] 该实施例中，通过对目标区域中的每个测量点进行校验，然后剔除其中的误差点，最后获取得到目标点。
[0076] 上述技术方案的有益效果是：通过遥感测量技术手段，获取边坡区域的面状形变特征，剔除测量误差点，提升形变测量精度，特别是能适用于难以布设地面传感器的高陡边
坡或危险边坡，最终为边坡形变监测和风险早期识别提供高效的技术支持。
[0077] 在一种可能实现的方式中，步骤1中，基于SAR数据，采用Offset Tracking技术计算目标区域的二维形变，包括：
[0078] 基于SAR数据，并利用影像互相关算法，对形变前后的SAR图像进行粗配准和空间重采样；
[0079] 选用滑动窗口，计算子像元间的相关性；
[0080] 提取对应方位向和距离向的偏移量，同时拟合并去除由两幅图像成像时间和空间不同所带来的系统误差偏移值；
[0081] 获得对应方位向和距离向的地表形变信息；
[0082] 其中，二维形变指的是方位向和距离向的地表形变信息。
[0083] 该实施例中，Offset Tracking技术利用影像互相关算法，通过追踪雷达幅度影像中特征目标的位置变化来监测地表形变。
[0084] 在一种可能实现的方式中，如图2所示，步骤2中，结合地形数据和卫星轨道数据，反演目标区域的准三维形变，包括：
[0085] 以目标区域中心点为坐标原点，建立空间直角坐标系，其中，x轴正方向指向正东，y轴正方向指向正北，z轴正方向指向天空；
[0086] 基于Offset Tracking技术的形变反演，获取目标区域中的目标点在距离向和方位向的二维形变，其中，距离向的形变向量dRa和方位向的形变向量dAz，在空间直角坐标系
中分别表达为：
[0087] dRa=[xRa,yRa,zRa]=|dRa|[uxRa,uyRa,uzRa]=|dRa|uRa
[0088] dAz=[xAz,yAz,zAz]=|dAz|[uxAz,uyAz,uzAz]=|dAz|uAz
[0089] 其中，uRa为距离向形变的单位向量，uAz为方位向形变的单位向量，|dRa|为距离向测得的形变值，|dAz|为方位向测得的形变值，uxRa为距离向形变单位向量uRa在x轴的分量，
uyRa为距离向形变单位向量uRa在y轴的分量，uzRa为距离向形变单位向量uRa在z轴的分量；xRa
为距离向形变向量dRa在x轴的分量，为距离向形变向量dRa在y轴的分量，zRa为距离向形变向
量dRa在z轴的分量；uxAa为方位向形变单位向量uAz在x轴的分量，uyAa为方位向形变单位向量
uAz在y轴的分量，uzAa为方位向形变单位向量uAz在z轴的分量；xAa为方位向形变向量dAz在x轴
的分量，yAa为方位向形变向量dAz在y轴的分量，zAa为方位向形变向量dAz在z轴的分量；
[0090] 基于地形数据，计算每个目标点位置处的坡度角θ和坡向角φ，并在每个目标点的临近区域建立一个滑动平面，其中，所述滑动平面的方程表示为：
[0091] f(θ,φ)x+g(θ,φ)y+z=0
[0092] 式中，f(θ,φ)和g(θ,φ)为滑动平面的方程系数，它们由坡度角θ和坡向角φ决定，根据空间几何知识，f(θ,φ)和g(θ,φ)通过以下代数方程求解：
[0093]
[0094] 基于上述理论，定义目标点的三维形变向量d=［x,y,z］，基于计算的二维形变以及目标点形变位于滑动平面的基本理论假设，建立以下方程组：
[0095]
[0096] 通过解算上述方程组，可获取目标点的准三维形变向量d，进而反演得到目标区域的准三维形变。
[0097] 该实施例中，结合滑坡的基本知识，目标点的形变向量是在该平面内的。
[0098] 上述技术方案的有益效果是：便于获取准三维形变，为边坡形变监测和风险早期识别提供高效的基础。
[0099] 在一种可能实现的方式中，步骤3：基于地形数据，校验每个测量点的准三维形变测量结果，剔除测量误差点，包括：
[0100] 基于三维形变反演结果进行第一校验；
[0101] 基于SAR几何照射关系进行第二校验；
[0102] 根据第一校验结果、第二校验结果，剔除测量误差点。
[0103] 在一种可能实现的方式中，基于三维形变反演结果进行第一校验，包括：
[0104] 获取目标点的坡度角θ和坡向角φ；
[0105] 计算所述目标点沿坡向下的梯度矢量S，其数学表达式为：
[0106] S=[tanθcosφ,tanθsinφ,‑1]
[0107] 所述目标点的准三维形变量d与梯度矢量S的夹角φ由下式计算：
[0108]
[0109] 将所述夹角φ的绝对值与夹角阈值φTh进行比较，根据比较结果对测量误差点进行第一剔除。
[0110] 在一种可能实现的方式中，基于SAR几何照射关系进行第二校验，包括：
[0111] 基于地形数据和卫星轨道数据，获取每个目标点的对应的雷达方向向量aLOS和本地法向量aS；
[0112] 确定雷达LOS方向向量aLOS和本地法向量aS的夹角关系，判断所述目标点是否处于阴影和叠掩状态；
[0113] 若是，对测量误差点进行第二剔除。
[0114] 在一种可能实现的方式中，确定雷达LOS方向向量aLOS和本地法向量aS的夹角关系，判断所述目标点是否处于阴影和叠掩状态，包括：
[0115] 计算雷达LOS方向向量aLOS和本地法向量aS的夹角;
[0116] 当夹角φ大于90°时，判定所述目标点在SAR图像中处于阴影状态；
[0117] 当夹角φ的绝对值小于φTh时，判定所述目标点在SAR图像中处于叠掩状态。
[0118] 该实施例中，在Offset Tracking的算法中，是基于主辅SAR图像窗口内数据的相关性来计算形变量。一般情况下，窗口的选取会影响算得的形变量。因此，Offset Tracking
技术直接生成的形变量可能在局部区域包含错误的测量结果。
[0119] 该实施例中，针对第一校验，一般来说，目标点主要的滑动方向是沿着滑动平面向下的。因此，每个目标点计算的准三维形变向量的方向与该点对应的地形梯度方向应基本
相似，可通过这个条件来剔除形变测量明显出错的目标点。
[0120] 通过设定合适的夹角阈值φTh，选出明显存在误差的形变测量结果（即｜φ｜＞φTh），进而将其选出并剔除。
[0121] 该实施例中，针对第二校验，考虑到处于阴影或叠掩状态下的目标信号质量较差，在后续分析过程中将剔除这些测量点。
[0122] 如图3所示，在判断目标点阴影状态的过程中，可先计算雷达LOS方向向量aLOS和本地法向量aS的夹角φ。当夹角φ大于90°时，可认为该目标点在SAR图像中处于阴影状态。
[0123] 如图4所示，在判断目标点叠掩状态的过程中，可先计算雷达LOS方向向量aLOS和本地法向量aS的夹角φ。当夹角φ的绝对值小于φTh时，可认为该目标点在SAR图像中处于叠
掩状态，一般来说，φTh可设置为5°
[0124] 上述技术方案的有益效果是：通过基于地形数据对形变监测结果进行校验，提升形变监测结果的可靠性。
[0125] 在一种可能实现的方式中，步骤4中，综合多时间段的形变测量结果，计算目标点的时序形变，包括：
[0126] 在目标区域内建立规则网格，并对网格上的点进行插值处理；
[0127] 根据插值处理结果，获取形变序列结果；
[0128] 形变序列结果，获取目标区域不同目标点的形变演化特征，并计算获取所述目标点的时序形变。
[0129] 上述技术方案的有益效果是：将不同时期的测量结果进行时序处理，获取目标区域不同位置的形变演化特征，更好地研究目标区域在长时间范围内的形变特征。
[0130] 针对上述技术方案，还具体通过下文对某边坡进行监测来进行具体说明：
[0131] SAR数据采用3m分辨率的COSMO‑SkyMed数据，其具体参数信息如表1所示。
[0132] 表1 COSMO‑SkyMed数据参数
[0133]
[0134] 步骤一：
[0135] 以20181110和20181122两期数据为例,首先采用OffsetTracking技术获得该时间间隔内的方位向和距离向的二维地表形变成果，基于此二维形变成果，结合准三维形变反
演模型得到准三维形变反演结果，结果如图5所示。由图5可知，初步的准三维形变成果存在
误差点。
[0136] 步骤二：
[0137] 基于本方法提出的模型，结合地形和卫星轨道数据可剔除掉形变测量误差点，进而能提取出目标区域准确的三维形变结果，结果如图6所示。时间上相邻的两期数据依次按
照上述流程处理，得到目标区域相邻时间间隔内的准三维形变信息。
[0138] 步骤三：
[0139] 基于上述得到的各期二维和准三维形变成果，考虑到不同时间段的形变测量结果是由不同的SAR数据对生成，其测量点的位置可能存在差异，先在目标区域内建立规则网
格，然后对网格上的点进行插值处理，最后获取形变序列结果。该滑坡在距离向和方位向的
累积形变分别如图7和图8所示。合成的形变量及形变方向示意图如图9所示，其最大累积形
变量达15m。在该滑坡中选择关键测量目标，其位置在图7中用蓝色五角星标识，该位置在距
离向的空间累积形变历史曲线如图10所示，该位置在方位向的空间累积形变历史曲线如图
11所示。
[0140] 显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围
之内，则本发明也意图包含这些改动和变型在内。
