CN116047519B - 一种基于合成孔径雷达干涉测量技术的选点方法 - Google Patents

一种基于合成孔径雷达干涉测量技术的选点方法 Download PDF

Info

Publication number
CN116047519B
CN116047519B CN202310322768.XA CN202310322768A CN116047519B CN 116047519 B CN116047519 B CN 116047519B CN 202310322768 A CN202310322768 A CN 202310322768A CN 116047519 B CN116047519 B CN 116047519B
Authority
CN
China
Prior art keywords
pixel
point
points
coherence
phase
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.)
Active
Application number
CN202310322768.XA
Other languages
English (en)
Other versions
CN116047519A (zh
Inventor
王志伟
李文慧
赵月
王翔
赵同龙
郭秋英
张国建
任晓芳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shandong Jianzhu University
Original Assignee
Shandong Jianzhu University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Shandong Jianzhu University filed Critical Shandong Jianzhu University
Priority to CN202310322768.XA priority Critical patent/CN116047519B/zh
Publication of CN116047519A publication Critical patent/CN116047519A/zh
Application granted granted Critical
Publication of CN116047519B publication Critical patent/CN116047519B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B15/00—Measuring arrangements characterised by the use of electromagnetic waves or particle radiation, e.g. by the use of microwaves, X-rays, gamma rays or electrons
    • G01B15/06—Measuring arrangements characterised by the use of electromagnetic waves or particle radiation, e.g. by the use of microwaves, X-rays, gamma rays or electrons for measuring the deformation in a solid

Abstract

本发明公开一种基于合成孔径雷达干涉测量技术的选点方法，属于采用无线电波测距或测速技术领域，用于矿区多源InSAR技术融合的高相干点选取，包括：SAR影像预处理、获取高相干点、获取慢失相关滤波相位像素目标点、识别同质像元集合，获取分布式目标点DS并优化DS相位，联合PS、SDFP与DS点进行时序处理，实现矿区地表形变监测。本发明将三种方法的优势互补，解决了采用常规InSAR技术进行地面沉降监测过程中测量点数有限，无法准确描述研究区的详细形变情况的缺点，提升了矿区长期植被覆盖、房屋、道路等人工建筑物上测量点的监测密度，进而对矿区形变进行精确识别和监测。

Description

一种基于合成孔径雷达干涉测量技术的选点方法
技术领域
本发明公开一种基于合成孔径雷达干涉测量技术的选点方法，属于采用无线电波测距或测速技术领域。
背景技术
现有技术的INSAR选点方法，例如DA方法，一般要求有25-30景以上SAR数据支持,由于SAR影像对于同一地物不同时间成像时的视角、轨道、大气状况等因素不同,会造成同一地物在不同SAR影像上的振幅不具可比性,因此必须先进行辐射校正。但该方法并未考虑目标本身的幅度值,只考虑了幅度的变化量,对于本身幅度值较小、变化量很小的物体（如水体）得到的DA值也很小,容易造成误判，难以有效地获取地表的完整形变信息。
发明内容
本发明的目的在于提供一种基于合成孔径雷达干涉测量技术的选点方法，以解决现有技术中，DA方法难以有效地获取地表的完整形变信息的问题。
一种基于合成孔径雷达干涉测量技术的选点方法，包括：将主影像分别与其他影像进行配准处理，获得干涉图，选取PS点、DS点和SDFP点，提取后向散射特性的分布式目标采用快速同质点选取方法进行同质像元的识别；通过时间相干系数阈值将PS、DS、SDFP点中相干系数低的点进行剔除，根据每个像元的信噪比计算每个像元所占的权重，按照加权平均的方式计算高相干点集的缠绕相位，并作为后续时序分析的对象，构成高相干分布式点目标的真实相位，将融合后得到的高相干分布式点目标进行常规时序InSAR技术处理。
配准处理前进行预处理，获取研究区域的SAR影像，将辅影像相对于主影像进行移动配准，对主影像和辅影像进行基线估计处理，通过计算卫星轨道参数获得基线，预估干涉像对的质量，获取干涉对组合方式；将主影像和辅影像间对应的像元共轭相乘，将相乘结果中的强度信息组成干涉强度图，生成干涉图；利用数据网站获取的精密轨道数据参数直接处理干涉图，去除平地效应。
预处理完成后，根据PS-InSAR技术中遵循时间基线、空间基线和多普勒质心频率差三者整体最小原则，选取最优的公共主影像，将公共主影像分别与其他影像进行配准处理，获取N幅干涉图，以干涉图像素点的幅度离差指数、相干系数值以及振幅强度作为衡量永久散射体目标的选取依据，提取强散射特性的点即PS点；
将满足振幅阈值和幅度离散指数阈值的点视为永久散射体的候选点，估计每个候选点的相位稳定性，基于时间相干系数阈值来挑选PS点；
采用相干系数法求解时序相干系数值，根据时间相干系数阈值，排除研究区范围内失相干严重的像元，得出初选目标点，求解相干系数的公式为：
Figure SMS_1
；
式中
Figure SMS_2
与/>
Figure SMS_3
表示构成干涉对的两幅图像，/>
Figure SMS_4
是共轭相乘，/>
Figure SMS_5
表示相干系数，
Figure SMS_6
表示像元坐标，m和n表示局部窗口总行数和总列数；
相干系数的取值范围为[0，1]，设定时间相干系数阈值为0.8来识别干涉点目标，选取相干系数时间序列值均大于阈值的像素点作为候选 PS 点；
将候选点对应的时序振幅均值从小到大依次排列，同时找出中间值作为精确选点的振幅阈值，从初选目标点中将其周围的误选点剔除，筛选出有效点，振幅阈值计算公式为：
Figure SMS_7
；
式中T为振幅阈值即所有影响振幅平均值的最小值，A为影像像元，像元坐标的振幅值且是经过辐射校正处理的，N是干涉对数量；
根据幅度平均值
Figure SMS_9
与标准差/>
Figure SMS_10
计算振幅离散指数/>
Figure SMS_11
，设定一个幅度离散指数阈值/>
Figure SMS_12
为0.4，比较/>
Figure SMS_13
与/>
Figure SMS_14
的大小，若/>
Figure SMS_15
小于/>
Figure SMS_8
时，该像元最终被去确定为PS点，否则为非PS点。
进行同质像元的识别包括：
快速同质点选取方法包括：基于SAR影像数据集任一像元时间维度上振幅的平均值服从高斯分布的假设，将假设检验方法转为置信区间估计，通过逻辑计算来判断两个像元是否服从相同函数分布，达到像元同质点识别的目的；
在N幅SAR影像组成的数据集中，任一像元L在时间维度上的时序振幅为：
Figure SMS_16
，平均振幅/>
Figure SMS_17
表示为
Figure SMS_18
；
随着样本数目N的增加，根据中心极限定理，
Figure SMS_19
逐步趋于高斯分布，当/>
Figure SMS_20
服从高斯分布时，/>
Figure SMS_21
的区间估计如下公式表示：
Figure SMS_22
；
其中
Figure SMS_23
表示概率，/>
Figure SMS_24
为标准正态分布中置信度为/>
Figure SMS_25
时对应的分为点，
Figure SMS_26
为像元 L 时间维度上幅度的真实方差；
在SAR影像数据服从正态分布的条件下，平均幅度的分布规律遵循瑞利分布，利用N幅时间序列的SAR影像计算的参考像元的平均幅度和待检测像元的平均幅度，根据标准正态分布的定义，平均幅度的置信区间为：
Figure SMS_27
；
其中L表示视数，N是干涉对数量，
Figure SMS_28
和/>
Figure SMS_29
分别表示参考像素和任意邻域像素的时间均值，/>
Figure SMS_30
表示标准正态分布中分位点，将/>
Figure SMS_31
当作像元L时间维度上强度均值的真值，即/>
Figure SMS_32
，则平均幅度的置信区间公式即为一个确定的区间，设置/>
Figure SMS_33
=0.5提炼中心像素均值，通过计算待估计像元时间维度上的平均强度值是否落入目标像元对应的区间，判断待估计像元与目标像元是否属于同质点；
依据同质像元集合估计相干矩阵并进行特征值分解，实现分布式目标相位的优化估计，求解最大特征值
Figure SMS_34
对应特征向量的相位分量/>
Figure SMS_35
作为优化相位估计值，求解最大特征值/>
Figure SMS_36
对应特征向量/>
Figure SMS_37
的表达式为：
Figure SMS_38
；
式中，argmax表示获取最大特征值对应的特征向量，
Figure SMS_39
为通过相干矩阵特征值分解方法相位优化方法从N(N-1)/2个多视干涉相位求解出的一组最优拟合相位估计值；同质像元集合对应的相干矩阵/>
Figure SMS_40
为在不同散射机制类型相互作用下得到的相干矩阵叠加总和；/>
Figure SMS_41
表示最大特征值对应的特征向量；/>
Figure SMS_42
表示共轭转置；
通过计算分布式散射体候选点时间相干性，对原始相干矩阵中自适应多视干涉相位与优化相位得到的干涉相位进行拟合比较，将拟合优度值高于0.7的像元最终确定为最终的分布式散射体目标点，分布式散射体候选点时间相干性计算公式为：
Figure SMS_43
；
式中，
Figure SMS_44
为分布式散射体的时间相干性，/>
Figure SMS_45
为第m幅影像与主影像分布式散射体的优化相位，/>
Figure SMS_46
为第n幅影像与主影像分布式散射体的优化相位，/>
Figure SMS_47
为第m幅和第n幅影像优化前差分干涉相位，N为SAR影像数量。
完成同质像元的识别后，根据SBAS-InSAR技术中短时空基线原则，将所有SAR影像分成多个短基线子集，每个子集有单独的主影像，影像之间进行差分干涉处理，获取干涉图，根据幅度法、信噪比法以及相干系数法获取慢失相关滤波相位像素目标点即SDFP点。
综合考虑PS、DS、SDFP点的信噪比，并通过时间相干系数阈值将PS、DS、SDFP点中相干系数低的点进行剔除，满足时间相干系数阈值要求的点构成高相干分布式点目标，根据每个像元的信噪比计算每个像元所占的权重，按照加权平均的方式计算高相干点集的缠绕相位，并作为后续时序分析的对象，构成高相干分布式点目标的真实相位；
根据SBAS干涉组合方式，重新计算PS、DS点的
Figure SMS_48
的残余相位变化，获取PS、DS从SBAS-InSAR干涉组合中的缠绕相位；
根据信噪比计算公式计算每个像元的信噪比，并通过时间相干系数阈值将PS、DS、SDFP点中相干系数低的点进行剔除，满足阈值要求的点构成StaMPS时序分析的高相干点集，信噪比计算公式为：
Figure SMS_49
；
Figure SMS_50
表示信噪比；/>
Figure SMS_51
为像元x处的相干值；
根据每个像元的信噪比计算每个像元所占的权重，按照加权平均的方式计算高相干点集的缠绕相位，并作为后续时序分析的对象：
Figure SMS_52
；
式中，W为每个像元的信噪比计算每个像元所占的权重；
计算高相干点集的缠绕相位:
Figure SMS_53
；
式中
Figure SMS_54
为重合点高相干像元最终缠绕相位，/>
Figure SMS_55
是PS点的信噪比，/>
Figure SMS_56
是SDFP点的信噪比，/>
Figure SMS_57
是DS点的信噪比，/>
Figure SMS_58
是PS点的缠绕相位，/>
Figure SMS_59
是SDFP点的缠绕相位，
Figure SMS_60
是DS点的缠绕相位。
高相干点集的缠绕相位后，将融合后得到的高相干分布式点目标进行常规时序InSAR技术处理，利用高相干点建立线性形变速率和高程误差的线性模型方程组，使用SVD算法求解模型参数获得线性形变相位；去除线性形变相位后，残余误差中包含大气相位、非线性形变和噪声，根据各分量的特征对残差进行时间域的高通滤波和空间域的低通滤波操作去除大气相位分量和噪声，获得非线性形变相位；使用SVD分解获得各相干点的形变速率并进行时间域积分获取累积形变值，最终得到矿区采空区的地表形变监测结果。
得到矿区采空区的地表形变监测结果后，根据采空区地表形变监测结果得到年平均沉降速率值及时序沉降值，识别形变区域，并结合矿区开采深度，开采时间进行预测分析。
相对比现有技术，本发明具有以下有益效果：在常规时序InSAR技术中融入慢失相关滤波相位像素目标（SDFP）与分布式目标点(DS)，解决了采用常规InSAR技术进行地面沉降监测过程中测量点数有限，无法准确描述研究区的详细形变情况的缺点，极大的提升了矿区长期植被覆盖、房屋、道路等人工建筑物上测量点的监测密度。在获取全面的地表形变结果的基础上，再结合矿区开采数据进行分析，可以得到生产活动对地表形变结果的影响情况，从而实现对矿区开采区可能存在的隐患进行及时的识别和判断，有效的对高风险区域进行重点排查。
附图说明
图1是本发明的技术流程图；
图2是高相干点数量图；
图3为PS点的高相干点点位分布图；
图4为SDFP点的高相干点点位分布图；
图5为DS点的高相干点点位分布图；
图6为PS+SDFP+DS点的高相干点点位分布图；
图7为研究区内获取5个监测点时间序列序沉降量对比图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
一种基于合成孔径雷达干涉测量技术的选点方法，包括：将主影像分别与其他影像进行配准处理，获得干涉图，选取PS点、DS点和SDFP点，提取后向散射特性的分布式目标采用快速同质点选取方法进行同质像元的识别；通过时间相干系数阈值将PS、DS、SDFP点中相干系数低的点进行剔除，根据每个像元的信噪比计算每个像元所占的权重，按照加权平均的方式计算高相干点集的缠绕相位，并作为后续时序分析的对象，构成高相干分布式点目标的真实相位，将融合后得到的高相干分布式点目标进行常规时序InSAR技术处理。
配准处理前进行预处理，获取研究区域的SAR影像，将辅影像相对于主影像进行移动配准，对主影像和辅影像进行基线估计处理，通过计算卫星轨道参数获得基线，预估干涉像对的质量，获取干涉对组合方式；将主影像和辅影像间对应的像元共轭相乘，将相乘结果中的强度信息组成干涉强度图，生成干涉图；利用数据网站获取的精密轨道数据参数直接处理干涉图，去除平地效应。
预处理完成后，根据PS-InSAR技术中遵循时间基线、空间基线和多普勒质心频率差三者整体最小原则，选取最优的公共主影像，将公共主影像分别与其他影像进行配准处理，获取N幅干涉图，以干涉图像素点的幅度离差指数、相干系数值以及振幅强度作为衡量永久散射体目标的选取依据，提取强散射特性的点即PS点；
将满足振幅阈值和幅度离散指数阈值的点视为永久散射体的候选点，估计每个候选点的相位稳定性，基于时间相干系数阈值来挑选PS点；
采用相干系数法求解时序相干系数值，根据时间相干系数阈值，排除研究区范围内失相干严重的像元，得出初选目标点，求解相干系数的公式为：
Figure SMS_61
；
式中
Figure SMS_62
与/>
Figure SMS_63
表示构成干涉对的两幅图像，/>
Figure SMS_64
是共轭相乘，/>
Figure SMS_65
表示相干系数，
Figure SMS_66
表示像元坐标，m和n表示局部窗口总行数和总列数；
相干系数的取值范围为[0，1]，设定时间相干系数阈值为0.8来识别干涉点目标，选取相干系数时间序列值均大于阈值的像素点作为候选 PS 点；
将候选点对应的时序振幅均值从小到大依次排列，同时找出中间值作为精确选点的振幅阈值，从初选目标点中将其周围的误选点剔除，筛选出有效点，振幅阈值计算公式为：
Figure SMS_67
；
式中T为振幅阈值即所有影响振幅平均值的最小值，A为影像像元，像元坐标的振幅值且是经过辐射校正处理的，N是干涉对数量；
根据幅度平均值
Figure SMS_69
与标准差/>
Figure SMS_70
计算振幅离散指数/>
Figure SMS_71
，设定一个幅度离散指数阈值/>
Figure SMS_72
为0.4，比较/>
Figure SMS_73
与/>
Figure SMS_74
的大小，若/>
Figure SMS_75
小于/>
Figure SMS_68
时，该像元最终被去确定为PS点，否则为非PS点。
进行同质像元的识别包括：
快速同质点选取方法包括：基于SAR影像数据集任一像元时间维度上振幅的平均值服从高斯分布的假设，将假设检验方法转为置信区间估计，通过逻辑计算来判断两个像元是否服从相同函数分布，达到像元同质点识别的目的；
在N幅SAR影像组成的数据集中，任一像元L在时间维度上的时序振幅为：
Figure SMS_76
，平均振幅/>
Figure SMS_77
表示为
Figure SMS_78
；
随着样本数目N的增加，根据中心极限定理，
Figure SMS_79
逐步趋于高斯分布，当/>
Figure SMS_80
服从高斯分布时，/>
Figure SMS_81
的区间估计如下公式表示：
Figure SMS_82
；
其中
Figure SMS_83
表示概率，/>
Figure SMS_84
为标准正态分布中置信度为/>
Figure SMS_85
时对应的分为点，
Figure SMS_86
为像元 L 时间维度上幅度的真实方差；
在SAR影像数据服从正态分布的条件下，平均幅度的分布规律遵循瑞利分布，利用N幅时间序列的SAR影像计算的参考像元的平均幅度和待检测像元的平均幅度，根据标准正态分布的定义，平均幅度的置信区间为：
Figure SMS_87
；
其中L表示视数，N是干涉对数量，
Figure SMS_88
和/>
Figure SMS_89
分别表示参考像素和任意邻域像素的时间均值，/>
Figure SMS_90
表示标准正态分布中分位点，将/>
Figure SMS_91
当作像元L时间维度上强度均值的真值，即/>
Figure SMS_92
，则平均幅度的置信区间公式即为一个确定的区间，设置/>
Figure SMS_93
=0.5提炼中心像素均值，通过计算待估计像元时间维度上的平均强度值是否落入目标像元对应的区间，判断待估计像元与目标像元是否属于同质点；
依据同质像元集合估计相干矩阵并进行特征值分解，实现分布式目标相位的优化估计，求解最大特征值
Figure SMS_94
对应特征向量的相位分量/>
Figure SMS_95
作为优化相位估计值，求解最大特征值/>
Figure SMS_96
对应特征向量/>
Figure SMS_97
的表达式为：
Figure SMS_98
；
式中，argmax表示获取最大特征值对应的特征向量，
Figure SMS_99
为通过相干矩阵特征值分解方法相位优化方法从N(N-1)/2个多视干涉相位求解出的一组最优拟合相位估计值；同质像元集合对应的相干矩阵/>
Figure SMS_100
为在不同散射机制类型相互作用下得到的相干矩阵叠加总和；/>
Figure SMS_101
表示最大特征值对应的特征向量；/>
Figure SMS_102
表示共轭转置；
通过计算分布式散射体候选点时间相干性，对原始相干矩阵中自适应多视干涉相位与优化相位得到的干涉相位进行拟合比较，将拟合优度值高于0.7的像元最终确定为最终的分布式散射体目标点，分布式散射体候选点时间相干性计算公式为：
Figure SMS_103
；
式中，
Figure SMS_104
为分布式散射体的时间相干性，/>
Figure SMS_105
为第m幅影像与主影像分布式散射体的优化相位，/>
Figure SMS_106
为第n幅影像与主影像分布式散射体的优化相位，/>
Figure SMS_107
为第m幅和第n幅影像优化前差分干涉相位，N为SAR影像数量。
完成同质像元的识别后，根据SBAS-InSAR技术中短时空基线原则，将所有SAR影像分成多个短基线子集，每个子集有单独的主影像，影像之间进行差分干涉处理，获取干涉图，根据幅度法、信噪比法以及相干系数法获取慢失相关滤波相位像素目标点即SDFP点。
综合考虑PS、DS、SDFP点的信噪比，并通过时间相干系数阈值将PS、DS、SDFP点中相干系数低的点进行剔除，满足时间相干系数阈值要求的点构成高相干分布式点目标，根据每个像元的信噪比计算每个像元所占的权重，按照加权平均的方式计算高相干点集的缠绕相位，并作为后续时序分析的对象，构成高相干分布式点目标的真实相位；
根据SBAS干涉组合方式，重新计算PS、DS点的
Figure SMS_108
的残余相位变化，获取PS、DS从SBAS-InSAR干涉组合中的缠绕相位；
根据信噪比计算公式计算每个像元的信噪比，并通过时间相干系数阈值将PS、DS、SDFP点中相干系数低的点进行剔除，满足阈值要求的点构成StaMPS时序分析的高相干点集，信噪比计算公式为：
Figure SMS_109
；
Figure SMS_110
表示信噪比；/>
Figure SMS_111
为像元x处的相干值；
根据每个像元的信噪比计算每个像元所占的权重，按照加权平均的方式计算高相干点集的缠绕相位，并作为后续时序分析的对象：
Figure SMS_112
；
式中，W为每个像元的信噪比计算每个像元所占的权重；
计算高相干点集的缠绕相位:
Figure SMS_113
；
式中
Figure SMS_114
为重合点高相干像元最终缠绕相位，/>
Figure SMS_115
是PS点的信噪比，/>
Figure SMS_116
是SDFP点的信噪比，/>
Figure SMS_117
是DS点的信噪比，/>
Figure SMS_118
是PS点的缠绕相位，/>
Figure SMS_119
是SDFP点的缠绕相位，
Figure SMS_120
是DS点的缠绕相位。
高相干点集的缠绕相位后，将融合后得到的高相干分布式点目标进行常规时序InSAR技术处理，利用高相干点建立线性形变速率和高程误差的线性模型方程组，使用SVD算法求解模型参数获得线性形变相位；去除线性形变相位后，残余误差中包含大气相位、非线性形变和噪声，根据各分量的特征对残差进行时间域的高通滤波和空间域的低通滤波操作去除大气相位分量和噪声，获得非线性形变相位；使用SVD分解获得各相干点的形变速率并进行时间域积分获取累积形变值，最终得到矿区采空区的地表形变监测结果。
得到矿区采空区的地表形变监测结果后，根据采空区地表形变监测结果得到年平均沉降速率值及时序沉降值，识别形变区域，并结合矿区开采深度，开采时间进行预测分析。
本发明中涉及的英文简写具体如下：PS(Persistent Scatterer)永久散射体；SDFP(slowly-decorrelating filtered phase,SDFP)滤波相位失相干缓慢目标点；DS(Distributed Scatterer) 分布式散射体；(Interferometric Synthetic ApertureRadar，In SAR) 合成孔径雷达干涉测量技术；合成孔径雷达(Synthetic Aperture Radar,SAR）；奇异值分解方法（Singular ValueDecomposition,SVD）。
本发明的技术流程图如图1所示，图2为研究区内不同方法获取监测点总数对比图，同一研究区域内，获取的PS点数量为21425，SDFP点数量为31539，DS点数量为40253，三种点融合之后点数量为92038；点的分布状况分别如图3、图4、图5、图6所示，其中PS点的高相干点点位分布图如图3所示，SDFP点的高相干点点位分布图如图4所示，DS点的高相干点点位分布图如图5所示，PS+SDFP+DS点的高相干点点位分布图如图6所示。
四种不同方法获取沉降速率值分布如表1所示。
表1 四种不同方法获取沉降速率值分布
Figure SMS_121
。
图7为研究区内获取5个监测点（C1、C2、C3、C4、C5）时间序列序沉降量，为更加详细分析研究区沉降情况，提取沉降区内沉降速率为-6.5mm/yr~-16.5mm/yr之间的五个点，获取时序累计沉降量，看出5个点表现为持续沉降，2021年04月24日的累计沉降量分别为-8mm、-13mm、-18mm、-19mm、-19mm。
以上实施例仅用于说明本发明的技术方案，而非对其限制，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换，而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (5)

1.一种基于合成孔径雷达干涉测量技术的选点方法，其特征在于，包括：将主影像分别与其他影像进行配准处理，获得干涉图，选取PS点、DS点和SDFP点，提取后向散射特性的分布式目标采用快速同质点选取方法进行同质像元的识别；通过时间相干系数阈值将PS、DS、SDFP点中相干系数低的点进行剔除，根据每个像元的信噪比计算每个像元所占的权重，按照加权平均的方式计算高相干点集的缠绕相位，并作为后续时序分析的对象，构成高相干分布式点目标的真实相位，将融合后得到的高相干分布式点目标进行常规时序InSAR技术处理；
配准处理前进行预处理，获取研究区域的SAR影像，将辅影像相对于主影像进行移动配准，对主影像和辅影像进行基线估计处理，通过计算卫星轨道参数获得基线，预估干涉像对的质量，获取干涉对组合方式；将主影像和辅影像间对应的像元共轭相乘，将相乘结果中的强度信息组成干涉强度图，生成干涉图；利用数据网站获取的精密轨道数据参数直接处理干涉图，去除平地效应；
预处理完成后，根据PS-InSAR技术中遵循时间基线、空间基线和多普勒质心频率差三者整体最小原则，选取最优的公共主影像，将公共主影像分别与其他影像进行配准处理，获取N幅干涉图，以干涉图像素点的幅度离差指数、相干系数值以及振幅强度作为衡量永久散射体目标的选取依据，提取强散射特性的点即PS点；
将满足振幅阈值和幅度离散指数阈值的点视为永久散射体的候选点，估计每个候选点的相位稳定性，基于时间相干系数阈值来挑选PS点；
采用相干系数法求解时序相干系数值，根据时间相干系数阈值，排除研究区范围内失相干严重的像元，得出初选目标点，求解相干系数的公式为：
Figure QLYQS_1
；
式中
Figure QLYQS_2
与/>
Figure QLYQS_3
表示构成干涉对的两幅图像，/>
Figure QLYQS_4
是共轭相乘，/>
Figure QLYQS_5
表示相干系数，/>
Figure QLYQS_6
表示像元坐标，m和n表示局部窗口总行数和总列数；
相干系数的取值范围为[0，1]，设定时间相干系数阈值为0.8来识别干涉点目标，选取相干系数时间序列值均大于阈值的像素点作为候选 PS 点；
将候选点对应的时序振幅均值从小到大依次排列，同时找出中间值作为精确选点的振幅阈值，从初选目标点中将其周围的误选点剔除，筛选出有效点，振幅阈值计算公式为：
Figure QLYQS_7
；
式中T为振幅阈值即所有影响振幅平均值的最小值，A为影像像元，像元坐标的振幅值且是经过辐射校正处理的，N是干涉对数量；
根据幅度平均值
Figure QLYQS_9
与标准差/>
Figure QLYQS_10
计算振幅离散指数/>
Figure QLYQS_11
，设定一个幅度离散指数阈值
Figure QLYQS_12
为0.4，比较/>
Figure QLYQS_13
与/>
Figure QLYQS_14
的大小，若/>
Figure QLYQS_15
小于/>
Figure QLYQS_8
时，该像元最终被去确定为PS点，否则为非PS点；
进行同质像元的识别包括：
快速同质点选取方法包括：基于SAR影像数据集任一像元时间维度上振幅的平均值服从高斯分布的假设，将假设检验方法转为置信区间估计，通过逻辑计算来判断两个像元是否服从相同函数分布，达到像元同质点识别的目的；
在N幅SAR影像组成的数据集中，任一像元L在时间维度上的时序振幅为：
Figure QLYQS_16
，平均振幅/>
Figure QLYQS_17
表示为
Figure QLYQS_18
；
随着样本数目N的增加，根据中心极限定理，
Figure QLYQS_19
逐步趋于高斯分布，当/>
Figure QLYQS_20
服从高斯分布时，/>
Figure QLYQS_21
的区间估计如下公式表示：
Figure QLYQS_22
；
其中
Figure QLYQS_23
表示概率，/>
Figure QLYQS_24
为标准正态分布中置信度为/>
Figure QLYQS_25
时对应的分为点，
Figure QLYQS_26
为像元 L 时间维度上幅度的真实方差；
在SAR影像数据服从正态分布的条件下，平均幅度的分布规律遵循瑞利分布，利用N幅时间序列的SAR影像计算的参考像元的平均幅度和待检测像元的平均幅度，根据标准正态分布的定义，平均幅度的置信区间为：
Figure QLYQS_27
；
其中L表示视数，N是干涉对数量，
Figure QLYQS_28
和/>
Figure QLYQS_29
分别表示参考像素和任意邻域像素的时间均值，/>
Figure QLYQS_30
表示标准正态分布中分位点，将/>
Figure QLYQS_31
当作像元L时间维度上强度均值的真值，即/>
Figure QLYQS_32
，则平均幅度的置信区间公式即为一个确定的区间，设置/>
Figure QLYQS_33
=0.5提炼中心像素均值，通过计算待估计像元时间维度上的平均强度值是否落入目标像元对应的区间，判断待估计像元与目标像元是否属于同质点；
依据同质像元集合估计相干矩阵并进行特征值分解，实现分布式目标相位的优化估计，求解最大特征值
Figure QLYQS_34
对应特征向量的相位分量/>
Figure QLYQS_35
作为优化相位估计值，求解最大特征值/>
Figure QLYQS_36
对应特征向量/>
Figure QLYQS_37
的表达式为：
Figure QLYQS_38
；
式中，argmax表示获取最大特征值对应的特征向量，
Figure QLYQS_39
为通过相干矩阵特征值分解方法相位优化方法从N(N-1)/2个多视干涉相位求解出的一组最优拟合相位估计值；同质像元集合对应的相干矩阵/>
Figure QLYQS_40
为在不同散射机制类型相互作用下得到的相干矩阵叠加总和；
Figure QLYQS_41
表示最大特征值对应的特征向量；/>
Figure QLYQS_42
表示共轭转置；
通过计算分布式散射体候选点时间相干性，对原始相干矩阵中自适应多视干涉相位与优化相位得到的干涉相位进行拟合比较，将拟合优度值高于0.7的像元最终确定为最终的分布式散射体目标点，分布式散射体候选点时间相干性计算公式为：
Figure QLYQS_43
；
式中，
Figure QLYQS_44
为分布式散射体的时间相干性，/>
Figure QLYQS_45
为第m幅影像与主影像分布式散射体的优化相位，/>
Figure QLYQS_46
为第n幅影像与主影像分布式散射体的优化相位，/>
Figure QLYQS_47
为第m幅和第n幅影像优化前差分干涉相位，N为SAR影像数量。
2.根据权利要求1所述的基于合成孔径雷达干涉测量技术的选点方法，其特征在于，完成同质像元的识别后，根据SBAS-InSAR技术中短时空基线原则，将所有SAR影像分成多个短基线子集，每个子集有单独的主影像，影像之间进行差分干涉处理，获取干涉图，根据幅度法、信噪比法以及相干系数法获取慢失相关滤波相位像素目标点即SDFP点。
3.根据权利要求2所述的基于合成孔径雷达干涉测量技术的选点方法，其特征在于，综合考虑PS、DS、SDFP点的信噪比，并通过时间相干系数阈值将PS、DS、SDFP点中相干系数低的点进行剔除，满足时间相干系数阈值要求的点构成高相干分布式点目标，根据每个像元的信噪比计算每个像元所占的权重，按照加权平均的方式计算高相干点集的缠绕相位，并作为后续时序分析的对象，构成高相干分布式点目标的真实相位；
根据SBAS干涉组合方式，重新计算PS、DS点的
Figure QLYQS_48
的残余相位变化，获取PS、DS从SBAS-InSAR干涉组合中的缠绕相位；
根据信噪比计算公式计算每个像元的信噪比，并通过时间相干系数阈值将PS、DS、SDFP点中相干系数低的点进行剔除，满足阈值要求的点构成StaMPS时序分析的高相干点集，信噪比计算公式为：
Figure QLYQS_49
；
Figure QLYQS_50
表示信噪比；/>
Figure QLYQS_51
为像元x处的相干值；
根据每个像元的信噪比计算每个像元所占的权重，按照加权平均的方式计算高相干点集的缠绕相位，并作为后续时序分析的对象：
Figure QLYQS_52
；
式中，W为每个像元的信噪比计算每个像元所占的权重；
计算高相干点集的缠绕相位:
Figure QLYQS_53
；
式中
Figure QLYQS_54
为重合点高相干像元最终缠绕相位，/>
Figure QLYQS_55
是PS点的信噪比，/>
Figure QLYQS_56
是SDFP点的信噪比，/>
Figure QLYQS_57
是DS点的信噪比，/>
Figure QLYQS_58
是PS点的缠绕相位，/>
Figure QLYQS_59
是SDFP点的缠绕相位，/>
Figure QLYQS_60
是DS点的缠绕相位。
4.根据权利要求3所述的基于合成孔径雷达干涉测量技术的选点方法，其特征在于，高相干点集的缠绕相位后，将融合后得到的高相干分布式点目标进行常规时序InSAR技术处理，利用高相干点建立线性形变速率和高程误差的线性模型方程组，使用SVD算法求解模型参数获得线性形变相位；去除线性形变相位后，残余误差中包含大气相位、非线性形变和噪声，根据各分量的特征对残差进行时间域的高通滤波和空间域的低通滤波操作去除大气相位分量和噪声，获得非线性形变相位；使用SVD分解获得各相干点的形变速率并进行时间域积分获取累积形变值，最终得到矿区采空区的地表形变监测结果。
5.根据权利要求4所述的基于合成孔径雷达干涉测量技术的选点方法，其特征在于，得到矿区采空区的地表形变监测结果后，根据采空区地表形变监测结果得到年平均沉降速率值及时序沉降值，识别形变区域，并结合矿区开采深度，开采时间进行预测分析。
CN202310322768.XA 2023-03-30 2023-03-30 一种基于合成孔径雷达干涉测量技术的选点方法 Active CN116047519B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310322768.XA CN116047519B (zh) 2023-03-30 2023-03-30 一种基于合成孔径雷达干涉测量技术的选点方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310322768.XA CN116047519B (zh) 2023-03-30 2023-03-30 一种基于合成孔径雷达干涉测量技术的选点方法

Publications (2)

Publication Number Publication Date
CN116047519A CN116047519A (zh) 2023-05-02
CN116047519B true CN116047519B (zh) 2023-06-16

Family

ID=86129845

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310322768.XA Active CN116047519B (zh) 2023-03-30 2023-03-30 一种基于合成孔径雷达干涉测量技术的选点方法

Country Status (1)

Country Link
CN (1) CN116047519B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116482684B (zh) * 2023-06-21 2023-08-22 深圳市城市公共安全技术研究院有限公司 区域五维成像方法、装置、设备及存储介质
CN116736305B (zh) * 2023-08-14 2023-10-27 北京观微科技有限公司 地物形变确定方法、装置、电子设备及存储介质
CN116908853B (zh) * 2023-09-13 2023-11-17 北京观微科技有限公司 高相干点选取方法、装置和设备

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114325706A (zh) * 2021-12-31 2022-04-12 中国科学院空天信息创新研究院 一种分布式散射体滤波方法

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6046695A (en) * 1996-07-11 2000-04-04 Science Application International Corporation Phase gradient auto-focus for SAR images
EP2304466B1 (en) * 2008-07-04 2015-10-28 Telespazio S.p.A. Identification and analysis of persistent scatterers in series of sar images
IT1394733B1 (it) * 2009-07-08 2012-07-13 Milano Politecnico Procedimento per il filtraggio di interferogrammi generati da immagini sar acquisite sulla stessa area.
CN104123464B (zh) * 2014-07-23 2017-02-22 中国国土资源航空物探遥感中心 一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法
CN106526590B (zh) * 2016-11-04 2018-09-25 山东科技大学 一种融合多源sar影像工矿区三维地表形变监测及解算方法
CN112014841A (zh) * 2020-08-31 2020-12-01 中国矿业大学 一种基于DS-InSAR技术监测油田区地表形变分析方法
CN113960595A (zh) * 2021-09-24 2022-01-21 中国科学院深圳先进技术研究院 一种地表形变监测方法及系统
CN115683094A (zh) * 2022-03-17 2023-02-03 山东建筑大学 一种复杂环境下车载双天线紧耦合定位方法及系统
CN115825955A (zh) * 2022-11-30 2023-03-21 中国矿业大学 一种基于相干矩阵自适应分解的极化时序InSAR方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114325706A (zh) * 2021-12-31 2022-04-12 中国科学院空天信息创新研究院 一种分布式散射体滤波方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
利用时序SAR影像集监测开采沉陷的试验研究;范洪冬;邓喀中;薛继群;祝传广;;煤矿安全(第02期);全文 *

Also Published As

Publication number Publication date
CN116047519A (zh) 2023-05-02

Similar Documents

Publication Publication Date Title
CN116047519B (zh) 一种基于合成孔径雷达干涉测量技术的选点方法
Almar et al. Wave-derived coastal bathymetry from satellite video imagery: A showcase with Pleiades persistent mode
Wang et al. Retrieval of phase history parameters from distributed scatterers in urban areas using very high resolution SAR data
CN113960595A (zh) 一种地表形变监测方法及系统
CN104123464A (zh) 一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法
CN109388887B (zh) 一种地面沉降影响因素定量分析方法及系统
Baier et al. A nonlocal InSAR filter for high-resolution DEM generation from TanDEM-X interferograms
CN106204539B (zh) 一种基于形态学梯度的反演城区建筑物沉降的方法
CN113281749A (zh) 一种顾及同质性的时序InSAR高相干点选取方法
Ma et al. A sequential approach for Sentinel-1 TOPS time-series co-registration over low coherence scenarios
Painam et al. A comprehensive review of SAR image filtering techniques: systematic survey and future directions
Bao et al. An improved distributed scatterers extraction algorithm for monitoring tattered ground surface subsidence with DSInSAR: A case study of loess landform in Tongren county
CN116908853B (zh) 高相干点选取方法、装置和设备
CN112269176B (zh) 一种矿山地表沉陷早期识别监测方法
Refice et al. On the use of anisotropic covariance models in estimating atmospheric DInSAR contributions
Huang et al. InSAR time-series analysis with a non-Gaussian detector for persistent scatterers
Sadeghi et al. Monitoring land subsidence in a rural area using a combination of ADInSAR and polarimetric coherence optimization
Azadnejad et al. Extending polarimetric optimization of multi-temporal InSAR techniques on dual polarized Sentinel-1 data
Liu et al. NL-MMSE: A hybrid phase optimization method in multimaster interferogram stack for DS-InSAR applications
Wang et al. Research on GBSAR deformation extraction method based on three-threshold
CN113625241A (zh) 差异沉降监控预警方法
Yue et al. A global seamless DEM based on multi-source data fusion (GSDEM-30): Product generation and evaluation
Li et al. Removing Milky Way from airglow images using principal component analysis
CN111983609B (zh) 基于雷达遥感影像的湿芦苇提取方法
Wu et al. Ground subsidence monitoring of high speed roads by using PS-InSAR method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant