地基干涉雷达实时处理软件设计与开发*

李一帆 杜柏林 付钰涵 杨红磊

（1.中国地质大学（北京）土地科学技术学院；2.中国矿业大学（北京）地球科学与测绘学院）

独特的地域和地质环境使我国成为世界上地质灾害最严重、受威胁人口最多的国家之一，其中滑坡是地质灾害中发生频率和危害最大的一种，特别是近年来由于人类活动及极端气候条件等因素，滑坡灾害发生的频率越来越高［1］。我国已将地质灾害防治视作生态文明和美丽中国建设的重要内容［2-3］，因此加强滑坡灾害的防灾减灾研究工作十分必要。

滑坡的发生、发展和演化过程伴随着大量的宏观可观测数据，如滑坡形变、物理场变化和地下水位变化等，在众多可观测的数据或物理量中，形变是表现滑坡变化最显著、最直接的变量，其形变趋势又与滑坡所处阶段存在良好的映射或函数关系，因此研究滑坡变形监测技术对于滑坡预报和边稳定性分析具有重要意义。现有的基于点位监测技术，对于高陡顺层边坡而言，存在监测点布设困难的问题，即使是布设了密集的监测设备也很难获取边坡的空间连续和关键部位全过程形变信息，而这些信息对于高风险的滑坡来讲又非常关键，近些年发展起来的地基雷达干涉测量技术（Ground-based Interferometry Synthetic Aperture Radar，GBInSAR）［4］不受雾、雨、雪、粉尘等影响，可以实现全天时、全天候、高精度和高时空分辨率的观测，是对局部区域形变非接触监测的一种新技术手段，获取的形变结果更有助于滑坡灾害力学模型的理解与稳定性评价。当前存在各种体制的地基干涉雷达设备，其中GAMMA Remote Sensing公司的研制的GPRI-II设备与同类设备相比，能够实现长距离360°监测，具有监测范围大、精度高的优势，同时该设备与专业数据处理软件GAMMA相适配，功能强大，市场占有率较高［5］。

与GPRI-II 设备配套的GAMMA 软件模块内的程序需在控制台上运行，无可视化界面，虽然程序的组合十分方便，但需要用户具有一定的编程基础来实现批处理，多采用事后的方式对数据进行处理，这种处理方式具有延后性，不利于滑坡灾害的监测预警。市面上具有可视化界面且能够实现实时处理的软件多于地基雷达硬件匹配，由于数据格式不兼容的原因，无法实现对GPRI-II 设备数据的处理。因此迫切需要在GAMMA 软件的基础上进行二次开发，编制一套具有可视化界面的地基干涉雷达实时处理软件。为解决上述问题，本文提出一种地基干涉雷达数据实时处理流程，在GAMMA 软件基础上，采用Python语言设计与开发地基干涉雷达实时处理软件，对于滑坡灾害监测具有重要的实际意义。

1 原 理

1.1 DInSAR技术

DInSAR 技术即差分干涉测量技术，与星载SAR不同，地基SAR观测方程如式（1）所示［5］：

式中，d 为干涉对时间间隔内的形变量；λ 为雷达波长；δφatmo是大气效应引起的相位延迟误差；δφnoise是噪声引起的相位误差；k为整周数；ΔφPP'为干涉相位，对应的值为相位主值，介于( -π,π ]，需要进行解缠才能获取到真实的形变相位。

地基干涉雷达大多采用Ku波段，为短波范畴，易受到大气延迟误差的影响，研究表明温度在20℃时，距离雷达1km 处，l%的相对湿度的变化可导致2mm的测量误差。在本文中基于大气同质（大气扰动引起的相位差与雷达和目标距离存在线性关系）的前提，构建二阶多项式拟合大气扰动特征［6］，如式（2）所示。

式中，r为目标点到雷达中心的距离；a、b、c为该方程的系数；φatm为该点解缠后的干涉相位。

利用稳定点扰动值和该点与雷达距离求解系数，最后根据求解的大气延迟相位模拟整个监测场景的大气延迟相位。采用DInSAR 技术获取监测目标形变的过程，包括影像配准、干涉图像生成、相位降噪滤波、相位解缠、去除大气延迟、形变计算、地理编码7 个步骤［7］。由于估算大气延迟相位时采用大气参数同质的假设较为理想，因此不能完整地消除影响，获得形变监测精度有限；对于采样间隔长的干涉，又易受到失相干因素的影响。因此，DInSAR技术不用于变形监测；但是其构建的函数模型可作为时序InSAR技术的基础模型。

1.2 时序InSAR技术

为解决大气延迟与时空失相干问题，在DInSAR技术基础上发展出时序InSAR技术，采用多时相SAR影像，依据每个像元在时间维的相位值质量的稳定性，识别出高质量的点，称之为PS 点，对离散的PS 点构网，建立相邻点监测相位增量，根据相邻点间的形变模式，建立合适的形变模型，进而获取形变速率和累计形变量。相应的技术有PSInSAR、IPTA 和SBAS等［8-9］，其中SBAS 技术以高时空相干性的原则，多参考影像的原则确定干涉对组合。对于地基干涉雷达，由于不存在空间基线，因此在干涉对组合时只考虑时间基线，短时间基线一方面可以减弱大气环境变化的影响，另一方面可以保持干涉对的相干性。大气延迟误差按照同质大气和非同质大气分两步进行剔除，第一步根据式（2）分离出同质大气的延迟相位；第二步根据初始获取的形变分布，确定出稳定的点，根据稳定点的相位，采用空间插值的方法拟合出形变的延迟相位；最终把两部分大气延迟相位累计相加，即为总体大气延迟相位。大气延迟相位剔除后，根据干涉对组合方式确定的观测方程系数举证，采用SVD分解方法计算出累计形变量。

1.3 形变预警原理

斜坡的稳定性与斜坡当前变形阶段联系紧密。松散土质斜坡和具有时效变形特点的岩质斜坡，变形往往呈现出蠕变的特点［10］，其变形过程一般可分为初始变形、等速变形和加速变形3 个阶段［11］，在加速变形阶段又包含初加速、中加速和临滑3 个亚阶段［12］。斜坡一旦进入加速阶段，迟早会发生滑坡，当处于临滑阶段时，斜坡整体失稳即将滑坡［13］。阈值预警法是目前国际上最常用的方法［14］，即通过统计分析等给定预警阈值，当监测数据达到阈值时发出警示信息。由于不同物质组成、不同规模、不同成因类型的滑坡变形阈值差异非常大，不存在统一的阈值［14］。本文使用速率曲线和加速度曲线的双重阈值法，达到阈值进行警报。

1.4 实时处理流程

本文基于DInSAR 技术与时序InSAR 技术提出一种实时处理的流程，如图1所示。基于地基SAR零空间基线和短时间基线的特点［15］，假设短时间内监测目标形变呈现线性趋势来对小基线集方法进行简化，每次选取连续的一定数量（20 景左右）的监测数据组合干涉对，根据时序干涉相干性的平均值和标注差双阈值确定高相干点。对每个干涉对进行大气延迟相位分离，剩下的相位作为观测值，采用最小二乘线性拟合［16］，即可得最终的形变量。

式中，A为时间间隔系数矩阵；b为形变值矩阵。

以一个像元点( i,j )为例，设有n副干涉影像参与解算，实际计算形式如式（4）所示，其中vˉ(i,j)是拟合出的该点的平均形变速率（一阶项），c(i,j)是拟合出的该点的常数项，d(i,j)n表示第n景干涉图像( i,j )点的形变值。

2 地基干涉雷达软件系统开发

地基干涉雷达实时处理软件是一个在GAMMA软件的基础上进行二次开发的软件。该软件使用并整合GAMMA 命令集，以InSAR 数据处理算法为核心，以合理的体系架构模块设计为基础，来为用户提供友好的界面和良好的使用体验。

2.1 开发环境与配置

软件开发采用Python 语言，使用Pyside2 库（Qt for Python）并结合如Numpy 等Python 众多第三方库在GAMMA软件的基础上进行二次开发。

Qt为Python编程环境提供支持，开发可在Python上使用的库PySide2，该库包含绝大多数Qt 的类与函数，用法与Qt 内的类与函数十分接近。安装PySide2及其它第三方库仅需通过python的包安装组件pip下载即可。

GAMMA 脚 本 文 件 使 用 如sh、csh、tcsh、bash、perl、python等多种脚本语言完成编写，要想成功运行GAMMA软件，除要含有GAMMA软件本体外，还应当含有对这些脚本语言的环境支持。由于本文软件开发使用Python语言，在文件打包时便已经包含Python语言环境支持，而且sh、csh、tcsh、bash 等都是Unix 或Linux 系统的命令交互语言，因此，在win 系统上需安装对Unix/Linux 命令及Perl 语言的支持，在Unix 和Linux仅需安装Perl语言支持即可。

同时，GAMMA 公司已经推出Python 模块文件py_gamma.py，在完成环境变量配置并正确导入该文件后，可实现在Python 中运行GAMMA 命令。原始的py_gamma.py文件需要存放在GAMMA 软件根目录下才能使用，这将在软件打包后产生找不到文件的错误。查看相关源代码后发现，该文件在原理上是通过添加环境变量的方式来查找同目录下GAMMA 软件的可执行文件与脚本文件，故可对该部分代码稍作修改，使得在打包后的软件根目录内也能向环境变量添加正确的文件路径。

2.2 软件组成

本文软件模块组成与关系如图2所示。

（1）格式转换模块。本软件不仅支持GPRI-II 设备，还可以将其他品牌的地基干涉雷达数据转换为GAMMA格式，如Fast GBSAR设备等。

（2）差分干涉处理模块。该模块对不同时间获取的两景SAR 影像进行处理，包括配准、干涉、滤波和相位解缠等功能，最终获得两景影像之间的形变图像。

（3）时序处理模块。该模块每次对连续的19 个干涉对间形变图进行处理。该模块流程含有生成累积形变图和形变速率2 个步骤。模块一次运行结束后获得19 个累积形变图和1 个形变速率图。其中19个累计形变图分别记录对应影像与第一景影像之间的累积形变，与之前数据有关联；那一个形变速率图记录该20景影像内的平均速率，与之前数据无关联。

（4）地理编码和成果输出模块。该模块将雷达坐标系下的成果数据转为地理坐标系下的成果，并转换输出可被如ArcGIS、QGIS 和Google Earth 等软件使用的文件。

（5）预警处理模块。该模块对所有累积形变图和形变速率图上的某一特定区域（人为设置的预警区域）进行处理。获得统计曲线并对标记形变速率和加速度高于预警阈值的时间段。

（6）项目文件与日志模块。该模块与软件设置和软件运行信息有关，可分为两部分：项目文件以及日志。项目文件储存软件相关的所有设置，能够保证软件在打开时恢复之前运行进度。日志部分将软件运行时的报错信息与运行信息输出到文件中，能够更好地发现错误的原因。

（7）显示模块。该模块将数据及相关信息显示在软件上，主要为影像的显示和统计曲线的显示两部分。影像的显示包含影像数据转为彩色图像、彩色图像的移动与缩放、彩色图颜色栏的显示等。统计曲线的显示包含曲线的展绘、曲线上某数据点信息的显示、曲线自适应标签等。此外，还有数据信息的显示如影像某一点的数据信息显示等。

（8） 自动监测与警报模块。在结果达到阈值时向外发送声音、短信、邮件等警报信息。

3 案例分析

3.1 数据来源

本文使用河北马兰庄矿区边坡监测数据进行分析，数据获取所使用的设备为中国地质大学（北京）与北京理工大学共同研制的地基MIMO 成像雷达系统，共计382景影像，采样间隔3 min。

3.2 结果分析

3.2.1 软件运行时间

测试软件的电脑配置为Intel Core i7-4810MQ CPU@2.80 GHz，内存32 GB。对382 景SLC 影像数据处理时间进行记录并成图，如图3 所示。可以看出，单张数据处理时间全都在25 s 以内，小于3min 采样间隔，速度较快。该软件处理效率与监测时刻大气环境的复杂性有关，对于多变的大气环境需要多次迭代，耗费时间较长，反之在20s内即可完成计算。

3.2.2 软件运行结果

图4 为获取的累计形变图，形变区域明显，深色部分相对于周围浅色部分有着明显不同，浅色部分数值表示形变值大概在0 mm 附近，而深色部分对应-30 mm 以下，能够直观地看出较大形变区域大概范围和所在的位置。

图5 位单点统计曲线。累积形变曲线前期变化较快，但变化率不断减小，后期趋近于直线，比较平缓；形变速率曲线能够直观地看出该区域前期有滑坡现象但后期逐渐稳定的态势，速率趋势与累积形变趋势对应，算法在时间线上无明显错误。

3.2.3 精度分析

在分析结果精度时，由于第三方同步监测数据的缺失，无法进行外符合精度评定。假设区域稳定点在监测时段内计算得到的累积形变量理想情况下应当为零。因此可以使用监测区域稳定点的累积形变曲线是否稳定，即残差的标准差是否较小，来判断处理结果精度的高低。残差为干涉值减去模型值的剩余部分，模型值由两部分构成，基于距离函数生成的大气模型与对形变进行时间拟合的形变模型。由于地基SAR 系统监测精度为mm 级，可以认为实时处理中数据标准差小于1 mm的点精度尚可。

在稳定点区域4个角落（图4中A、B、C和D）选取5×5大小共100个点进行数据残差的标准差计算，如图6所示。

由图6 可知，残差的标准差小于1 mm 的较高精度的点所占百分比即为77%。与后处理相比，实时处理不能对整体数据在处理时进行参数的调节，较高精度稳定点占比未能达到接近100% 的水准，精度有所损失，但在该精度下能够满足对累积形变曲线走势和滑坡预警的判断。

4 结 论

针对当前GPRI-II 设备配套的GAMMA 软件不具备实时处理及可视化的问题，在滑坡灾害监测预警方面应用存在缺陷，提出一种地基干涉雷达影像数据的实时处理流程，并在GAMMA 软件的基础上，结合监测预警的实际需要，采用Python语言进行设计和二次开发，编制了地基干涉雷达实时处理软件，该软件具有良好的人机交互界面，可以实现地基干涉雷达数据预处理到形变监测产品制作全过程。以河北马兰庄矿区边坡监测数据为例进行分析，实时处理时单景数据处理时间优于25s，小于数据采集时间；监测结果精度优于1mm，满足形变监测预警的需求。该软件的研制能够弥补GPRI-II 在监测应用中的不足，同时为拓展地基干涉雷达技术提供了技术支持。