基于arima与lstm的海岸带地面沉降预测方法
技术领域
1.本发明涉及深度学习领域，更具体的说是涉及一种基于arima与lstm 的海岸带地面沉降预测方法。
背景技术：2.海陆过渡带是陆地与海洋之间的交界地带，其中包含了复杂的物理过程、化学过程、生物过程以及地质过程等，因而它成为了地球表面动态环境最复杂多变的区域之一。目前海陆过渡带区域有数十亿人生活着，是世界上最繁荣发达的区域之一，其演变过程与人类的可持续发展直接相关联，因而该区域也受到全球各个沿海国家的高度关注。
3.地面沉降是指地下支撑物的移动导致地面标高损失的一种地质灾害，它是自然因素与人为因素共同影响的地质过程。其中自然因素包括构造下沉、地震及火山活动、气候变化等，而人为因素包括开采地下水、油气资源、矿产资源等。海陆过渡带区域的长期地面沉降的累积，会导致洪涝灾害的加剧、建筑物寿命的缩短以及交通网络的损失等，这些都会严重威胁到人们的生命财产安全。因而能否有效监测地面沉降并给出相应的预警信息，是十分必要的。
4.地面沉降监测方法主要有水准测量、gps(global positioning system,全球定位系统)测量以及insar(interferometric synthetic aperture radar,合成孔径干涉雷达)测量等。其中insar技术具有全天时、全天候、高精度、实时性、连续性等优点，广泛应用于陆地地面沉降的监测。
5.而对于地面沉降预测方法，主要有基于物理机制的方法、基于数理统计的方法以及基于机器学习的方法。其中基于物理机制的方法需要大量的实地监测数据，提出基于岩石、水文特性演化的物理模型，并对大量复杂参数进行求解，并且受到时效性的严重影响；基于数理统计的方法则是根据大量历史测量数据得到统计规律，并进一步进行分析；基于机器学习的方法则是从更高的数据维度提取现有历史数据的特征，进而得到对现有数据的反馈，不会受限于物理参数的影响。
6.上述预测方法，其本质是对于时间序列的预测问题。而在时序预测方面，主要的模型包括：传统的arima模型(auto regressive integrated movingaverage model，差分整合移动平均自回归模型)，机器学习的支持向量机、贝叶斯网络、矩阵分解、高斯过程、深度学习等。
7.目前，在对insar数据得到的地面沉降时间序列经过单一预测模型处理后，得到的预测效果并不理想，因此，如何解决这一问题是本领域技术人员亟需研究的。
技术实现要素：8.有鉴于此，本发明提供了一种混合了arima和深度学习中的lstm的预测方法，即对insar得到的形变量时间序列与arima方法得到的预测时间序列作差，然后利用lstm对得到的时间序列差值进行学习，并把lstm 模型预测的结果加在arima预测数据上，就得到了原
始insar形变量时间序列的预测结果。
9.为了实现上述目的，本发明采用如下技术方案：
10.一种基于arima与lstm的海岸带地面沉降预测方法，包括以下步骤：
11.insar的slc影像数据经过预处理得到干涉图；
12.根据干涉图利用stamps得到第一时间序列；
13.基于第一时间序列，利用arima模型进行初步的预测，得到预测时间序列；
14.所述预测时间序列与所述第一时间序列做差，得到第二时间序列；
15.利用lstm模型对第二时间序列进行学习，得到第三时间序列；
16.将第三时间序列与预测时间序列相加，得到原始insar形变量时间序列的预测结果。
17.可选的，干涉图的相位组成如式(1)所示：
[0018][0019]
其中，代表了地形变化相位，代表了大气延迟相位，代表了残余轨道误差相位，代表了残余视向角误差相位，代表了噪声相位。
[0020]
可选的，根据干涉图利用stamps得到第一时间序列的具体步骤如下：
[0021]
stamps利用三维解缠算法，在时间上计算每个ps点的相位差异；
[0022]
设置参考点用最小二乘法在空间上进行解缠，得到第一时间序列。
[0023]
可选的，arima模型具体为：
[0024][0025]
其中，x
t
是t时刻的随机变量，αi以及βi是权重系数，p以及q表示时间窗长度，ε
t
则代表了t时刻的白噪声。
[0026]
可选的，所述lstm模型包括多个细胞单元，
[0027]
其中，x
t
表示第t时刻该细胞单元的输入，c
t
表示t时刻该细胞单元的状态，h
t
表示该细胞单元的隐状态输出；
[0028]
每个细胞单元中含有三个控制门，分别是遗忘门f、输入门i以及输出门 o；遗忘门决定上一时刻的细胞单元状态c
t-1
有多少成分保存到c
t
中；输入门决定当前时刻输入x
t
有多少成分保存到c
t
之中；输出门决定c
t
中有多少成分被传递到h
t
之中；故其前向计算方法如下所示：
[0029]ft
＝σ(wf·
[h
t-1
,x
t
]+bf)#(5)
[0030]it
＝σ(wi·
[h
t-1
,x
t
]+bi)#(6)
[0031]ct
＝f
tct-1
+i
t
tanh(wc·
[h
t-1
,x
t
]+bc)#(7)
[0032]ot
＝σ(wo·
[h
t-1
,x
t
]+bo)#(8)
[0033]ht
＝o
t
tanh(c
t
)#(9)
[0034]
式中，w和b分别为对应的权重系数矩阵和偏置项；σ和tanh分别为 sigmoid函数和双曲正切激活函数。
[0035]
可选的，所述第二时间序列求解过程如下：
[0036]
利用stamps方法得到若干个ps点，每个ps点对应的形变量时间序列记作{y
t
}，其长度为t。
[0037]
由最大似然估计得到arima(p,d,q)的最优参数，进而得到了预测序列那么将两者作差得到序列：
[0038][0039]
可选的，所述lstm模型的损失函数求解过程如下：
[0040]
采用adam优化算法，设定损失函数为绝对误差：
[0041][0042]
其中，测试集的长度为t
test
，训练用时间序列长度l，标准化的残差序列为为测试集序列。
[0043]
经由上述的技术方案可知，与现有技术相比，本发明公开提供了一种基于arima与lstm的海岸带地面沉降预测方法，相对单一的预测方法得到的精度指标有所提升。
附图说明
[0044]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
[0045]
图1为本发明的整体流程示意图；
[0046]
图2为本发明的lstm隐藏层细胞结构图；
[0047]
图3为本发明的主副影像时间基线与空间垂直基线关系图；
[0048]
图4(a)-图4(b)为本发明年平均形变速率；
[0049]
图5(a)为总形变量图；图5(b)年平均形变量分布直方图；
[0050]
图6(a)-图6(b)选取的ps点；
[0051]
图7(a)-图7(b)为时序预测结果对比图。
具体实施方式
[0052]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0053]
本发明实施例公开了一种基于arima与lstm的海岸带地面沉降预测方法，是一种混合了arima和深度学习中的lstm(long short-term memory, 长短期记忆单元)的预测方法，即对insar得到的形变量时间序列与arima 方法得到的预测时间序列作差，然后利用lstm对得到的时间序列差值进行学习，并把lstm模型预测的结果加在arima预测数据上，就得到了原始 insar形变量时间序列的预测结果。具体如图1所示，包括以下步骤：
[0054]
insar的slc影像数据经过预处理得到干涉图；
[0055]
根据干涉图利用stamps得到第一时间序列；
[0056]
基于第一时间序列，利用arima模型进行初步的预测，得到预测时间序列
[0057]
预测时间序列与第一时间序列做差，得到第二时间序列；
[0058]
利用lstm模型对第二时间序列进行学习，得到第三时间序列；
[0059]
将第三时间序列与预测时间序列相加，得到原始insar形变量时间序列的预测结果。
[0060]
其中，insar的slc(singlelookcomplex,单视复数)影像数据经过欧洲航天局的snap(sentinelapplicationplatform)软件进行预处理得到干涉图，从而进一步利用stamps(stanfordmethodforpersistentscatterers)得到特定点的时间序列。后续则利用arima方法进行初步的预测，与原始序列作差，并利用深度学习对该差值序列进行学习，改进了arima方法的误差项，在此基础上进行预测。
[0061]
stamps方法是由hooper发展而来的一种从合成孔径雷达采集的时间序列中提取地面形变量的算法软件包。相比于传统的ps-insar(persistentscattererinterferometricsyntheticapertureradar,永久散射体合成孔径干涉雷达)选取永久散射体点的方法，stamps方法的选取规则可以识别散射强度较低但稳定性较高的目标，因而可以改善在非城镇地区传统ps-insar方法选取ps点不足的情况，更加适用于海陆过渡带区域。
[0062]
干涉相位图的相位组成如式(1)所示：
[0063][0064]
其中，代表了地形变化相位，代表了大气延迟相位，代表了残余轨道误差相位，代表了残余视向角误差相位，代表了噪声相位。所有这些相位纠缠在一起构成了干涉相位图的相位信息。
[0065]
stamps利用了三维解缠算法，在时间上计算每个ps点的相位差异，然后设置参考点用最小二乘法在空间上进行解缠
[26]
，相比于传统的方法精度更高。
[0066]
解缠后的相位通过时间上的高通滤波，空间上的低通滤波，得到大气延迟相位并减去，从而矫正了大气误差。最后经过相位反演得到形变信息。
[0067]
经典的时间序列预测模型包括自回归模型(autoregressive,ar)，移动平均模型(movingaverage,ma)，以及两者相结合的自回归移动平均模型(autoregressivemovingaverage,arma)
[8]
，进一步利用差分运算，就得到了arima模型。ar(p)模型如式(2)所示，ma(q)模型如式(3)所示，arma(p,q)模型如式(4)所示：
[0068][0069]
[0070][0071]
其中x
t
是t时刻的随机变量，αi以及βi是权重系数，p以及q表示时间窗长度，ε
t
则代表了t时刻的白噪声。
[0072]
上述三种模型是弱平稳的，适用于平稳时间序列预测，但涉及到非平稳时间序列，则需要结合差分运算得到平稳序列后再进行预测，即arima模型。经过实验比对，实际上arima模型的预测误差稳定性较好，因而可以用于后面与insar数据作差，进而可以使利用lstm进行残差序列学习结果的收敛性得到保障。
[0073]
循环神经网络(recurrent neural network,rnn)是一种能够将当前的数据与之前的数据建立联系的一种深度神经网络。但由于非线性函数的存在，其中离现在数据较远的数据会受到梯度衰减的影响，进而造成梯度消失问题；但另一方面，如果离当前数据较近的数据记录过多的情况下，也会造成梯度爆炸问题。
[0074]
为解决rnn的潜在问题，lstm模型被提出。它继承了rnn的记忆性，并增加了门控制单元对信息进行取舍，然后经过反向传播的梯度调整，从而完成模型的建立。其中以隐藏层最为重要，它由多个细胞单元组成，其细胞单元结构如图2所示。
[0075]
其中x
t
表示第t时刻该细胞单元的输入，c
t
表示t时刻该细胞单元的状态， h
t
表示该细胞单元的隐状态输出。
[0076]
每个细胞单元中含有三个控制门，分别是遗忘门f、输入门i以及输出门 o。其中，遗忘门决定了上一时刻的细胞单元状态c
t-1
有多少成分保存到c
t
中；输入门决定了当前时刻输入x
t
有多少成分保存到c
t
之中；输出门决定了c
t
中有多少成分被传递到h
t
之中。故其前向计算方法如下所示：
[0077]ft
＝σ(wf·
[h
t-1
,x
t
]+bf)#(5)
[0078]it
＝σ(wi·
[h
t-1
,x
t
]+bi)#(6)
[0079]ct
＝f
tct-1
+i
t
tanh(wc·
[h
t-1
,x
t
]+bc)#(7)
[0080]ot
＝σ(wo·
[h
t-1
,x
t
]+bo)#(8)
[0081]ht
＝o
t
tanh(c
t
)#(9)
[0082]
式中，w和b分别为对应的权重系数矩阵和偏置项；σ和tanh分别为 sigmoid函数和双曲正切激活函数。
[0083]
lstm的模型训练过程与循环神经网络一致，误差项δ的反向传播也是包括两个方向：一个是沿时间的反向传播，计算每个时刻的误差项；一个是将误差项向上一层神经元传播。然后根据误差项计算每个权重的梯度，进而利用梯度优化算法更新权重完成训练。
[0084]
利用lstm对arima的残差序列进行学习的预测方法：假设利用stamps 方法得到若干个ps点，每个ps点对应的形变量时间序列记作{y
t
}，其长度为t。
[0085]
由最大似然估计得到arima(p,d,q)的最优参数，进而得到了预测序列
[0086]
那么将两者作差得到序列：
[0087][0088]
由于两者作差得到的时间序列为平稳的，因而该序列不必在lstm训练前做差分处
理，因而只需做标准化处理即可。设标准化的残差序列为那么，
[0089][0090]
然后利用lstm模型对上述标准化残差序列进行训练。设测试集的长度为t
test
，训练用时间序列长度l。实践表明，单步预测的误差是最小的，因而对时间序列进行标注，使其成为监督数据时，只需要最后一个作为标签即可。故训练用的样本输入为(l-1)个。
[0091]
lstm的主要参数包括训练用时间序列长度l，神经网络层数k，隐藏层神经元个数n，训练次数epochs以及学习率η。可以采用网格搜索法得到最优参数，同时采用随训练次数增加逐步衰减的学习率。
[0092]
确定lstm参数后，就可以开始对其进行训练了。采用adam优化算法
[29]
，设定损失函数为绝对误差：
[0093][0094]
对数据进行反向标准化，得到同时lstm模型的权重与偏置项构成了一个训练好的网络。
[0095]
将测试集中的时序样本输入到上述网络中，就得到了预测样本然后将其加到上，就得到了{y
t
}的一个预测。
[0096]
最后，通过计算网络训练及网络预测过程中的输出结果与残差序列样本的偏差，定量地给出模型训练和预测的精度
[7]
：
[0097][0098][0099][0100]
[0101]
式中，rmse(root mean squard error)表示均方根误差，mae(meanabsolute error)表示平均绝对误差，map(mean absolute practice)表示平均实际绝对量，μ表示平均预测精度。
[0102]
本实施例选取杭州湾部分地区(30
°
n-30
°
30
′
n,120
°
45
′
e-120
°
55
′
e)为研究区域。由于水准数据采用的是南边的绍兴市上虞区数据，故通过相对稳定的嘉绍大桥沟通南北两岸，以验证北边的形变数据的准确度。
[0103]
以2017年1月至2019年12月获取的sentinel-1a影像的上升轨道，vv 极化方式的slc影像为遥感数据源，利用ps-insar技术选取2018年7月16 日的影像作为主影像，共得到配准的主副影像干涉图83幅，其主副影像的时间基线以及空间垂直基线的关系如图3所示。
[0104]
然后利用stamps方法筛选一共得到274926个高相干点(永久散射体)，并进一步得到沿雷达视线方向的地形变化时间序列，其中每个点拥有83个形变记录，其年平均形变速率及其叠加在卫星图上的情况如图4(a)-图4(b)所示。从图中可以看到，沿杭州湾两岸处地形整体呈现相对下降的趋势，而处在更靠近陆地的地方则在整体上呈现相对上升的趋势。
[0105]
对所有ps点的形变量数据、年平均形变量数据做统计分析，3年时间内，该地区的形变量在-71.98mm～90.25mm范围内，年平均形变量为
ꢀ‑
23.65mm/year～10.93mm/year。图5(a)代表了所有ps点的形变量数据分布情况，可以看到其中的绝大多数形变量数据集中在-20mm～20mm之间。图5(b) 为年平均形变量直方图分布，可以看到平均形变速率大多集中在
ꢀ‑
10mm/year～5mm/year的范围内。
[0106]
针对地面沉降漏斗区域的监测情况，选取两处地面沉降较为明显的ps点，分别利用arima、lstm以及arima-lstm三种方法进行进一步预测以及对比。用深蓝色框标注了选取的两个点，并分别命名为a和b如图6(a)-图6(b) 所示。图7则分别表明了两个点各自对应的预测方法对比结果。其rmse、 mae以及μ指标结果如表1所示。
[0107]
这里采用的是最大似然估计得到arima的最优参数，网格搜索法得到 lstm的最优参数。arima方法的(p,d,q)参数最终选择为(5,1,0)；lstm方法则按照时间先后顺序划分训练集以及测试集，2017年和2018年的数据作为训练集，2019年作为测试集，epochs次数为100，隐藏神经元个数为10个； arima-lstm方法的以2019年的前2/3年数据作为训练集，后1/3年为测试集，epochs次数为50，隐藏层神经元个数为12。
[0108]
从图7(a)-图7(b)中可以看到arima和lstm方法均会出现预测的结果偏后出现的情况，这源自于两种方法对过去信息的一种记忆方式。假设某时刻的实际值为y
t
,前一个时刻的实际值为y
t-1
,arima方法在对t时刻的值进行预测的时候，往往与t时刻的前几个时刻的值有一定的关系，当信息数量不是很大的时候，前一个时刻的实际值y
t-1
所对应的权重较大，因而其预测结果会更偏向于y
t-1
，反映在图中的结果就会显示为滞后。同样的道理也适用于 lstm方法，因为本质上两者都是利用函数模型去逼近。而arima和lstm 相结合的方法则会改善这种情况，并且在某些情况下会得到比较接近insar 原始形变量时间序列的预测点，例如b点在2019年的最后预测值。
[0109]
表1 三种方法的精度指标对比
[0110][0111]
从上表中可以看到，arima-lstm方法相比于单一的arima方法或者 lstm方法有着明显的提升。rmse至少减小了48.75％，mae至少减小了 24.84％，μ至少提升了15.19％。借助于预测结果的可接受性，可以将lstm 训练的模型合理外推至没有相应insar数据的未来年份，并及时对地面沉降结果进行适当预警。
[0112]
本发明通过建立lstm模型，对insar形变量时间序列与arima预测形变量时间序列之间的差值序列进行预测，得到的均方根误差、平均绝对误差均比单一的lstm、arima模型预测的小，且平均预测精度更高。在数值上具体体现为rmse至少减小了48.75％，mae至少减小了24.84％，μ至少提升了15.19％。
[0113]
而对于海陆过渡带地面沉降监测而言，及时发现出现的严重沉降区域是必要的。并且如果能够在未来一段时间内没有或者难以获取insar数据监测的情况下，合理的预测算法将具有重大的意义，它为潜在的地质灾害的预警提供相应的帮助。就目前的情况来看，通过arima预测只能针对某个ps点。对大量的ps点而言，通过lstm学习arima与原始时间序列差值的方法高度依赖于arima模型的建立，因而两者结合的方法有一定局限性。后续的工作可以寻找大范围预测中能够生成较为稳定的差值序列的方法，进而用lstm 模型进行拟合，得到更加精确的预测结果。
[0114]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0115]
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。
技术特征：1.一种基于arima与lstm的海岸带地面沉降预测方法，其特征在于，包括以下步骤：insar的slc影像数据经过预处理得到干涉图；根据干涉图利用stamps得到第一时间序列；基于第一时间序列，利用arima模型进行初步的预测，得到预测时间序列；所述预测时间序列与所述第一时间序列做差，得到第二时间序列；利用lstm模型对第二时间序列进行学习，得到第三时间序列；将第三时间序列与预测时间序列相加，得到原始insar形变量时间序列的预测结果。2.根据权利要求1所述的一种基于arima与lstm的海岸带地面沉降预测方法，其特征在于，干涉图的相位组成如式(1)所示：其中，代表了地形变化相位，代表了大气延迟相位，代表了残余轨道误差相位，代表了残余视向角误差相位，代表了噪声相位。3.根据权利要求1所述的一种基于arima与lstm的海岸带地面沉降预测方法，其特征在于，根据干涉图利用stamps得到第一时间序列的具体步骤如下：stamps利用三维解缠算法，在时间上计算每个ps点的相位差异；设置参考点用最小二乘法在空间上进行解缠，得到第一时间序列。4.根据权利要求1所述的一种基于arima与lstm的海岸带地面沉降预测方法，其特征在于，arima模型具体为：其中，x
t
是t时刻的随机变量，α
i
以及β
i
是权重系数，p以及q表示时间窗长度，ε
t
则代表了t时刻的白噪声。5.根据权利要求1所述的一种基于arima与lstm的海岸带地面沉降预测方法，其特征在于，所述lstm模型包括多个细胞单元，其中，x
t
表示第t时刻该细胞单元的输入，c
t
表示t时刻该细胞单元的状态，h
t
表示该细胞单元的隐状态输出；每个细胞单元中含有三个控制门，分别是遗忘门f、输入门i以及输出门o；遗忘门决定上一时刻的细胞单元状态c
t-1
有多少成分保存到c
t
中；输入门决定当前时刻输入x
t
有多少成分保存到c
t
之中；输出门决定c
t
中有多少成分被传递到h
t
之中；故其前向计算方法如下所示：f
t
＝σ(w
f
·
[h
t-1
，x
t
]+b
f
)#(5)i
t
＝σ(w
i
·
[h
t-1
，x
t
]+b
i
)#(6)c
t
＝f
t
c
t-1
+i
t
tanh(w
c
·
[h
t-1
，x
t
]+b
c
)#(7)o
t
＝σ(w
o
·
[h
t-1
，x
t
]+b
o
)#(8)h
t
＝o
t
tanh(c
t
)#(9)式中，w和b分别为对应的权重系数矩阵和偏置项；σ和tanh分别为sigmoid函数和双曲正切激活函数。
6.根据权利要求1所述的一种基于arima与lstm的海岸带地面沉降预测方法，其特征在于，所述第二时间序列求解过程如下：利用stamps方法得到若干个ps点，每个ps点对应的形变量时间序列记作{y
t
}，其长度为t。由最大似然估计得到arima(p，d，q)的最优参数，进而得到了预测序列那么将两者作差得到序列：7.根据权利要求1所述的一种基于arima与lstm的海岸带地面沉降预测方法，其特征在于，所述lstm模型的损失函数求解过程如下：采用adam优化算法，设定损失函数为绝对误差：其中，测试集的长度为t
test
，训练用时间序列长度l，标准化的残差序列为为测试集序列。
技术总结本发明公开了一种基于ARIMA与LSTM的海岸带地面沉降预测方法，涉及深度学习领域。本发明包括以下步骤：InSAR的SLC影像数据经过预处理得到干涉图；根据干涉图利用StaMPS得到第一时间序列；基于第一时间序列，利用ARIMA模型进行初步的预测，得到预测时间序列；所述预测时间序列与所述第一时间序列做差，得到第二时间序列；利用LSTM模型对第二时间序列进行学习，得到第三时间序列；将第三时间序列与预测时间序列相加，得到原始InSAR形变量时间序列的预测结果。本发明相对单一的预测方法得到的精度指标有所提升。指标有所提升。指标有所提升。
技术研发人员：吴自银 罗孝文 朱宝 尚继宏 赵荻能 陈小伦
受保护的技术使用者：自然资源部第二海洋研究所
技术研发日：2022.07.09
技术公布日：2022/11/1
