CN109613531B - 一种微变感知预警雷达的多阈值优化形变反演方法及系统 - Google Patents

一种微变感知预警雷达的多阈值优化形变反演方法及系统 Download PDF

Info

Publication number
CN109613531B
CN109613531B CN201910011969.1A CN201910011969A CN109613531B CN 109613531 B CN109613531 B CN 109613531B CN 201910011969 A CN201910011969 A CN 201910011969A CN 109613531 B CN109613531 B CN 109613531B
Authority
CN
China
Prior art keywords
time sequence
pixel
threshold value
deformation quantity
calculating
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
CN201910011969.1A
Other languages
English (en)
Other versions
CN109613531A (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.)
North China University of Technology
Original Assignee
North China University of Technology
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 North China University of Technology filed Critical North China University of Technology
Priority to CN201910011969.1A priority Critical patent/CN109613531B/zh
Publication of CN109613531A publication Critical patent/CN109613531A/zh
Application granted granted Critical
Publication of CN109613531B publication Critical patent/CN109613531B/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/886—Radar or analogous systems specially adapted for specific applications for alarm systems
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
    • G01B7/16—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Image Processing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开一种微变感知预警雷达的多阈值优化形变反演方法及系统，属于雷达探测技术领域，能够基于获取的雷达时序影像图自动的识别出最优阈值，提高了形变反演结果的精确度。该方法包括：S1，基于微变感知预警雷达依次获取N+1幅时序影像；S2，采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，将其定义为稳定点后计算形变量；S3，逐步放宽阈值条件，对应从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后计算形变量；S4，判断当前动态点形变量与稳定点形变量是否发生突变，若没发生突变则执行步骤S3，若发生突变则将当前阈值条件输出为最优阈值。

Description

一种微变感知预警雷达的多阈值优化形变反演方法及系统
技术领域
本发明涉及雷达探测技术领域，尤其涉及一种微变感知预警雷达的多阈值优化方法及系统。
背景技术
我国是世界上地质灾害发生最为频繁的国家之一，在我国常见的各类地址灾害中，滑坡数量占地质灾害总量的70％以上，是发生在山区最主要的地质灾害类型。对滑坡体展开形变监测可以客观真实的记录坡体形变的发展演变过程，对了解坡体演化规律及准确预测坡体发展趋势具有重要意义。微变感知预警雷达是地基雷达的一种，可用于边坡形变监测，具有非接触、高精度、大区域、全天时全天候、24小时连续监测的技术优势，是垮塌灾害监测预警的重要技术手段，可以极大的减少或避免因边坡位移形变引起的滑坡灾害对国家和人民生命财产造成的损失。
微波遥感虽然能穿透大气进行成像，但是大气会对形变结果的最终精度造成影响。永久散射体(permanent scatterers,PS)干涉测量技术是地基SAR抑制大气、水汽等扰动因素获取高精度形变信息的有效手段；在探测地形形变的过程中，PS点的选择对地形形变探测结果的精度影响尤为关键，现有的PS点的选择方法可分为单阈值筛选法和多阈值信息筛选法，其中，单阈值筛选法如相干系数法，由于只是基于PS点的某一特性进行筛选，容易造成错判或漏判，因此存在得到的地形形变探测结果精度差的缺陷，现有的多阈值信息筛选法，如相位误差指数-振幅指数结合法，虽然能够有效解决PS点筛选过程中存在的错判或漏判问题，但由于现有的多阈值信息筛选法需要在形变监测开始前通过人工反复选择或调试阈值，若人工设置的阈值不合理则会导致PS选取的质量和数量出现偏差，进而影响形变反演结果的精确度。
发明内容
本发明的目的在于提供一种微变感知预警雷达的多阈值优化方法及系统，能够基于获取的雷达时序影像图自动的识别出最优阈值，提高了形变反演结果的精确度。
为了实现上述目的，本发明的一方面提供一种微变感知预警雷达的多阈值优化形变反演方法，包括：
S1，基于微变感知预警雷达依次获取N+1幅时序影像，计算时序影像中各像素的平均相干系数、时间序列幅度离差指数以及相位误差指数；
S2，采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，将其定义为稳定点后计算形变量；
S3，逐步放宽阈值条件，对应从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后计算形变量；
S4，判断当前动态点形变量与稳定点形变量是否发生突变，若没发生突变则执行步骤S3，若发生突变则将当前阈值条件输出为最优阈值，执行步骤S5；
S5，依次对相邻时序影像中的各像素共轭相乘得到N幅干涉图，滤波处理后从中提取符合最优阈值的PS点相位差信息；
S6，根据各PS点的相位差信息计算对应的形变量，得到形变反演结果。
优选地，所述步骤S1中，计算时序影像中各像素的平均相干系数、时间序列幅度离差指数以及相位误差指数的方法包括：
获取时序影像中的各像素坐标，通过公式
Figure GDA0002559534490000021
依次计算第x个像素在相邻的时序影像M和时序影像S的相干系数，并对获得的N个γx值求均值得到第x个像素的平均相干系数，其中，γx表示第x个像素的相干系数，M(i,j)和S(i,j)表示相邻时序影像中第x个像素的坐标，m和n为滑动窗口大小；
采用公式
Figure GDA0002559534490000031
计算每一像素的时间序列幅度离差指数，其中，δA表示所述像素在N+1幅时序影像中时间序列的标准差，mA表示所述像素在N+1幅时序影像中时间序列的均值；
采用公式
Figure GDA0002559534490000032
计算每一像素的噪声相位，其中，
Figure GDA0002559534490000033
为噪声相位，
Figure GDA0002559534490000034
为干涉相位，Filter表示均值滤波。
较佳地，所述步骤S2中，采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，将其定义为稳定点后计算形变量的方法包括：
所述门限值条件包括相干系数T1、时间序列幅度离差指数T2、相位误差指数T3；
从时序影像中选取大于相干系数T1，小于时间序列幅度离差指数T2且小于相位误差指数T3的像素，将其定义为稳定点；
采用形变量公式
Figure GDA0002559534490000035
计算所述稳定点的形变量，其中，λ表示雷达电磁波波长，
Figure GDA0002559534490000036
为大气相位。
进一步地，所述步骤S3中，逐步放宽阈值条件，对应从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后计算形变量的方法包括：
阈值条件的放宽方法为将阈值条件相干系数T1逐步减小和/或时间序列幅度离差指数T2逐步增大；
从时序影像中选取大于放宽后相干系数T1、小于放宽后时间序列幅度离差指数T2且小于相位误差指数T3的像素，将其定义为动态点；
通过形变量公式计算所述动态点的形变量。
优选地，所述步骤S4中，判断当前动态点形变量与稳定点形变量是否发生突变的方法包括：
计算当前动态点形变量与稳定点形变量的差值Δt；
将差值Δt与预设的突变阈值做大小比较，当差值Δt大于预设的突变阈值时认为发生突变，当差值Δt小于或等于预设的突变阈值时认为未发生突变。
优选地，所述步骤S5中，依次对相邻时序影像中的各像素共轭相乘得到N幅干涉图，滤波处理后从中提取符合最优阈值的PS点相位差信息的方法包括：
定义相邻的时序影像M和时序影像S中像素的读写格式，
Figure GDA0002559534490000041
Figure GDA0002559534490000042
且S＝M+1，S≤N+1；
对应像素采用共轭相乘公式
Figure GDA0002559534490000043
计算第M个干涉图，其中，干涉相位
Figure GDA0002559534490000044
重复上述步骤依次对各相邻时序影像中的像素进行共轭相乘计算，对应得到N幅干涉图；
采用均值滤波法分别对N幅干涉图进行滤波处理，从中筛选出符合最优阈值条件的全部PS点；
依次将每个PS点对应的N个干涉相位
Figure GDA0002559534490000045
累加，得到各PS点的相位差信息。
较佳地，采用均值滤波法对干涉图进行滤波处理的方法包括：
将干涉图中的全部像素按照a*a的像素模板模块化划分；
分别计算各模块像素中的像素坐标均值，并将模块的像素坐标均值视为新像素的坐标实现滤波处理。
与现有技术相比，本发明提供的微变感知预警雷达的多阈值优化形变反演方法具有以下有益效果：
本发明提供的微变感知预警雷达的多阈值优化形变反演方法中，首先通过微变感知预警雷达依次获取N+1幅时序影像，并基于N+1幅时序影像计算其中每个相应像素的平均相干系数、时间序列幅度离差指数和相位误差指数，然后基于预先设定的门限值条件采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，并将筛选出的像素定义为稳定点后计算各自的形变量，为了得到最优的阈值条件，本发明采用逐步放宽门限值条件的方法重新从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后对应计算形变量，此时，为了验证当前放宽的门限值条件是否达到最优阈值条件，需要判断当前动态点形变量与稳定点形变量是否发生突变，若没发生突变则说明还未达到最优阈值条件，需要重返步骤S3继续放宽门限值条件，若发生了突变则说明当前阈值条件已被优化为最优阈值，之后，再依次对相邻时序影像中的各像素共轭相乘得到N幅干涉图，滤波处理后从中提取符合最优阈值的PS点相位差信息，进而得到形变反演结果。
可见，本发明采用严苛阈值筛选和逐步放宽阈值条件相结合的方式能够自动识别出最优阈值，相比较于现有技术中人工调试阈值条件来说具有高效率和高精度的特点，通过最优阈值能够从时序影像中获取尽可能多的正确PS点，提高了反演结果的精确度。
本发明的另一方面提供一种微变感知预警雷达的多阈值优化形变反演系统，应用于上述技术方案所述的微变感知预警雷达的多阈值优化形变反演方法中，所述系统包括：
处理单元，用于基于微变感知预警雷达依次获取N+1幅时序影像，计算时序影像中各像素的平均相干系数、时间序列幅度离差指数以及相位误差指数；
严苛筛选单元，用于采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，将其定义为稳定点后计算形变量；
阈值宽限单元，用于逐步放宽阈值条件，对应从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后计算形变量；
判断单元，用于判断当前动态点形变量与稳定点形变量是否发生突变，若没发生突变则继续启动阈值宽限单元继续阈值宽限单元，若发生突变则将当前阈值条件输出为最优阈值；
滤波单元，用于依次对相邻时序影像中的各像素共轭相乘得到N幅干涉图，滤波处理后从中提取符合最优阈值的PS点相位差信息；
结果输出单元，用于根据各PS点的相位差信息计算对应的形变量，得到形变反演结果。
优选地，所述处理单元包括：
平均相干系数计算模块，用于获取时序影像中的各像素坐标，通过公式
Figure GDA0002559534490000061
依次计算第x个像素在相邻的时序影像M和时序影像S的相干系数，并对获得的N个γx值求均值得到第x个像素的平均相干系数，其中，γx表示第x个像素的相干系数，M(i,j)和S(i,j)表示相邻时序影像中第x个像素的坐标，m和n为滑动窗口大小；
时间序列标准差计算模块，用于采用公式
Figure GDA0002559534490000062
计算每一像素的时间序列幅度离差指数，其中，δA表示所述像素在N+1幅时序影像中时间序列的标准差，mA表示所述像素在N+1幅时序影像中时间序列的均值；
相位误差计算模块，用于采用公式
Figure GDA0002559534490000063
计算每一像素的噪声相位，其中，
Figure GDA0002559534490000064
为噪声相位，
Figure GDA0002559534490000065
为干涉相位，Filter表示均值滤波。
优选地，所述阈值宽限单元包括：
阈值条件设置模块，设置的阈值条件放宽方法为将阈值条件相干系数T1逐步减小和/或时间序列幅度离差指数T2逐步增大；
动态点筛选模块，用于从时序影像中选取大于放宽后相干系数T1、小于放宽后时间序列幅度离差指数T2且小于相位误差指数T3的像素，将其定义为动态点；
动态点形变量计算模块，用于通过形变量公式计算所述动态点的形变量。
与现有技术相比，本发明提供的微变感知预警雷达的多阈值优化形变反演系统的有益效果与上述技术方案提供的微变感知预警雷达的多阈值优化形变反演方法的有益效果相同，在此不做赘述。
附图说明
此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1为本发明实施例一中微变感知预警雷达的多阈值优化形变反演方法的流程示意图；
图2为本发明实施例一中基于多幅时序影像处理得到干涉图的示例图；
图3为本发明实施例二中微变感知预警雷达的多阈值优化形变反演系统的结构框图。
附图标记：
1-处理单元， 2-严苛筛选单元；
3-阈值宽限单元， 4-判断单元；
5-滤波单元， 6-结果输出单元。
具体实施方式
为使本发明的上述目的、特征和优点能够更加明显易懂，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其它实施例，均属于本发明保护的范围。
实施例一
图1为本发明实施例一中微变感知预警雷达的多阈值优化形变反演方法流程示意图。请参阅图1，本实施例提供一种微变感知预警雷达的多阈值优化形变反演方法，包括：
S1，基于微变感知预警雷达依次获取N+1幅时序影像，计算时序影像中各像素的平均相干系数、时间序列幅度离差指数以及相位误差指数；S2，采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，将其定义为稳定点后计算形变量；S3，逐步放宽阈值条件，对应从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后计算形变量；S4，判断当前动态点形变量与稳定点形变量是否发生突变，若没发生突变则执行步骤S3，若发生突变则将当前阈值条件输出为最优阈值，执行步骤S5；S5，依次对相邻时序影像中的各像素共轭相乘得到N幅干涉图，滤波处理后从中提取符合最优阈值的PS点相位差信息；S6，根据各PS点的相位差信息计算对应的形变量，得到形变反演结果。
本实施例提供的微变感知预警雷达的多阈值优化形变反演方法中，首先通过微变感知预警雷达依次获取N+1幅时序影像，并基于N+1幅时序影像计算其中每个相应像素的平均相干系数、时间序列幅度离差指数和相位误差指数，然后基于预先设定的门限值条件采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，并将筛选出的像素定义为稳定点后计算各自的形变量，为了得到最优的阈值条件，本实施例采用逐步放宽门限值条件的方法重新从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后对应计算形变量，此时，为了验证当前放宽的门限值条件是否达到最优阈值条件，需要判断当前动态点形变量与稳定点形变量是否发生突变，若没发生突变则说明还未达到最优阈值条件，需要重返步骤S3继续放宽门限值条件，若发生了突变则说明当前阈值条件已被优化为最优阈值，之后，再依次对相邻时序影像中的各像素共轭相乘得到N幅干涉图，滤波处理后从中提取符合最优阈值的PS点相位差信息，进而得到形变反演结果。
可见，本实施例采用严苛阈值筛选和逐步放宽阈值条件相结合的方式能够自动识别出最优阈值，相比较于现有技术中人工调试阈值条件来说具有高效率和高精度的特点，通过最优阈值能够从时序影像中获取尽可能多的正确PS点，提高了反演结果的精确度。
具体地，上述实施例步骤S1中，计算时序影像中各像素的平均相干系数、时间序列幅度离差指数以及相位误差指数的方法包括：
获取时序影像中的各像素坐标，通过公式
Figure GDA0002559534490000081
依次计算第x个像素在相邻的时序影像M和时序影像S的相干系数，并对获得的N个γx值求均值得到第x个像素的平均相干系数，其中，γx表示第x个像素的相干系数，M(i,j)和S(i,j)表示相邻时序影像中第x个像素的坐标，m和n为滑动窗口大小；采用公式
Figure GDA0002559534490000091
计算每一像素的时间序列幅度离差指数，其中，δA表示像素在N+1幅时序影像中时间序列的标准差，mA表示像素在N+1幅时序影像中时间序列的均值；采用公式
Figure GDA0002559534490000092
计算每一像素的噪声相位，其中，
Figure GDA0002559534490000093
为噪声相位，
Figure GDA0002559534490000094
为干涉相位，Filter表示均值滤波。通过上述计算后，可分别得到各像素在N+1幅时序影像中对应的平均相干系数，时间序列幅度离差指数和相位误差指数。
上述实施例步骤S2中，采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，将其定义为稳定点后计算形变量的方法包括：
门限值条件包括相干系数T1、时间序列幅度离差指数T2、相位误差指数T3；从时序影像中选取大于相干系数T1，小于时间序列幅度离差指数T2且小于相位误差指数T3的像素，将其定义为稳定点；采用形变量公式
Figure GDA0002559534490000095
Figure GDA0002559534490000096
计算稳定点的形变量，其中，λ表示雷达电磁波波长，
Figure GDA0002559534490000097
为大气相位。
具体实施时，相干系数T1的设置需参考上述实施例中各像素的平均相干系数，时间序列幅度离差指数T2的设置需参考上述实施例中各像素的时间序列幅度离差指数，相位误差指数T3的设置需参考上述实施例中各像素的相位误差指数，通常来说在设置初始门限值条件时，相干系数T1的取值应当远大于各像素平均相干系数的均值，时间序列幅度离差指数T2的取值应当远大于各像素时间序列幅度离差指数的均值，相位误差指数T3的取值应当远大于各像素相位误差指数的均值，示例性地，相干系数T1取值为0.9，时间序列幅度离差指数T2为取值0.1，相位误差指数T3为取值0.2；此外，对于形变量公式
Figure GDA0002559534490000098
中，由于大气相位
Figure GDA0002559534490000099
的取值影响很小，因此在实际计算过程中可将
Figure GDA00025595344900000910
的取值视为零忽略不计。
进一步地，上述实施例步骤S3中，逐步放宽阈值条件，对应从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后计算形变量的方法包括：
阈值条件的放宽方法为将阈值条件相干系数T1逐步减小和/或时间序列幅度离差指数T2逐步增大；从时序影像中选取大于放宽后相干系数T1、小于放宽后时间序列幅度离差指数T2且小于相位误差指数T3的像素，将其定义为动态点；通过形变量公式计算动态点的形变量。
具体实施时，由于相位误差指数T3对动态点的选择影响不大，因此在放宽阈值条件时，可只对相干系数T1逐步减小和/或时间序列幅度离差指数T2逐步增大，在每次放宽阈值条件后，均会从时序影像中选取大于放宽后相干系数T1、小于放宽后时间序列幅度离差指数T2且小于相位误差指数T3的像素，将其定义为动态点，同时通过上述实施例中的形变量公式计算这些动态点对应的形变量；需要说明的是，阈值条件的每次放宽梯度越小越好，这样得到的最优阈值的精度对应也就会越高。
具体地，上述实施例步骤S4中，判断当前动态点形变量与稳定点形变量是否发生突变的方法包括：
计算当前动态点形变量与稳定点形变量的差值Δt；将差值Δt与预设的突变阈值做大小比较，当差值Δt大于预设的突变阈值时认为发生突变，当差值Δt小于或等于预设的突变阈值时认为未发生突变。
上述实施例步骤S5中，依次对相邻时序影像中的各像素共轭相乘得到N幅干涉图，滤波处理后从中提取符合最优阈值的PS点相位差信息的方法包括：
定义相邻的时序影像M和时序影像S中像素的读写格式，
Figure GDA0002559534490000101
Figure GDA0002559534490000102
且S＝M+1，S≤N+1；对应像素采用共轭相乘公式
Figure GDA0002559534490000103
计算第M个干涉图，其中，干涉相位
Figure GDA0002559534490000104
重复上述步骤依次对各相邻时序影像中的像素进行共轭相乘计算，对应得到N幅干涉图；采用均值滤波法分别对N幅干涉图进行滤波处理，从中筛选出符合最优阈值条件的全部PS点；依次将每个PS点对应的N个干涉相位
Figure GDA0002559534490000105
累加，得到各PS点的相位差信息，然后再基于各PS点的相位差信息利用形变量公式
Figure GDA0002559534490000111
计算各PS点的形变量，最终得到形变反演结果。
具体实施时，如图2所示，通过将相邻时序影像中的各对应像素共轭相乘计算后依次会得到第1幅至第N幅干涉图；另外，采用均值滤波法对干涉图进行滤波处理的方法包括：将干涉图中的全部像素按照a*a的像素模板模块化划分；分别计算各模块像素中的像素坐标均值，并将模块的像素坐标均值视为新像素的坐标实现滤波处理，示例性的a的取值为3、5、7等奇数。
可以理解的是，根据微变感知预警雷达特点，若在不良天气采集数据则需进行大气去除以提高形变反演精度，例如，在PS点筛选后选用经典方法进行大气去除即可。
实施例二
请参阅图1和图3，本实施例提供一种微变感知预警雷达的多阈值优化形变反演系统，包括：
处理单元1，用于基于微变感知预警雷达依次获取N+1幅时序影像，计算时序影像中各像素的平均相干系数、时间序列幅度离差指数以及相位误差指数；
严苛筛选单元2，用于采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，将其定义为稳定点后计算形变量；
阈值宽限单元3，用于逐步放宽阈值条件，对应从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后计算形变量；
判断单元4，用于判断当前动态点形变量与稳定点形变量是否发生突变，若没发生突变则继续启动阈值宽限单元继续阈值宽限单元，若发生突变则将当前阈值条件输出为最优阈值；
滤波单元5，用于依次对相邻时序影像中的各像素共轭相乘得到N幅干涉图，滤波处理后从中提取符合最优阈值的PS点相位差信息；
结果输出单元6，用于根据各PS点的相位差信息计算对应的形变量，得到形变反演结果。
具体地，处理单元1包括：
平均相干系数计算模块，用于获取时序影像中的各像素坐标，通过公式
Figure GDA0002559534490000121
依次计算第x个像素在相邻的时序影像M和时序影像S的相干系数，并对获得的N个γx值求均值得到第x个像素的平均相干系数，其中，γx表示第x个像素的相干系数，M(i,j)和S(i,j)表示相邻时序影像中第x个像素的坐标，m和n为滑动窗口大小；
时间序列标准差计算模块，用于采用公式
Figure GDA0002559534490000122
计算每一像素的时间序列幅度离差指数，其中，δA表示像素在N+1幅时序影像中时间序列的标准差，mA表示像素在N+1幅时序影像中时间序列的均值；
相位误差计算模块，用于采用公式
Figure GDA0002559534490000123
计算每一像素的噪声相位，其中，
Figure GDA0002559534490000124
为噪声相位，
Figure GDA0002559534490000125
为干涉相位，Filter表示均值滤波。
具体地，阈值宽限单元3包括：
阈值条件设置模块，设置的阈值条件放宽方法为将阈值条件相干系数T1逐步减小和/或时间序列幅度离差指数T2逐步增大；
动态点筛选模块，用于从时序影像中选取大于放宽后相干系数T1、小于放宽后时间序列幅度离差指数T2且小于相位误差指数T3的像素，将其定义为动态点；
动态点形变量计算模块，用于通过形变量公式计算所述动态点的形变量。
与现有技术相比，本发明实施例提供的微变感知预警雷达的多阈值优化形变反演系统的有益效果与上述实施例一提供的微变感知预警雷达的多阈值优化形变反演方法的有益效果相同，在此不做赘述。
本领域普通技术人员可以理解，实现上述发明方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，上述程序可以存储于计算机可读取存储介质中，该程序在执行时，包括上述实施例方法的各步骤，而所述的存储介质可以是：ROM/RAM、磁碟、光盘、存储卡等。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (7)

1.一种微变感知预警雷达的多阈值优化形变反演方法，其特征在于，包括：
S1，基于微变感知预警雷达依次获取N+1幅时序影像，计算时序影像中各像素的平均相干系数、时间序列幅度离差指数以及相位误差指数；
S2，采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，将其定义为稳定点后计算形变量；
S3，逐步放宽阈值条件，对应从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后计算形变量；
S4，判断当前动态点形变量与稳定点形变量是否发生突变，若没发生突变则执行步骤S3，若发生突变则将当前阈值条件输出为最优阈值，执行步骤S5；
S5，依次对相邻时序影像中的各像素共轭相乘得到N幅干涉图，滤波处理后从中提取符合最优阈值的PS点相位差信息；
S6，根据各PS点的相位差信息计算对应的形变量，得到形变反演结果；
其中，所述步骤S2中，采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，将其定义为稳定点后计算形变量的方法包括：
所述门限值条件包括相干系数T1、时间序列幅度离差指数T2、相位误差指数T3；
从时序影像中选取大于相干系数T1，小于时间序列幅度离差指数T2且小于相位误差指数T3的像素，将其定义为稳定点；
采用形变量公式
Figure FDA0002559534480000011
计算所述稳定点的形变量，其中，λ表示雷达电磁波波长，
Figure FDA0002559534480000012
为大气相位；
所述步骤S3中，逐步放宽阈值条件，对应从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后计算形变量的方法包括：
阈值条件的放宽方法为将阈值条件相干系数T1逐步减小和/或时间序列幅度离差指数T2逐步增大；
从时序影像中选取大于放宽后相干系数T1、小于放宽后时间序列幅度离差指数T2且小于相位误差指数T3的像素，将其定义为动态点；
通过形变量公式计算所述动态点的形变量。
2.根据权利要求1所述的方法，其特征在于，所述步骤S1中，计算时序影像中各像素的平均相干系数、时间序列幅度离差指数以及相位误差指数的方法包括：
获取时序影像中的各像素坐标，通过公式
Figure FDA0002559534480000021
依次计算第x个像素在相邻的时序影像M和时序影像S的相干系数，并对获得的N个γx求均值得到第x个像素的平均相干系数，其中，γx表示第x个像素的相干系数，M(i,j)和S(i,j)表示相邻时序影像中第x个像素的坐标，m和n为滑动窗口大小；
采用公式
Figure FDA0002559534480000022
计算每一像素的时间序列幅度离差指数，其中，δA表示所述像素在N+1幅时序影像中时间序列的标准差，mA表示所述像素在N+1幅时序影像中时间序列的均值；
采用公式
Figure FDA0002559534480000023
计算每一像素的噪声相位，其中，
Figure FDA0002559534480000024
为噪声相位，
Figure FDA0002559534480000025
为干涉相位，Filter表示均值滤波。
3.根据权利要求1所述的方法，其特征在于，所述步骤S4中，判断当前动态点形变量与稳定点形变量是否发生突变的方法包括：
计算当前动态点形变量与稳定点形变量的差值Δt；
将差值Δt与预设的突变阈值做大小比较，当差值Δt大于预设的突变阈值时认为发生突变，当差值Δt小于或等于预设的突变阈值时认为未发生突变。
4.根据权利要求1所述的方法，其特征在于，所述步骤S5中，依次对相邻时序影像中的各像素共轭相乘得到N幅干涉图，滤波处理后从中提取符合最优阈值的PS点相位差信息的方法包括：
定义相邻的时序影像M和时序影像S中像素的读写格式，
Figure FDA0002559534480000031
Figure FDA0002559534480000032
且S＝M+1，S≤N+1；
对应像素采用共轭相乘公式
Figure FDA0002559534480000033
计算第M个干涉图，其中，干涉相位
Figure FDA0002559534480000034
重复上述步骤依次对各相邻时序影像中的像素进行共轭相乘计算，对应得到N幅干涉图；
采用均值滤波法分别对N幅干涉图进行滤波处理，从中筛选出符合最优阈值条件的全部PS点；
依次将每个PS点对应的N个干涉相位
Figure FDA0002559534480000035
累加，得到各PS点的相位差信息。
5.根据权利要求4所述的方法，其特征在于，采用均值滤波法对干涉图进行滤波处理的方法包括：
将干涉图中的全部像素按照a*a的像素模板模块化划分；
分别计算各模块像素中的像素坐标均值，并将模块的像素坐标均值视为新像素的坐标实现滤波处理。
6.一种微变感知预警雷达的多阈值优化形变反演系统，其特征在于，包括：
处理单元，用于基于微变感知预警雷达依次获取N+1幅时序影像，计算时序影像中各像素的平均相干系数、时间序列幅度离差指数以及相位误差指数；
严苛筛选单元，用于采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，将其定义为稳定点后计算形变量；
阈值宽限单元，用于逐步放宽阈值条件，对应从时序影像中筛选出符合放宽门限值条件的像素，将其定义为动态点后计算形变量；
判断单元，用于判断当前动态点形变量与稳定点形变量是否发生突变，若没发生突变则继续启动阈值宽限单元继续阈值宽限单元，若发生突变则将当前阈值条件输出为最优阈值；
滤波单元，用于依次对相邻时序影像中的各像素共轭相乘得到N幅干涉图，滤波处理后从中提取符合最优阈值的PS点相位差信息；
结果输出单元，用于根据各PS点的相位差信息计算对应的形变量，得到形变反演结果；
其中，所述严苛筛选单元采用严苛阈值筛选从时序影像中筛选出符合门限值条件的像素，将其定义为稳定点后计算形变量的方法包括：
所述门限值条件包括相干系数T1、时间序列幅度离差指数T2、相位误差指数T3；
从时序影像中选取大于相干系数T1，小于时间序列幅度离差指数T2且小于相位误差指数T3的像素，将其定义为稳定点；
采用形变量公式
Figure FDA0002559534480000041
计算所述稳定点的形变量，其中，λ表示雷达电磁波波长，
Figure FDA0002559534480000043
为大气相位；
所述阈值宽限单元包括：
阈值条件设置模块，设置的阈值条件放宽方法为将阈值条件相干系数T1逐步减小和/或时间序列幅度离差指数T2逐步增大；
动态点筛选模块，用于从时序影像中选取大于放宽后相干系数T1、小于放宽后时间序列幅度离差指数T2且小于相位误差指数T3的像素，将其定义为动态点；
动态点形变量计算模块，用于通过形变量公式计算所述动态点的形变量。
7.根据权利要求6所述的系统，其特征在于，所述处理单元包括：
平均相干系数计算模块，用于获取时序影像中的各像素坐标，通过公式
Figure FDA0002559534480000042
依次计算第x个像素在相邻的时序影像M和时序影像S的相干系数，并对获得的N个γx求均值得到第x个像素的平均相干系数，其中，γx表示第x个像素的相干系数，M(i,j)和S(i,j)表示相邻时序影像中第x个像素的坐标，m和n为滑动窗口大小；
时间序列标准差计算模块，用于采用公式
Figure FDA0002559534480000051
计算每一像素的时间序列幅度离差指数，其中，δA表示所述像素在N+1幅时序影像中时间序列的标准差，mA表示所述像素在N+1幅时序影像中时间序列的均值；
相位误差计算模块，用于采用公式
Figure FDA0002559534480000052
计算每一像素的噪声相位，其中，
Figure FDA0002559534480000053
为噪声相位，
Figure FDA0002559534480000054
为干涉相位，Filter表示均值滤波。
CN201910011969.1A 2019-01-07 2019-01-07 一种微变感知预警雷达的多阈值优化形变反演方法及系统 Active CN109613531B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910011969.1A CN109613531B (zh) 2019-01-07 2019-01-07 一种微变感知预警雷达的多阈值优化形变反演方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910011969.1A CN109613531B (zh) 2019-01-07 2019-01-07 一种微变感知预警雷达的多阈值优化形变反演方法及系统

Publications (2)

Publication Number Publication Date
CN109613531A CN109613531A (zh) 2019-04-12
CN109613531B true CN109613531B (zh) 2020-10-02

Family

ID=66016745

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910011969.1A Active CN109613531B (zh) 2019-01-07 2019-01-07 一种微变感知预警雷达的多阈值优化形变反演方法及系统

Country Status (1)

Country Link
CN (1) CN109613531B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2021056008A (ja) * 2019-09-27 2021-04-08 株式会社パスコ 地すべり領域検出装置及びプログラム
CN111950361B (zh) * 2020-07-07 2022-09-20 内蒙古农业大学 基于单时序ndvi的甜菜识别方法、系统、设备、介质及终端
CN113537240B (zh) * 2021-07-09 2023-09-05 北方工业大学 一种基于雷达序列图像的形变区智能提取方法及系统

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003500658A (ja) * 1999-05-25 2003-01-07 ポリテクニコ ディ ミラノ 市街地領域及び地滑り地帯の運動に関するレーダー測定のための手順
CN104111456A (zh) * 2014-07-23 2014-10-22 中国国土资源航空物探遥感中心 一种高速铁路沿线地表形变高分辨率InSAR监测方法
CN105678716A (zh) * 2016-02-25 2016-06-15 内蒙古工业大学 一种地基sar大气干扰相位校正方法及装置
CN105866777A (zh) * 2016-03-29 2016-08-17 北京理工大学 多角度多时段导航卫星双基地PS-InSAR三维形变反演方法
CN106932777A (zh) * 2017-04-12 2017-07-07 西南石油大学 基于温度基线的合成孔径雷达干涉对优化选取方法
CN107516318A (zh) * 2017-08-25 2017-12-26 四川长虹电器股份有限公司 基于模式搜索算法与萤火虫算法的图像多阈值分割方法
CN108663017A (zh) * 2018-08-13 2018-10-16 伟志股份公司 一种监测城市地铁沿线地表沉降的方法
CN108802727A (zh) * 2018-04-13 2018-11-13 长沙理工大学 一种顾及流变参数的时序InSAR公路形变监测模型及解算方法
CN108983233A (zh) * 2018-06-13 2018-12-11 四川大学 GB-InSAR数据处理中的PS点组合选取方法
CN109031300A (zh) * 2018-09-03 2018-12-18 中科卫星应用德清研究院 合成孔径雷达监测危岩体变形方法及系统
CN109116354A (zh) * 2018-09-03 2019-01-01 北京市测绘设计研究院 一种基于信杂比加权的振幅离差ps点选取方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003500658A (ja) * 1999-05-25 2003-01-07 ポリテクニコ ディ ミラノ 市街地領域及び地滑り地帯の運動に関するレーダー測定のための手順
CN104111456A (zh) * 2014-07-23 2014-10-22 中国国土资源航空物探遥感中心 一种高速铁路沿线地表形变高分辨率InSAR监测方法
CN105678716A (zh) * 2016-02-25 2016-06-15 内蒙古工业大学 一种地基sar大气干扰相位校正方法及装置
CN105866777A (zh) * 2016-03-29 2016-08-17 北京理工大学 多角度多时段导航卫星双基地PS-InSAR三维形变反演方法
CN106932777A (zh) * 2017-04-12 2017-07-07 西南石油大学 基于温度基线的合成孔径雷达干涉对优化选取方法
CN107516318A (zh) * 2017-08-25 2017-12-26 四川长虹电器股份有限公司 基于模式搜索算法与萤火虫算法的图像多阈值分割方法
CN108802727A (zh) * 2018-04-13 2018-11-13 长沙理工大学 一种顾及流变参数的时序InSAR公路形变监测模型及解算方法
CN108983233A (zh) * 2018-06-13 2018-12-11 四川大学 GB-InSAR数据处理中的PS点组合选取方法
CN108663017A (zh) * 2018-08-13 2018-10-16 伟志股份公司 一种监测城市地铁沿线地表沉降的方法
CN109031300A (zh) * 2018-09-03 2018-12-18 中科卫星应用德清研究院 合成孔径雷达监测危岩体变形方法及系统
CN109116354A (zh) * 2018-09-03 2019-01-01 北京市测绘设计研究院 一种基于信杂比加权的振幅离差ps点选取方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
InSAR若干关键算法及其在地表沉降监测中的应用研究;范洪冬;《中国博士学位论文全文数据库 基础科学辑》;20110415;全文 *
Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry;Ferretti A等;《IEEE Transactions on Geoscience ＆ Ｒemote Sensing》;20001231;第38卷(第5期);全文 *
Permanent scatterers in SAR interferometry;Ferretti A等;《IEEE Transactions on Geoscience ＆ Ｒemote Sensing》;20011231;第39卷(第1期);全文 *
地基SAR形变监测误差分析与实验;曲世勃等;《电子与信息学报》;20110131;第33卷(第1期);全文 *
基于动态PS的地基合成孔径雷达高精度形变测量技术研究;朱茂;《中国博士学位论文全文数据库 信息科技辑》;20160715;全文 *
小数据集PS-DInSAR的PS点探测方法;韩洁等;《信号处理》;20150630;第31卷(第6期);全文 *

Also Published As

Publication number Publication date
CN109613531A (zh) 2019-04-12

Similar Documents

Publication Publication Date Title
CN109613531B (zh) 一种微变感知预警雷达的多阈值优化形变反演方法及系统
CN101727662B (zh) Sar图像非局部均值去斑方法
CN103645476B (zh) 一种合成孔径雷达差分干涉图序列的时空同质滤波方法
CN111596366B (zh) 一种基于地震信号优化处理的波阻抗反演方法
Tang et al. Downscaling remotely sensed imagery using area-to-point cokriging and multiple-point geostatistical simulation
CN103902802B (zh) 一种顾及空间信息的植被指数时间序列数据重建方法
CN101980293A (zh) 一种基于刃边图像的高光谱遥感系统mtf检测方法
CN104200471A (zh) 基于自适应权值图像融合的sar图像变化检测方法
CN115512231A (zh) 适用于国土空间生态修复的遥感解译方法
CN115797335B (zh) 用于桥梁振动测量的欧拉运动放大效果评估及优化方法
CN109146797A (zh) 一种基于Lp伪范数与交叠组稀疏的脉冲噪声古籍图像修复方法
CN103208101A (zh) 一种基于局部信噪比的干涉图滤波方法
CN112347992B (zh) 一种荒漠地区时序agb遥感估算方法
CN112051572A (zh) 一种融合多源sar数据三维地表形变监测方法
CN110033434A (zh) 一种基于纹理显著性的表面缺陷检测方法
Warner et al. Pine Island Glacier (Antarctica) velocities from Landsat7 images between 2001 and 2011: FFT-based image correlation for images with data gaps
Zhang et al. A study on coastline extraction and its trend based on remote sensing image data mining
CN111426616A (zh) 碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质
CN108230365A (zh) 基于多源差异图像内容融合的sar图像变化检测方法
CN107944497A (zh) 基于主成分分析的图像块相似性度量方法
CN110532969B (zh) 基于多尺度图像分割的斜坡单元划分方法
CN103955943A (zh) 基于融合变化检测算子与尺度驱动的无监督变化检测方法
Zhang et al. Segmented noise reduction based on Brillouin-spectrum-partition in Brillouin optical time domain sensors
Silvetti et al. Quadratic self-correlation: An improved method for computing local fractal dimension in remote sensing imagery
CN114779249A (zh) 一种gb-sar干涉影像滤波方法及系统

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