InSAR和GNSS融合的露天矿边坡形变测量方法

摘要

本发明涉及一种InSAR与GNSS融合的露天矿边坡形变测量方法，包括InSAR和GNSS两种技术手段，其特征在于，包括以下步骤：1、利用InSAR监测方法获取边坡的初步形变监测结果；2、在在过形变中心的剖面位置布设GNSS在线监测点；3、利用最小二乘迭代法获取修正后的瞬时形变场；4、利用自适应滤波方法、卡尔曼滤波器方程组和克里金插值方法计算，获取时空域全覆盖的形变监测结果。本发明的优点是：通过融合InSAR和GNSS在空间分辨率和时间分辨率方面的优势进行强强互补，既可实现露天矿边坡重点区域的精确点式实时监测，同时也可实现全矿区的面式综合监测，为露天矿边坡形变实时测量预警提供支撑。

说明书

技术领域

本发明属于边坡工程变形监测技术领域，具体涉及一种InSAR与GNSS融合的露天矿边坡形变测量方法。

背景技术

随着我国经济的快速发展，露天矿的开采范围和开采深度在不断的增加，从而形成了众多大型的高陡边坡。一方面，高陡边坡极大地改变了局部区域的应力条件，进而影响露天矿周边的建构筑物；另一方面，由于高陡边坡很容易发生滑坡等事故，对矿区安全生产造成了很大的威胁。目前，传统的露天矿边坡变形监测技术采用的是GNSS监测手段。由于GNSS定位技术具有精度高、无需通视等特点，在露天矿山边坡稳定性监测工作中可综合运用GNSS定位系统与数据无线远程传输、监测数据自动分析处理等技术，构建边坡稳定性GNSS自动化监测系统。虽然GNSS监测手段和方法已被广泛应用，但其在应用过程中暴露出的不足显而易见，归结为：1)观测装置及供电设施的保存与维护比较困难；2)获取数据为离散点形变信息，空间域覆盖度严重不足，难以反映连续形变规律。

空间遥感技术即合成孔径雷达干涉测量(Synthetic Aperture RadarInterferometry)，简称InSAR，通过对覆盖同一地区的一定时间范围内的多幅雷达图像的联合处理可以提取出相位差图，进而可提取该地区相应时间段内的地表形变信息，观测精度可以达到cm级甚至mm级。它具有比GNSS高得多的采样密度(约1-20m)、空间覆盖度好、非接触性和无需建立地面接收站等优点，是一种极具潜力的空间对地观测新技术。因此，利用InSAR技术开展露天矿形变监测具有重要的应用价值。但InSAR技术使用的数据目前主要来源于商业星载雷达卫星，而雷达卫星运行的重访周期从11天到46天左右，较低的时间域覆盖成为制约InSAR应用的主要因素。

InSAR与GNSS技术在地表形变测量应用方面各有优势和不足，具体情况归纳如下：1)两种技术的空间分辨率：在空间范围上，GNSS的监测范围仅仅局限于一定区域的点式测量，而利用InSAR可以监测大范围的滑坡变形，采用时序InSAR技术对地观测的精度可达到mm级，且雷达差分干涉测量所得图像是连续覆盖的。由此得到的地面形变也是连续覆盖的，这对分析露天矿边坡滑坡形变分布及发展规律是非常有用的。而GNSS采集数据的空间分辨率则远不如遥感，GNSS连续运行监测站点之间的间隔较远且不易维护，而且需要事先建立监测网，会受到地理环境和运作成本等因素的限制。2)两种技术的时间分辨率：GNSS可以在很短的时间间隔采集数据(例如1-30秒)，如果建立了GNSS连续运行监测网，更可以提供连续的、区域性的大气层数据。GNSS网已经达到了非常高的时间分辨率。而InSAR目前主要的商业数据还是来源于星载SAR，雷达卫星运行的重访周期从11天到46天左右。因此，InSAR与GNSS技术在地表形变监测应用方面具有较强的互补性，InSAR和GNSS的技术互补性包括以下几个方面：

1)GNSS定位精度高，但是其空间覆盖度较低，GNSS监测点之间相隔很远，而InSAR提供的是整个区域面上的连续信息；

2)GNSS可提供时间分辨率很高的观测数据(采样率为30秒乃至20Hz)，GNSS允许长时间连续观测，而SAR卫星通常重访周期较长(11-46天)，很难提供高时间分辨率的监测结果。

3)InSAR对高程变化信息特别敏感，利用时序InSAR进行形变监测的精度可达到厘米级甚至毫米级，GNSS对高程变化信息不敏感，单历元GNSS获取的高程变化远达不到这一精度；

综上所述，如何利用InSAR技术与GNSS技术的互补性，将二者紧密的结合在一起，实现优势互补，用于满足我国露天矿边坡监测的需求，已经构成了一个重要的研究方向，具有重要的实际意义。许多学者开展了这方面的研究，如独知行等[4]将GPS和InSAR结果综合应用于地下开采导致的沉降监测中；剧成宇[5]等将GPS应用修正InSAR监测中的大气延迟误差以提高InSAR监测结果的精度等。现有研究以GPS监测结果验证InSAR监测结果或利用GPS修正InSAR监测误差居多，而考虑露天矿高陡边坡特性，基于InSAR和GNSS进行时间域与空间域融合的研究较少。本专利提出面向高陡边坡的GNSS高时间分辨率和高平面位置精度与InSAR高空间分辨率和高垂向监测精度的融合方法，实现露天矿边坡时空域全覆盖式监测，对推动我国边坡监测技术发展、防止露天矿高陡边坡地质灾害发生、保证矿山生产安全，具有重要的理论与实际意义。

发明内容

为克服现有技术存在的缺陷，本发明的目的在于提供一种InSAR与GNSS融合的露天矿边坡测量方法，通过对露天矿边坡形变采用InSAR融合GNSS手段进行监测，并对所获监测数据进行优势互补的联合处理和分析，实现了露天矿边坡时空域全覆盖式监测，监测结果更加清晰准确，对保障矿区安全具有重大意义。

本发明的目的是通过以下技术方案实现的：

本发明是一种InSAR与GNSS融合的露天矿边坡形变测量方法，包括InSAR和GNSS两种技术手段，其特征在于，包括以下步骤：

步骤1、根据边坡的坡度和坡向信息，选择具有合适入射角的升降轨SAR数据，利用InSAR监测方法获取边坡的初步形变监测结果，所述的入射角为即升轨观测时地形坡度小于雷达入射角或降轨观测时地形坡度在雷达入射的背坡面的角度小于雷达入射角的余角；

步骤2、以InSAR监测结果为依据，在过形变中心的剖面位置布设GNSS在线监测点，开展GNSS在线监测；

步骤3、以InSAR监测的瞬时形变场为基础，联合GNSS监测结果利用最小二乘迭代法获取修正后的瞬时形变场；

步骤4、以GNSS观测获取的高时间分辨率形变时间序列数据为基础，利用自适应滤波方法，对每一GNSS测站建立InSAR视线方向形变量的动态模型，依据动态模型和两种数据联合处理得到的修正后的瞬时形变场，应用卡尔曼滤波器方程组和克里金插值方法计算InSAR图像上其他像元的形变时间序列，获取时空域全覆盖的形变监测结果。

所述的最小二乘迭代法为基于最小二乘方法，将所有的SAR影像组合成若干个集合，基于最小二乘方法得到每个小集合的地表形变序列，利用奇异值分解即SVD方法将多个小基线联合迭代求解，获取修正后的瞬时形变场。

所述的卡尔曼滤波器方程组为引入地面控制点，递归算法推进，对卫星轨道、卫星姿态以及影像中心点的大地定位参数进行估值，从而实现大地精校正处理。

所述的克里金插值方法为根据未知形变点的有限领域内的若干已知形变点数据，考虑点位的形状、大小和空间方位，与未知形变点的相互空间关系，以及变异函数提供的结构信息之后，对未知形变点进行的一种线性无偏最优估计。

下面对InSAR原理和InSAR变形测量方法做一简要介绍。

InSAR原理

InSAR技术是以同一地物的两幅SAR图像为基本处理数据，求得两幅SAR图像相位差，进而获取干涉图像，经相位解缠，最终从干涉条纹中获取高程数据或形变信息。SAR影像既记录了地面目标的辐射强度信息，也记录了与斜距有关的相位信息。由于雷达卫星两次过境时轨道不可能完全重复，因此导致雷达卫星两次获取SAR数据的轨道位置和视角存在差别。图1给出了雷达卫星以两次稍有差异的视角对地表进行观测时卫星轨道位置与地面目标的相对几何关系。

其中，S

两幅影像共轭相乘后：

式中，A

其中，△R为两次观测雷达波单程传播距离差。

式中，φ表示实际的干涉相位，通常利用主影像和从影像进行复共轭相乘得到一幅干涉复图像，再由该复图像的实部和虚部之比的反正切得到实际干涉相位在[-π,π)的相位主值，即缠绕相位，而为了还原实际的干涉相位还必须对缠绕相位进行相位解缠。

InSAR变形监测方法

InSAR变形监测方法包括以下步骤：

1、影像配准和重采样

所述的遥感影像配准是将同一目标区域的两幅来自不同传感器或不同时相的影像在空间位置上最佳的套和起来，是制约DInSAR数据处理自动化的关键步骤。

InSAR影像配准采用影像自动配准，首先在参考影像和待配准影像之间找到足够的同名点，由这些同名点作为控制点确定影像之间的函数转换模型及其转换参数，然后根据确定的函数转换模型将待配准影像应用最邻近内插法进行重采样。

所述的重采样是根据目标像元与源像元的比值，取源像元相对位置的像素点作为目标像素点的值。

2、干涉图和相干系数图生成

对雷达影像进行精确配准后，就可计算出每一同名点上的相位差，生成干涉条纹图。

主副影像配准后，将主副影像对应像元的复数进行共轭相乘即可得到干涉图；干涉图本身也是复数数据，其复角即为干涉相位：

相干系数是衡量干涉测量精度高低的一个重要标准，其计算公式为：

式中：C

相干系数越大，信噪比越高，干涉测量的精度越高。

3、形变相位分离

形变相位包括平地相位和地形相位，干涉相位还可以表示为：

φ＝φ

3.1平地相位和地形相位的去除

在干涉相位中，参考椭球面的相位贡献表现为平行的条纹，垂直基线越长，干涉条纹越密集；而地形坡度越大，地形引起的干涉条纹越密集，地形变化越大，干涉条纹曲率变化越明显。

3.2其他噪声的去除

理想的干涉条纹图应该具有连续性和周期性等特点，但实际上真是的干涉图往往出现周期性不明显，连续性中断、干涉条纹不明晰等现象，主要是由于大气延迟噪声、算法模型局限性等的影响。因此，要想得到比较精确的地表形变相位必须削弱噪声相位的影响。目前，比较常用的方法有滤波法，雷达干涉数据处理可以根据噪声的特征采用相应的滤波方法。

4、相位解缠

在干涉处理中智能获得[-π,π]的缠绕相位，为了得到与地形高程直接关联的真实相位，必须通过解缠算法来恢复实际相位，相位解缠的准确性直接影响地形监测结果。

φ＝φ

式中，φ

目前，所采用的相位解缠算法一般是基于路径跟踪和最小二乘原理。所谓路径跟踪是通过像元到像元的局部运算来解缠，最小二乘是通过解缠前后的相位梯度差整体最小进行求解，实际上是基于同一个基本的数学模型。即L

式中，

5、地理编码

地理编码是将前期处理得到的雷达坐标系下的结果转化到地理坐标系下的过程。

地球方程模型如下：

式中R

进行地理编码的主要原因是：

①纠正由地形起伏造成的SAR影像几何变形；

②与已有的地理坐标系下的其他资料进行比较和融合

与现有技术相比，本发明的优点是：

1、利用雷达遥感的全天候作业以及穿透性特征，实现露天矿复杂环境下的边坡形变全覆盖监测。

2、通过InSAR和GNSS在时间分辨率和空间分辨率方面的强强互补，实现露天矿边坡不同形变速率(包括缓慢形变和快速形变)的联合有效监测。

3、通过轨道精炼和重去平方法去除监测区的大气噪声影响，提高InSAR的监测精度和可靠性。

4、通过组合不同时空基线的像对进行联合处理以提高影像的利用率，实现高精度的全时序形变监测。

通过融合InSAR和GNSS两种方法，既可实现露天矿边坡重点区域的精确点式实时测量，同时也可实现全矿区的面式综合监测，为露天矿边坡形变实时监测预警提供支撑。

附图说明

图1为本发明的InSAR测量原理几何示意图。

图2鞍千矿区雷达遥感影像。

图3 Sentinel干涉像对时空基线连接图。

图4配准后的影像强度图。

图5去平并滤波后的干涉图。

图6矿区InSAR与GNSS结果对比。

具体实施方式

下面结合附图对本发明作进一步说明。

如图1，本发明的一种InSAR和GNSS融合的露天矿边坡形变测量方法，包括InSAR和GNSS两种技术手段，其特征在于，包括以下步骤：

步骤1、根据边坡的坡度和坡向信息，选择具有合适入射角的升降轨SAR数据，利用InSAR监测方法获取边坡的初步形变监测结果，所述的入射角为即升轨观测时地形坡度小于雷达入射角或降轨观测时地形坡度在雷达入射的背坡面的角度小于雷达入射角的余角；

步骤2、以InSAR监测结果为依据，在过形变中心的剖面位置布设GNSS在线监测点，开展GNSS在线监测；

步骤3、以InSAR监测的瞬时形变场为基础，联合GNSS监测结果利用最小二乘迭代法获取修正后的瞬时形变场；

步骤4、以GNSS观测获取的高时间分辨率形变时间序列数据为基础，对每一GNSS测站建立InSAR视线方向形变量的动态模型，依据动态模型和两种数据联合处理得到的修正后的瞬时形变场，应用卡尔曼滤波器方程组和克里金插值方法计算InSAR图像上其他像元的形变时间序列，获取时空域全覆盖的形变监测结果。

具体详细过程如下：

1)InSAR数据及处理

如图2-图6所示，持续、稳定、可靠的合成孔径雷达(SAR)影像数据源是进行连续监测的前提。目前主要能够提供SAR数据的雷达传感器如下，可提供用于InSAR监测的单视复数(SLC)影像数据。可用的几种新型SAR传感器相继发射升空，包括RadarSAT-2、Cosmo-SkyMed、TerraSAR-X、ALOS-2、Sentinel-1A/1B以及PAZ。其中，采用C波段的Sentinel-1A/1B继承了ERS的优点，由于传感器具有非常高的定位精度，从而使得在较长时间内都能保持很好的相干性。

本次研究采用的是从2017年11月8日至2018年10月22日之间由欧空局的Sentinel卫星获取的26景雷达影像，其中2018年1月31日获取的鞍千矿区的雷达强度影像如图2所示。Sentinel-1卫星在近极地太阳同步轨道上运行，轨道高度约700km，重访周期为12d。为获取较好的干涉产品，Sentinel-1采用了严格的轨道控制技术，沿既定轨道运行时的卫星位置非常精确。Sentinel-1A于2014年4月发射升空，可持续提供重访周期为12天的覆盖全球的开放数据。其中，TOPS成像技术获得的IW(Interferometric Wide Swath)模式影像覆盖范围为250km×170km，距离向分辨率为2.3m，方位向分辨率为14.0m，空前的成图能力足以满足大区域多时相InSAR处理和应用分析。Sentine-1携带一个单独的C波段垂直极化辐射计，波长在3.8-7.5cm之间，中心频率为5.405GHz，带宽为0.100MHz，脉冲宽度为5～10μs，接收机噪声系数-22dB，脉冲重复频率1000～3000Hz，数据存储容量为1410GB，X波段下行链路容量为520MB/s，入射角在20°～45°之间。C波段天线尺寸为12.3m×0.821m，质量为800kg(占发射质量的40％)，支持双极化(HH，VV-VH)运行。数据压缩采用熵约束分块自适应量化计算法和弹性动力分块自适应量化算法。

2)影像对组合

对鞍千矿区2017年11月至2018年10月间所获取的全部Sentinel影像进行像对组合，针对Sentinel影像数据的特点，以及本项目的实际需求，选取20180123获取的影像作为影像序列的主影像，以短基线集的原则为基础，选定空间垂直基线阈值为临界基线的15％，时间基线阈值为72天，在这个过程中，可以对基线组合进行修改编辑。时空基线连接图如图3所示，构成干涉像对的主辅影像和空间基线与时间基线数值见表1。

表1干涉像对及基线值部分列表

3)影像配准和重采样

在主影像被精确校正配准到DEM的情况下，我们在主影像上选取合适的地面控制点(GCP)形成一个xml文件，同时利用抚顺地区的DEM和基线组合生成的辅助文件作为输入文件。首先是配准，分为三步进行，粗配准时我们采用窗口大小距离向和方位向为1024*4096，距离向和视向窗口数分别为6和9，在此基础上计算影像的交互相干性。象元级配准，我们采用窗口大小距离向和方位向为256*512，粗配准相干性小于0.25的部分无法进行象元级配准。最后是亚象元级配准，这一步可以实现0.1像素的精度。距离向和方位向窗口数为20和40，窗口大小为32*32，过采样参数为16，相干性过采样参数为2，信噪比阈值默认为3.2，计算相关系数，计算相关参数时要充分考虑各方面影响参数计算的因素，因此，平时要经常积累积算这些系数所涉及的内容，以及准确计算这些参数的方法，大量的数据处理以及反复的比较验证等方法对配准效果是非常有用的，争取做到能进行精密配准。配准后的强度如图4所示。

4)干涉图生成与去平

在影像重采样完成后，两幅影像共轭生成干涉图，设置多视数方位向与距离向分别为1和4，得到部分干涉图如下。很显然干涉条纹比较复杂，需要进行去平地操作。

去平地时，采用的是基于DEM和主影像上控制点的去除参考相位的方法，首先将抚顺地区DEM数据转换到SAR影像坐标系并模拟出DEM对应的雷达幅度图；然后对幅度图和高程数据进行克里金插值，使其与干涉图具有相同的分辨率；最后将模拟的雷达幅度图与主影像配准，根据配准参数对SAR坐标系下的DEM数据进行重采样，将高程信息转换成地形相位，影像去除地形相位得到去平后的干涉图。

本项目选取达尔曼滤波器方程组，设定距离向和视向的窗口大小均为5，同时对光谱频移进行滤波，计算得到去平后和滤波后的_fint文件。去平并滤波后的干涉图如图5所示。

在干涉图生成的同时，为了评价干涉图的质量，得到了相干系数图。针对本项目的目标以及影像和测区的实际情况，依据如下原则确定相干性临界值：

(1)每一景影像的相干性分布；

(2)每一个像素影像序列相干性统计；

(3)参考国内外沉降监测的成果，经多次反复试验。

最终相干性临界值设为0.18，相干系数图反映了干涉图质量的文件_cc文件。

5)相位解缠

在得到的去平滤波干涉图_fint文件和相干系数图_cc文件基础上进行相位解缠，可选用的方法包括枝切法、最小二乘法和最小费用流法等多种方法，主要是为了恢复实际相位。本项目采用的解缠方法是最小费用流法，该算法先以轨道参数计算各点位的大约位置，再配合点位的相干性，以相干性高的点逐步推求周围的点，将错误点的影响降到最低程度，且可以自动计算完成整副影像的相位值，相位解缠的象元相干阈值设为0.18，结果得到解缠后的相位。

6)轨道精炼与重去平

精炼与重去平，对卫星轨道和位相偏移进行纠正，进行轨道精炼和位相偏移的计算，消除可能的斜坡位相。利用从主影像上选取的含GCP信息的文件和解缠后相位文件进行精炼，由于控制点数量足够，我们选择在控制点基础上进行精炼。重去平时控制点的高程和干涉相位上的高程差的平方根是被限定的，这样可以保证重去平的质量。重去平后还需两个重要的步骤：

第一步：利用重去平后的干涉图、辅助文件，对形变和高程变化相关信息进行估计，选取三次方程模型对残余高程，位移速率，加速度和加速度的速率进行估计，并移除这些低通部分，这些部分在干涉图相位解缠前移除。得到结果如下：高程和精确高程文件和斜距的DEM文件、速率文件和平均速率文件、精确速率文件和加速度相关信息的文件。

第二步：移除高通变量，主要对大气影响部分进行估计与消除，大气中的低通部分与空间有关，用一个正方形窗口加以滤除。大气中的高通部分与时间有关，用一个时间窗口滤除，得到与时间相对应的沉降情况有关的文件。该步骤类似于常规D-InSAR中的将相位转换到形变的步骤。利用时间序列分析工具可以对提取出的形变信息进行分析。经上述处理后得到的未地理编码前沉降图。

6)地理编码

地理编码，将雷达坐标系转换到用户自定义坐标系，即将雷达视线向的形变量转换到规定的坐标系中，可以转换到垂直方向的形变量，也可以转换到带有坡度的形变情况。我们设定x和y方向上的栅格大小为20*20，这一步将之前的速率文件、高程文件、加速度文件和形变文件转换到相应坐标系中，从而得到对应的一系列编码后的结果。

具体露天矿InSAR与GNSS监测分析结果为：由于目前InSAR和GNSS监测数据在时间上没有重叠，所以暂时无法进行定量对比，但从已有的GNSS结果和InSAR形变结果来看，在变形区域上具有非常好的吻合性。InSAR监测结果和GNSS的监测结果均表明测区存在较大变形，而且InSAR监测结果表明变形的范围不仅限于边邦陡坡，在其后面仍然存在一个较大范围的变形区，建议进一步加强GNSS在线监测以及后续InSAR监测。在边坡容易形变位置合理布设GNSS在线监测点，开展GNSS在线监测，以InSAR监测的瞬时形变场为基础，结合GNSS监测结果利用最小二乘迭代法获取修正后的瞬时形变场以GNSS观测获取的高时间分辨率形变时间序列数据为基础，对每一GNSS测站建立InSAR视线方向形变量的动态模型，依据动态模型和两种数据联合处理得到的修正后的瞬时形变场，应用卡尔曼滤波器方程组和克里金插值方法计算InSAR图像上其他像元的形变时间序列，获取时空域全覆盖的形变测量结果。