CN118191841A - 一种基于相关性分析的地表沉降形变测量校正方法、装置、设备及介质 - Google Patents
一种基于相关性分析的地表沉降形变测量校正方法、装置、设备及介质 Download PDFInfo
- Publication number
- CN118191841A CN118191841A CN202410338689.2A CN202410338689A CN118191841A CN 118191841 A CN118191841 A CN 118191841A CN 202410338689 A CN202410338689 A CN 202410338689A CN 118191841 A CN118191841 A CN 118191841A
- Authority
- CN
- China
- Prior art keywords
- phase
- deformation
- aperture
- view
- sub
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 87
- 238000010219 correlation analysis Methods 0.000 title claims abstract description 50
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 53
- 239000005436 troposphere Substances 0.000 claims abstract description 38
- 238000012937 correction Methods 0.000 claims abstract description 34
- 238000005259 measurement Methods 0.000 claims abstract description 32
- 238000012876 topography Methods 0.000 claims abstract description 31
- 238000012545 processing Methods 0.000 claims abstract description 18
- 238000011160 research Methods 0.000 claims abstract description 8
- 230000008569 process Effects 0.000 claims description 13
- 238000004062 sedimentation Methods 0.000 claims description 12
- 239000011159 matrix material Substances 0.000 claims description 11
- 238000007781 pre-processing Methods 0.000 claims description 9
- 230000011218 segmentation Effects 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 8
- 238000001914 filtration Methods 0.000 claims description 7
- 238000012952 Resampling Methods 0.000 claims description 6
- 230000008859 change Effects 0.000 claims description 6
- 230000006870 function Effects 0.000 claims description 6
- 238000004458 analytical method Methods 0.000 claims description 5
- 238000004590 computer program Methods 0.000 claims description 5
- 238000000605 extraction Methods 0.000 claims description 4
- 238000004091 panning Methods 0.000 claims description 3
- 238000000926 separation method Methods 0.000 claims description 3
- 101000581402 Homo sapiens Melanin-concentrating hormone receptor 1 Proteins 0.000 claims description 2
- 102000037055 SLC1 Human genes 0.000 claims description 2
- 238000001228 spectrum Methods 0.000 claims description 2
- 230000000694 effects Effects 0.000 description 9
- 238000011439 discrete element method Methods 0.000 description 7
- 238000012544 monitoring process Methods 0.000 description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000005305 interferometry Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 241001420622 Meris Species 0.000 description 1
- 208000004350 Strabismus Diseases 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000005388 cross polarization Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001934 delay Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000010587 phase diagram Methods 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 238000012502 risk assessment Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Landscapes
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种基于相关性分析的地表沉降形变测量校正方法、装置、设备及介质，方法：获取研究区域的按时间序列分布的多景单基线全孔径SAR影像，构建多组全孔径影像对及对应的前视、后视子孔径影像对，并干涉处理和去除平地相位、“裸地表”对应的地形相位、轨道误差相位，获得差分干涉相位；利用小波分解对前视、后视的差分干涉相位进行相关性分析，去除地形残差和噪声相位；再对时序相邻且无时间交叉的两个差分干涉图进行相关性分析，提取相邻时序干涉图的共同形变相位，进而分离对流层延迟相位和形变相位；最终进行堆栈法反演，完成地表形变测量的校正。本发明在未引入外部数据的情况下，提供一种稳健可行的地表沉降形变测量校正方法。
Description
技术领域
本发明属于InSAR影像数据处理领域，具体涉及一种基于相关性分析的地表沉降形变测量校正方法、装置、设备及介质。
背景技术
地震、滑坡等地质构造活动或地下水开采、工程建设等人类活动都会引起地表变形现象，地表变形又会加剧这些现象的发生，因此，对地表变形进行监测具有重要的研究意义。地表形变测量是指通过现代化的技术手段对地表进行监测和测量，以获取地表形变信息，可用于地质灾害风险评估、资源开发与环境保护、基础设施安全、气候变化研究以及地表科学研究。相较传统测绘、GPS、光学遥感等测量技术，合成孔径雷达干涉(Interferometry Synthetic Aperture Radar,InSAR)技术全天候、全天时、高精度、低成本、大范围观测的优势，使其在地表形变测量领域得到了广泛应用。但受时空失相干、轨道误差、对流层延迟等因素影响，单基线InSAR难以获取高质量数据，从而影响InSAR形变监测的精度。其中，对流层延迟是重轨InSAR测量工作的主要误差源之一，因其分布随机且随时空变化差异较大，且与形变信息混杂在一起，一直是地表形变测量数据处理难点之一。
现有大气误差校正方法主要有堆叠法、高程相关分析法、时序InSAR校正法和借助外部数据的方法。现有方法面临的问题主要有以下几个方面：
(1)基于外部数据的InSAR校正方法对外部数据有依赖性：
目前，借助外部数据校正InSAR对流层延迟的方法是常用手段，尤其是GACOS的出现，为这一方法提供了更大便利。然而，不可否认的是，此类方法对外部数据可靠性的依赖从未减小。此类方法的局限主要有以下几点：①外部数据的获取：GPS校正方法要求实验区有密集GPS站布网，GPS数据获取过程中的云雨、光照环境；仪器设备是否正常运作以及传播路径延迟等都会影响数据质量；②外部数据的使用场景限制：MERIS、MODIS等大气水气数据无法在多云层、夜间天气使用；③最重要的是需要的数据量较大，且外部数据的分辨率会影响校正效果；
(2)基于时空滤波的InSAR大气校正方法具有局限性：
首先，基于时空滤波的高程相关分析法受限于垂直分层原理，对于平坦地区的大气效应的去除效果不佳。其次，基于时空滤波的时序InSAR大气校正方法在时序InSAR技术的使用上会影响对流层延迟误差的去除效果。从MT-InSAR的角度讲，这类技术首先要求大量的SAR影像数据集，提高了实验数据的要求和成本，也增加了研究工作量；此外，不同的时序InSAR技术有不同的局限，例如PS-InSAR要求实验区稳定可靠。另外，在InSAR形变监测应用中，形变信号和大气信息可能存在重叠，无法通过时空滤波方法提取大气相位信息；
(3)缺少利用单基线InSAR数据去除InSAR地形测量大气误差的简便方法：
由(1)(2)两点可知，现有的方法中结合常规InSAR技术的对流层延迟方法存在所需数据量大的不足。虽然也出现了基于“极化不变性”的利用双极化InSAR数据提取大气相位分量的方法，但在森林等复杂区域并不存在满足相干性要求的交叉极化数据。
然而，这些校正方法需要大量数据集支持，工作量大，且易受云雨雾天气影响，难以支持复杂试验区的InSAR变形监测大气校正工作。基于此，研究需要解决数据量需求大、依赖外部水汽数据资料的问题，针对长重访周期的InSAR平台(例如ALOS、TerraSAR等卫星)，提出一种无需依赖外部水汽数据资料，仅通过少量InSAR数据实现InSAR沉降形变测量大气误差校正的简便方法。
发明内容
本发明提供一种基于相关性分析的地表沉降形变测量校正方法、装置、设备及介质，在不借助外部水汽数据的前提下，实现复杂区域形变测量的对流层延迟校正工作，从而获得符合精度要求的形变数据。
为实现上述技术目的，本发明采用如下技术方案：
一种基于相关性分析的地表沉降形变测量校正方法，包括：
步骤1，获取研究区域的按时间序列分布的多景单基线全孔径SAR影像，并进行统一配准、重采样，将每两个时间相邻的影像组对得到多组全孔径影像对，再对各全孔径影像对均沿方位向进行子孔径分解处理，得到对应的前视、后视子孔径影像对；
步骤2，对每一组全孔径影像对、前视子孔径影像对、后视子孔径影像对分别进行干涉处理，获得干涉图；
步骤3，对步骤2获得的各干涉图分别进行预处理，去除其中的平地相位、“裸地表”对应的地形相位、轨道误差相位，获得对应的差分干涉相位；
步骤4，根据地形残差、对流层延迟、形变相位分量的频率特征，利用小波分解对前视、后视子孔径影像对的差分干涉相位进行相关性分析，去除其中的地形残差和噪声相位，得到相似信息，包括对流层延迟相位和形变相位；
步骤5，根据时序上对流层延迟相位、形变相位的分布特性，对时序相邻且无时间交叉的两个干涉图所得到的相似信息进行批量相关性分析，提取相邻时序干涉图的共同形变相位，进而分离每个干涉图的对流层延迟相位和形变相位；
步骤6，对时序干涉图所得到的时序形变相位进行堆栈法反演，完成地表形变测量的校正。
进一步地，将n景单基线全孔径SAR影像按时序记为SLCm，m＝1,2,3...n，以第1景的单基线全孔径SAR影像SLC1为参考，对其余时间上的单基线全孔径SAR影像SLC2～SLCn进行统一配准。
进一步地，步骤2对配准后的各影像SLCm沿方位向进行子孔径分割处理的方法为：首先将配准后的影像SLCm重采样至主影像坐标框架，然后分别沿方位向进行一维傅里叶变换；再对傅里叶变化后的方位向频谱进行带通滤波，截取方位向上正值域进行傅里叶逆变换得到前视子孔径影像SLCm.F，截取方位向上负值域进行傅里叶逆变换得到后视子孔径影像SLCm.B。
进一步地，设全孔径影像SLCm和SLC(m+1)得到的干涉图为Inf_m.m+1，其干涉相位表示为：
式中，为由影像对得到的干涉相位，/>分别为平地相位、地形相位、轨道误差相位、对流层延迟相位、地形形变相位、噪声相位；
全孔径干涉图为Inf_m.m+1的前视、后视子孔径影像的干涉相位表示为：
式中，F、B分别代表前视、后视两个子孔径，为前视子孔径干涉图中的地形相位和噪声相位，/>为后视子孔径干涉图中的地形相位和噪声相位；
步骤3中预处理时，基于以下平地相位和地形相位的计算式，计算并去除平地相位和“裸地表”对应的地形相位为：
式中，λ为雷达微波的波长，R1、R2为别是不同时间下雷达与散射目标的距离；为由影像SLCm和SLC(m+1)得到的地形相位，B为总基线长，B⊥为垂直基线，θ为入射角，α为基线倾角；h为“裸地表”高度，Δhm.m+1为散射目标高度；
采用分块的方式，对轨道误差相位进行多项式拟合，并去除该轨道误差相位，得到全孔径、前视、后视差分干涉相位分别表示为：
式中，m为单基线全孔径SAR影像的时序编号，为全孔径、前视子孔径和后视子也孔径干涉图中去除“裸地表”对应地形相位后剩余的地形相位，即地形残差相位。
进一步地，步骤4具体过程为：
首先，将利用小波分解对尺寸大小为p×q的前视、后视差分干涉图进行的多分辨率分解，表示为：
式中，Φ和Ψ分别是平移函数和母小波函数；J是小波分解尺度的数量，j表示第j层尺度；和/>分别代表多分辨率分析后前视、后视差分干涉图的低频和高频系数；参数ε＝1,2,3分别表示水平、垂直和对角方向；p，q为前视、后视差分干涉图的尺寸，x,y是第j层尺度下对应的分解矩阵(高频或低频)尺寸，ix和iy分别表示分解矩阵的第ix行和第iy列；；
然后计算每层尺度j下前视、后视差分干涉图在ε方向的高频系数之间的相关系数
式中，df、dmax分别第f个高频系数点到拟合直线的垂直距离和所有高频系数点到拟合直线的最大垂直距离；λi表示对应的行权重值；为分解尺度j下的前视、后视信息对应的高频系数矩阵；
最终通过相关性分析得到j尺度下ε方向高频系数中的对流层延迟、形变信号：
式中，右下角标中的i表示当前分解尺度j下对应的矩阵，表示j分解尺度下ε方向提取到的对应对流层延迟、形变信号；
再将前视干涉图在J尺度下低频信息和式(11)或(12)提取到的每一层分解尺度下的高频信息进行小波重构，得到相似信息。
进一步地，小波分解尺度的数量J，由均方根误差法确定；其中，均方根误差计算式为：
式中，为前视、后视干涉图的原始信号，/>为前视、后视干涉图J级分解重构信号；
计算相邻两级的信号均方根误差比rJ+1为：
将rJ+1达到预设值时的分解尺度作为最优的小波分解尺度的数量。
进一步地，步骤5具体过程为：
将相邻全孔径影像SLCm和SLC(m+1)的干涉图Inf_m.m+1经步骤4后得到的相似信息记为其由流层延迟相位和形变相位构成，表示为：
则干涉图Inf_m.m+1时序相邻且无时间交叉的另一个干涉图Inf_m+2.m+3，其相位构成表示为：
基于干涉图Inf_m.m+1与干涉图Inf_m+2.m+3两者无时间交叉得到与/>相互独立，而形变相位/>与/>包含共同的形变相位，即/> 为Inf_m.m+1与Inf_m+2.m+3两组干涉图在时间间隔内发生的沉降形变；由此将式(18)变形为：
此时，通过相关性分析提取式(17)(19)中的形变相位完成对形变相位和对流层延迟相位/>的分离；
按同样的方法分离所有干涉图的对流层延迟相位和形变相位。
一种基于相关性分析的地表沉降形变测量校正装置，包括：
子孔径分割模块，用于：对每两个时间相邻的影像组对得到多组全孔径影像对，均沿方位向进行子孔径分解处理，得到对应的前视、后视子孔径影像对；
数据预处理模块，用于：对每一组全孔径影像对、前视子孔径影像对、后视子孔径影像对分别进行干涉处理，并去除其中的平地相位、“裸地表”对应的地形相位、轨道误差相位，获得对应的差分干涉相位；
子孔径相似相位提取模块，用于：根据地形残差、对流层延迟、形变相位分量的频率特征，利用小波分解对前视、后视子孔径影像对的差分干涉相位进行相关性分析，去除其中的地形残差和噪声相位，得到相似信息，包括对流层延迟相位和形变相位；
时序相似相位提取模块，用于：根据时序上对流层延迟相位、形变相位的分布特性，对时序相邻且无时间交叉的两个干涉图所得到的相似信息进行相关性分析，提取相邻时序干涉图的共同形变相位，进而分离干涉图的对流层延迟相位和形变相位；
形变信息反演模块，用于：对时序干涉图所得到的时序形变相位进行堆栈法反演，完成地表形变测量的校正。
一种电子设备，包括存储器及处理器，所述存储器中存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器实现上述任一项所述的基于相关性分析的地表沉降形变测量校正方法。
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任一项所述的基于相关性分析的地表沉降形变测量校正方法。
有益效果
本发明通过对感兴趣实验区时序SAR影像进行子孔径分割处理，获得前视、后视雷达影像；对雷达影像进行预处理：通过引入两种不同的外部DEM，分别对时序全孔径SAR影像对和其对应的前视、后视影像对构建差分干涉相位公式，实现平地相位、“裸地表”对应的地形相位的去除；再通过多项式拟合模拟轨道误差；使用小波分解对前视、后视干涉图进行相关性分析，对地形残差相位进行二次预处理，成功提取时序干涉对的对流层延迟、形变等相似信息；利用时序干涉对对流层延迟互相独立，沉降形变含有相似信息的特征，再次使用小波分解对相邻时序干涉对进行相关性分析，实现时序干涉对的对流层延迟误差校正；最后，使用InSAR堆栈法反演时序InSAR的沉降形变信息。
本发明通过子孔径分割和相关性分析实现沉降形变测量误差校正，使形变测量精度更高；而且本发明仅使用自身SAR数据，未引入气象资料等外部辅助数据，避免外部误差引入的同时，具有良好的稳健性；除此之外，本发明可以同时校正垂直分层相关的大气延迟和湍流大气延迟，适用性更广。为长周期地表沉降形变测量校正提供了一种新的思路与方法，可用于复杂区域的地表沉降形变测量作业等方面。
附图说明
图1为本发明实施例所述方法的流程图。
图2为本发明实施例所述方子孔径InSAR相关性分析流程图。
图3为本发明所选择美国的南加州作为验证区域的位置及SAR影像的覆盖范围。
图4为实验区地形反演结果示意图。其中，数字序号(1-6)分别表示干涉图DInf_1.2、DInf_3.4、DInf_5.6、DInf_7.8和DInf_2.3、DInf_4.5、DInf_6.7对应的实验结果图；图(1.a)和图(1.a.2)分别表示DInf_1.2对应的前视、后视干涉图，图(1.b)为图(1.a)和图(1.a.2)进行相关性分析得到的相似信息，即DInf_1.2的对流层延迟相位和形变相位；图(2.a)和图(2.a.2)分别表示DInf_3.4对应的前视、后视干涉图，图(2.b)为图(1.b)和图(1.b.2)进行相关性分析得到的相似信息，即DInf_3.4的对流层延迟相位和形变相位；图(1.b.2)和图(2.b.2)为模拟沉降形变；(1-6，c)则依次为添加模拟形变后得到的差分干涉图Dinf_1.2、Dinf_3.4、Dinf_5.6、Dinf_7.8和Dinf_2.3、Dinf_4.5、Dinf_6.7；(1.d)为(1.c)和(2.c)通过相关性分析得到的相似形变信息，即干涉图DInf_1.2的形变相位；同样地，(2,d)为Dinf_3.4的形变相位，(3,d)和(4,d)为Dinf_5.6的形变相位，(5,d)为Dinf_2.3的形变相位，(6,d)和(7,d)为Dinf_4.5的形变相位。
具体实施方式
为了更好的说明本发明的方法和步骤，以位于美国西海岸的南加州实验区的ALOSPALSAR数据，对本发明进行进一步详细说明，利用瑞士GAMMA公司发布的GAMMA软件辅助数据预处理；注意，此处所描述的具体实施仅用以解释本发明，并不用于限定本发明。
本实施例提供一种基于相关性分析的地表沉降形变测量校正方法，如图1、2所示，包括以下步骤：
步骤1，根据所选试验区域的坐标范围，获取L波段ALOS卫星PALSAR数据，和对应外部DEM：30m分辨率的SRTM DEM和AW3D DEM。对获取的SAR影像进行配准、重采样和子孔径分割处理。
本实施实例选择了如图3所示位于美国西海岸的南加州实验区2008年08月20日到2010年10.月11日的ALOS PALSAR数据，其具体信息如表1所示。对获取的SLC影像进行配准、重采样，借助Gamma处理软件完成子孔径分割，获得前视、后视影像对。
表1
步骤2，将表1按照表2的方式组成时序干涉图，通过GAMMA软件对每一组原始全孔径影像对、前视子孔径影像对、后视子孔径影像对进行干涉处理，获得干涉图。分析干涉图的相位组成及其是否受子孔径分割的影响，构建全孔径、前视、后视干涉相位公式；
全孔径干涉图Inf_m.m+1(m＝1,2,3,4,5,6,7)的干涉相位表示为：
式中，为SLCm和SLCm+1组成的干涉相位，其中m和m+1分别表示相邻的全孔径SLC影像；/>为地形相位，/>是轨道误差相位，/>为对流层延迟相位，/>为形变相位，/>代表噪声相位；
Inf_m.m+1对应的前视、后视子孔径影像的干涉相位可表示为：
F、B分别代表前视、后视两个子孔径，“裸地表”高度对应的部分地形相位表现相同，但地面散射目标高度相关的部分地形相位会受前视、后视影响；轨道误差由轨道参数不准确引起的，是一种具有弱二次曲面特征的系统误差，在前视、后视干涉图中表现相同；对于子孔径分解后的前视、后视干涉图而言，虽然其频率、观测视角发生了变化，但相对卫星到散射目标的距离而言，极短时间内所穿越的大气变化不显著，因此可认为对流层相位误差独立于斜视角的变化；/>同样也不受子孔径分解的影响；/>属于随机误差；
表2
步骤3，对步骤2得到的每一组原始全孔径干涉图、前视干涉图、后视干涉图进行实验预处理，实现平地相位和“裸地表”对应的地形相位及轨道误差相位的去除；
将表1按照表2的方式组成时序干涉图，以GAMMA软件作为数据处理平台，分别引入SRTM、AW3D作为外部DEM，对这些时序影像对及其对应的前视、后视影像对进行差分干涉处理，对得到的差分干涉图滤波、相位解缠，生成前视、后视解缠差分干涉图，差分干涉处理的具体参数设置如表3所示。
平地相位GAMMA软件基于SLC基线参数，可以将其去除；
地形相位为了使前视、后视干涉相位残余地形误差的差异更明显，利用后续相关性分析提取相似的对流层延迟相位和形变相位，对前视、后视干涉图分别引入两种不同外部DEM去除地形相位分量，通过GAMMA软件模拟地形相位，去除其与“裸地表”对应的部分地形相位，剩余全孔径、前视、后视干涉图的地形残差相位/>
轨道误差具有差弱二次曲面的特征，以MATLAB为平台，编写基于多项式拟合代码，实现轨道误差的去除；
完成平地相位、“裸地表”对应的地形相位及轨道误差的校正，全孔径、前视、后视解缠差分干涉相位表示为：
表3差分干涉数据处理步骤
步骤4，根据地形残差、对流层延迟、形变等干涉相位分量的频率特征，通过小波分解相关性分析法，实现地形残差的第二次预处理，提取前视、后视差分干涉图相同的对流层延迟相位和形变相位；
地形残余误差相位和噪声相位主要由高频信号组成，大气误差相位和形变相位则混杂在高频、低频信息中，简单的滤波方法很难区分这些信号。由于前视、后视干涉图中的 相同，其对应的小波系数也呈现相似性。为了评估这种相似性，可以计算每层尺度j下高频系数之间的相关系数。已有研究通过小波分解对InSAR距离变化图和DEM进行相关性分析，提取与地形相关的对流层延迟相位。多分辨率分析(Multi-resolutionAnalysis，MRA)是小波变换的重要应用之一，通过这种变换，可以把前视、后视干涉图由粗到细逐渐分解成不同尺度下具有不同频率特征的多个分量，计算其小波系数的相关性，评估干涉相位在不同尺度下的空间分布，将相关系数接近1的小波系数归结为大气和形变信号，成功提取对流层延迟相位和形变相位。
去除地形残差、噪声等信息，实现地形残差的第二次预处理：
对尺寸大小为p×q的前视、后视差分干涉图的多分辨率分析可表示为：
式中，Φ和Ψ分别是平移函数和母小波函数；J是小波分解尺度的数量，j表示第j层尺度；和/>分别代表多分辨率分析后前视、后视差分干涉图的低频和高频系数；参数ε＝1,2,3分别表示水平、垂直和对角方向。p，q为前视、后视差分干涉图的尺寸，x,y是第j层尺度下对应的分解矩阵(高频或低频)尺寸，ix和iy分别表示分解矩阵的第ix行和第iy列。由于前视、后视干涉图中的/>相同，其对应的小波系数也呈现相似性。为了评估这种相似性，可以计算每层尺度j下高频系数之间的相关系数。对前视、后视两幅干涉图做小波分解，分解后得到两个x×y大小的高频系数矩阵，将两个高频矩阵对应行列号的值作为ix和iy的坐标，形成一个高频点，根据这些点可以拟合得到拟合直线。然后计算每个高频系数点到这条直线的距离，这个距离和最远垂直距离的比可以表示该点对应高频信息的相关性大小。为增强相关性分析对信号随机性的敏感度，将各个高频系数点到拟合直线的垂直距离di和最大垂直距离dmax作为权重参考，相关系数/>为：
式中，df、dmax分别第f个高频系数点到拟合直线的垂直距离和所有高频系数点到拟合直线的最大垂直距离；λi表示对应的行权重值。为分解尺度j下的前视、后视信息对应的高频系数矩阵。
通过相关性分析得j尺度下ε方向高频系数中的对流层延迟、形变信号
将J尺度下低频信息和式(12)提取到的每一层分解尺度下的高频信息进行小波重构，得到相似信息，即
小波分解尺度J由均方根误差法确定，具体过程为：
式中，为前视、后视干涉图的原始信号，/>为前视、后视干涉图的J级分解重构信号。
相邻两级的信号均方根误差比为：
将rJ+1接近1时对应的分解尺度作为分解高频、低频信号的最优分解尺度。此时分解得到J+1尺度和J尺度下的高频小波系数最相似，可认为地形残差、噪声等信号已基本去除；
对所有干涉图做上述相关性分析，得到其对应的对流层延迟相位和形变相位ATP_DEF_m.m+1(m＝1,2,3,4,5,6,7)，如图4(1,b)、(2,b)分别为干涉图DInf_1.2、DInf_3.4提取到的ATP_DEF_1.2、ATP_DEF_3.4。
步骤5，通过分析对流层延迟、形变相位的时空特征，对相邻的全孔径干涉图进行小波分解相关性分析，实现对流层延迟和形变信息的分离：
为使实验效果展示更明显，对步骤4获得的ATP_DEF_m.m+1(m＝1,2,3,4,5,6,7)分别加入逐渐变化的沉降形变信息，得到图4(1-6.c)。如图4(1,b,2)和(2,b,2)即为图(1,b)、(2,b)包含的模拟形变。
以按时间顺序排列且没有时间交叉的干涉图DInf_1.2、DInf_3.4为例， 去除后，其相位组成表示为：
式中，分别为DInf_1.2、DInf_3.4的对流层延迟，DInf_1.2、DInf_3.4两者无时间交叉，所以/>相互相互独立；/>为DInf_1.2、DInf_3.4共同包含的形变信息，/>为DInf_1.2、DInf_3.4两组干涉图时间间隔内发生的沉降形变；
因此，分别对DInf_1.2、DInf_3.4，DInf_3.4、DInf_5.6，DInf_5.6、DInf_7.8做小波分解相关性分析，得到图4(1,d)、(2,d)、(3,d)所示的形变相位Def_1.2、Def_3.4、Def5.6；对DInf_2.3、DInf_4.5，DInf_4.5、DInf_6.7做小波分解相关性分析，得到图4(5,d)、(6,d)所示的形变相位Def_2.3、Def_4.5。
步骤6，进一步地，步骤5完成后获得的时序差分干涉图进行线性堆栈法解算，获得沉降形变信息。
通过GAMMA软件使用堆栈法(Stacking)，对所有形变相位使用Stacking法，计算研究区域的年平均形变速率和累计形变量；分别对相邻两幅形变相位图使用Stacking法，统计观察沉降形变的时序变化，完成地表形变测量的校正。
对计算校正前后的地表沉降形变均方根误差，得到表4，说明该发明对地表沉降形变测量的误差校正效果良好。
表4均方根误差
以上实施例为本申请的优选实施例，本领域的普通技术人员还可以在此基础上进行各种变换或改进，在不脱离本申请总的构思的前提下，这些变换或改进都应当属于本申请要求保护的范围之内。
本发明由湖南省自然科学基金优秀青年基金2023JJ20061提供支持。
Claims (10)
1.一种基于相关性分析的地表沉降形变测量校正方法，其特征在于，包括：
步骤1，获取研究区域的按时间序列分布的多景单基线全孔径SAR影像，并进行统一配准、重采样，将每两个时间相邻的影像组对得到多组全孔径影像对，再对各全孔径影像对均沿方位向进行子孔径分解处理，得到对应的前视、后视子孔径影像对；
步骤2，对每一组全孔径影像对、前视子孔径影像对、后视子孔径影像对分别进行干涉处理，获得干涉图；
步骤3，对步骤2获得的各干涉图分别进行预处理，去除其中的平地相位、“裸地表”对应的地形相位、轨道误差相位，获得对应的差分干涉相位；
步骤4，根据地形残差、对流层延迟、形变相位分量的频率特征，利用小波分解对前视、后视子孔径影像对的差分干涉相位进行相关性分析，去除其中的地形残差和噪声相位，得到相似信息，包括对流层延迟相位和形变相位；
步骤5，根据时序上对流层延迟相位、形变相位的分布特性，对时序相邻且无时间交叉的两个干涉图所得到的相似信息进行批量相关性分析，提取相邻时序干涉图的共同形变相位，进而分离每个干涉图的对流层延迟相位和形变相位；
步骤6，对时序干涉图所得到的时序形变相位进行堆栈法反演，完成地表形变测量的校正。
2.根据权利要求1所述的基于相关性分析的地表沉降形变测量校正方法，其特征在于，将n景单基线全孔径SAR影像按时序记为SLCm，m＝1,2,3...n，以第1景的单基线全孔径SAR影像SLC1为参考，对其余时间上的单基线全孔径SAR影像SLC2～SLCn进行统一配准。
3.根据权利要求1所述的基于相关性分析的地表沉降形变测量校正方法，其特征在于，步骤2对配准后的各影像SLCm沿方位向进行子孔径分割处理的方法为：首先将配准后的影像SLCm重采样至主影像坐标框架，然后分别沿方位向进行一维傅里叶变换；再对傅里叶变化后的方位向频谱进行带通滤波，截取方位向上正值域进行傅里叶逆变换得到前视子孔径影像SLCm.F，截取方位向上负值域进行傅里叶逆变换得到后视子孔径影像SLCm.B。
4.根据权利要求1所述的基于相关性分析的地表沉降形变测量校正方法，其特征在于，设全孔径影像SLCm和SLC(m+1)得到的干涉图为Inf_m.m+1，其干涉相位表示为：
式中，为由影像对得到的干涉相位，/>分别为平地相位、地形相位、轨道误差相位、对流层延迟相位、地形形变相位、噪声相位；
全孔径干涉图为Inf_m.m+1的前视、后视子孔径影像的干涉相位表示为：
式中，F、B分别代表前视、后视两个子孔径，为前视子孔径干涉图中的地形相位和噪声相位，/>为后视子孔径干涉图中的地形相位和噪声相位；
步骤3中预处理时，基于以下平地相位和地形相位的计算式，计算并去除平地相位和“裸地表”对应的地形相位为：
式中，λ为雷达微波的波长，R1、R2为别是不同时间下雷达与散射目标的距离；为由影像SLCm和SLC(m+1)得到的地形相位，B为总基线长，B⊥为垂直基线，θ为入射角，α为基线倾角；h为“裸地表”高度，Δhm.m+1为散射目标高度；
采用分块的方式，对轨道误差相位进行多项式拟合，并去除该轨道误差相位，得到全孔径、前视、后视差分干涉相位分别表示为：
式中，m为单基线全孔径SAR影像的时序编号，为全孔径、前视子孔径和后视子也孔径干涉图中去除“裸地表”对应地形相位后剩余的地形相位，即地形残差相位。
5.根据权利要求1所述的基于相关性分析的地表沉降形变测量校正方法，其特征在于，步骤4具体过程为：
首先，将利用小波分解对尺寸大小为p×q的前视、后视差分干涉图进行的多分辨率分解，表示为：
式中，Φ和Ψ分别是平移函数和母小波函数；J是小波分解尺度的数量，j表示第j层尺度；和/>分别代表多分辨率分析后前视、后视差分干涉图的低频和高频系数；参数ε＝1,2,3分别表示水平、垂直和对角方向；p，q为前视、后视差分干涉图的尺寸，x,y是第j层尺度下对应的分解矩阵尺寸，ix和iy分别表示分解矩阵的第ix行和第iy列；
然后计算每层尺度j下前视、后视差分干涉图在ε方向的高频系数之间的相关系数
式中，df、dmax分别第f个高频系数点到拟合直线的垂直距离和所有高频系数点到拟合直线的最大垂直距离；λi表示对应的行权重值；为分解尺度j下的前视、后视信息对应的高频系数矩阵；
最终通过相关性分析得到j尺度下ε方向高频系数中的对流层延迟、形变信号：
式中，右下角标中的i表示当前分解尺度j下对应的矩阵，表示j分解尺度下ε方向提取到的对应对流层延迟、形变信号；
再将前视干涉图在J尺度下低频信息和式(11)或(12)提取到的每一层分解尺度下的高频信息进行小波重构，得到相似信息。
6.根据权利要求5所述的基于相关性分析的地表沉降形变测量校正方法，其特征在于，小波分解尺度的数量J，由均方根误差法确定；其中，均方根误差计算式为：
式中，为前视、后视干涉图的原始信号，/>为前视、后视干涉图J级分解重构信号；
计算相邻两级的信号均方根误差比rJ+1为：
将rJ+1达到预设值时的分解尺度作为最优的小波分解尺度的数量。
7.根据权利要求1所述的基于相关性分析的地表沉降形变测量校正方法，其特征在于，步骤5具体过程为：
将相邻全孔径影像SLCm和SLC(m+1)的干涉图Inf_m.m+1经步骤4后得到的相似信息记为其由流层延迟相位和形变相位构成，表示为：
则干涉图Inf_m.m+1时序相邻且无时间交叉的另一个干涉图Inf_m+2.m+3，其相位构成表示为：
基于干涉图Inf_m.m+1与干涉图Inf_m+2.m+3两者无时间交叉得到与/>相互独立，而形变相位/>与/>包含共同的形变相位，即/> 为Inf_m.m+1与Inf_m+2.m+3两组干涉图在时间间隔内发生的沉降形变；由此将式(18)变形为：
此时，通过相关性分析提取式(17)(19)中的形变相位完成对形变相位/>和对流层延迟相位/>的分离；
按同样的方法分离所有干涉图的对流层延迟相位和形变相位。
8.一种基于相关性分析的地表沉降形变测量校正装置，其特征在于，包括：
子孔径分割模块，用于：对每两个时间相邻的影像组对得到多组全孔径影像对，均沿方位向进行子孔径分解处理，得到对应的前视、后视子孔径影像对；
数据预处理模块，用于：对每一组全孔径影像对、前视子孔径影像对、后视子孔径影像对分别进行干涉处理，并去除其中的平地相位、“裸地表”对应的地形相位、轨道误差相位，获得对应的差分干涉相位；
子孔径相似相位提取模块，用于：根据地形残差、对流层延迟、形变相位分量的频率特征，利用小波分解对前视、后视子孔径影像对的差分干涉相位进行相关性分析，去除其中的地形残差和噪声相位，得到相似信息，包括对流层延迟相位和形变相位；
时序相似相位提取模块，用于：根据时序上对流层延迟相位、形变相位的分布特性，对时序相邻且无时间交叉的两个干涉图所得到的相似信息进行相关性分析，提取相邻时序干涉图的共同形变相位，进而分离干涉图的对流层延迟相位和形变相位；
形变信息反演模块，用于：对时序干涉图所得到的时序形变相位进行堆栈法反演，完成地表形变测量的校正。
9.一种电子设备，包括存储器及处理器，所述存储器中存储有计算机程序，其特征在于，所述计算机程序被所述处理器执行时，使得所述处理器实现如权利要求1～7中任一项所述的方法。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1～7中任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410338689.2A CN118191841A (zh) | 2024-03-25 | 2024-03-25 | 一种基于相关性分析的地表沉降形变测量校正方法、装置、设备及介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410338689.2A CN118191841A (zh) | 2024-03-25 | 2024-03-25 | 一种基于相关性分析的地表沉降形变测量校正方法、装置、设备及介质 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN118191841A true CN118191841A (zh) | 2024-06-14 |
Family
ID=91408365
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202410338689.2A Pending CN118191841A (zh) | 2024-03-25 | 2024-03-25 | 一种基于相关性分析的地表沉降形变测量校正方法、装置、设备及介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN118191841A (zh) |
-
2024
- 2024-03-25 CN CN202410338689.2A patent/CN118191841A/zh active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Braun | Retrieval of digital elevation models from Sentinel-1 radar data–open applications, techniques, and limitations | |
| US11269072B2 (en) | Land deformation measurement | |
| CN111474544B (zh) | 一种基于sar数据的滑坡形变监测及预警方法 | |
| CN106772342B (zh) | 一种适用于大梯度地表沉降监测的时序差分雷达干涉方法 | |
| Catalão et al. | Merging GPS and atmospherically corrected InSAR data to map 3-D terrain displacement velocity | |
| CN113624122A (zh) | 融合GNSS数据与InSAR技术的桥梁变形监测方法 | |
| CN106526590A (zh) | 一种融合多源sar影像工矿区三维地表形变监测及解算方法 | |
| CN113866764B (zh) | 基于InSAR和LR-IOE模型的滑坡易发性改进评估方法 | |
| CN113340191B (zh) | 时间序列干涉sar的形变量测量方法及sar系统 | |
| CN109100719B (zh) | 基于星载sar影像与光学影像的地形图联合测图方法 | |
| CN112284332B (zh) | 基于高分辨率insar的高层建筑沉降监测结果三维定位方法 | |
| CN104316920A (zh) | 一种雷达高度计小入射角干涉的海面高度高精度提取方法 | |
| CN115201825B (zh) | 一种InSAR震间形变监测中的大气延迟校正方法 | |
| CN112051572A (zh) | 一种融合多源sar数据三维地表形变监测方法 | |
| Aobpaet et al. | Land subsidence evaluation using InSAR time series analysis in Bangkok metropolitan area | |
| CN112685819A (zh) | 一种用于大坝及滑坡变形gb-sar监测的数据后处理方法及系统 | |
| CN109118520A (zh) | 一种采动区电网杆塔位移监测方法及系统 | |
| CN113238228B (zh) | 基于水准约束的三维地表形变获取方法、系统及装置 | |
| CN115079172A (zh) | 一种MTInSAR滑坡监测方法、设备及存储介质 | |
| Bryksin et al. | Using of SAR data and DInSAR-PSInSAR technique for monitoring Western Siberia and Arctic. | |
| CN114415178A (zh) | 一种适用于变形地质灾害识别的InSAR快速处理方法—GHR-InSAR | |
| Mao et al. | Estimation and compensation of ionospheric phase delay for multi-aperture InSAR: An azimuth split-spectrum interferometry approach | |
| CN118191841A (zh) | 一种基于相关性分析的地表沉降形变测量校正方法、装置、设备及介质 | |
| CN113281748B (zh) | 一种地表形变监测方法 | |
| Osmanoglu | Applications and development of new algorithms for displacement analysis using InSAR time series |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication |