一种用于大地测量速度场的改进双边滤波应变率计算方法

本发明涉及地质灾害监测，更具体地说，本发明涉及一种用于大地测量速度场的改进双边滤波应变率计算方法。

背景技术：

1、随着合成孔径雷达干涉测量技术和地面卫星定位观测技术的发展，基于升轨雷达影像序列和降轨雷达影像序列获取地表形变信息，并进一步形成区域形变速率场，已成为地壳形变监测和地震危险性评估中的常见技术路线。为了获得能够反映区域构造变形状态的应变率空间分布，现有技术通常先获取视线向形变结果或目标向形变速率结果，再基于形变速率场开展后续空间分析。

2、例如，现有公开号为：cn101770027a的专利申请公开了一种基于insar与gps数据融合的地表三维形变监测方法，通过gps点解算、sar图像处理、轨道误差改正、大气延迟反演、相位解缠以及insar与gps形变信息融合，获得地表三维形变结果。该类方案能够在一定程度上提高地表形变监测结果的精度和可靠性，但其主要侧重于多源观测信息的融合监测和形变结果获取。当后续处理对象进一步转为形变速率场时，前述融合过程并不直接解决速率场中的异常噪声抑制和边界保持问题。

3、又如，现有公开号为：cn112050725a的专利申请公开了一种融合insar和gps的地表形变监测方法，同样围绕sar图像预处理、轨道误差改正、大气延迟改正、相位解缠以及融合监测结果展开。该类方案说明现有技术已经认识到多源观测结果融合对于区域形变监测的重要性，但其重点仍在于形变监测结果本身的获取。当地表形变结果进一步转换为形变速率场并进入后续滤波处理阶段时，速率场中的局部离群噪声和断层附近速度不连续边界仍然可能同时存在。

4、又如，现有公告号为：cn115469308b的发明专利公告了一种多轨道insar震间形变速率场拼接方法、装置、设备及介质，通过获取gnss三维形变速率和sar数据，将gnss三维形变速率投影至目标向，对各轨道目标向形变速率进行参考基准统一、入射角校正和拼接，以获取大范围目标向形变速率场。该类方案能够改善多轨道速率场之间的统一性和覆盖范围，说明现有技术已经认识到区域速率场构建过程中多源数据协同和基准一致性的重要性。但是，该类方案主要解决的是不同轨道速率场之间的统一和拼接问题。当形变速率场进入后续滤波和平滑处理阶段时，若邻域像元加权运算不能同时兼顾异常噪声抑制和断层边界保持，则断层一侧的速度信息仍可能通过邻域权重传递到另一侧，形成跨断层信息混合。

5、再如，现有公告号为：cn109035152b的发明专利公告了一种合成孔径雷达图像非局部均值滤波方法，其通过改进权值测量方式和自适应退化参数函数，在滤除相干斑的同时保持纹理边缘信息。该类方案说明现有技术已经注意到滤波处理中噪声抑制与边缘保持之间的平衡问题。但该类方案主要针对sar图像本身的相干斑噪声抑制和纹理边缘保持。当处理对象转为大地测量形变速率场时，问题不再仅仅是一般意义上的图像边缘模糊，而在于断层两侧速度不连续边界在后续滤波加权过程中是否发生跨边界权重传递。若仍采用一般性的平滑方式对速率场进行处理，则虽然能够在一定程度上压制局部异常点，但也容易在断层邻域引入跨断层权重泄漏，使滤波结果出现边界模糊或者局部失真。

6、进一步地，应变率空间分布的计算通常建立在形变速率场空间变化分析的基础上。若输入的形变速率场在滤波阶段已经发生跨断层信息混合，则该混合效应会继续传递至后续应变率计算过程中，进而影响断层邻域应变集中区域的位置表达和幅值表达，降低结果对构造变形分析和地震危险性评估的支撑能力。结合上述现有技术可以看出，现有方案分别侧重于形变信息融合、速率场拼接或者图像滤波，但在形变速率场滤波处理中，对于如何通过加权运算同时兼顾离群噪声抑制和断层边界保持，尚未形成与后续应变率分析直接对应的处理方式。

7、因此，现有技术存在的主要问题在于：在形变速率场滤波处理中，现有加权运算方式难以同时兼顾离群噪声抑制和断层两侧速度不连续边界保持，容易产生跨断层权重泄漏，导致滤波后的形变速率场出现边界模糊或者局部失真，并进一步影响后续应变率空间分布结果的准确性。

技术实现思路

1、为了克服现有技术中在形变速率场滤波处理中难以同时兼顾离群噪声抑制和断层两侧速度不连续边界保持、易产生跨断层权重泄漏，进而影响滤波后形变速率场及后续应变率空间分布结果准确性的缺陷，本发明提供如下技术方案：

2、对采集的升轨雷达影像序列与降轨雷达影像序列分别进行视线向时序反演，得到升轨视线速率场与降轨视线速率场，并对升轨视线速率场与降轨视线速率场进行缺测识别，生成缺测掩膜图；

3、基于采集的gnss站点速度场观测数据，在缺测掩膜图中对升轨视线速率场和降轨视线速率场进行联合反演，得到东向速率场和北向速率场；

4、基于缺测掩膜图生成缺测边界权重图；基于东向速率场和北向速率场生成跨断层屏蔽权重；

5、基于缺测边界权重图和跨断层屏蔽权重，对东向速率场进行嵌套式双边滤波处理，获得东向滤波速率场；

6、在缺测掩膜图中，对东向滤波速率场和北向速率场分别进行局部趋势拟合求导，获得应变率空间分布结果。

7、进一步地，基于缺测掩膜图生成缺测边界权重图的方法包括：

8、缺测掩膜图包括有效像元和缺测边界线，有效像元为缺测掩膜图中未被标记为缺测的像元，缺测边界线为有效像元与被标记为缺测的像元在网格上的相邻分界线；

9、计算各有效像元到缺测边界线的最短距离，并对各有效像元对应的最短距离进行距离归一化处理，得到缺测边界权重；将各有效像元对应的缺测边界权重按缺测掩膜图中的像元位置写入，生成缺测边界权重图。

10、进一步地，基于东向速率场和北向速率场生成跨断层屏蔽权重的方法包括：

11、基于东向速率场和北向速率场生成速度变化强度图；

12、基于速度变化强度图生成断层边界图；

13、根据断层边界图，对未被标记为断层边界的有效像元进行连通区域划分，得到连通区域；对未被标记为断层边界且已经划入连通区域的有效像元赋予所在连通区域对应的分区标签，得到断层分区图；

14、分别读取中心像元和邻域像元在断层分区图中的分区标签，并判断中心像元和邻域像元对应的分区标签是否一致；

15、当中心像元和邻域像元对应的分区标签一致时，将中心像元与邻域像元之间的跨断层屏蔽权重记为1；

16、当中心像元和邻域像元对应的分区标签不一致时，将中心像元与邻域像元之间的跨断层屏蔽权重记为0。

17、进一步地，基于东向速率场和北向速率场生成速度变化强度图的方法包括：

18、以各有效像元作为中心像元，读取中心像元的东向速率值和北向速率值；

19、在预设邻域内读取各邻域有效像元的东向速率值和北向速率值，并分别计算各邻域有效像元与中心像元之间的东向速率差值和北向速率差值；

20、对各东向速率差值取绝对值并选取最大值作为东向变化量，对各北向速率差值取绝对值并选取最大值作为北向变化量；

21、对东向变化量和北向变化量进行平方求和并开方，得到速度变化强度，并按各中心像元对应的速度变化强度生成速度变化强度图。

22、进一步地，基于速度变化强度图生成断层边界图的方法包括：

23、逐一读取速度变化强度图中各有效像元对应的速度变化强度，并将速度变化强度大于或者等于预设断层边界候选阈值的有效像元标记为候选像元；

24、判断各候选像元之间是否在空间位置上相互邻接，并将相互邻接的候选像元归并为候选连通域；

25、对各候选连通域计算连通长度，删除连通长度小于预设连通长度阈值的候选连通域；

26、对保留的候选连通域执行迭代细化处理，将删除后仍保持候选连通域连续且与两个以上相邻候选像元邻接的候选像元，确定为本轮可删除候选像元；逐轮删除本轮可删除候选像元，直至不存在本轮可删除候选像元，得到经细化后的候选连通域；

27、将经细化后的候选连通域中的剩余候选像元标记为断层边界，生成断层边界图。

28、进一步地，获得东向滤波速率场的方法包括：

29、以未被标记为断层边界且已经写入分区标签的有效像元作为中心像元，针对每个中心像元预设参考估计窗口和主滤波窗口；

30、在参考估计窗口内对东向速率场进行局部趋势拟合，确定中心位置的参考速度；

31、在主滤波窗口内，根据邻域有效像元与中心像元之间的位置差确定各邻域有效像元对应的空间域权重，根据邻域有效像元的东向速率值与参考速度之间的差值，计算各邻域有效像元对应的速度域权重；

32、根据空间域权重、速度域权重、跨断层屏蔽权重和缺测边界权重，计算各邻域有效像元对应的组合权重；

33、根据各邻域有效像元对应的组合权重，对主滤波窗口内的东向速率值执行归一化加权平均，得到中心像元对应的东向滤波速率值；

34、将各中心像元对应的东向滤波速率值按网格位置写入，得到东向滤波速率场。

35、进一步地，确定中心位置的参考速度的方法包括：

36、在参考估计窗口内读取有效像元的平面位置坐标和东向速率值；

37、以平面位置坐标为输入、以东向速率值为输出，进行局部趋势拟合，得到东向速率关于平面位置坐标的局部平面表达式；

38、将中心像元的平面位置坐标代入局部平面表达式，得到中心位置对应的拟合东向速率值，并将拟合东向速率值确定为中心位置的参考速度。

39、进一步地，计算各邻域有效像元对应的速度域权重的方法包括：

40、将参考估计窗口内各有效像元的平面位置坐标分别代入局部平面表达式，得到各有效像元对应的拟合东向速率值；

41、计算各有效像元的东向速率值与拟合东向速率值之间差值的绝对值，得到拟合残差；

42、对各拟合残差取绝对值，得到拟合残差绝对值序列；读取东向速率场的数据记录精度作为速率记录分辨量，从拟合残差绝对值序列中提取中位数，并将中位数与速率记录分辨量相加，得到噪声尺度；并将噪声尺度与预设比例系数相乘，得到速度相似性阈值；

43、计算主滤波窗口内邻域有效像元的东向速率值与参考速度之间差值的绝对值，得到速度差异值，并将速度差异值和速度相似性阈值代入预设的速度域权重核函数，得到各邻域有效像元对应的速度域权重。

44、进一步地，计算各邻域有效像元对应的组合权重的方法包括：

45、在主滤波窗口内，读取各邻域有效像元的跨断层屏蔽权重和缺测边界权重；

46、计算邻域有效像元与中心像元之间的位置差，并将位置差代入预设空间域权重核函数，得到各邻域有效像元对应的空间域权重；

47、将各邻域有效像元对应的空间域权重、速度域权重、跨断层屏蔽权重和缺测边界权重相乘，得到各邻域有效像元对应的组合权重。

48、进一步地，得到中心像元对应的东向滤波速率值的方法包括：

49、统计主滤波窗口内组合权重大于0的邻域有效像元数量占主滤波窗口内邻域有效像元数量的比例，得到滤波有效像元比例；

50、当滤波有效像元比例大于或等于预设的滤波有效像元比例阈值时，基于各邻域有效像元对应的组合权重对主滤波窗口内的东向速率值执行归一化加权平均，得到中心像元对应的东向滤波速率值；

51、当滤波有效像元比例小于滤波有效像元比例阈值时，使中心像元对应的东向滤波速率值取参考速度。

52、与现有技术相比，本发明的一种用于大地测量速度场的改进双边滤波应变率计算方法的技术效果和优点为：

53、本发明通过对升轨视线速率场和降轨视线速率场进行缺测识别，并在缺测掩膜图限定的有效像元范围内结合gnss站点速度场观测数据进行联合反演，得到东向速率场和北向速率场，使东向速率场和北向速率场按照缺测掩膜图中的有效像元范围输出，并为后续缺测边界权重图生成、跨断层屏蔽权重生成和局部趋势拟合求导提供统一网格位置基础。

54、本发明基于缺测掩膜图生成缺测边界权重图，并基于东向速率场和北向速率场生成跨断层屏蔽权重，使缺测边界影响和断层两侧速度不连续边界影响不再以隐含方式混入后续滤波运算，而是分别转化为可直接参与加权计算的边界约束信息，能够减轻缺测边界附近样本分布不完整对邻域加权结果的拉动作用，并降低断层一侧速度信息通过邻域权重向另一侧传递的风险，从而为后续形变速率场滤波提供更贴合边界条件的约束基础。

55、本发明对东向速率场执行嵌套式双边滤波处理，在参考估计窗口内确定中心位置的参考速度，在主滤波窗口内联合空间域权重、速度域权重、跨断层屏蔽权重和缺测边界权重计算组合权重，并基于组合权重对邻域有效像元的东向速率值执行归一化加权平均；其中，分区标签不一致的邻域有效像元对应的跨断层屏蔽权重记为0，分区标签一致的邻域有效像元继续结合空间域权重、速度域权重和缺测边界权重参与组合权重计算。

56、本发明在主滤波窗口内进一步结合有效像元占比对归一化加权平均结果进行约束，当有效像元支撑不足时使中心像元对应的东向滤波速率值取参考速度，能够避免在样本不足或者边界邻接关系复杂条件下继续使用不稳定的加权结果，从而减轻缺测区域附近由于参与样本不完整导致的局部波动和异常偏移，提高滤波结果的稳定性。

57、本发明在缺测掩膜图中对东向滤波速率场和北向速率场分别进行局部趋势拟合求导，获得应变率空间分布结果，使应变率计算不再直接建立在未经边界约束和噪声抑制的速率场基础上，而是建立在已经完成边界控制和滤波处理的速率场基础上，能够减轻离群噪声、缺测边界扰动以及跨断层信息混合在求导阶段的继续放大，减少断层邻域应变集中区域在位置表达和幅值表达上的失真，从而提高应变率空间分布结果的稳定性和可用性。

58、综上，本发明通过围绕有效像元限定、边界约束构建、嵌套式双边滤波加权和局部趋势拟合求导的连续处理链，针对形变速率场滤波处理中现有加权运算方式难以同时兼顾离群噪声抑制和断层边界保持、易产生跨断层权重泄漏的问题，提高了滤波后形变速率场的边界表达能力、稳定性以及后续应变率空间分布结果的准确性和可靠性。