一种表征地面沉降时序演变的预测算法的制作方法

1.本发明涉及地质灾害预防技术领域，特别是涉及一种表征地面沉降时序演变的预测算法。


背景技术：

2.地面沉降是指在自然因素或人为因素的影响下，由于地下松散地层固结压缩，导致表面标高降低的地质现象。当地面沉降发生严重时，会诱发一系列的地质灾害，如地面塌陷、地裂缝、海水倒灌等，形成灾害链。区域不均匀沉降威胁了城市居住安全，会造成地下管道破裂、楼房开裂等问题。因此高效快速的预测地面沉降有助于对地面沉降进行精准防控。
3.目前，地面沉降预测模型主要包括物理模型、数据驱动模型和geo
‑
ai智能模型等。其中，物理模型需要详尽的基础地质、水文地质参数，准确刻画复杂的地层结构，才能具有较高的模拟精度。这类先验知识和参数信息难以获取且数据密度较低，计算效率低下，实际应用难度较高。数据驱动模型(如灰色模型(grey model,gm)及改进型gm等)从时序数据出发，挖掘数据特征进而实现短期预测，具有较好的结果，但难以进行长期预测达到防控的目的。
4.针对上述问题，本发明提出一种基于麻雀搜索算法改进的bp神经网络算法具有寻优能力强、收敛速度快的优点，弥补了原始bp神经网络收敛速度慢等缺点。地面沉降是多因素叠加作用形成的一种地质现象，各因素之间为非线性关系。人工智能模型能良好的解决非线性问题，因此，近五年运用机器学习开展地面沉降预测成为了前沿研究。本发明创新地将ssa
‑
bp算法引用到地面沉降预测研究中。首先运用空间分析、数理统计等技术方法，联合研究区地质条件和城市功能区划分等，厘定了研究沉降的最优格网尺度；其次，优化集成google earth engine(gee)、机器学习等新技术，联合北京通州地下水水位、可压缩层厚度、viirs夜间灯光、osm道路等多源空间数据，获得多因素驱动下的ssa
‑
bp沉降预测模型。能够为地面沉降灾害的防治提供依据，对合理规划城市发展、预防和控制城市地面沉降具有重要的应用价值。


技术实现要素：

5.本发明的目的是为了获取地面沉降在时间序列上的演变信息，将演变规律进行定量的描述，从而提出了一种表征地面沉降时序演变的预测算法，进行高效的地面沉降防治。
6.为了实现上述目的，本发明采用了如下技术方案：
7.一种表征地面沉降时序演变的预测算法，包括有以下步骤；
8.s1、获取覆盖研究区范围的长时序合成孔径雷达影像数据，利用时序合成孔径雷达干涉测量(ps
‑
insar)方法，获取地表视线向的形变信息，通过实测gps、水准点数据对insar监测结果进行校准和精度验证，将视线向的值转化为垂向，获取ps点的高精度垂向形变数据，通过克里金插值/反距离权重等插值法获取年均沉降速率；
9.s2、构建麻雀优化器，将麻雀优化器代入到bp神经网络代码中，得到基于麻雀搜索
算法优化的bp神经网络模型(ssa
‑
bp模型)；
10.s3、将多源地学数据通过空间分析、数理统计等技术方法进行处理，得到权重因子；
11.s4、将s1中计算所得的年均沉降速率和s3中得到的权重因子代入到s2中提到的基于麻雀搜索算法优化的bp神经网络模型中，构建基于麻雀搜索算法优化的bp神经网络的地面沉降预测模型(ssa
‑
bp地面沉降预测模型)；
12.s5、选取rmse、mae、r2作为精度参数，对s4中的得到的基于麻雀搜索算法优化的bp神经网络的地面沉降预测模型(ssa
‑
bp地面沉降预测模型)进行精度分析，验证模型的可行性，并将验证结果输出。
13.优选地，所述s2中提到的构建麻雀优化器，具体包括以下步骤：
14.a1、将n只麻雀组成的种群表示为：
[0015][0016]
式中，r表示待优化问题变量的维数，n为麻雀的数量；初始化种群，迭代次数，初始化捕食者和加入者比例；
[0017]
a2、所有麻雀的适应度值表示为：
[0018][0019]
式中，f表示适应度值；利用式(2)计算种群的适应度值，并对其进行排序；
[0020]
a3、依据ssa算法规则，每次迭代时发现者位置更新为：
[0021][0022]
式中，k代表当前迭代数；j＝1，2，3，
…
，d；iter
max
为常数，表示最大迭代次数；x
i,j
为第i只麻雀在第j维中的位置信息；α∈(0,1]为随机数；r2和st分别表示预警值和安全值；q服从正态分布的随机数；l表示一个1xd的矩阵，其中该矩阵内每个元素全部为1；利用式(3)更新捕食者位置；
[0023]
a4、加入者的位置更新描述为：
[0024]
[0025]
式中，x
p
为目前发现者最优位置，x
worst
为目前全局最差位置，a为一个1xd矩阵，其中每个元素为随机值1或
‑
1，有a

＝a
t
(aa
t
)
‑1；利用式(4)更新加入者位置；
[0026]
a5、假设意识到危险的麻雀占总数量的10％～20％，则有：
[0027][0028]
式中，x
best
为全局最优位置，β为控制参数，服从均值为0方差为1的正态分布随机数；k∈[
‑
1,1]的随机数，f
i
为目前麻雀个体的适应度值；f
g
和f
w
分别为目前全局最好和最差的适应度值；ε为最小的常数，避免分母为0；利用式(5)更新警戒者位置；
[0029]
a6、利用a2中提到的式(2)计算适应度值，并实时更新麻雀位置；
[0030]
a7、根据a6中所得的计算结果以及位置信息判断是否满足停止条件，满足则退出，输出结果，否则，重复执行a2
‑
a6所述操作。
[0031]
优选地，所述s1在处理雷达观测影像数据时，需用到雷达影像处理软件，所述雷达影像处理软件为stamps、sarproz和gamma。
[0032]
优选地，所述式(1)
‑
(5)的计算可以通过matlab、python、c 等计算获得，所述s3中提到的空间分析技术和数理统计技术中的计算通过arcgis、matlab、excel软件计算所得。
[0033]
优选地，所述s3中权重因子由基础地质、遥感数据处理得到，其具体处理方法包括有：
[0034]
b1、地下水水位数据处理方法：
[0035]
b1.1、将地下水等值线构建不规则三角网(tin)，具体包括：
[0036]
b1.1.1、等高线作为特征线建立tin：将等高线作为断裂线或结构线；使用等高线上的特征点，并将等高线段作为约束线段处理；
[0037]
b1.1.2、自动增加特征点及优化tin的方法：该方法实质仍将等高线离散化建立tin，但采用增加特征点的方式来消除tin中的平三角形，并使用优化tin的当时来消除不合理的三角化；特征点的增加需要利用一定的算法自动提取，这些算法原理大都基于原始等高线的拓扑关系；
[0038]
b1.1.3、等高线离散点直接生成tin的方法：直接将等高线离散化，然后利用常用tin的生成算法；
[0039]
b1.2、转化为栅格数据：通过插值法将不规则三角网(tin)转换为栅格；为输出中的每个像元指定高度或nodata值，具体取决于像元中心是否落在tin的插值区内；该方法可以采用linear或natural_neighbors；线性插值法可将tin三角形显示为平面；通过查找落在2d空间中的三角形并计算像元中心相对于三角形平面的位置来为每个输出像元指定高度(自然邻域法插值可产生比线性插值更平滑的结果；它对在每个输出像元中心周围所有方向上找到的最近tin结点使用基于区域的权重方案；tin中的硬隔断线将影响结果，因为自然邻域表面在跨越该硬隔断线时不是连续平滑的)；采样距离变量表示输出像元大小，可以采用两种方式之一进行定义；observations关键字后面是指示输出中所需列数的整数；基于输出范围和列数确定像元大小；默认情况下，使用observations 250；cellsize关键字
提供明确的像元大小，后面紧跟可以浮点数形式给出的所需大小；z因子变量用于转换输出栅格的z单位；输出高度将乘以该值；当且仅当x、y和z单位在输入tin的空间参考中定义，而且都是标准投影单位(例如，米和英尺)时，才会计算默认z因子变量以将z单位转换为与x和y单位相同的测量单位；如果x、y和z单位相同，则使用值1.0；如果未定义任何单位，则使用值1.0；
[0040]
b1.3、将30m分辨率的栅格数据转为点图层，在输出要素类中创建一个点；这些点将定位于它们所代表的像元的中心；nodata像元不会转换为点；输入栅格的像元大小可以任意，并且可以属于任何有效的栅格数据集；求交集取平均，获取各个格网的地下水水位数据信息；
[0041]
b2、可压缩层数据处理方法：
[0042]
将可压缩层厚度等值线数据转化为点数据，并进行空间idw插值，使用一组采样点的线性权重组合来确定像元值，权重是一种反距离函数，进行插值处理的表面具有局部因变量；主要依赖于反距离的幂值；幂参数可基于距输出点的距离来控制已知点对内插值的影响；幂参数是一个正实数，默认值为2；通过定义更高的幂值，可进一步强调最近点；因此，邻近数据将受到最大影响，表面会变得更加详细(更不平滑)；随着幂数的增大，内插值将逐渐接近最近采样点的值；指定较小的幂值将对距离较远的周围点产生更大影响，从而导致更加平滑的表面；由于反距离权重公式与任何实际物理过程都不关联，因此无法确定特定幂值是否过大；作为常规准则，认为值为30的幂是超大幂，因此不建议使用；如果距离或幂值较大，可能生成错误结果；可将所产生的最小平均绝对误差最低的幂值视为最佳幂值；从而获取研究区可压缩层组厚度栅格数据；
[0043]
b3、夜间灯光数据处理方法：
[0044]
将研究区的夜间灯光数据进行裁剪、下载并导出，导出数据为栅格图像，空间分辨率为30m；联合空间分析和数理统计技术，将栅格图像转化为点图层，取交集求平均，获取研究区各格网的灯光像元累计值；
[0045]
b4、道路数据处理方法：
[0046]
按照既有道路数据的分类结合中国道路的不同平均速度，得到不同道路的权重。
[0047]
与现有技术相比，本发明提供了一种表征地面沉降时序演变的预测算法，具备以下有益效果：
[0048]
本发明利用雷达影像数据，利用ps
‑
insar方法，获取区域地面沉降信息。根据时序沉降信息，与其它基础地质、遥感数据相结合，将区域地面上的沉降时序演变特征进行模拟，建立ssa
‑
bp地面沉降预测模型，更加准确的预测区域地面沉降，在应对地面沉降灾害时能够提出更加高效的防治措施，为城市灾害防治节约时间和经济成本，同时也可以作为城市规划的依据。
附图说明
[0049]
图1为本发明提出的一种表征地面沉降时序演变的预测算法的流程示意图；
[0050]
图2为本发明提出的一种表征地面沉降时序演变的预测算法的ps
‑
insar流程示意图。
具体实施方式
[0051]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。
[0052]
实施例1：
[0053]
请参阅图1，一种表征地面沉降时序演变的预测算法，包括有以下步骤；
[0054]
s1、获取覆盖研究区范围的长时序合成孔径雷达观测影像数据，假设影像获取时间为t0，t1,...t
n
。利用永久散射体合成孔径雷达干涉测量(ps
‑
insar)方法获取地表每一个ps点的形变信息，通过实测gps、水准点数据对形变数据进行校准和精度验证，获取稳定ps点的高精度视线向形变数据，其中包括每个ps点的年平均沉降速率v和观测时间的累积沉降量s0，s1，...s
n
，通过插值法计算获得年均沉降速率；
[0055]
s1在获取雷达观测影像数据时，需用到雷达影像处理软件，雷达影像处理软件为stamps、sarproz和gamma。
[0056]
s2、构建麻雀优化器，将麻雀优化器代入到bp神经网络代码中，得到基于麻雀搜索算法优化的bp神经网络模型(ssa
‑
bp模型)；
[0057]
s2中提到的构建麻雀优化器，具体包括以下步骤：
[0058]
a1、将n只麻雀组成的种群表示为：
[0059][0060]
式中，r表示待优化问题变量的维数，n为麻雀的数量；初始化种群，迭代次数，初始化捕食者和加入者比例；
[0061]
a2、所有麻雀的适应度值表示为：
[0062][0063]
式中，f表示适应度值；利用式(2)计算种群的适应度值，并对其进行排序；
[0064]
a3、依据ssa算法规则，每次迭代时发现者位置更新为：
[0065][0066]
式中，k代表当前迭代数；j＝1，2，3，
…
，d；iter
max
为常数，表示最大迭代次数；x
i,j
为第i只麻雀在第j维中的位置信息；α∈(0,1]为随机数；r2和st分别表示预警值和安全值；q服从正态分布的随机数；l表示一个1xd的矩阵，其中该矩阵内每个元素全部为1；利用式(3)更新捕食者位置；
[0067]
a4、加入者的位置更新描述为：
[0068][0069]
式中，x
p
为目前发现者最优位置，x
worst
为目前全局最差位置，a为一个1xd矩阵，其中每个元素为随机值1或
‑
1，有a

＝a
t
(aa
t
)
‑1；利用式(4)更新加入者位置；
[0070]
a5、假设意识到危险的麻雀占总数量的10％～20％，则有：
[0071][0072]
式中，x
best
为全局最优位置，β为控制参数，服从均值为0方差为1的正态分布随机数；k∈[
‑
1,1]的随机数，f
i
为目前麻雀个体的适应度值；f
g
和f
w
分别为目前全局最好和最差的适应度值；ε为最小的常数，避免分母为0；利用式(5)更新警戒者位置；
[0073]
a6、利用a2中提到的式(2)计算适应度值，并实时更新麻雀位置；
[0074]
a7、根据a6中所得的计算结果以及位置信息判断是否满足停止条件，满足则退出，输出结果，否则，重复执行a2
‑
a6操作。
[0075]
式(1)
‑
(5)的计算可以通过matlab、python、c 等计算获得，s3中提到的空间分析、数理统计技术中的计算通过arcgis、matlab、excel软件计算所得。
[0076]
s3、将多源地学数据通过空间分析、数理统计等技术方法进行处理，得到权重因子；
[0077]
s4、将s1中计算所得的年均沉降速率和s3中得到的权重因子代入到s2中提到的基于麻雀搜索算法优化的bp神经网络模型中，构建基于麻雀搜索算法优化的bp神经网络的地面沉降预测模型；
[0078]
s5、选取rmse、mae、r2作为精度参数，对s4中的得到的基于麻雀搜索算法优化的bp神经网络的地面沉降预测模型进行精度分析，证明模型的可行性，并将证明结果输出。
[0079]
本发明利用雷达影像数据，利用ps
‑
insar方法，获取区域地面沉降信息。依据时序沉降信息，与其它基础地质、遥感数据相结合，将区域地面沉降时序演变特征进行模拟，建立ssa
‑
bp地面沉降预测模型，更加准确的预测区域地面沉降，在应对地面沉降灾害时能够提出更加高效的防治措施，为城市灾害防治节约时间和经济成本，同时也可以作为城市规划的依据。
[0080]
实施例2：
[0081]
基于实施例1但有所不同之余在于，本发明先通过gamma软件处理2010
‑
2020年的radarsat
‑
2雷达遥感影像，利用水准点数据验证ps
‑
insar监测结果的精度，获取高精度视线向地面沉降信息，并将其转化为垂向位移信息。将北京市通州区部分区域作为研究区，利用监测方法获得地面沉降平均速率的变化，分析2010
‑
2020年研究区地面沉降的时序演变规律，并进行预测。主要步骤如下：
[0082]
步骤一：将2010
‑
2019年覆盖北京平原区的81景radarsat
‑
2数据利用gamma软件进行处理，获取北京市通州区2010
‑
2019年的视线向地面沉降信息，并将其转化为垂向。利用研究区内水准点监测数据对ps
‑
insar结果进行精度验证，获取高精度的地面沉降信息。
[0083]
步骤二：
[0084]
确定研究区范围，依据北京市通州行政区划、地下水位及地下水含水层边界，通过arcgis平台进行叠加分析，确定研究区范围。其次，综合考虑通州新城裂缝分布、地质构造单元划分、控制性详细规划等因素，将研究区划分不同分辨率的网格。
[0085]
步骤三：
[0086]
将ps
‑
insar中的累计沉降量进行反距离权重插值，获得2010年11月
‑
2020年1月的累计沉降量，通过计算，获得每年的平均沉降量。
[0087]
将可压缩层数据、地下水位数据、夜间灯光数据、道路osm数据通过gee平台及arcgis软件处理，得到权重因子。
[0088]
步骤四：
[0089]
将沉降速率和权重因子带入已建立的ssa
‑
bp沉降预测模型，通过下述参数，运行代码，得到预测数据。
[0090]
具体参数如表所示：
[0091][0092]
步骤五：
[0093]
对代码结果进行精度分析(rmse、mae、r2)，判断预测数据的准确性。
[0094]
以上，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。