地质灾害风险调查中遥感解译的应用分析*

乔天荣 崔 剑 马涛峰

（河南省地质研究院）

为深入贯彻落实习近平总书记对地质灾害防治工作的重要指示、批示，切实做好2021年度河南省地质灾害防治工作，河南省地质研究院开展了河南省罗山县1∶5 万地质灾害风险调查评价遥感技术服务工作。主要任务是摸清罗山县域自然灾害风险隐患和抗灾能力，客观认识区域内的自然灾害综合风险水平，为地方政府有效开展自然灾害防治工作、切实保障经济社会可持续发展提供权威的灾害风险信息和科学决策依据［1］。遥感解译是本次工作的主要技术手段。工作内容主要包括收集相关资料，通过开展遥感解译工作，了解该区域地质环境背景条件、灾害点特征、地质灾害分布状况、承灾体类型等；开展地质灾害与孕灾地质条件、承灾体调查；判别地质灾害隐患；总结地质灾害分布规律和分析其成灾模式；开展地质灾害易发性、危险性和风险评价，编制地质灾害风险调查评价相关图件；建立地质灾害风险调查空间数据库；提出地质灾害风险管控对策建议，为防灾减灾管理、国土空间规划和用途管制等提供基础依据。

1 遥感解译的质量控制

本次研究采用的遥感影像数据源为美国陆地资源卫星Landsat-8、GF-1、GF-2、GF-7等数据，结合收集地理图、地质图、地形图等图件及文字资料，采用人机交互的形式，使用遥感信息处理软件ENVI 和地理信息系统软件MapGIS、ArcGIS进行影像信息提取，工作过程中，按照遥感环境地质调查规范制定的解译要求，严格进行质量控制，满足全区1∶1 万遥感解译需求。

遥感解译作为环境地质工作的一种手段，目的是配合面上调查。主要解译地质环境背景，如色彩反映较大的地层、规模较大的地质构造体、灾害体，不同土地类型、水体等；质量控制主要进行数据源的选择、数据预处理、初步判读、遥感解译、目视修改、野外验证等几个关键步骤。

数据源质量控制：在选择原始影像的时候，尽量选择覆盖整个工作区范围、分辨率高、色彩信息丰富、含云量较少的影像，各个时相的季节基本一致。其中选用的多光谱卫星遥感影像，时相均为秋冬季，无云，包括高、中空间分辨率，包含多个可见光—近红外波段，波段数量较多，影像信息量丰富。雷达影像数据选择卫星升轨状态的数据，时相周期为12 d，保证了影像之间的相关性。

（1）数据预处理质量控制。①多光谱遥感数据几何纠正选择多项式纠正法，以1∶5 万的地形图为参照；②控制点尽量选择不易随时间变化、固定地物的交叉点；③拟合误差山区控制在2 个像元，平原控制在1 个像元内；④辐射校正按照影像获取时间；⑤模拟大气状态，保证了象元反射率的可靠性；⑥雷达影像的预处理，以1∶5 万地形图为参考，选择一景影像与其配准，之后以该影像为基准，每一景影像按基准影像进行配准，保证像元误差控制在0.5个像元以内。

（2）初步判读质量控制。初步判读阶段，参考目标地物的影像特征，并结合工作区的实际情况和项目的具体需求，确定合适的地物类别。通过影像初步判读提取典型地物影像图，设计合适的路线进行野外踏勘，踏勘范围包含各类目标信息的特征，以获取整个工作区的地形地貌特征和不同地物的图像，踏勘样点记录坐标信息、地物图片和文字描述资料。

（3）遥感解译质量控制。遥感解译阶段，为了突出相关的专题信息，在信息提取之前对遥感影像进行特定的图像增强处理，主要进行波谱信息增强和空间信息增强。进行地质构造解译时，在整理分析野外踏勘结果的基础上结合影像特征，选取尽可能多而且能代表本区目标信息特征的直接解译标志和间接解译标志，把握目标信息的宏观特征及微观现象，对地质构造信息进行细致的分析及提取。在进行地表形变反演时，对各个影像之间的相关性进行评价，剔除解缠效果较差的影像，选择山脚、高层建筑等相关性高的点作为控制点，控制点拟合误差小于1。

（4）目视修改质量控制。目视修改阶段，以项目成果精度要求为依据，以原始影像为基础，选择合适的窗口大小。对细碎图斑进行合并，并对错分、漏分的图斑进行修改，保证结果与原始数据的一致性。

（5）野外验证质量控制。野外验证阶段，在室内要根据验证目标地物类型选择具有代表性的验证点，设计好合适的路线，野外验证过程中，详细记录验证点地物特征，尽可能在实地进行初步修改。

2 遥感解译的方法

本次工作采用SBAS-InSAR 技术开展罗山县地表形变监测。具体流程见图1。

SBAS（Small Baseline Subsets）是短基线集差分干涉测量技术，主要用于提取大规模、形变速率缓慢的地表形变信息。该技术是将获取的SAR 影像按照子集内的干涉对空间基线距相对紧凑、子集间干涉对空间基线距相对较大的原则进行自由组合，形成若干短基线集合，建立干涉像对连接图，并利用空间基线阀值法选取短基线干涉对，削弱空间失相关的影响［2］。进行地表反演时引入奇异值分解（Singular Value Decomposition，SVD），联合多个小基线集求解，得到形变参数、高程误差在最小范数意义上的最小二乘解，最后估计出非线性形变和大气相位。该方法不仅解决了长空间基线造成的时间不连续问题，还增加了监测的时间采样率，又结合稳定的干涉相位，提高了监测的空间分辨率［3］。

将N+1 景覆盖同一区域的SAR 影像序列，影像获取时间依次为(t0,…tn)，根据空间基线和时间基线的条件组合干涉对，得到M幅差分干涉图，M的数量与影像数量满足以下关系：

经过2个时刻tA和tB（tA＞tB）获取的影像生成的第j幅差分干涉图中像元（r，x）的相位值可以表示为

式中，λ为入射的雷达波长，φ(tA,r,x)和φ(tB,r,x)分别为生成干涉相位的2 个影像相位，d(tA,r,x)和d(tB,r,x)分别为像元于tA，和tB时刻相对于初始t0时刻沿雷达视线方向的形变。

此相位值表示除去DEM 引入的残余高程误差相位、大气延迟相位和噪声误差等因素的形变相位。设d(t0,r,x)=0，则干涉图的相位时间序列可以表示为

SBAS 进行逐像元计算，那么对于SAR 影像中的某个像元来讲，该像元点的形变量可用向量表示为

解缠差分干涉图中相位组成的向量值为

式中，δφi(i=1,···,M)为相对于解缠参考点的相位值。

设参与干涉的主、辅影像对应的时间分别为

将主、辅影像按时间排列，即IMj＞ISj(j=1,…,M)，则差分干涉图中的相位值如下：

可以简化为δφ=Aφ。

矩阵A[M×N]的每一行对应1 个干涉像对，即有[j,IMj]=1 和[j,ISj]=-1(j=1,…,M)，矩阵中其他值为0。假设δφ1=φ4-φ2和δφ2=φ3-φ0，则矩阵为

若所有干涉对属于同一个子基线集，则矩阵为1个N阶矩阵（M≥N），当M＞N时，该方程是1个超定解方程，其最小二乘解可以表示为

当基线集中含多个子集时，矩阵A秩亏，ATA为奇异矩阵。如果有L个不同的子基线集，矩阵A的秩就为N-L+1，δφ=Aφ会产生多组解。利用奇异值分解（SVD）法能简便地解出式δφ=Aφ的唯一解，对矩阵A进行奇异值分解可得

式（10）中U为M×N的正交矩阵；∑的对角元素为奇异值δi(i=1,···,N)；V为M×N正交矩阵。

δφ=Aφ的最小二乘解为

其中，∑-1=diag(1∕σ1,···,1∕σN-L+1,0,···,0)。为了得到符合物理意义的解，将对相位的求解转化为对相位变化速度的求解问题，则相位变化速度如下。

则式δφ=Aφ可以转化为

可简化为

式 中，B为M×N矩阵，矩阵元素B[i,j]=tj+1-tj(ISj+1 ≤j≤IMi,∀i=1,···,M)，其它元素值都为0。

对B进行奇异值分解，就能求出各时间段中的平均速度V，在时间域上通过积分运算即可得该像元的形变信息［4］。

综上所述，SBAS 技术就是利用奇异值分解来解决因时间基线较长而失相干的问题，根据各个集合间基线较大、集合中基线较小的原则，组成若干短基线集合，最后用最小二乘法解出集合内地表形变的方法。

3 关键技术处理分析

3.1 数据导入

数据导入采用精密轨道文件，将使数据的定位更准确，同时能减少处理时由轨道误差引入的相位误差。Sentinel-1的精密轨道文件有2种，一种是在接收到GNSS 数据的21 d 后产生的，精度高于5 cm，覆盖1 d 加1 d 的前1 h 和后1 h，共26 h，为最精确的精密定轨星历数据文件；另一种是在接收到GNSS 数据的3 h 内就可以使用，精度高于10 cm，覆盖1 个卫星的轨道，比较精确的回归轨道数据文件。项目为了得到更精确的监测结果，在Sentinel-1数据导入时，采用精度高于5 cm 的精密定轨星历数据文件，罗山县位置跨两景遥感数据，同时整景影像的处理所占空间较高，为了减少数据量，节省数据的处理时间，需要对遥感数据集进行裁剪合并。

3.2 生成图像连接图

根据集合间SAR 影像的基线较大、集合内基线较小的原则，生成若干小基线集，产生相对于主影像（2017-11-09）的位置图和时间—空间基线图。经过不同参数调试后，确定时间基线为120 d，空间基线阈值为极限基线的2%。46 景数据共生成像对130 对，每个图像的平均连接数5.65，最大连接数6，最小连接数3；平均时间基线绝对值为73.015 4 d，最大时间基线为120 d，最小时间基线值为36 d。平均空间基线绝对值为54.580 6 m，最大绝对基线为169.598 m，最小绝对基线为1.252 93 m。

3.3 滤波与干涉处理

主辅影像基于相位的精确配准是生成干涉图的基本条件，采用相干系数法基于Sentinel-1 精密轨道数据和30 m 的ASTERDEM 数据对所有像素进行逐一配准，为了达到方位向和距离向一致的分辨率，对距离向做4 倍多视处理。通过2 个相位共轭相乘得到具有平地效应的初始干涉图。从初始干涉图中去除平地效应，只保留了因地形和高程起伏所导致的干涉条纹。采用Goldstein 滤波方法对干涉图进行滤波处理，可提高相干性，减少2点噪声，使干涉条纹更加平滑［5］。查看去平和滤波后的干涉图，针对低相干区域较多，增加滤波强度，根据工作区的相干情况，Goldstein Max Alpha的参数设为2.5，Goldstein Min Alpha 为0.3。相位解缠采用MinimumCostFlow 方法，解缠分解等级设置为1，适当减少低相干性区域的解缠错误，提高解缠质量和效率，解缠相干系数阈值为0.15，低于此阈值的区域不进行解缠。经过上述一系列步骤对所有配对的像对进行了干涉处理，通过查看去平干涉图、相干性图和解缠图的效果（图2），剔除相干性较低的像对，为下一步轨道精炼和重去平做好数据准备。

3.4 轨道精炼和去平

轨道精炼和去平是估算和去除残余相位与解缠后还存在的相位斜坡，重点是控制点的选择，要求控制点相关性高且稳定。本次工作采用PS-InSAR技术获取区域内形变量为0，且相关性高的点作为控制点。

4 地质灾害遥感解译

本次遥感初步解译结合任务要求，充分利用2021年第一季度镶嵌融合影像，开展了1∶5万地质灾害遥感综合调查工作。本次遥感初步解译在较短的时间内，准确而迅速地查明了该区地质灾害的类型、数量、分布范围等，编制了1∶5 万地质灾害遥感解译图，为罗山县地质灾害风险调查评价工作提供了基础资料和科学依据，同时也提高了地质灾害风险初步评价的质量和精度。其取得的初步遥感解译认识如下。

（1）本次遥感解译遵循从已知到未知、从区域到局部、从总体到个别、从定性到定量，按先易后难、循序渐进、不断反馈和逐步深化的原则进行，对罗山县遥感影像进行解译，得到罗山县地质灾害分布发育规律草图。经过解译和野外验证，解译出了滑坡、崩塌和泥石流等地灾隐患54处。

（2）根据解译成果可见，罗山县以崩塌为主，遥感解译发现崩塌38 处；其次为滑坡，有14 处；随后为泥石流，有2处，无地面塌陷；另外在解译过程和实际解译过程中，对于变形特征只有5～10 cm 的小型地质灾害，其判译结果较为困难；对于滑坡等大范围内的地质灾害，其判译特征较明显，解译结果较好。

（3）通过遥感解译，罗山县地质灾害分布表现出大分散小集中的特点，集中分布在定远乡中部及南部、周党镇中部及南部、灵山镇中部及东部。这3 个乡镇崩塌的形成与地貌、构造、岩性、植被和人为活动有关，具有突发性强、速度快、分布范围广和具有一定隐蔽性等特点。根据解译成果，罗山县解译出的38处崩塌，主要分布于南部山区定远乡和周党镇。

5 结论和建议

通过野外排查对罗山县地质灾害和孕灾地质条件建立解译标志，根据解译标志及专家经验，对获取的遥感影像进行人机交互解译，并对解译结果进行定量分析，本次遥感解译的具体成果如下。

（1）建立了地质灾害解译标志和孕灾地质条件解译标志。构建直接解译法、对比法、逻辑推理法、图像处理法相结合的方法进行遥感解译。根据任务要求，选择空间分辨率优于1 m 的GF-2、GF-7 卫星为主、空间分辨率优于2 m 的GF-1 作为补充，以人机交互方式开展解译工作，解译数据可靠，该解译方法可为河南省同类项目提供技术方法支持。

（2）基于INSAR 技术提取了河南省罗山县地表形变信息，初步圈定形变量较大20 处，为地质灾害（隐患）早期识别提供了依据。

（3）通过遥感技术手段，初步查明了工区地形地貌、地质构造、植被覆盖度、工程地质等因子，为后续易发性评价和风险评价提供数据支持。

（4）掌握了针对河南省罗山县1∶5万地质灾害风险调查评价遥感解译技术方法体系，建立了罗山县最新地质灾害本底数据，为今后开展多时相地质灾害遥感动态监测奠定了基础。

（5）从遥感影像图上看到的地质特征和野外核查人员实地看到的眼前现象存在着一定差距，因此野外验证与遥感影像特征描述和判断存在一定偏差。

（6）地面塌陷在影像上一般没有直接解译标志，通过INSAR 技术提取的地表形变信息，一定程度上反映地面沉降特征，但解译地面塌陷单点位置较为困难，建议根据圈定的地面塌陷（隐患区）进行现场排查复核。

（7）由于遥感影像为正射影像，一些小型崩塌不易发现。建议增加野外补充。