融合像元偏移跟踪和短基线集的矿区地表形变解算方法
【专利摘要】一种融合像元偏移跟踪和短基线集的矿区地表形变解算方法。适用于地表矿区地表形变研究使用。其步骤为：对矿区地表形变区域进行预估；选取时间序列上相干图中的高相干点；利用时序偏移量跟踪算法恢复大形变区域高相干点的相位整周数N；大形变区域和小形变区域的高相干点联合建立短基线集解算模型；矿区地表高程误差及形变速率的解算。其监测精度高、范围大，实现操作过程简单，费用低，克服了传统时序InSAR方法无法正确获取矿区大形变梯度下的地表沉降问题，也解决了基于SAR幅度信息的时序像元偏移跟踪算法和基于差分相位的短基线集技术难以联合解算地表形变速率和高程误差的问题，具有广泛的实用性。
【专利说明】
融合像元偏移跟踪和短基线集的矿区地表形变解算方法
技术领域
[0001] 本发明设及一种矿区地表形变解算方法，尤其适用于一种融合像元偏移跟踪和短 基线集的矿区地表形变解算方法。 技术背景
[0002] 合成孔径雷达干设测量技术（InSAR，Interferometric Synthetic Aperture Radar)具有全天候、覆盖范围大、精度高、具有穿透力等优点，在地形测绘、形变监测、地表 参数反演等领域的应用越来越广。DInSAR(DifferentialInSAR)主要用于地表形变监测， 在火山监测、地面沉降、地震形变场获取、滑坡监测等领域具有极大的应用潜能。但该技术 受到时间基线、空间基线、大气延迟等因素的影响，仅对部分高相干点或永久散射点进行处 理分析的时序111541?方法应运而生。该类方法主要包括:？5-1]1541?(化1'1]1日]1日]118。日1:1日的1·- InSAR)、SBAS(Small Baseline Subsets)、IPTA(Interferometric Point Target Analysis)等。
[0003] 时序InSAR处理技术中，因主影像不唯一，时空基线较短，因此相干性较高的短基 线集技术(SBAS)应用较广。但该方法在解算时首先要对差分相位进行解缠，而矿区开采沉 陷速度快、形变量大，相位解缠只能解算相邻点形变梯度在(-η，π]间的沉降值，大形变区域 的地表下沉难W通过相位解缠方法正确得到，导致运些区域监测量级较小，难W满足矿区 大形变监测要求。同时，像元偏移跟踪算法能够解算地表大形变，目前主要应用于冰川移 动、地震形变等，但利用SAR幅度信息得到的地表形变精度较低，并且非常耗时，解算结果中 的地形误差与基于差分相位建立的短基线集技术中的地形误差并不等价，无法将其带入短 基线集模型进行直接解算。因此，目前没有将像元偏移跟踪和短基线集融合解算地表形变 和局程误差的方法。
【发明内容】
[0004] 针对上述技术问题，提供一种克服了时序InSAR方法缺点，计算量小，精度高的融 合像元偏移跟踪和短基线集的矿区地表形变解算方法。
[0005] 为实现上述技术目的，本发明的融合像元偏移跟踪和短基线集的矿区地表形变解 算方法步骤如下：
[0006] 步骤1:预估矿区地表形变区域：
[0007] 选取被测因煤炭开采导致的矿区，利用卫星雷达获取选取矿区的SAR影像，利用 DInSAR或开采沉陷预计得到的地表沉降结果，将SAR影像区域中相邻像元形变差大于SAR影 像波长一半的下沉盆地初步划分成大形变区域，将像元形变量小于10mm设定为无形变区 域，其余区域设定为小形变区域；
[0008] 步骤2:选取时间序列上相干图中的高相干点：
[0009] 设定SAR影像的时间和空间基线阔值，选取小于阔值的SAR影像进行干设处理，再 依据相干性获取整个或大部分时间序列上具有高相干性的点位；
[0010] 步骤3:利用时序偏移量跟踪算法恢复大形变区域高相干点的相位整周数N:
[0011] 在划分出的无形变区域内选取相干性大于0.3的高相干点计算整体偏移量，利用 最小二乘原理建立偏移量的双线性多项式内插函数，将形变区域内相干点的影像坐标带入 内插函数计算形变区域内相干点的整体偏移量；
[0012] 利用小窗口最大频谱配准对大形变区域的高相干点进行精确配准，获取局部区域 偏移量；
[0013] 从局部精确偏移量中去除整体偏移量，并去除未下沉区域因地形起伏引起的残 差，从而获取最终的视线向像兀偏移量off set;
[0014] 利用视线向像元偏移量offset得到大形变区域高相干点的相位整周数N，
[0015] 利用相位整周数N和大形变区域内高相干点的原差分缠绕相位联合恢复大形变区 域高相干点的真实相位差S Ob;
[0016] 步骤4:利用大形变区域和小形变区域的高相干点联合建立短基线集解算模型：
[0017] 在随机时间点ti、随机时间点t2分别采集了两幅被测矿区的SAR影像信息，判断步 骤2选取的相干点位于SAR影像的小形变和无形变区域时，通过传统最小费用流方法解缠差 分干设图中高相干点（1，m)的真实相位差δ Φ S;
[0018] 判断高相干点（l，m)位于大形变区域时，不能采用传统最小费用流方法解缠差分 干设图中高相干点（l，m)真实相位差，此时则由步骤3获取高相干点（l，m)的真实相位差δ 巫b ;
[0019]最终，利用δΦ3、δ(1Π 联合构建短基线集的解算模型；
[0020] 步骤5,解算矿区地表高程误差及形变速率：
[0021] 利用短基线集处理模型，采用最小二乘解算矿区高相干点处地表高程误差及形变 速率，在获取矿区地表形变后，通过反演矿区概率积分法预计参数，最终获取矿区地表沉陷 预计及分析模型，为矿区沉陷工程应用提供基础数据。
[002^ 所述未下沉区域因地形起伏引起的残差利用公式:峨妊,。,。=。。巧+。2(//-豆V 计算得到；式中，offsettopo为地形起伏引起的偏移量误差；曰日，曰1，曰2为最小二乘的拟合系 数;Η为对应位置的高程;忌为研究区域的平均高程；
[0023] 利用公式:N = ceil(2offset · Δ/λ)，计算得到大形变区域高相干点的相位整周 数Ν，式中，ceil()为向下取整运算，Δ为斜距向像元尺寸，λ为SAR影像波长;从而克服像元 偏移跟踪算法求解的偏移量仅用到SAR影像幅度信息，不能同传统短基线技术得到的差分 相位进行最小二乘建模的问题。
[0024] 利用公式：S〇b = 2炯+δφ<Η??，计算大形变区域高相干点的真实相位差，式中，δ Φ diff为缠绕的差分相位；
[0025] 利用真实相位差δ Ob建立短基线集处理模型A方法为:首先利用真实相位差δ Ob构建解 缠相位矩阵δ Φ = [ δ Φ sS Φ b]，再利用公式：
建立短基线集处理模型;式中：A hi,m为相干点（l，m)的高程误差;V为相干点的下沉速率;W 为残余相位;λ为卫星雷达波长;B丄,i,m为两幅SAR影像相干点处的垂直基线；0(l，m)为相干 点处的雷达入射角;ri,m为卫星距离点目标(1，m)的斜距；
[0026] 局部区域偏移量还可W通过最小二乘匹配方法对大形变区域的高相干点进行精 确配准获取；
[0027] 所述矿区高相干点处地表高程误差及形变速率利用公式：
计算，式中P为权阵，可取单位阵。
[0028] 有益效果:本发明将采用幅度获取地表大形变的时序偏移量跟踪算法和利用相位 获取地表微小形变的时序短基线算法有机融合，将两者纳入到一个模型中统一计算，利用 两种方法的互补性，通过利用真实相位差建立短基线集处理模型解决了两者无法统一建模 计算和传统时序InSAR方法无法正确获取矿区大形变梯度下的地表沉降问题，改变了原有 短基线集算法中只能解算地表微小形变的方式，现有时序偏移量方法不能有效解算地形误 差的问题也得到了的解决。方法简单，监测精度高、范围大，费用低，对指导矿区生产、预警 地质灾害、治理生态环境、维护矿群关系等具有重要的实际意义和应用价值。
【附图说明】
[0029] 图1为本发明实施例一种融合像元偏移跟踪和短基线集的矿区地表形变解算方法 流程图
[0030] 图2为本发明实施例采用相干性阔值选取的高相干点分布图 [0031 ]图3为本发明实施例解算的相干点处地表年形变速率图
[0032] 图4为本发明实施例估计的相干点处DEM误差分布图
【具体实施方式】
[0033] 下面将结合图和具体实施过程对本发明做进一步详细说明：
[0034] 如图1所示，本发明的融合像元偏移跟踪和短基线集的矿区地表形变解算方法，包 括如下步骤：
[0035] 步骤1:矿区地表形变区域预估：
[0036] 本实例选取被测因煤炭开采导致的矿区，采用的是裁切后的利用卫星获取选取矿 区的11景Im分辨率的Terrasar-X卫星影像，范围为2000 X 1600像元。利用DInSAR或开采沉 陷预计得到的地表沉降结果，将SAR影像区域中相邻像元形变差大于SAR影像波长一半的下 沉盆地初步划分成大形变区域，将像元形变量小于10mm设定为无形变区域，其余区域设定 为小形变区域；
[0037] 步骤2:选取时间序列上相干图中的高相干点：
[003引设定SAR影像的时间和空间基线阔值，选取小于阔值的SAR影像进行干设处理，再 依据相干性获取整个或大部分时间序列上具有高相干性的点位;设定时间和空间基线阔值 分别为11天和160m，选取小于阔值的SAR影像进行干设处理，W相干性阔值为0.3选取时间 序列上的高相干性点位，本发明实施例采用相干性阔值选取的高相干点分布图如图2所示， 共选出885997个点，中部大形变区域内因失相干影响，只选出了部分点位；
[0039] 步骤3:利用时序偏移量跟踪算法恢复大形变区域高相干点的相位整周数N:
[0040] 在划分出的无形变区域内选取相干性大于0.3的高相干点计算整体偏移量，基于 最小二乘原理建立偏移量的双线性多项式内插函数，将形变区域内相干点的影像坐标带入 内插函数计算形变区域内相干点的整体偏移量；
[0041] 利用小窗口最大频谱配准对大形变区域的高相干点进行精确配准，获取局部区域偏移 量;从局部精确偏移量中去除整体偏移量，并去除利用公式:<斯'啤I aW W I 3:W //)2 得到的未下沉区域因地形起伏引起的残差，从而获取最终的视线向像元偏移量offset,式 中，offsettopo为地形起伏引起的偏移量误差;曰日，曰1，恥为最小二乘的拟合系数;Η为对应位 置的高程;77为研究区域的平均高程；
[0042] 由于利用上述像元偏移跟踪算法求解的偏移量仅用到SAR影像幅度信息，不能同 传统短基线技术得到的差分相位进行最小二乘建模。利用公式:N = ceil(2offset . Δ/λ) 得到大形变区域高相干点的相位整周数Ν，式中，ceil()为向下取整运算，Δ为斜距向像元 尺寸，λ为SAR影像波长。由相位整周数Ν和大形变区域内高相干点的原差分缠绕相位联合恢 复大运些点的真实相位差δΦ6 = 2侃+Sd)diff，式中，Sd)diff为缠绕的差分相位；
[0043] 步骤4:大形变区域和小形变区域的高相干点联合建立短基线集解算模型：
[0044] 在随机两个时间点ti，t2时刻分别获取了两幅被测矿区的SAR影像，判读相干点位 于SAR影像的小形变和无形变区域时，通过传统最小费用流方法解缠差分干设图中高相干 点（l，m)的相位δΦ3相位;判读相干点位于大形变区域时，相位解缠结果则不可靠，需利用 步骤3得到的5〇b构建解缠相位矩阵：δΦ = [S0sS0b]，然后直接利用公式：
i立模型，改变了原有短基线集算 法中只能解算地表微小形变的方式；
[0045] 式中：Ahi,m为相干点（l，m)的高程误差;V为相干点的下沉速率;W为残余相位;λ为 雷达波长;为两幅SAR影像相干点处的垂直基线；0(l，m)为相干点处的雷达入射角； ri,m为卫星距离点目标（1，m)的斜距；
[0046] 步骤5,矿区地表高程误差及形变速率的解算
[0047] 利用步骤4中的短基线集处理模型，采用最小二乘解算矿区高相干点处地表高程 误差及形变速率
其中P为权阵，可取单位阵。
[004引本发明实施例解算的相干点处地表年形变速率如图3所示，从中可W看出，每年下 沉大于1000mm的部分地表点已经可W得到，解决了短基线集技术无法监测地表大变形的问 题。图4为本发明实施例估计的相干点处DEM误差分布图。因此，本发明解决了基于SAR幅度 信息的时序像元偏移跟踪算法和基于差分相位的短基线集技术难W联合解算地表形变速 率和高程误差的问题，为矿区地表沉降时序监测提供了新的处理方法。
【主权项】
1. 一种融合像元偏移跟踪和短基线集的矿区地表形变解算方法，其特征在于步骤如 下： 步骤1:预估矿区地表形变区域： 选取被测因煤炭开采导致的矿区，利用卫星雷达获取选取矿区的SAR影像，利用DInSAR 或开采沉陷预计得到的地表沉降结果，将SAR影像区域中相邻像元形变差大于SAR影像波长 一半的下沉盆地初步划分成大形变区域，将像元形变量小于l〇mm设定为无形变区域，其余 区域设定为小形变区域； 步骤2:选取时间序列上相干图中的高相干点： 设定SAR影像的时间和空间基线阈值，选取小于阈值的SAR影像进行干涉处理，再依据 相干性获取整个或大部分时间序列上具有高相干性的点位； 步骤3:利用时序偏移量跟踪算法恢复大形变区域高相干点的相位整周数N: 在划分出的无形变区域内选取相干性大于0.3的高相干点计算整体偏移量，利用最小 二乘原理建立偏移量的双线性多项式内插函数，将形变区域内相干点的影像坐标带入内插 函数计算形变区域内相干点的整体偏移量； 利用小窗口最大频谱配准对大形变区域的高相干点进行精确配准，获取局部区域偏移 量； 从局部精确偏移量中去除整体偏移量，并去除未下沉区域因地形起伏引起的残差，从 而获取最终的视线向像元偏移量offset; 利用视线向像元偏移量off set得到大形变区域高相干点的相位整周数N， 利用相位整周数N和大形变区域内高相干点的原差分缠绕相位联合恢复大形变区域高 相干点的真实相位差 步骤4:利用大形变区域和小形变区域的高相干点联合建立短基线集解算模型： 在随机时间点t、随机时间点t2分别采集了两幅被测矿区的SAR影像信息，判断步骤2选 取的相干点位于SAR影像的小形变和无形变区域时，通过传统最小费用流方法解缠差分干 涉图中高相干点（1，m)的真实相位差δ Φ s; 判断高相干点（l，m)位于大形变区域时，不能采用传统最小费用流方法解缠差分干涉 图中高相干点（1，m)真实相位差，此时则由步骤3获取高相干点（1，m)的真实相位差δ Φ b; 最终，利用S Φ s、S Φ b联合构建短基线集的解算模型； 步骤5,解算矿区地表高程误差及形变速率： 利用短基线集处理模型，采用最小二乘解算矿区高相干点处地表高程误差及形变速 率，在获取矿区地表形变后，通过反演矿区概率积分法预计参数，最终获取矿区地表沉陷预 计及分析模型，为矿区沉陷工程应用提供基础数据。2. 根据权利要求1所述的融合像元偏移跟踪和短基线集的矿区地表形变解算方法，其特征 在于:所述未下沉区域因地形起伏引起的残差利用公式:喻… 计算得到；式中，off sett%。为地形起伏引起的偏移量误差;a〇，ai，a2为最小二乘的拟合系 数;Η为对应位置的高程;互为研究区域的平均高程。3. 根据权利要求1所述的融合像元偏移跟踪和短基线集的矿区地表形变解算方法，其 特征在于:利用公式:N=ceil(2offset · △ /λ)，计算得到大形变区域高相干点的相位整周 数N，式中，ceil〇为向下取整运算，Λ为斜距向像元尺寸，λ为SAR影像波长;从而克服像元 偏移跟踪算法求解的偏移量仅用到SAR影像幅度信息，不能同传统短基线技术得到的差分 相位进行最小二乘建模的问题。4. 根据权利要求1所述的融合像元偏移跟踪和短基线集的矿区地表形变解算方法，其 特征在于：利用公式：δ Φb = 2Νπ+δ φ diff，计算大形变区域高相干点的真实相位差，式中，δ Φ diff为缠绕的差分相位。5. 根据权利要求1所述的融合像元偏移跟踪和短基线集的矿区地表形变解算方法，其特征在 于:利用真实相位差S Φ b建立短基线集处理模型A方法为:首先利用真实相位差δ Φ b构建解缠相 位矩阵δΦ = [δΦ3 δΦ」，再利用公式：建立短基线集处理模型;式中：Ain,?为相干点（l，m)的高程误差;V为相干点的下沉速率;W 为残余相位;λ为卫星雷达波长;为两幅SAR影像相干点处的垂直基线；0(l，m)为相干 点处的雷达入射角;n, m为卫星距离点目标(1，m)的斜距。6. 根据权利要求1所述的融合像元偏移跟踪和短基线集的矿区地表形变解算方法，其 特征在于:局部区域偏移量还可以通过最小二乘匹配方法对大形变区域的高相干点进行精 确配准获取。7. 根据权利要求1或6所述的融合像元偏移跟踪和短基线集的矿区地表形变解算方法，其特 - Λ - 征在于:所述矿区高相干点处地表高程误差及形变速率利用公式： _ V _ ; 计算，式中Ρ为权阵，可取单位阵。
【文档编号】G01S13/90GK106066478SQ201610362009
【公开日】2016年11月2日
【申请日】2016年5月27日 公开号201610362009.6, CN 106066478 A, CN 106066478A, CN 201610362009, CN-A-106066478, CN106066478 A, CN106066478A, CN201610362009, CN201610362009.6
【发明人】范洪冬, 杜森, 黄继磊, 闫世勇, 邓喀中, 汪云甲
【申请人】中国矿业大学