一种城市岩溶致陷因子定量识别和灾害风险评估方法与流程

1.本发明涉及地质灾害监测、机理分析和灾害评估技术领域，尤其涉及一种城市岩溶致陷因子定量识别和灾害风险评估方法。


背景技术：

2.岩溶地表塌陷是指覆盖在溶蚀洞穴之上的松散土体，在自然因素或人为因素作用下突发性下沉而形成圆锥形塌陷坑的一种地质灾害。岩溶地表塌陷事件不仅会直接造成经济损失，还会严重危害居民区、交通、工程建设和供水的安全。
3.目前的岩溶地表监测技术方法主要是围绕岩溶发育程度(溶洞、裂隙、溶槽的发育和分布)、上覆土结构和厚度、水文动力(地下水水位、水/气压力、流速/向、坡降、ph值等)以及地表位移来开展岩溶塌陷调查和监测的。监测技术方法主要包括接触式原位监测传感器、地质雷达、光纤技术、时域反射技术、综合物探法、gps观测等。这些方法直接面向对象监测，监测精度高，但由于密集的城市建设和自然侵蚀过程导致频繁的地貌变化，使得上述物原位接触性(甚至破坏性)的现场监测或调查方法在广域岩溶地表监测上存在工作效率低下、密度不足的问题。鉴于岩溶地表塌陷影响因子层的结构和层次非常明晰，目前我国对岩溶地表塌陷灾害的风险评估主要是基于层次分析法(ahp)法进行的。然而这种定性方法过于依赖专家对各因子层的打分，且没有顾及与塌陷坑相伴相生的地表形变因素的影响，对缺少先验信息的地区，这类方法的盲目性较大。


技术实现要素：

4.针对传统方法难以对广域岩溶地表塌陷的潜在影响因子进行定量识别以及灾害风险等级分区的问题，本发明提供一种城市岩溶致陷因子定量识别和灾害风险评估方法，充分利用高时空分辨率、长时序的地表形变结果，实现了不同水系对岩溶地表塌陷的影响程度的定量识别，以及所有潜在影响因子对岩溶地表塌陷的影响程度的定量识别，在此基础上，通过建立二元逻辑回归模型和基于加权角值变形法(weighted angular distortion，wad)分别实现了对研究区的岩溶地表塌陷的风险分区。
5.本发明提供一种城市岩溶致陷因子定量识别和灾害风险评估方法，包括以下步骤：
6.步骤1：根据研究区的长时序terrasar-x高分辨sar影像，利用stamps-sbas方法提取所述研究区的地表形变速率和地表位移时序；
7.步骤2：建立所述地表形变速率与逐点的地表位移时序、地下水位时序和长江水位时序之间的交叉小波变换模型，从特殊点位尺度定量提取时滞效应，进而识别出不同水系对岩溶地表塌陷的影响程度；
8.步骤3：对研究区的岩溶地表塌陷的潜在影响因子进行分级，并基于历史记录的岩溶地表塌陷事故确定所有潜在影响因子各分级的确定性系数；
9.步骤4：建立所述地表形变速率和所有潜在影响因子各分级的确定性系数之间的
空间地理探测模型，根据所述空间地理探测模型计算得到的空间分异程度从空间尺度上确定所有潜在影响因子对岩溶地表塌陷的影响程度；
10.步骤5：联合步骤2的逐点多时序分析结果和步骤4的空间分异因子探测结果选取影响程度较大的前n个潜在影响因子作为岩溶地表塌陷的重要影响因子；
11.步骤6：建立岩溶地表塌陷事件与地表形变速率、各个重要影响因子的确定性系数之间的二元逻辑回归模型，实现对研究区的岩溶地表塌陷进行风险分区。
12.进一步地，所述方法还包括：
13.步骤7：根据步骤1中计算得到的所述地表形变速率，提取研究区的沉降水平梯度；
14.步骤8：根据所述沉降水平梯度、地表形变速率和人类工程建设空间分布，采用加权角值变形法计算研究区的岩溶地表塌陷的风险值；
15.步骤9：根据所述风险值，基于gis中自然断点算法实现对研究区的岩溶地表塌陷进行风险分区；
16.步骤10：联合步骤2中的逐点多时序分析结果和步骤6中的二元逻辑回归模型的风险分区结果和步骤9中的加权角值变形法风险分区结果，最终得到研究区的岩溶地表塌陷风险评估结果。
17.进一步地，所述潜在影响因子包括：地表形变速率、地层岩性、岩溶发育程度、上覆土结构与厚度、距离四级以上水系的距离、第四系沉积物富水性、距离地铁和大型工地的距离以及城市法定图则。
18.进一步地，步骤2具体包括：
19.步骤2.1：确定研究区内的历史岩溶地表塌陷点位，去除所述历史岩溶地表塌陷点位的地表位移时序、地下水位时序和长江水位时序的趋势项；
20.步骤2.2：利用交叉小波变换分别提取去除趋势项后的地表位移时序、地下水水位时序和长江水位时序的通用功率谱和相角，将所述相角转换为时间延迟；
21.步骤2.3：基于所述通用功率谱和所述时间延迟，定量识别不同水系对岩溶地表塌陷的影响程度。
22.进一步地，步骤4中，所述地表形变速率与所有潜在影响因子各分级的确定性系数之间的空间地理探测模型具体为：以所有潜在影响因子各分级的确定性系数作为自变量，以所述地表形变速率作为属性变量，得到单个潜在影响因子和多个潜在因子的耦合作用对岩溶地表形变空间分布的空间分异程度。
23.进一步地，步骤6中，按照公式(1)建立岩溶地表塌陷事件与地表形变速率、各个重要影响因子的确定性系数之间的二元逻辑回归模型：
[0024][0025]
其中，p是发生岩溶地表塌陷的概率；b0是截距；b1,
…bn
是n个重要影响因子的系数；x1,
…
xn为n个重要影响因子的确定性系数，其中，n个所述重要影响因子中包括地表形变速率。
[0026]
进一步地，步骤8中，按照公式(4)计算研究区的岩溶地表塌陷的风险值risk：
[0027]
risk＝((3
×
shg) v
sub
)
×
mc
den
ꢀꢀ
(4)
[0028]
其中，shg表示沉降水平梯度，v
sub
表示地表形变速率，mc
den
表示人类工程建设空间
分布对岩溶地表沉降的贡献程度。
[0029]
本发明的有益效果：
[0030]
相较于传统的岩溶地表致陷因子识别方法(主要是对岩溶分布与沉降特征在空间上进行定性分析；或者是基于少数单点比较的位移时序与地下水变化的定量分析，而没有从空间尺度定量比较单个潜在影响因子或潜在影响因子之间的耦合作用对岩溶地表塌陷的贡献程度以及因子显著性)，本发明提出的一种城市岩溶致陷因子定量识别和灾害风险评估方法，充分顾及了基于多时相insar获取的高精度广域地表形变速率，将其作为岩溶地表塌陷的一种直观指标，建立其与岩溶地表塌陷潜在影响因子的空间对应关系，从而能够定量客观的反映不同潜在影响因子对岩溶地表沉降分布的贡献程度，进而有针对性的对岩溶地表塌陷灾害进行风险评估建模，避免了因子的冗余问题，提高了模型的精度。另外，直接基于insar实测的岩溶地表形变速率和重要影响因子的加权角值变形法也能够提供研究区最客观的风险情况，从而用于指导那些缺乏先验信息的岩溶地区进行影响因子标度确定和分级，进而反过来提高二元逻辑回归模型的建模精度。
附图说明
[0031]
图1为本发明实施例提供的一种城市岩溶致陷因子定量识别和灾害风险评估方法的流程示意图之一；
[0032]
图2为本发明实施例提供的一种城市岩溶致陷因子定量识别和灾害风险评估方法的流程示意图之二；
[0033]
图3为本发明实施例提供的地表位移时序、地下水水位时序和长江水位时序的通用功率谱和相角图；
[0034]
图4为本发明实施例提供的纳入地表形变速率因子和未纳入地表形变速率因子的岩溶地表塌陷的风险分区及统计结果；
[0035]
图5为本发明实施例提供的基于wad法的岩溶地表塌陷的风险分区图；
具体实施方式
[0036]
为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0037]
实施例1
[0038]
如图1所示，本发明实施例提供一种城市岩溶地表塌陷灾害风险评估方法，包括以下步骤：
[0039]
s101：根据研究区的长时序terrasar-x高分辨sar影像，利用stamps-sbas方法提取所述研究区的地表形变速率和地表位移时序；
[0040]
s102：建立所述地表形变速率与逐点的地表位移时序、地下水位时序和长江水位时序之间的交叉小波变换模型，从特殊点位尺度定量提取时滞效应，进而识别出不同水系对岩溶地表塌陷的影响程度；
[0041]
s103：对研究区的岩溶地表塌陷的潜在影响因子进行分级，并基于历史记录的岩
溶地表塌陷事故确定所有潜在影响因子各分级的确定性系数；
[0042]
s104：建立所述地表形变速率和所有潜在影响因子各分级的确定性系数之间的空间地理探测模型，根据所述空间地理探测模型计算得到的空间分异程度从空间尺度上确定所有潜在影响因子对岩溶地表塌陷的影响程度；
[0043]
s105：联合步骤2的逐点多时序分析结果和步骤4的空间分异因子探测结果选取影响程度较大的前n个潜在影响因子作为岩溶地表塌陷的重要影响因子；
[0044]
s106：建立岩溶地表塌陷事件与地表形变速率、各个重要影响因子的确定性系数之间的二元逻辑回归模型，实现对研究区的岩溶地表塌陷进行风险分区。
[0045]
本发明实施例提供的一种城市岩溶地表塌陷灾害风险评估方法，充分考虑了基于多时相insar的高精度广域地表形变速率，将其作为岩溶地表塌陷的一种直观指标，建立其与岩溶地表塌陷潜在影响因子之间的空间对应关系，从而实现定量客观的反映不同潜在影响因子对岩溶地表沉降分布的贡献程度，进而有针对性的对岩溶地表塌陷灾害进行风险评估建模，避免了影响因子的冗余问题，提高了风险评估模型的精度。
[0046]
实施例2
[0047]
在上述实施例的基础上，本发明实施例还提供一种城市岩溶致陷因子定量识别和灾害风险评估方法，与上述实施例的不同之处在于，本发明实施例还包括以下步骤：
[0048]
s107：根据步骤s101中计算得到的所述地表形变速率，提取研究区的沉降水平梯度；
[0049]
s108：根据所述沉降水平梯度、地表形变速率和人类工程建设空间分布，采用加权角值变形法计算研究区的岩溶地表塌陷的风险值；
[0050]
s109：根据所述风险值，基于gis中自然断点算法实现对研究区的岩溶地表塌陷进行风险分区；
[0051]
s110：联合步骤s102中的逐点多时序分析结果和步骤s106中的二元逻辑回归模型的风险分区结果和步骤s109中的加权角值变形法风险分区结果，最终得到研究区的岩溶地表塌陷风险评估结果。
[0052]
本发明实施例提供的一种城市岩溶地表塌陷灾害风险评估方法，直接基于多时相insar实测的岩溶地表形变速率和重要影响因子的加权角值变形法能提供研究区最客观的风险情况，也可指导那些缺乏先验信息的岩溶地区进行影响因子标度确定和分级，反过来提高二元逻辑回归模型的建模精度。
[0053]
实施例3
[0054]
在上述各实施例的基础上，以武汉地区为例，如图2所示，本发明实施例提供一种城市岩溶致陷因子定量识别和灾害风险评估方法，包括以下步骤：
[0055]
s201：根据研究区的长时序terrasar-x高分辨sar影像，利用stamps-sbas方法提取所述研究区的地表形变速率和地表位移时序；
[0056]
具体地，收集覆盖研究区的terrasar-x高时空分辨率sar影像并对主影像进行配准重采样。按照一定的时空基线阈值对sar时序进行组合并差分干涉、缓慢失相干滤波点(sdfp)选取、相位解缠，然后利用insar通用大气改正在线产品(gacos)迭代进行大气估计、相位解缠结果的精炼，最后对正确解缠得到的m个解缠图利用最小二乘法求取武汉地区的地表平均形变速率和地表位移时序。
[0057]
本实施例中，根据几何失相干理论和数据处理经验确定空间基线阈值为250m，时间基线阈值为90天。
[0058]
s202：建立所述地表形变速率与逐点的地表位移时序、地下水位时序和长江水位时序之间的交叉小波变换模型，从特殊点位尺度定量提取时滞效应，进而识别出不同水系对岩溶地表塌陷的影响程度；
[0059]
具体地，本步骤主要包括以下子步骤：
[0060]
s2021：对步骤s201得到的岩溶地表形变速率进行90m*90m采样处理，对采样数据中的空值采取邻近内插的方式进行补空。
[0061]
s2022：确定研究区内的历史岩溶塌陷点位，去除所述历史岩溶塌陷点位的地表位移时序、地下水位时序和长江水位时序的趋势项；
[0062]
s2023：利用交叉小波变换分别提取去除趋势项后的地表位移时序、地下水水位时序和长江水位时序的通用功率谱和相角，将所述相角转换为时间延迟(如图3)；
[0063]
s2024：基于所述通用功率谱和所述时间延迟，定量识别不同水系对岩溶地表塌陷的影响程度。
[0064]
s203：对研究区的岩溶地表塌陷的潜在影响因子进行分级，并基于历史记录的岩溶地表塌陷事故确定所有潜在影响因子各分级的确定性系数；
[0065]
具体地，根据自然断点算法和《岩溶塌陷调查规范(1:50000)》对岩溶地表塌陷的潜在影响因子进行分级；依据落入各级别中的天坑数量和该级别的面积，利用确定性系数法(cf)确定该分级的确定性系数cf
ij
{i＝1,2,3,4,5；j＝1,2,
…
,4}。
[0066][0067]
上式中，cf
ij
表示第i个影响因子的第j个分级的确定性系数；ppa是某个潜在影响因子的第j个等级中天坑的个数与该等级面积的比值；pps是研究区的总天坑个数与武汉地区总岩溶面积的比值。
[0068]
本实施例中，所述潜在影响因子包括：地表形变速率、地层岩性、岩溶发育程度、上覆土结构与厚度、距离四级以上水系的距离、第四系沉积物富水性、距离地铁和大型工地距离以及城市法定图则。
[0069]
从地表形变速率来看，使用gis中自然间断算法将武汉地区内地表形变速率按-89.7～5.8mm/yr，-5.7～1.3mm/yr，-1.2～2.3mm/yr和2.4～29mm/yr进行分类。
[0070]
从地层岩性来看，研究区以覆盖型岩溶、埋藏型岩溶以及非碳酸盐岩区为主，它们共同构成了武汉岩溶地表塌陷的高、中、低敏感区。
[0071]
从岩溶发育程度来看，根据《岩溶塌陷调查规范(1:50000)》，将岩溶发育程度划分为发育良好(κ》10％)、中等发育(10％≥κ≥3％)、轻度发育(κ《3％)和非岩溶区域，其中κ为钻孔岩溶遇洞率。
[0072]
从上覆土结构与厚度来看，根据历史记录的岩溶塌陷事件，大多数陷坑发生在典型二元结构覆盖区，对应的土壤厚度一般《15m；接着是厚度为15-30m的多层软土结构区，再来就是覆盖厚度为30-40m的埋藏型岩溶和单层土壤结构区。此外，在土壤厚度超过40m的地
区也记录到了一些岩溶塌陷。
[0073]
从距离四级以上水系的距离与第四系沉积物富水性来看，按照0-1000m，1000-3000m，3000-5000m和》5000m对距离四级以上水系的距离进行多元缓冲区提取。根据361个钻井数据，对第四系沉积物富水性分类为》1000m3/d，100-1000m3/d和《100m3/d。
[0074]
从距离地铁和大型工地距离与城市法定图则来看，对所有地铁线路及大型建筑工地按照《500m，500-1000m，1000-2000m和》2000m进行多元缓冲区分析。对城市规划图则按制造业(m)-交通(t)-市政公用事业(u)-仓库(w)用地、住宅(r)-商业(c)用地、绿地和农业(g)用地以及生态控制用地和水域进行分类。
[0075]
s204：建立所述地表形变速率和所有潜在影响因子各分级的确定性系数之间的空间地理探测模型，根据所述空间地理探测模型计算得到的空间分异程度从空间尺度上确定所有潜在影响因子对岩溶地表塌陷的影响程度；
[0076]
具体地，本步骤中，所述地表形变速率与所有潜在影响因子各分级的确定性系数之间的空间地理探测模型具体为：以所有潜在影响因子各分级的确定性系数作为自变量，以所述地表形变速率作为属性变量，得到单个潜在影响因子和多个潜在因子的耦合作用对岩溶地表形变空间分布的空间分异程度。
[0077]
本实施例中，通过步骤s203和s204可以实现从特殊点位到区域对岩溶地表塌陷影响因子的定量识别。
[0078]
s205：联合步骤s202的逐点多时序分析结果和步骤s204的空间分异因子探测结果选取影响程度较大的前n个潜在影响因子作为岩溶地表塌陷的重要影响因子，
[0079]
具体地，本步骤实现了从逐点尺度和空间尺度上选取影响程度较大的前n个潜在影响因子作为重要影响因子。
[0080]
s206：建立岩溶塌陷事件与地表形变速率、各个重要影响因子的确定性系数之间的二元逻辑回归模型(lr模型)，实现对研究区的岩溶地表塌陷进行风险分区；
[0081]
具体地，按照公式(1)建立岩溶塌陷事件与地表形变速率、各个重要影响因子的确定性系数之间的二元逻辑回归模型：
[0082][0083]
其中，p是发生岩溶地表塌陷的概率；b0是截距；b1,
…bn
是n个重要影响因子的系数；x1,
…
xn为n个重要影响因子的确定性系数，其中，n个所述重要影响因子中包括地表形变速率。
[0084]
接着，根据计算得到的各采样单元的logit(p)后，代入公式(2)求得各个采样单元内发生岩溶塌陷事故的概率：
[0085][0086]
本实施例中，分别纳入和排除地表形变速率作为重要影响因子，按照gis中的自然间断算法对基于二元逻辑回归计算出的风险概率值p进行分类，获得岩溶地表塌陷的风险分区图(如图4所示)。比较在考虑地表形变速率和不考虑地表形变速率两种情况下的二元逻辑回归模型对岩溶地表塌陷的风险分区结果，可以看到地表形变速率对岩溶地表塌陷的风险评估结果的影响(如图4所示)。
[0087]
本实施例中，评价指标为：
[0088]
1)灾害点落在高易发区的百分比应该最大和低易发区占整个研究区面积的百分比应该最大；
[0089]
2)灾害点落在高易发区的百分比/各等级区的面积占整个研究区的总面积的百分比应该最大；
[0090]
从图4中可以看出：考虑地表形变速率和不考虑地表形变速率两种情况下的二元逻辑回归模型对岩溶地表塌陷的风险分区结果差异显著。如，高、中、低风险类别的面积分别减小3.9％、减小36.7％和增大37.1％。排除地表形变速率的二元逻辑回归模型的风险分区导致落入高风险等级的岩溶塌陷事故预测降低了13％，同时增加8％的随机生成的非岩溶塌陷点落入岩溶区。可见，本发明实施例中，纳入基于多时相insar的地表形变速率因子的二元逻辑回归模型能够提高对岩溶地表塌陷风险分区的合理性和准确度。
[0091]
s207：根据步骤s201中计算得到的所述地表形变速率，提取研究区的沉降水平梯度；
[0092]
具体地，利用角值变形法(ad)提取相邻sdfp点的沉降水平梯度，具体为，对地理坐标系下的岩溶地表形变投影后，按照公式(3)计算得到研究区的沉降水平梯度：
[0093]
沉降水平梯度＝(相邻sdfp点的地表形变速率差/相邻sdfp点的水平距离)
[0094]
×
100％(3)
[0095]
沉降水平梯度可以反映岩溶地表塌陷对地表结构的撕裂破坏程度。
[0096]
s208：根据所述沉降水平梯度、地表形变速率v
sub
和人类工程建设空间分布，采用加权角值变形法(wad)计算研究区的岩溶地表塌陷的风险值；
[0097]
本发明实施例中，距离地铁和大型工地的距离对岩溶地表形变影响较大，将距离地铁和大型工地的距离点线密度作为人类工程建设空间分布对岩溶地表沉降的贡献程度mc
den
。
[0098]
因此，按照公式(4)计算得到研究区的岩溶地表塌陷的风险值：
[0099]
risk＝((3
×
shg) v
sub
)
×
mc
den
ꢀꢀ
(4)
[0100]
s209：根据所述风险值，基于gis中自然断点算法实现对研究区的岩溶地表塌陷进行风险分区；
[0101]
s210：联合步骤s202中的逐点多时序分析结果和步骤s206中的二元逻辑回归模型的风险分区结果和步骤s209中的加权角值变形法风险分区结果，最终得到研究区的岩溶地表塌陷风险评估结果。
[0102]
对比图4和图5，被wad法识别到拥有较大的shg值的区域，被二元逻辑回归模型识别为中高风险等级的可能性也越大。例如，后湖大道和和平大道(图4和图5中的
①
标记)，墨水湖北路、汉阳鹦鹉-建港路(图4和图5中的
②
标记)和烽火村地区(图4和图5中的
③
标记)，这些区域均被两种方法识别为中高风险区。
[0103]
另外，相较于传统的基于层次分析法(ahp)，本发明对岩溶地表塌陷风险预测的精度更高。基于二元逻辑回归模型获得和较未引入基于多时相insar的地表形变速率因子的ahp法分别提高了12.5％和18.9％。本发明中的二元逻辑回归模型和传统ahp法
得到的受试者工作特征曲线的面积(auc)分别为0.911和0.812，也有所提高。
[0104]
本发明能够完成形变监测、诱因识别、灾害评估全过程的定量信息获取，为城市岩溶地表塌陷预警和治理提供技术支撑。
[0105]
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。