CN113866764B - 基于InSAR和LR-IOE模型的滑坡易发性改进评估方法 - Google Patents

基于InSAR和LR-IOE模型的滑坡易发性改进评估方法 Download PDF

Info

Publication number
CN113866764B
CN113866764B CN202110928751.XA CN202110928751A CN113866764B CN 113866764 B CN113866764 B CN 113866764B CN 202110928751 A CN202110928751 A CN 202110928751A CN 113866764 B CN113866764 B CN 113866764B
Authority
CN
China
Prior art keywords
evaluation
landslide
deformation
slope
factor
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110928751.XA
Other languages
English (en)
Other versions
CN113866764A (zh
Inventor
邓辉
张文江
李雪聪
薛东剑
刘严松
王战卫
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chengdu Univeristy of Technology
Original Assignee
Chengdu Univeristy of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chengdu Univeristy of Technology filed Critical Chengdu Univeristy of Technology
Priority to CN202110928751.XA priority Critical patent/CN113866764B/zh
Publication of CN113866764A publication Critical patent/CN113866764A/zh
Application granted granted Critical
Publication of CN113866764B publication Critical patent/CN113866764B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88Radar or analogous systems specially adapted for specific applications
    • G01S13/89Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021SAR image post-processing techniques
    • G01S13/9023SAR image post-processing techniques combined with interferometric techniques
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/18Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06QINFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00Administration; Management
    • G06Q10/06Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063Operations research, analysis or management
    • G06Q10/0639Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G06Q10/06393Score-carding, benchmarking or key performance indicator [KPI] analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06QINFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/10Services
    • G06Q50/26Government or public services
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Human Resources & Organizations (AREA)
  • Remote Sensing (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Tourism & Hospitality (AREA)
  • Development Economics (AREA)
  • Educational Administration (AREA)
  • Economics (AREA)
  • Strategic Management (AREA)
  • Databases & Information Systems (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • General Business, Economics & Management (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Marketing (AREA)
  • Evolutionary Biology (AREA)
  • Electromagnetism (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Game Theory and Decision Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Quality & Reliability (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computing Systems (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Primary Health Care (AREA)

Abstract

基于InSAR和LR‑IOE模型的滑坡易发性改进评估方法,包括如下步骤:S1:收集评价区地质地理数据,提取出评价因子,建立滑坡易发性评价指标体系;S2:利用SBAS‑InSAR技术获取评价区地表沿雷达视线向的形变速率,筛选出可靠形变点,将形变速率由雷达视线方向转换为沿最大坡度方向,作为滑坡易发性评价因子。S3:构建熵指数‑逻辑回归耦合模型,进行滑坡易发性评价。本发明利用了SBAS‑InSAR技术采集长时序地表的变形信息,将沿最大坡度方向形变速率作为滑坡易发性评价因子,基于熵指数‑逻辑回归耦合模型对滑坡易发性进行评价进行优化,使模型对滑坡易发性的预测精度显著提升。

Description

基于InSAR和LR-IOE模型的滑坡易发性改进评估方法
技术领域
本发明涉及地质灾害评价技术领域,具体涉及一种基于InSAR和LR-IOE模型的滑坡易发性改进评估方法。
背景技术
我国由于山区较多,地形复杂,构造发育,地质灾害分布广泛。频发的地质灾害不仅对生态环境造成巨大的破坏,还严重威胁着当地居民的生命财产安全。地质灾害的防治关系到人类的安全与发展问题,成为涉及国家安全与社会稳定的重大问题。对滑坡灾害易发性进行评价分区能够为滑坡地质灾害的防治提供重要的基础资料,因此具有重要的意义。
合成孔径雷达(InSAR),是新近发展起来的空间对地观测技术,是传统的SAR遥感技术与射电天文干涉技术相结合的产物。它利用雷达向目标区域发射微波,然后接收目标反射的回波,得到同一目标区域成像的SAR复图像对,若复图像对之间存在相干条件,SAR复图像对共轭相乘可以得到干涉图,根据干涉图的相位值,得出两次成像中微波的路程差,从而计算出目标地区的地形、地貌以及表面的微小变化,可用于数字高程模型建立、地壳形变探测等。
进入21世纪后地质灾害易发性进入快速发展的阶段,不同种类的数学方法模型在分析地质灾害易发性中得到实践应用。这些数学模型方法主要分为定性和定量方法,定性方法主要为AHP法、多标准分析等,定量方法主要为统计方法和机器学习方法,这些方法被广泛应用于滑坡易发性评价中。在进行滑坡易发性评价中,对评价指标的研究也主要是集中在地形地貌、地层岩性、水系、地质构造、工程活动等方面,未充分考虑InSAR获取的地表形变速率信息对评价模型精度的影响。
滑坡易发性评价是研究滑坡地质灾害发生的可能性,能够为区内的滑坡地质灾害防治与管理提供科学的依据。滑坡的形变与时间息息相关,是一个变化的过程,因此滑坡的易发性评价也需要是动态的。传统的研究一般选用地质、地理等因子构建评价指标体系,未考虑滑坡是一个动态变化过程。
发明内容
为克服现有技术存在的技术缺陷,本发明公开了一种基于InSAR和LR-IOE模型的滑坡易发性改进评估方法。
本发明所述基于InSAR和LR-IOE模型的滑坡易发性改进评估方法,其特征在于,包括如下步骤:
S1:收集评价区地质地理数据,根据滑坡清单数据统计分析各类滑坡隐患分布特征,提取出滑坡易发性评价因子,建立滑坡易发性评价指标体系;
S2:利用SBAS-InSAR技术进行形变信息提取,获取评价区地表沿雷达视线向的形变速率,筛选出可靠形变点;
将形变速率由雷达视线方向转换为沿最大坡度方向,对沿最大坡度方向形变的速率进行分级,作为滑坡易发性评价因子之一;
S3:利用熵指数模型计算各评价因子分级比系数,将其作为逻辑回归模型中的输入数据,构建熵指数-逻辑回归耦合模型,计算逻辑回归系数,并进行滑坡易发性概率计算。
具体的,所述S1步骤中滑坡易发性评价指标体系建立具体为:
S101:将提取出的评价因子,通过GIS平台统一为同一分辨率和坐标系统,并以栅格单元形式存储因子图层,完成因子数据集的量化与属性赋值;
S102:根据评价区已有滑坡点信息清单,对选择的评价因子进行多重共线性诊断,对评价因子的相关性进行评价,去除存在高相关性的评价因子,剩余评价因子构建滑坡易发性评价指标体系;
所述多重共线性分析采用TOL和VIF进行评价。
优选的,所述步骤S102中TOL和VIF的计算公式如下所示:
Figure BDA0003210234240000031
TOL=1/VIF;
Ri为自变量Xi对其余自变量作回归分析的负相关系数;
去除的标准为设置TOL下限或VIF上限,低于TOL下限或高于VIF上限的评价因子被去除。
具体的,所述步骤S2中获取评价区地表沿雷达视线向的形变速率具体方法为:
第j幅干涉图的形变相位值:
Figure BDA0003210234240000041
其中k为干涉图获取时间点,(Tk—Tk-1)为相邻两个获取时间点的时间差;vk为k时刻的形变速率;
将所有的形变相位值组合成矩阵形式:
Δφ=Bv
采用奇异值分解方法求得矩阵Bv的广义逆矩阵,计算得到最终的形变速率。
优选的,所述步骤S2将形变速率由雷达视线方向转换为沿最大坡度方向具体步骤为:
Vslope=VLOS/cosβ
cosβ=nLOS×nslope
nLOS=(-sinθ×cosαs,sinθ×sinαs,cosθ)
nslope=(-sinα×cosφ,-cosα×cosφ,sinφ);
其中Vslope表示沿最大斜坡方向的年均形变速率值;VLOS表示沿雷达视线方向的年均形变速率结果值;α是该点斜坡坡向;φ是该点斜坡的最大坡度;θ为雷达入射角;αs为正北方向和卫星轨道方向的夹角。
优选的,设定当0<|cosβ|<0.3时,取cosβ=0.3。
优选的,所述步骤S3中,评价因子的分级比系数由以下公式计算得到:
Figure BDA0003210234240000042
式中,Pij为各评价因子分级比系数;a为对应分级的面积百分比;b为相应的滑坡百分比。
优选的,所述步骤S3中滑坡易发性概率计算公式为:
Figure BDA0003210234240000051
其中Z为滑坡易发性概率,x1,x2,...,xn代表各滑坡易发性评价因子,β12,...,βn代表逻辑回归系数;α为基础常数。
本发明将熵指数(index of entropy,IOE)模型与逻辑回归(logisticregression,LR)模型相耦合,对滑坡易发性进行评价进行优化,通过IOE模型计算具有相同维度的评价因子分级比系数,可以避免出现滑坡灾害点与评价因子间出现线性相关的问题以及在建模过程减少噪声的出现,因此耦合模型既可以避免评价因子在选取过程中受人为因素的干扰,同时也可以避免机器学习过程中的过拟合问题的出现。
而以往滑坡易发性评价只采用静态因子,本发明利用SBAS-InSAR技术提取形变作为动态因子,表征滑坡的动态变化;而由于滑坡的形变与时间息息相关,是一个变化的过程,而SBAS-InSAR技术可以长时间监测大范围滑坡失稳的迹象,表征滑坡位移在时间序列上的运动特征,因此用地表动态形变情况来更新斜坡的滑坡易发性,能够使模型对滑坡易发性的预测精度显著提升。
附图说明
图1为本发明所述评估方法的一个具体实施方式流程示意图;
图2为本发明所述时间基线图的一个具体实施方式流程示意图;
图2中横坐标为获取数据时间,最小刻度的单位为月,纵坐标为空间基线长度,最小刻度的单位为10米;
图3为本发明所述具体实施例的优化前后评价模型的AUC值(模型评估指标,areaunder the curve)示意图;图3中横坐标为特异性,纵坐标为敏感度,均无单位。
具体实施方式
下面结合附图,对本发明的具体实施方式作进一步的详细说明。
本发明所述基于InSAR和LR-IOE模型的滑坡易发性改进评估方法包括以下步骤:
S1:收集评价区地质地理数据,根据滑坡清单数据统计分析各类滑坡隐患分布特征,提取出滑坡易发性评价因子,建立滑坡易发性评价指标体系;
S2:利用SBAS-InSAR技术进行形变信息提取,获取评价区地表沿雷达视线向的形变速率,筛选出可靠形变点;
将形变速率由雷达视线方向转换为沿最大坡度方向,对沿最大坡度方向形变的速率进行分级,作为滑坡易发性评价因子之一;
S3:利用熵指数模型计算各评价因子分级比系数,将其作为逻辑回归模型中的输入数据,构建熵指数-逻辑回归耦合模型,计算逻辑回归系数,并代入进行滑坡易发性概率计算。
进一步地,所述步骤S1中,根据评价区历史滑坡清单数据,结合区域地质、地理环境的特点选择高程、坡度、坡向、剖面曲率、起伏度、平面曲率、距断层距离、年均降水量、TWI(地形湿度)指数、距水系距离、地层岩性、工程岩组12个评价因子,以空间格网作为基本评价单元,滑坡易发性评价指标体系建立包括以下步骤:
S101:将选择的高程、坡度、坡向等12个评价因子,通过GIS(地理信息系统)平台统一为同一分辨率和坐标系统,并以空间栅格单元形式存储因子图层,完成因子数据集的量化与属性赋值。
S102:根据评价区已有滑坡点信息清单,对选择的高程、坡度、坡向等12个评价因子进行多重线性分析,对评价因子的相关性进行评价,去除存在高相关性的评价因子。
多重线性分析采用容忍度TOL(Tolerance)和方差膨胀因子VIF(VarianceInflation Factor)进行评价,二者互为倒数。
进一步地,步骤S102中计算VIF和TOL的公式如下所示:
Figure BDA0003210234240000071
TOL=1/VIF(2);
式(1)和(2)中,Ri为自变量Xi对其余自变量作回归分析的负相关系数。
方差膨胀因子VIF越大,说明自变量之间存在共线性的可能性越大。一般来讲,如果方差膨胀因子超过10,则回归模型存在严重的多重共线性。
可以设置当自变量的容忍度大于0.1,方差膨胀系数小于10的范围是可以接受的,表明自变量之间没有共线性问题存在。
进一步地,所述步骤S2中,利用SBAS-InSAR技术获取评价区地表形变速率信息,SBAS-InSAR技术原理和SBAS-InSAR技术形变信息提取方法步骤具体说明如下:
小基线合成孔径雷达干涉测量技术(SBAS-InSAR),利用相干系数法对目标存在的高相干点进行识别,得到长时间序列的地表形变。SBAS-InSAR通过设置时间、空间基线的阈值,将SAR数据集分为多个小集合,很好地解决了原影像集间基线过长导致的基线失相干的问题,然后采用最小二乘法(LS)对相位图进行解缠求解,可以有效地减小计算过程中的误差。
原始的干涉相位主要包括以下组成部分:形变相位φdem、地形相位φtopo、参考椭球相位φflat、大气相位φatm和随机噪声相位φnoise,可表示为:
φ=φdemtopoflatnoiseatm (3)
初始相位中的地形相位、大气相位、噪音相位等都需要基于一定的理论方法进行移除。假设获取了覆盖特定研究区域的N+1幅SAR影像,将其他的影像进行配准至主影像。将时间基线和空间基线设置特定值进行多景SAR影像干涉对组合,运用小基线算法对差分干涉图相位解缠。对于从影像Ta和主影像Tb生成的第j幅差分干涉图(Tb>Ta),其任意点的干涉相位Δφj可以表示为:
Figure BDA0003210234240000081
式(4)中:λ是雷达波长,j∈[1,...,M],M为满足时间基线和空间基线阈值约束的多视差分干涉对数量;φTa和φTb是时间Ta和Tb引起的相位变化,dTa和dTb是在时间Ta、Tb相对于初始时刻(dT0=0)的视线向的累计形变量。φtopo和φatm分别为地形相位和大气延迟相位,φnoise为随机噪声相位。
为方便分析,在式(4)中忽略去地形相位、大气延迟相位和随机噪声相位相关;则公式可以简化为:
Figure BDA0003210234240000091
假设Ta和Tb时间间隔内的形变是线性的,即在整个时间段内地表变形是分段线性的,那么第j景干涉图的形变相位值可以写成:
Figure BDA0003210234240000092
其中k为干涉图获取时间点,(Tk—Tk-1)为相邻两个获取时间点的时间差;vk为k时刻的形变速率;
将所有的解缠差分干涉相位组合成矩阵形式:
Δφ=Bv (7)
Bv是M×N的矩阵,采用奇异值分解方法可以求得矩阵的广义逆矩阵,计算得到最终的形变速率或累计位移时间序列,得到沿最大斜坡方向的年均形变速率值Vslope
形变速率分级通常可以按照以下步骤进行:根据评价区形变速率平均值以及历史滑坡的分布,将平均速率作为一类评价指标进行重分类分析,以设置一定的速率间隔将Vslope分为五类,对Vslope进行共线性分析得出该因子的TOL与VIF值分别为0.804与1.243,通过共线性检验,可以参与评价。
SBAS-InSAR技术形变信息提取方法包括以下步骤:
1.短基线(SBAS)技术的数据处理
1.1连接图生成
连接图生成是对输入的雷达数据进行干涉像对配对,按照公式(N×(N-1))/2计算得到最多配对数,N为输入的SAR影像幅数。为了减少运行过程中的冗余量过高以及提高结果精度,设置合理的空间基线阀值和时间基线阀值,进行3D解缠,将获得的干涉像对进行干涉工作流处理,然后用于SBAS反演。
1.2干涉工作流
干涉处理过程将所有影像视为一个集合,对其设置时间和空间基线阈值,将小于两阈值条件的雷达影像进行干涉组合。SBAS技术进行差分干涉处理时,为方便监测长时间的间隔形变信息,通常先对较短时间跨度内的进行处理,然后把相邻时间段的干涉纹图相加,进而获取相应的形变数据。在干涉处理过程中,为了削弱斑点噪声的影响,对图像进行了滤波处理,生成了去平和滤波后的干涉图和相干系数图。
1.3轨道精炼和重去平
根据生成干涉像对滤波后的干涉图、解缠图和相干性图结合强度影像图等去除掉像对质量不好的连接像对,对剩余的干涉像对中依然存在的地形相位和相位跃变进行估算和去除,最主要的是利用GCP点进行相位去除。
1.4两次反演
该步骤是SBAS反演的核心步骤,第一次反演估算了地表形变速率和残余地形;本次评价选择健壮性最好的线性模型(linear)公式如下:
Phase=(Hres×K1)+(V×T1×4π/λ) (8)
Phase是输入的一系列解缠后的相位,Hres是估算残余地形的高程,K1是空间基线和残余地形的比例,V是估算出的平均速率,T1是时间基线,λ为波长。
第二次反演是在第一次反演的基础上去除大气滤波,要将大气相位从残余相位中分离出,首先时间序列上的高通滤波将非线性形变相位φNL去除,剩余了大气相位φatm和失相关噪声φnoi;然后高通滤波基础上,空间域上的低通滤波滤除失相关噪声φnoi,剩余相位即为大气相位φatm
1.5地理编码
为了更好的对比SAR图像几何和辐射特征,将反演结果转化成统一的地理坐标,也就是将SAR数据从斜距坐标系转换到地理坐标系。
2.地表形变速率信息提取
2.1基于RI指数的可视性分析
SBAS-InSAR获取的形变速率方向为沿雷达视线向,由于视线向的几何特性会影响形变速率的贡献值,因此可视性分析是获取研究区准确形变结果的基础,RI指数进一步提升了地形可视化分析的实用意义,R-Index指数的定义如下:
R-Index=sin[θ-β×sin(A)] (9)
其中θ为视线入射角,β为坡度,A为坡向校正系数
当RI=1,这时坡面和雷达波束平行,比较适合获取山区地面信息,可视性最好;当1>RI>sinθ时,该区域坡面背向传感器,斜距方向的分辨率得到很大的提高,为好可视区域;当RI=sinθ时,代表的是坡度角β等于0的平坦区域;当0<RI<sinθ,这些区域发生了透视收缩;当-1<RI≤0时,为发生叠掩和阴影的区域,可视性很差。
2.2可靠形变点获取
除了地形和卫星自身因素对形变点精度有影响外,研究区植被随季节变化也比较大,植被地区的地貌特征变化的也比较快,即便SAR数据时间间隔短,也会导致失相干,可根据评价区土地利用类型分布图,剔除不可靠区域的形变点,又例如根据图像成像质量和可能发生滑坡的地形,例如河流处的河岸线经常发生变化,但通常不发生滑坡,可以剔除部分区域的形变,仅保留可能发生滑坡的地段作为可靠区域形变点。
2.3形变点转换
包括滑坡在内的斜坡地质灾害多沿斜坡滑动面滑动,这就使得沿LOS方向的地表形变速率与实际斜坡上的滑坡形变速率信息并不一致,因此需要将将年均形变速率结果值由雷达视线方向转换为沿最大坡度方向,其理论依据是数字高程模型与二者间的几何关系,其变换关系如下所示:
Vslope=VLOS/cosβ (10)
cosβ=nLOS×nslope (11)
nLOS=(-sinθ×cosαs,sinθ×sinαs,cosθ) (12)
nslope=(-sinα×cosφ,-cosα×cosφ,sinφ) (13)
其中Vslope表示沿最大斜坡方向的年均形变速率值;VLOS表示沿雷达视线方向的年均形变速率结果值;α是该点斜坡坡向;φ是该点斜坡的最大坡度;θ为雷达入射角;αs为正北方向和卫星轨道方向的夹角,即雷达卫星飞行方向,升轨数据与降轨数据方向相反,升轨为负,降轨为正。
为了避免Vslope绝对值被夸大,避免转换系数影响投影测量的精度,设定当0<|cosβ|<0.3时,取cosβ=0.3,即Vslope不能大于VLOS的3.33倍,这样就可以将可接受的投影速度值与投影值不可靠的数据区分开来。且由于滑坡的滑向与斜坡坡向相同,此时Vslope<0,因此只需保留Vslope<0的形变点。
进一步地,所述步骤S3中,构建熵指数-逻辑回归耦合模型(LR-IOE模型)进行滑坡易发性评价,主要包括步骤:
在滑坡易发性分析评价中,不同孕灾因子对滑坡发生的影响贡献程度由不同的熵值代表。熵指数值越大,代表评价因子的影响作用越大,反之亦然,各评价因子的权重值由以下公式计算得到:
Figure BDA0003210234240000131
Figure BDA0003210234240000132
Figure BDA0003210234240000133
Himax=log2 S (17)
Figure BDA0003210234240000134
Figure BDA0003210234240000141
Wi=Ii×Pi (20)
其中,Pij为各评价因子分级比系数;a为对应分级的面积百分比;b为相应的滑坡百分比;(Pij)代表概率密度;Hi、Himax表示熵值;S为评价因子分级数;Ii为评价因子的信息率;Wi为评价因子所占的权重。
逻辑回归是一种用于二分类或多分类的分类预测模型,在滑坡易发性评价中,各评价因子作为模型的输入自变量,滑坡的发生与否是因变量的分类依据,滑坡发生与未发生分别用1和0代表。
设滑坡以概率P发生,以1-P为滑坡不发生的概率,P设置取值范围为[0,1],P的表达式为:
Figure BDA0003210234240000142
式中,x1,x2,...,xn代表各评价因子,β12,...,βn代表逻辑回归系数;α为基础常数,此项为滑坡发生与不发生概率之比的对数值,假设环境为没有任何影响滑坡发生的因素存在。逻辑回归方程为对P/(1-P)取自然对数,可得到:
Figure BDA0003210234240000143
其中Z为滑坡易发性概率。
逻辑回归模型得到的Pij具有相同的维度,可以避免出现滑坡灾害点与评价因子间出现线性相关的问题以及在建模过程减少噪声的出现,将各评价因子按照对应的分级比系数Pij进行分类,并将Pij作为逻辑回归模型中的输入数据,构建熵指数-逻辑回归耦合模型,计算得到混合模型下各指标的β值,对滑坡易发性进行评价。
进一步地,通过ROC曲线(受试者工作特征曲线,receiver operatingcharacteristic curve)对滑坡易发性评价结果精度进行验证。
具体实施例
参考图1,本发明的实施例提供了一种基于InSAR和熵指数-逻辑回归耦合模型的滑坡易发性改进评估方法,包括以下步骤:
S1:收集整理评价区地质、地理环境等多源数据,选择了高程、坡度、坡向、剖面曲率、起伏度、平面曲率、距断层距离、年均降水量、TWI指数、距水系距离、地层岩性、工程岩组等12个滑坡致灾因子,将数据统一到同一投影坐标系和空间分辨率,本次具体实施例选取的评价区为中国青海省西宁市湟中县,该区域内共有历史滑坡地质灾害点147个。
S102:根据评价区已有滑坡点的信息清单,对选择的高程、坡度、坡向等12个评价因子进行多重共线性诊断,多重线性分析采用TOL(Tolerance,容忍度)和VIF(VarianceInflation Factor,方差膨胀因子)进行评价,二者互为倒数,评价区滑坡易发性评价指标共线性特征如表1所示。
表1评价区滑坡易发性评价指标共线性诊断
Figure BDA0003210234240000161
由表1分析得出,上述12个致灾因子变量的VIF均小于5,TOL均大于0.2,这能保证各评价因子之间的独立性,所以都可参与到滑坡易发性模型计算中。
短基线(SBAS)技术的数据处理详细过程如下:
(1)连接图生成
本次采用32幅Sentinel-1A(哨兵1A卫星)数据,地面分辨率为5m×20m,工作模式为IW(干涉宽视场模式),极化方式为VV(垂直同向极化),数据产品为SLC(单视复数影像),间跨度为2018年5月至2020年5月。
根据数据的实际环境情况,将空间基线阈值最大百分比设置为临界基线的3%(5162.12×0.03=154.86m),时间基线阈值设置为180天,最后得到136个干涉对,最大多普勒中心频率基线为291.89Hz;超级主影像为2018年10月2日,其余影像为辅影像,得到的时间基线图如图2所示。
(2)滤波与干涉处理
主辅影像基于相位的精确配准是生成干涉图的基本条件,采用相干系数法基于Sentinel-1(哨兵1号卫星)精密轨道数据和30米的ALOS DEM(先进对地观测卫星的数字高程模型)数据对所有像素进行逐一配准。本次评价滤波方法采用Goldstein法,相位解缠采用Minimum Cost Flow方法,为更好地避免解缠误差,提升处理过程效率,将解缠阈值设为0.2。
(3)轨道精炼和重去平
根据生成干涉像对滤波后的干涉图、解缠图和相干性图结合强度影像图等去除掉像对质量不好的连接像对,这一步是对剩余的干涉像对中依然存在的地形相位和相位跃变进行估算和去除,主要是利用地面控制点(GCP)进行相位去除。
在对评价区地理环境进行评价的基础上,选择道路、建筑等远离残余相位的地面控制点46个,本研究使用健壮性强的线性优化(Polynomial Refinement)进行轨道优化,在基线距较小的条件下提高结果精度。
(4)两次反演
残余地形的去除和地表形变速率的估算均在该步骤第一次反演中操作,采用的方法为奇异值分解(SVD)法。同时,为进行下一步处理,还需优化处理干涉图,需要进行二次解缠。
第二次反演是以第一次反演得到的形变速率为基础,从而反映出形变位移在时间序列上的变化,本次反演根据步骤(3)中选择的地面控制点再次进行轨道精炼和重去平,大气相位的估算和去除利用特定的大气滤波方法,最终结果即为地表形变在时间序列上的位移。
(5)地理编码
地理编码即为将数据从斜距坐标系转换为地理坐标系,本次研究经过地理编码生成视线入射角(ILOS)与视线方位角(ALOS),分别为35.5°与83.1°。编码结果经过重新投影得到了研究区视距传播方向(LOS)方向的地表平均形变速率,结果为-149.83~56.19mm/yr,时间跨度从2018年5月至2020年5月,正、负值表分别代表形变沿朝向与远离雷达视线向位移。速率结果直方图中峰值出现在0左右,形变大多发生在上五庄镇的高山峡谷小范围内,高程均超过3000m,西堡镇,总寨镇等人口聚居区形变速率几乎为0,并未出现大面积沉降或抬升,说明湟中县大部分区域一直保持稳定状态,对形变点进一步筛选与验证以确保SBAS结果的可靠度。
评价区地表形变速率信息提取详细过程
(1)基于R-Index指数的可视性分析
SBAS方法获取的形变速率方向为沿雷达视线向,由于视线向的几何特性会影响形变速率的贡献值,因此可视性分析是获取研究区准确形变结果的基础。R-Index指数的定义如下:R-Index=sin[θ-β×sin(A)],其中θ为视线入射角,β为坡度,A为坡向校正系数,由于本研究数据为升轨,所以A=α+ε+180,其中α为坡向,ε为卫星飞行方向与正北方向的夹角。本次评价以sinθ=0.58为分界值,以区分透视收缩区域,统计R-Index值结果。本次实验保留1>RI>0.58的可视区域,面积为1389.7km2,超过总面积占比一半以上。
(2)可靠形变点筛选
筛选评价区2019年的Landsat8 OLI(美国陆地卫星计划第八颗卫星的陆地成像仪)影像。采用指数综合和CART(分类与回归树)结合的方法获取研究区土地利用数据,其中指数包括归一化水指数(NDWI)、土壤调节植被指数(SAVI)、归一化建筑指数(NDBI)等五类,均为与研究区土地利用类型密切相关的指数,样本点的选择依据Google Earth Pro中的历史高分辨率影像,分类结果表示林地面积为:1202.32km2,耕地面积为1072.08km2,人类工程面积为268.99km2,水域面积为6.17km2,草地面积为325.79km2,林地和水域总面积占研究区总面积的42.02%。受C波段的穿透能力制约,植被高覆盖区域地表形变点的可靠性不高,河流湖泊区域也可排除滑坡的发生,本评价依据土地利用分类结果剔除林地区域及水域的形变点。
(3)形变点转换与可靠性验证
为了避免Vslope绝对值被夸大,避免转换系数影响投影测量的精度,设定当0<|cosβ|<0.3时,取cosβ=0.3,即Vslope不能大于VLOS的3.33倍,这样就可以将可接受的投影速度值与投影值不可靠的数据区分开来。且由于滑坡的滑向与斜坡坡向相同,此时Vslope<0,因此只需保留Vslope<0的形变点,根据以上条件,综合RI指数、土地利用类型,以及滑坡通常不会在坡度低于5°的区域分布的特点,对地表形变点进行筛选,删除研究区边缘处过于分散的点以避免大气误差对精度的影响,最终保留了265622个点。区内沉降速率最大达202.06mm/yr(毫米/年),为上五庄镇北部高山地带,由于物源条件差人迹罕至,因此未对人类造成真正危害,区内总寨镇、土门关乡、田家寨镇形变点速率大多不超过5mm/yr,其他平原人类聚居区分布少量形变点,因此判断最大坡度向形变速率空间分布较符合实际情况。
(4)基于IOE-LR耦合模型滑坡易发评价
评价区的样本数据选取了历史滑坡灾害点与非滑坡灾害点,二者选取的数量相同,其中非滑坡点为利用ArcGIS软件随机生成。训练样本集为从总训练样本集中随机选取,数量为滑坡点与非滑坡点总数的70%,共计205个,剩余的88个则作为测试样本集。IOE-LR耦合模型以本评价选取的12个评价指标的熵值作为自变量,因变量则以滑坡是否发生为依据进行滑坡易发性分析。利用逻辑回归模型对训练样本集进行学习,得到所有评价指标的逻辑回归系数β12,...,βn
(5)InSAR和熵指数-逻辑回归耦合模型的滑坡易发性改进评估方法
以InSAR形变速率作为一类评价因子来进行易发性评价优化,结论得出Vslope作为衡量滑坡变化的物理量可以反映出滑坡的形变特征,所以Vslope值可以表示动态变化下的滑坡易发性大小。因此在易发性IOE-LR耦合模型评价基础上,利用SBAS-InSAR监测成果Vslope作为新增的评价因子进行易发性评价优化,以期能更精确地预测滑坡灾害发生概率。根据评价区形变速率平均值为-6.6mm/yr以及历史滑坡的分布,将平均速率作为一类评价指标进行重分类分析,以-6mm/yr为间隔将Vslope分为五类,对Vslope进行共线性分析得出该因子TOL与VIF分别为0.804与1.243,通过共线性检验,可以参与评价。
(6)滑坡易发性评价结果及检验
本次评价采用IOE-LR耦合模型,分别对选用高程、坡度、坡向等12个评价因子和添加SBAS-InSAR技术提取的地表形变因子的滑坡易发性进行评价。对IOE-LR耦合模型与优化后模型ROC曲线(受试者工作特征曲线,receiver operating characteristic curve)进行分析,如图3所示,左边a部分为优化前,右边b部分为优化后,优化前后评价模型的AUC值分别为0.898,0.907,说明优化模型更能准确反映湟中县的滑坡灾害分布规律。
将历史滑坡点与滑坡易发性分区图进行分析得出,处于极高易发区和高易发区的所占比例由优化前的66.7%提高到73.5%,将优化后的滑坡易发性评价图与研究区265622个Vslope矢量点进行叠加分析,可见大部分形变点均落在高或极高易发性区域上,平原地区及高山森林中的易发性等级与形变点分布相符,聚集效应明显,综合以上分析得出优化后的滑坡易发性评价结果更准确,符合实际情况。
前文所述的为本发明的各个优选实施例,各个优选实施例中的优选实施方式如果不是明显自相矛盾或以某一优选实施方式为前提,各个优选实施方式都可以任意叠加组合使用,所述实施例以及实施例中的具体参数仅是为了清楚表述发明人的发明验证过程,并非用以限制本发明的专利保护范围,本发明说明书括号中的英文为括号前中文在计算机语言或函数中的英文名称,不在括号中的英文为没有严格中文译文的数学函数或数据库等的命名,本发明的专利保护范围仍然以其权利要求书为准,凡是运用本发明的说明书及附图内容所作的等同结构变化,同理均应包含在本发明的保护范围内。

Claims (4)

1.基于InSAR和LR-IOE模型的滑坡易发性改进评估方法,其特征在于,包括如下步骤:
S1:收集评价区地质地理数据,根据滑坡清单数据统计分析各类滑坡隐患分布特征,提取出滑坡易发性评价因子,建立滑坡易发性评价指标体系;
S2:利用SBAS-InSAR技术进行形变信息提取,获取评价区地表沿雷达视线向的形变速率,筛选出可靠形变点;
将形变速率由雷达视线方向转换为沿最大坡度方向,对沿最大坡度方向形变的速率进行分级,作为滑坡易发性评价因子之一;
S3:利用熵指数模型计算各评价因子分级比系数,将其作为逻辑回归模型中的输入数据,构建熵指数-逻辑回归耦合模型,计算逻辑回归系数,并进行滑坡易发性概率计算;
所述步骤S2中获取评价区地表沿雷达视线向的形变速率具体方法为:
第j幅干涉图的形变相位值:
Figure FDA0004156204340000011
其中k为干涉图获取时间点,(Tk-Tk-1)为相邻两个获取时间点的时间差;vk为k时刻的形变速率;
将所有的形变相位值组合成矩阵形式:
Δφ=Bv
采用奇异值分解方法求得矩阵Bv的广义逆矩阵,计算得到最终的形变速率;
所述步骤S2将形变速率由雷达视线方向转换为沿最大坡度方向具体步骤为:
Vslope=VLOS/cosβ
cosβ=nLOS×nslope
nLOS=(-sinθ×cosαs,sinθ×sinαs,cosθ)
nslope=(-sinα×cosφ,-cosα×cosφ,sinφ);
其中Vslope表示沿最大斜坡方向的年均形变速率值;VLOS表示沿雷达视线方向的年均形变速率结果值;α是该点斜坡坡向;φ是该点斜坡的最大坡度;θ为雷达入射角;αs为正北方向和卫星轨道方向的夹角;
滑坡易发性概率
Figure FDA0004156204340000021
其中
Figure FDA0004156204340000022
x1,x2,...,xn代表各评价因子,β12,...,βn代表逻辑回归系数;α为基础常数;
各评价因子的权重值由以下公式计算得到:
Figure FDA0004156204340000023
Figure FDA0004156204340000024
/>
Figure FDA0004156204340000025
Himax=log2S (17)
Figure FDA0004156204340000031
Figure FDA0004156204340000032
Wi=Ii×Pi (20)
其中,Pij为各评价因子分级比系数;a为对应分级的面积百分比;b为相应的滑坡百分比;(Pij)代表概率密度;Hi、Himax表示熵值;S为评价因子分级数;Ii为评价因子的信息率;Wi为评价因子所占的权重。
2.如权利要求1所述的评估方法,其特征在于,所述S1步骤中滑坡易发性评价指标体系建立具体为:
S101:将提取出的评价因子,通过GIS平台统一为同一分辨率和坐标系统,并以栅格单元形式存储因子图层,完成因子数据集的量化与属性赋值;
S102:根据评价区已有滑坡点信息清单,对选择的评价因子进行多重共线性诊断,对评价因子的相关性进行评价,去除存在高相关性的评价因子,剩余评价因子构建滑坡易发性评价指标体系;
所述多重共线性分析采用TOL和VIF进行评价。
3.如权利要求2所述的评估方法,其特征在于,所述步骤S102中TOL和VIF的计算公式如下所示:
Figure FDA0004156204340000033
TOL=1/VIF;
Ri为自变量Xi对其余自变量作回归分析的负相关系数;
去除的标准为设置TOL下限或VIF上限,低于TOL下限或高于VIF上限的评价因子被去除。
4.如权利要求1所述的评估方法,其特征在于,设定当0<|cosβ|<0.3时,取cosβ=0.3。
CN202110928751.XA 2021-08-13 2021-08-13 基于InSAR和LR-IOE模型的滑坡易发性改进评估方法 Active CN113866764B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110928751.XA CN113866764B (zh) 2021-08-13 2021-08-13 基于InSAR和LR-IOE模型的滑坡易发性改进评估方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110928751.XA CN113866764B (zh) 2021-08-13 2021-08-13 基于InSAR和LR-IOE模型的滑坡易发性改进评估方法

Publications (2)

Publication Number Publication Date
CN113866764A CN113866764A (zh) 2021-12-31
CN113866764B true CN113866764B (zh) 2023-05-26

Family

ID=78990482

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110928751.XA Active CN113866764B (zh) 2021-08-13 2021-08-13 基于InSAR和LR-IOE模型的滑坡易发性改进评估方法

Country Status (1)

Country Link
CN (1) CN113866764B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114185048B (zh) * 2022-02-15 2022-05-31 湖南吉赫信息科技有限公司 地基InSAR提取滑坡位移向量的方法、系统及存储介质
CN114299402B (zh) * 2022-03-07 2022-05-20 成都理工大学 隐患点自动识别方法、电子设备及计算机可读存储介质
CN114818547A (zh) * 2022-06-10 2022-07-29 重庆地质矿产研究院 一种基于数据模型的浅层滑坡易感性评估方法
CN115345511B (zh) * 2022-08-29 2023-06-06 中咨数据有限公司 一种公路走廊滑坡危险性动态评价方法、评价系统及设备
CN115114807B (zh) * 2022-08-29 2022-12-20 成都理工大学 一种水库库岸滑坡易发性评价方法
CN115512531B (zh) * 2022-09-28 2023-06-23 重庆地质矿产研究院 一种基于形变有序性的滑坡灾害多监测点融合预警方法
CN116050120B (zh) * 2023-01-06 2023-09-01 中国自然资源航空物探遥感中心 一种滑坡隐患活动性遥感评价建模方法、系统和存储介质
CN116561536B (zh) * 2023-07-11 2023-11-21 中南大学 一种滑坡隐患的识别方法、终端设备及介质

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109165424A (zh) * 2018-08-03 2019-01-08 四川理工学院 一种基于国产gf-1卫星数据的滑坡易发性评估方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104715255B (zh) * 2015-04-01 2017-11-21 电子科技大学 一种基于sar图像的滑坡信息提取方法
CN110569554B (zh) * 2019-08-13 2020-11-10 成都垣景科技有限公司 一种基于空间逻辑回归与地理探测器的滑坡易发性评价方法
CN111798135B (zh) * 2020-07-06 2022-03-22 天津城建大学 一种基于多源数据集成的高铁沉降危害性评估方法
AU2020103096A4 (en) * 2020-10-29 2021-01-07 Bajrang Tapase, Anand DR Movements/shifts/displacements monitoring SMART box of Earth Retaining Structures in Landslides Mitigation
CN113158570A (zh) * 2021-04-26 2021-07-23 电子科技大学 一种融合多源卫星遥感的全天候地表温度近实时反演方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109165424A (zh) * 2018-08-03 2019-01-08 四川理工学院 一种基于国产gf-1卫星数据的滑坡易发性评估方法

Also Published As

Publication number Publication date
CN113866764A (zh) 2021-12-31

Similar Documents

Publication Publication Date Title
CN113866764B (zh) 基于InSAR和LR-IOE模型的滑坡易发性改进评估方法
Braun Retrieval of digital elevation models from Sentinel-1 radar data–open applications, techniques, and limitations
US11269072B2 (en) Land deformation measurement
Jebur et al. Using ALOS PALSAR derived high-resolution DInSAR to detect slow-moving landslides in tropical forest: Cameron Highlands, Malaysia
CN109388887B (zh) 一种地面沉降影响因素定量分析方法及系统
CN110427857A (zh) 一种基于遥感数据融合的输电线路地质灾害分析方法
CN109782282A (zh) 一种集成对流层大气延迟改正的时间序列InSAR分析方法
Shamshiri et al. Improving tropospheric corrections on large-scale Sentinel-1 interferograms using a machine learning approach for integration with GNSS-derived zenith total delay (ZTD)
Kwak et al. Near real-time flood volume estimation from MODIS time-series imagery in the Indus River Basin
Zhang et al. Reconstruction of GF-1 soil moisture observation based on satellite and in situ sensor collaboration under full cloud contamination
Kouhartsiouk et al. The application of DInSAR and Bayesian statistics for the assessment of landslide susceptibility
Zhang et al. Deformations monitoring in complicated-surface areas by adaptive distributed Scatterer InSAR combined with land cover: Taking the Jiaju landslide in Danba, China as an example
Tang et al. Changes of Chinese coastal regions induced by land reclamation as revealed through TanDEM-X DEM and InSAR analyses
Goyal et al. Comparison and validation of satellite-derived digital surface/elevation models over India
Zhang et al. Reduction of atmospheric effects on InSAR observations through incorporation of GACOS and PCA into small baseline subset InSAR
Dominici et al. Multiscale documentation and monitoring of L’Aquila historical centre using UAV photogrammetry
Castaneda et al. Dedicated SAR interferometric analysis to detect subtle deformation in evaporite areas around Zaragoza, NE Spain
Miyamoto et al. Using multimodal learning model for earthquake damage detection based on optical satellite imagery and structural attributes
Luo et al. Ice flow velocity mapping in East Antarctica using historical images from 1960s to 1980s: recent progress
Wang et al. Mapping three-dimensional urban structure by fusing landsat and global elevation data
Yang et al. Monitoring, prediction, and evaluation of mountain geological hazards based on InSAR technology
Aydöner et al. The role of the integration of remote sensing and GIS in land use/land cover analysis after an earthquake
CN114091274A (zh) 一种滑坡易发性评价方法及系统
CN115079172A (zh) 一种MTInSAR滑坡监测方法、设备及存储介质
Karami et al. Monitoring of land surface displacement based on SBAS-InSAR time-series and GIS techniques: A case study over the Shiraz Metropolis, Iran

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant