一种北斗与InSAR融合的地表形变监测方法与流程

本发明涉及大地测量与地表形变监测，具体为一种北斗与insar融合的地表形变监测方法。

背景技术：

1、围海造陆是沿海地区拓展土地空间的重要工程途径。围海造陆区域普遍分布厚层海相淤泥质软土，上层叠加新吹填的松散土层，在土体自重固结与长期蠕变作用下，会产生持续时间长、空间差异显著的地面沉降，对区域工程建设安全、地质灾害防控存在不利影响，因此开展高精度、高时空分辨率的地表形变监测具有重要的现实意义。

2、当前地表形变监测的常用技术手段包括三类：第一类为水准测量，其垂直形变测量精度较高，但存在监测点布设稀疏、外业作业成本高、观测周期长的局限，难以满足大范围面状监测需求；第二类为北斗卫星导航系统连续观测，可获取离散基准站的高精度三维位移时间序列，但空间分辨率受测站布设密度限制，无法得到全覆盖的面状形变结果；第三类为时序合成孔径雷达干涉测量技术，可实现大范围、高空间分辨率的面状形变观测，但观测值为沿雷达视线方向的一维形变，且易受大气延迟、轨道偏差等系统误差影响；围海造陆区临近海面，大气水汽时空变化剧烈，大气延迟误差的干扰尤为突出。

3、为获取三维地表形变场，现有技术中已出现合成孔径雷达干涉测量与卫星导航定位数据融合的监测方案，此类方案结合通用应力应变物理模型，采用方差分量估计迭代确定各类观测值的权重，通过最小二乘平差求解地表三维形变场。但此类方案应用于围海造陆场景时存在三方面不足：一是所采用的应力应变模型为通用地球物理模型，依赖杨氏模量、泊松比等弹性力学参数，无法适配围海造陆区软土固结蠕变时间效应显著、不同吹填分区土体物理力学性质差异大的特点，模型适用性不足；二是围海造陆区的专属工程地质信息未被利用于约束形变场解算，场景专属数据的价值未得到有效利用；三是现有围海造陆沉降物理机制的研究成果多停留在定性规律总结层面，未转化为可与大地测量平差模型耦合的定量约束矩阵。

技术实现思路

1、为了解决上述背景技术中存在的问题，本发明提出一种北斗与insar融合的地表形变监测方法，以解决围海造陆区地表形变监测场景下通用应力应变模型对软土固结蠕变特性适配性不足、专属工程信息未转化为可计算平差约束、物理约束机制未充分利用的问题，实现多源形变观测数据与围海造陆专属沉降物理规律的耦合联合平差，达成高精度、高空间分辨率、结果稳定可靠的围海造陆区三维地表形变监测。

2、为实现上述目的，本发明采用以下技术方案：

3、一种北斗与insar融合的地表形变监测方法，包括以下步骤：

4、s1、获取覆盖围海造陆区的多源形变监测数据，所述多源形变监测数据包括北斗卫星导航系统基准站观测数据、多轨道合成孔径雷达影像数据以及水准测量数据；

5、s2、对所述多源形变监测数据分别进行预处理，对应得到北斗三维形变观测量、各成像几何下的合成孔径雷达干涉测量视线向形变观测量以及水准垂直形变观测量；

6、s3、对各观测量进行时空基准统一，根据合成孔径雷达影像的成像几何参数建立合成孔径雷达干涉测量视线向形变与地表三维形变分量之间的几何投影关系；分别构建北斗观测方程、各轨道合成孔径雷达干涉测量观测方程以及水准观测方程；

7、s4、构造围海造陆区的三类物理机制约束方程，三类物理机制约束方程包括固结-时间约束、地层压缩分层约束以及各向异性空间梯度约束；将三类物理机制约束方程与所述北斗观测方程、各轨道合成孔径雷达干涉测量观测方程、水准观测方程采用统一残差形式联立，共同构建联合平差模型；

8、s5、在所述联合平差模型中，将北斗三维形变观测量、不同卫星和不同轨道方向的合成孔径雷达干涉测量视线向形变观测量、水准垂直形变观测量分别设置为独立观测组，将三类物理机制约束方程对应设置为独立约束组；

9、s6、构建分组随机模型，采用方差分量估计方法迭代求解所述联合平差模型，迭代过程中同步更新各观测组与各约束组的权重，直至满足收敛条件；

10、s7、根据收敛后的解算结果输出三维形变最优估值，生成围海造陆区三维地表形变场。

11、具体地，步骤s2中对所述北斗卫星导航系统基准站观测数据的预处理过程包括：

12、对北斗卫星导航系统基准站观测数据进行基线解算与网平差处理，得到国际地球参考框架下的三维坐标时间序列；

13、对所述三维坐标时间序列依次执行粗差剔除、阶跃修正以及周期项分离处理，得到北斗三维形变观测量。

14、具体地，步骤s2中对所述多轨道合成孔径雷达影像数据的预处理过程包括：

15、采用时序合成孔径雷达干涉测量方法对多轨道合成孔径雷达影像进行处理，所述时序合成孔径雷达干涉测量方法包括永久散射体法、小基线集法；处理过程依次包括干涉处理、相位解缠、大气相位改正以及形变量解算，得到各成像几何下的合成孔径雷达干涉测量视线向形变观测量。

16、具体地，步骤s4中所述固结-时间约束的构造过程为：

17、按照吹填完成时间将围海造陆区划分为若干个成陆分区，获取各成陆分区的成陆距今时间，约束第个成陆分区的垂直沉降速率，约束公式为：

18、；

19、其中，为第个成陆分区的垂直沉降速率，为成陆距今时间为零时的初始参考沉降速率，为随成陆距今时间增大而衰减的时间衰减函数。

20、具体地，所述时间衰减函数为分段函数，具体形式满足：

21、当时：

22、；

23、当时：

24、；

25、当时：

26、；

27、其中，分界点处函数值满足；为主固结阶段时间阈值，为次固结阶段时间阈值，为主固结阶段结束时的衰减系数，为指数衰减时间常数，为双曲线衰减系数。

28、具体地，步骤s4中所述地层压缩分层约束的构造过程为：

29、根据钻孔柱状图资料与土工试验资料，确定各监测单元下方各压缩土层的层厚参数与体积压缩系数，结合各压缩土层的有效应力增量，计算第个压缩土层的沉降贡献权重，计算公式为：

30、；

31、其中，为第个压缩土层的层厚，为第个压缩土层的体积压缩系数，为第个压缩土层的有效应力增量；

32、以所述沉降贡献权重作为地层压缩分层约束的系数矩阵元素，构造关于监测单元垂直形变参数的约束残差方程，使监测单元垂直形变在各压缩土层间的沉降贡献分配比例与相一致；所述约束残差方程作为地层压缩分层约束方程纳入所述联合平差模型。

33、具体地，步骤s4中所述各向异性空间梯度约束的构造过程为：

34、针对相邻的第个监测单元与第个监测单元，根据二者的空间距离、成陆年龄差、地层相似度构造各向异性邻接权重，权重计算公式为：

35、；

36、其中，为第个监测单元与第个监测单元之间的空间距离，为距离衰减特征尺度参数，为成陆时序方向权重系数，为第个监测单元的成陆年龄，为第个监测单元的成陆年龄，为第个监测单元与第个监测单元之间的成陆年龄差值，为研究区内最大成陆年龄差值，为符号函数，为第个监测单元与第个监测单元之间的地层相似度因子；

37、基于所述各向异性邻接权重构造差分约束算子，并将所述差分约束算子作为各向异性空间梯度约束方程并纳入所述联合平差模型；沿成陆时序方向相邻的监测单元对应设置梯度允许规则，跨越施工分区边界或地质分界的相邻监测单元对应降低邻接权重取值。

38、具体地，步骤s6中所述分组随机模型的方差-协方差矩阵形式为：

39、；

40、其中，为观测误差向量的方差-协方差矩阵，为独立分组的总数量，所述独立分组包括所有观测组与所有约束组；为第个独立分组的单位权方差分量，为第个独立分组对应的协因数矩阵；

41、采用helmert型方差分量估计方法迭代求解，方差分量估计方程的形式为：

42、；

43、其中，为方差分量系数矩阵，矩阵元素由各独立分组的协因数矩阵与残差投影矩阵做迹运算得到；为各独立分组单位权方差分量的估值向量；为残差二次型向量，向量元素为各独立分组对应的残差二次型；所述残差投影矩阵由平差设计矩阵与当前迭代的权矩阵计算得到。

44、具体地，步骤s3中所述北斗三维形变与视线向形变的几何投影关系满足：

45、；

46、其中，为合成孔径雷达干涉测量视线向形变量，为东西方向形变投影系数，为南北方向形变投影系数，为垂直方向形变投影系数，三个投影系数由合成孔径雷达影像的航向角与入射角计算得到；为东西方向形变分量，为南北方向形变分量，为垂直方向形变分量。

47、具体地，步骤s6中所述收敛条件为：相邻两次迭代得到的所有独立分组单位权方差分量估值的相对变化最大值小于预设收敛阈值，或迭代次数达到预设最大迭代次数；

48、步骤s7的解算结果输出的最优估值包括各监测点的三维形变速率与合成孔径雷达系统误差改正参数。

49、综上所述，本发明的有益技术效果为：

50、1.场景适配性更强：本发明构造固结-时间约束、地层压缩分层约束、各向异性空间梯度约束三类围海造陆专属物理机制约束，替代通用应力应变模型，可适配软土固结蠕变的沉降演化规律，充分利用吹填时序、地质勘察等场景专属工程信息。

51、2.约束定权更合理：本发明将三类物理约束作为独立约束组，与各类观测组同步纳入方差分量估计迭代流程，约束强度由观测数据自适应确定，可降低人工设定约束系数的主观性，平衡物理规律与实测数据的可信度。

52、3.解算稳定性更好：本发明将不同卫星、不同轨道的合成孔径雷达观测分设为独立观测组，联合北斗三维观测与水准垂直观测，并将雷达系统误差参数纳入平差解算，可增强三维形变分解能力，减弱近海大气延迟、轨道偏差等系统误差的干扰。

53、4.监测覆盖度更高：本发明融合高空间分辨率的合成孔径雷达面状观测与北斗、水准高精度点式观测，可生成高密度全覆盖三维地表形变场，可减少围海造陆区监测盲区，在保证监测精度的同时降低外业作业成本。