CN114779249A - 一种gb-sar干涉影像滤波方法及系统 - Google Patents

一种gb-sar干涉影像滤波方法及系统 Download PDF

Info

Publication number
CN114779249A
CN114779249A CN202210407680.3A CN202210407680A CN114779249A CN 114779249 A CN114779249 A CN 114779249A CN 202210407680 A CN202210407680 A CN 202210407680A CN 114779249 A CN114779249 A CN 114779249A
Authority
CN
China
Prior art keywords
sliding window
noise
phase
points
standard deviation
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
Application number
CN202210407680.3A
Other languages
English (en)
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.)
Northeastern University China
Original Assignee
Northeastern University China
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 Northeastern University China filed Critical Northeastern University China
Priority to CN202210407680.3A priority Critical patent/CN114779249A/zh
Publication of CN114779249A publication Critical patent/CN114779249A/zh
Pending legal-status Critical Current

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
    • 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

Landscapes

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

Abstract

本发明涉及一种GB‑SAR干涉影像滤波方法及系统，包括：获取GB‑SAR监测的露天矿边坡的时间序列干涉影像并进行粗滤波和形变分离，得到粗滤波后的影像中的稳定点；建立稳定点相位数据集，并基于最小二乘算法确定稳定点相位数据集的高斯函数半宽度信息，进而确定标准差最优补偿系数和标准差限差；构建噪声识别规则；基于当前滑动窗口结合八邻域规则和噪声识别规则得到真实噪声点；基于反距离和相位偏差双加权算法对真实噪声点进行相位重构；移动滑动窗口直至遍历所有的干涉影像数据。不仅能够精确地识别GB‑SAR干涉影像中的噪声点，还能够准确地重构噪声点的干涉相位。

Description

一种GB-SAR干涉影像滤波方法及系统
技术领域
本发明涉及边坡工程变形监测技术领域，特别是涉及一种基于八邻域离群值判别的GB-SAR干涉影像滤波方法及系统。
背景技术
露天矿随开采度的不断加深，边坡长度不断加长，其稳定性将会越来越差，近年来国内外诸多露天矿边坡滑坡事故时有发生，其给矿山安全生产带来了极大威胁。为解决露天矿边坡安全问题，国内外学者多年来一直致力于边坡失稳分析及预测预报工作的研究。目前，GB-SAR边坡监测技术近期发展而来的无接触式地面新型监测设备，相比传统D-InSAR、PS-InSAR、GPS静态测量等监测技术，它拥有更高的分辨率、监测精度和影像采集频率，能够实现全天时、短周期、连续地高精度监测。然而其监测精度受噪声、大气延迟误差、重复轨道误差等因素的影响较大，目前众多学者仅针对大气延迟误差补偿问题进行了大量的研究，提出了多种方法可实现对该误差的有效补偿。而重复轨道误差经过理论推导在方位向具有较好的回归特性，因此该误差可以利用数学模型准确修正。但针对GB-SAR干涉影像中噪声的去除研究较少，缺乏关于GB-SAR干涉影像噪声的特征分析，且噪声点具有误差传递性，严重影响形变解算的精度，此外GB-SAR位移监测数据的准确性是边坡失稳分析及预警的关键，因而GB-SAR噪声的精准滤波问题亟待解决。理论上GB-SAR干涉影像的去噪问题可以借鉴星载InSAR干涉影像滤波的经典goldstein自适应滤波和Baran滤波等方法，但考虑到上述方法针对具有干涉条纹的干涉影像效果明显，而GB-SAR成像原理具有距离短、周期小、无空间基线影响等特征，在短周期连续监测中干涉影像一般无明显干涉条纹，因此上述方法的适用性受到了限制。通常情况下针对GB-SAR干涉影像噪声可通过中值滤波和小波去噪等方法进行剔除，上述方法虽然能取得一定的去噪效果，但也存在着一定的缺陷与不足，其一，中值滤波和小波去噪均会引起形变区域边缘的相位模糊问题，且在滤波过程中干涉相位有效信息丢失较多，最终不同程度地削减形变区域位移量的峰值，严重时将会引起GB-SAR预警失效的问题。其二，两种方法均对干涉影像执行全局滤波，忽略了非边坡点的冗余相位对噪声点判断和估值的影响，且无法解决形变信息的损失问题。因此，如何实现GB-SAR影像中相位噪声的准确去除和重构是本领域技术人员期望克服的。
发明内容
本发明的目的是提供一种GB-SAR干涉影像滤波方法，引入格拉布斯准则和八邻域规则以及基于反距离和相位偏差双加权算法，不仅能够精确地识别GB-SAR干涉影像中的噪声点，还能够准确地重构噪声点的干涉相位，有效的实现了干涉影像滤波。
为实现上述目的，本发明提供了如下方案：
一种GB-SAR干涉影像滤波方法，包括：
获取GB-SAR监测的露天矿边坡的时间序列干涉影像；
对所述时间序列干涉影像进行粗滤波，并对粗滤波后的影像进行形变分离，得到所述粗滤波后的影像中的稳定点；
对分离出的各所述稳定点建立稳定点相位数据集，根据所述稳定点相位数据集应用最小二乘算法确定所述稳定点相位数据集的高斯函数半宽度信息；
根据所述高斯函数半宽度信息确定标准差最优补偿系数和标准差限差；
基于所述标准差最优补偿系数、所述标准差限差、永久散射体点的相位均值、永久散射体点的相位标准差以及格拉布斯准则构建噪声识别规则；
确定所述时间序列干涉影像的滑动窗口，基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点；
基于反距离和相位偏差双加权算法对当前所述滑动窗口中的所述真实噪声点进行相位重构，实现GB-SAR干涉影像滤波；
移动所述滑动窗口至下一所述滑动窗口，令下一所述滑动窗口为当前所述滑动窗口，并返回步骤“基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点”，直至遍历所述时间序列干涉影像的全部影像数据。
一种GB-SAR干涉影像滤波系统，包括：
数据获取模块，用于获取GB-SAR监测的露天矿边坡的时间序列干涉影像；
粗滤波及形变分离模块，用于对所述时间序列干涉影像进行粗滤波，并对粗滤波后的影像进行形变分离，得到所述粗滤波后的影像中的稳定点；
第一计算模块，用于对分离出的各所述稳定点建立稳定点相位数据集，根据所述稳定点相位数据集应用最小二乘算法确定所述稳定点相位数据集的高斯函数半宽度信息；
第二计算模块，用于根据所述高斯函数半宽度信息确定标准差最优补偿系数和标准差限差；
噪声识别规则构建模块，用于基于所述标准差最优补偿系数、所述标准差限差、永久散射体点的相位均值、永久散射体点的相位标准差以及格拉布斯准则构建噪声识别规则；
噪声识别模块，用于确定所述时间序列干涉影像的滑动窗口，基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点；
噪声重构模块，用于基于反距离和相位偏差双加权算法对当前所述滑动窗口中的所述真实噪声点进行相位重构。
滑动窗口移动模块，用于移动所述滑动窗口至下一所述滑动窗口，令下一所述滑动窗口为当前所述滑动窗口，并返回执行噪声识别模块中的“基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点”，直至遍历所述时间序列干涉影像的全部影像数据。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明涉及一种GB-SAR干涉影像滤波方法及系统，包括：获取GB-SAR监测的露天矿边坡的时间序列干涉影像并进行粗滤波和形变分离，得到粗滤波后的影像中的稳定点；建立稳定点相位数据集，并基于最小二乘算法确定稳定点相位数据集的高斯函数半宽度信息，进而确定标准差最优补偿系数和标准差限差；基于标准差最优补偿系数、标准差限差、永久散射体点的相位均值、永久散射体点的相位标准差以及格拉布斯准则构建噪声识别规则；确定时间序列干涉影像的滑动窗口，基于当前滑动窗口结合八邻域规则和所述噪声识别规则对当前滑动窗口中的永久散射体点进行离群值判别，得到当前滑动窗口中的真实噪声点；基于反距离和相位偏差双加权算法对当前滑动窗口中的所述真实噪声点进行相位重构，实现GB-SAR干涉影像滤波；移动滑动窗口直至遍历所有的干涉影像数据。通过构建噪声识别规则以及应用八邻域规则精确地识别GB-SAR干涉影像中的噪声点，同时还应用基于反距离和相位偏差双加权算法进行噪声点的干涉相位精确重构，大大提高了干涉影像滤波的效果，提高GB-SAR的边坡监测精度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例1提供的一种GB-SAR干涉影像滤波方法流程图；
图2为本发明实施例1提供的八邻域去噪示意图；
图3为本发明实施例1提供的GB-SAR干涉影像去噪前后对比图；
图4为本发明实施例2提供的一种GB-SAR干涉影像滤波系统框图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种GB-SAR干涉影像滤波方法及系统，引入格拉布斯准则和八邻域规则以及基于反距离和相位偏差双加权算法，不仅能够精确地识别GB-SAR干涉影像中的噪声点，而且还能够准确地重构噪声点的干涉相位，有效的实现了干涉影像滤波。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
实施例1
如图1所示，本实施例提供一种GB-SAR干涉影像滤波方法，包括：
步骤S1：获取GB-SAR监测的露天矿边坡的时间序列干涉影像；
在干涉计算过程中只需将两个回波矩阵共轭相乘，即为两幅影像的干涉图，设回波矩阵中任意信号为I1，对应相邻影像中相同位置的回波信号为I2，则干涉信号I12为：
Figure BDA0003602435970000051
步骤S2：对所述时间序列干涉影像进行粗滤波，并对粗滤波后的影像进行形变分离，得到所述粗滤波后的影像中的稳定点。
GB-SAR影像粗滤波过程旨在去除非边坡回波信号对有效边坡回波信号的干扰，提高影像的质量。理论上边坡回波信号较非边坡回波信号强度较大，且信号相干性较好，后向散射强度在时间序列上较为稳定，因此定义了时间序列平均幅度、时间序列相干系数和幅度离差三个指标对GB-SAR影像进行粗滤波，准确提取边坡有效回波信号。
步骤S2具体包括：
根据所述时间序列干涉影像确定影像的时间序列平均幅度均值、时间序列幅度离差、时间序列相干系数、时间序列极相干系数；
具体的，计算所述时间序列幅度均值的公式为：
Figure BDA0003602435970000052
式中m，n为干涉影像距离向和方位向的尺寸，
Figure BDA0003602435970000053
为像元(i，j)的时间序列后向散射强度均值；
计算所述时间序列幅度离差的公式为：
Figure BDA0003602435970000054
式中，dt为像元(i，j)在t时刻的后向散射强度；
Figure BDA0003602435970000055
为像元(i，j)的时间序列后向散射强度均值；k为时间序列干涉影像数目。
计算时间序列相干系数的公式为：
Figure BDA0003602435970000061
式中：k为时间序列影像数目，It(i，j)和It+1(i，j)为相邻时刻两幅干涉影像中像元(i，j)的散射复数信号，m，n为影像距离向和方位向的尺寸；
计算所述极相干系数的公式为：
Figure BDA0003602435970000062
式中：γ为选取的影像数目，k为总影像数目，It，1和It，k+1-γ为干涉影像像元复信号。
根据所述时间序列平均幅度均值、所述时间序列幅度离差和所述时间序列相干系数对所述时间序列干涉影像进行粗滤波，并根据粗滤波后的影像结合所述时间序列极相干系数分离出所述粗滤波后的影像中的稳定点，具体包括：
(1)令所述时间序列平均幅度均值为时间序列幅度均值阈值，选择时间序列内平均幅度小于所述时间序列幅度阈值的像元作为一次粗滤波干涉影像；
(2)设置时间序列幅度离差阈值，并基于所述一次粗滤波干涉影像，选择所述时间序列内幅度离差小于所述时间序列幅度离差阈值的像元作为二次粗滤波干涉影像。时间序列幅度离差阈值可以设为0.5，也可以根据需求设定其他数值。
(3)设置时间序列相干系数阈值，并基于二次粗滤波干涉影像选择所述时间序列内平均相干系数大于所述时间序列相干系数阈值的像元作为最终粗滤波干涉影像，所述最终粗滤波干涉影像即为所述粗滤波后的影像。时间序列相干系数阈值可以设为0.85，也可以根据需求设定其他数值。
(4)设置极相干系数阈值；基于所述粗滤波后的影像，提取所述极相干系数大于所述极相干系数阈值的永久散射点(PS点)作为稳定点。极相干系数阈值可以设为0.8，也可以根据需求设定其他数值。
步骤S3：对分离出的各所述稳定点建立稳定点相位数据集，根据所述稳定点相位数据集应用最小二乘算法确定所述稳定点相位数据集的高斯函数半宽度信息；
在曲线峰值特征方面，基于样本方差和均值的正态分布与数据样本概率密度分布曲线差异较大，其原因是由于稳定区域PS点干涉相位得噪声绝对值多集中在区间[1,2]，无噪声相位信息多集中在区间[-0.3,0.3]，因而样本方差估算偏大，导致正态分布曲线分布较为平缓。显然该正态分布不是样本的最优概率密度分布曲线。为此需要基于最小二乘算法估计了与样本数据最优拟合的高斯概率分布函数，计算滤波所用的必要参数提升滤波效果：干涉影像稳定点相位数据集的高斯函数峰值、高斯函数估计均值、高斯函数半宽度信息、最优补偿系数和标准差限差。本实施例干涉影像滤波应用了高斯函数半宽度信息、最优补偿系数和标准差限差3个参数。
具体的，步骤S3具体包括：
(1)设定所述稳定点相位数据集为(xf，f(xf))，f＝1,2,3……，M0；xf为第f个稳定点的相位；f(xf)为第f个稳定点的相位的概率；f(xf)为高斯函数，
Figure BDA0003602435970000071
a为高斯函数峰值，μg为高斯函数估计均值，S为高斯函数半宽度信息；
(2)对所述高斯函数f(xf)取对数得到
Figure BDA0003602435970000072
(3)令zf＝lnf(xf),
Figure BDA0003602435970000073
得出所述稳定点相位数据集的矩阵形式
Figure BDA0003602435970000074
(4)将所述稳定点相位数据集的矩阵形式简化为Z＝XB，并基于最小二乘原理，得到矩阵B的广义最小二乘解；
(5)基于所述矩阵B的广义最小二乘解得到所述稳定点相位数据集的高斯函数半宽度信息。基于所述矩阵B的广义最小二乘解也可得出稳定点相位数据集的高斯函数峰值a、高斯函数估计均值μg。
步骤S4：根据所述高斯函数半宽度信息确定标准差最优补偿系数和标准差限差，具体包括：
(1)基于所述高斯函数半宽度信息和基于极大似然估计的稳定点干涉相位样本标准差确定标准差最优补偿系数；
标准差最优补偿系数的表达式为：
Figure BDA0003602435970000081
S为高斯函数的半宽度信息；σ为基于极大似然估计的稳定点干涉相位样本标准差，
Figure BDA0003602435970000082
(2)根据所述标准差最优补偿系数和所述高斯半宽度信息确定标准差限差。
标准差限差的表达式为：
Figure BDA0003602435970000083
步骤S5：基于所述标准差最优补偿系数、所述标准差限差、永久散射体点的相位均值、永久散射体点的相位标准差以及格拉布斯准则构建噪声识别规则；
考虑PS点数目的不确定性，引入格拉布斯离群值判别准则，同时考虑了样本数据的最优高斯拟合方差与样本极大似然估计方差的差异性，需要对样本方差进行最优高斯估计补偿。考虑PS点数目的不确定性，引入格拉布斯离群值判别准则，该准则以正态分布为前提，针对不同数量的数据集具有不同的Grubbs判别临界值，满足PS点数目不确定的要求。若单个PS点相位的残余误差绝对值大于理论阈值，则该点被识别为初始噪声点。同时考虑了样本数据的最优高斯拟合方差与样本极大似然估计方差的差异性，需要对样本方差进行最优高斯估计补偿。则所述噪声识别规则包括：
(1)当滑动窗口中的永久散射体点的个数大于第一预设值时，且所述滑动窗口内的所述永久散射体点的相位标准差小于所述标准差限差时，采用W1公式进行噪声点估计；
W1＝(xps-μ)＞g(s,α)·λ1σs，其中，
Figure BDA0003602435970000084
λ1为标准差最优补偿系数；S为高斯函数的半宽度信息；σ为基于极大似然估计的稳定点干涉相位样本标准差，；xps为永久散射体点相位；μ为滑动窗内永久散射体点相位均值；g为Grubbs临界值，s为永久散射体点数目；α为显著性水平，可取值为0.05；σs为滑动窗口内PS点相位标准差；M0为稳定点个数；
(2)当滑动窗口中的永久散射体点的个数大于第一预设值时，且所述滑动窗口内的所述永久散射体点的相位标准差大于所述标准差限差时，采用W2公式进行噪声点估计；
Figure BDA0003602435970000091
其中，σg为标准差限差；
(3)当滑动窗口中的永久散射体点的个数小于第一预设值时，采用W3公式进行噪声点估计；
W3＝(xps-μ)＞g(s,α)·λ1σg。
除上述情况外，还可以针对干涉影像中边坡边缘的滤波问题，若滑动窗口位于边坡边缘缺少特定方向的邻域，则噪声点识别中仅考虑存在邻域的噪声交集进行判别。
步骤S6：确定所述时间序列干涉影像的滑动窗口，基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点；
步骤S6具体包括：
(1)依据八邻域判别规则，将当前所述滑动窗口的邻域按形变区域分布的空间位置分割为八种类型，得到八种邻域分割区域；
(2)依据所述噪声识别规则分别对当前所述滑动窗口内的永久散射体点和每一所述邻域分割区域内的永久散射体点进行噪声识别，得到当前所述滑动窗口内的初始噪声点集和每一所述邻域分割区域内的初始噪声点集；
(3)对所述滑动窗口内的初始噪声点集和每一所述邻域分割区域内的初始噪声点集取交集，得到当前所述滑动窗口内的真实噪声点集。
在边坡稳定的前提下，基于上述离群值检测方法可以准确识别滑动窗内相位误差过大的噪声点，然而一旦边坡部分区域发生位移，滑动窗口区域与边坡稳定区域和位移区域的临界线相交，则会引起少量位移PS点被误识别噪声点的问题，最终导致位移区域边缘形变信息丢失，严重降低GB-SAR的监测精度。为此发明提出了基于八邻域噪声点探测的综合判别方法，该方法在识别滑动窗口内噪声点的基础上，考虑了噪声点在滑动窗口空间邻域内具有同样的性质，理论上在滑动窗口被识别的噪声点在其邻域范围内也一定被识别。假设边坡少量形变点位于滑动窗口内部，虽然利用上述方法在识别噪声点的同时会将少量形变点错误识别为噪声点，但是必有特定邻域分布在形变区域内，且该邻域内形变点不被误识别为噪声点，而噪声点则会被准确识别。因此在识别出滑动窗口A中初始噪声点集NA的基础上，将滑动窗口的邻域按形变区域可能分布的空间位置分割为八种类型，如图2所示，将分割区域所在邻域称之为A1，A2，A3，A4，A5，A6，A7，A8。再基于Grubbs准则分别对八个邻域进行离群值判别，识别出各个邻域内的噪声点集N1，N2，N3，N4，N5，N6，N7，N8，最终将各个邻域与滑动窗口公共区域的噪声点集交集识别为真实噪声点集。
步骤S7：基于反距离和相位偏差双加权算法对所述真实噪声点进行相位重构，实现GB-SAR干涉影像滤波。
识别出噪声点位后，需要对其相位进行重构以提升GB-SAR监测精度。为此发明提出了基于反距离和相位偏差双加权算法对噪声点干涉相位进行相位重构。则步骤S7具体包括：
(1)基于相位偏差最小补偿原则，对当前所述滑动窗口中的所述真实噪声点选取预设个数的参考点；
(2)计算对每一所述真实噪声点与每一所述参考点之间的距离参数和相位偏差；
(3)对每一所述真实噪声点，基于对应的归一化后的所述距离参数和归一化后的所述相位偏差计算权重值；
权值公式为
Figure BDA0003602435970000101
噪声点数目为sn；参考点数目为E。若s-sn<10，E＝3；当s-sn≥10时，E＝[(s-sn)/5]+1。
(4)将所述权重值带入到噪声补偿公式中得到每一所述真实噪声点的估计相位。
噪声补偿公式为：
Figure BDA0003602435970000111
其中，φnoise为噪声点相位的估计值，φe为第e个参考点的相位值。
步骤S8：移动所述滑动窗口至下一所述滑动窗口，令下一所述滑动窗口为当前所述滑动窗口，并返回步骤S6中的“基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点”直至遍历所述时间序列干涉影像的全部影像数据。
移动前和移动后的两个滑动窗口是不存在重叠区域的，移动前和移动后的两个滑动窗口首尾顺次相接。
下面以唐山市马兰庄露天铁矿采场GB-SAR边坡监测为例，该边坡为岩质边坡。首先对露天矿GB-SAR边坡监测数据进行了干涉和粗滤波，并基于本发明八邻域离群值判别去噪方法进行了精滤波，最终解算了边坡监测结果，以此来验证本发明的效果。
1)GB-SAR监测数据采集
首先在鞍山市大孤山露天铁矿采场的稳定区域安置GB-SAR，调整好其相对被监测边坡的最佳观测姿态。利用GB-SAR采集为期14天的时间序列原始影像数据，并对影像数据进行聚焦处理生成时间序列单视复影像。
2)计算时间序列幅度均值
通过对GB-SAR单视复影像中的复数据求取绝对值，以获取边坡后向散射信号的幅值，并在时间序列上将对应像素的幅值利用文中计算公式计算时间序列幅度均值影像。
3)计算幅度离差
幅度离差指数D是像元在时间序列中幅度标准差与幅度均值的比值，其用于评价一个像元在时间序列上后向散射强度的稳定性。利用文中计算公式计算幅度离差。
4)计算相干系数
根据用于干涉的两幅单视复影像可以计算相干系数，计算相干系数时滑动窗口设置为4×12(方位向×距离向)，利用文中计算公式计算时间序列相干系数。参数计算时设置的滑动窗口根据参数计算要求设定，与噪声识别时的滑动窗口的设置是完全独立的过程。
5)计算极相干系数
极相干系数可衡量稳定点和形变点在长时间序列中的相干性，根据首末两幅或多幅单视复影像利用文中计算公式可以计算极相干系数，计算相干系数时滑动窗口设置为3×9(方位向×距离向)。
6)粗滤波和分离形变点
上述四种指标计算完成后，执行以下步骤进行粗滤波和形变分离：
步骤1、设置时间序列幅度均值阈值M2＝M，并选择在时间序列内平均幅度Mi，j>M2的像元作为一次粗滤波结果。
步骤2、设置时间序列幅度离差阈值D2＝0.5，并基于一次粗滤波结果的基础上选择时间序列幅度离差Di，j<D2的像元作为二次粗滤波结果。
步骤3、设置时间序列相干系数阈值L2＝0.85，并基于二次粗滤波结果的基础上选择时间序列内平均相干系数Li，j>L2的像元作为最终粗滤波的结果。
步骤4、在粗滤波后影像的基础上，设置极相干系数阈值Sl＝0.8，提取S>S1的PS点作为稳定点，S<Sl的PS点作为形变点实现两者分离。
7)精滤波参数的计算
计算干涉影像稳定点相位数据集的最优补偿系数和标准差限差等参数，用于精滤波的最优分布估计和样本的异常方差限制。
8)基于八邻域离群值判别的精滤波
设置精滑动窗口为4×16(方位向×距离向)，执行步骤S5和S6可基于八邻域判别规则依据格拉布斯准则实现对GB-SAR干涉影像的噪声点识别，获取噪声点集。
9)噪声点相位重构
基于反距离和相位偏差双加权算法对噪声点干涉相位进行相位重构，具体计算方式可根据步骤S7提出的计算方法进行计算。
本发明选取了典型的露天GB-SAR干涉影像进行说明，该干涉影像中包含大量的噪声误差，去噪前后效果如图3所示。
实施例2
如图4所示，本实施例提供一种GB-SAR干涉影像滤波系统，包括：
数据获取模块T1，用于获取GB-SAR监测的露天矿边坡的时间序列干涉影像；
粗滤波及形变分离模块T2，用于对所述时间序列干涉影像进行粗滤波，并对粗滤波后的影像进行形变分离，得到所述粗滤波后的影像中的稳定点；
第一计算模块T3，用于对分离出的各所述稳定点建立稳定点相位数据集，根据所述稳定点相位数据集应用最小二乘算法确定所述稳定点相位数据集的高斯函数半宽度信息；
第二计算模块T4，用于根据所述高斯函数半宽度信息确定标准差最优补偿系数和标准差限差；
噪声识别规则构建模块T5，用于基于所述标准差最优补偿系数、所述标准差限差、永久散射体点的相位均值、永久散射体点的相位标准差以及格拉布斯准则构建噪声识别规则；
噪声识别模块T6，用于确定所述时间序列干涉影像的滑动窗口，基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点；
噪声重构模块T7，用于基于反距离和相位偏差双加权算法对当前所述滑动窗口中的所述真实噪声点进行相位重构，实现GB-SAR干涉影像滤波。
滑动窗口移动模块T8，用于移动所述滑动窗口至下一所述滑动窗口，令下一所述滑动窗口为当前所述滑动窗口，并返回执行噪声识别模块中的“基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点”，直至遍历所述时间序列干涉影像的全部影像数据。
对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种GB-SAR干涉影像滤波方法，其特征在于，包括：
获取GB-SAR监测的露天矿边坡的时间序列干涉影像；
对所述时间序列干涉影像进行粗滤波，并对粗滤波后的影像进行形变分离，得到所述粗滤波后的影像中的稳定点；
对分离出的各所述稳定点建立稳定点相位数据集，根据所述稳定点相位数据集应用最小二乘算法确定所述稳定点相位数据集的高斯函数半宽度信息；
根据所述高斯函数半宽度信息确定标准差最优补偿系数和标准差限差；
基于所述标准差最优补偿系数、所述标准差限差、永久散射体点的相位均值、永久散射体点的相位标准差以及格拉布斯准则构建噪声识别规则；
确定所述时间序列干涉影像的滑动窗口，基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点；
基于反距离和相位偏差双加权算法对当前所述滑动窗口中的所述真实噪声点进行相位重构；
移动所述滑动窗口至下一所述滑动窗口，令下一所述滑动窗口为当前所述滑动窗口，并返回步骤“基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点”，直至遍历所述时间序列干涉影像的全部影像数据。
2.根据权利要求1所述的方法，其特征在于，所述对所述时间序列干涉影像进行粗滤波，并对粗滤波后的影像进行形变分离，得到所述粗滤波后的影像中的稳定点，具体包括：
根据所述时间序列干涉影像确定所述时间序列干涉影像的时间序列平均幅度均值、时间序列幅度离差、时间序列相干系数、时间序列极相干系数；
根据所述时间序列平均幅度均值、所述时间序列幅度离差和所述时间序列相干系数对所述时间序列干涉影像进行粗滤波，并根据粗滤波后的影像结合所述时间序列极相干系数分离出所述粗滤波后的影像中的稳定点。
3.根据权利要求2所述的方法，其特征在于，所述根据所述时间序列平均幅度均值、所述时间序列幅度离差和所述时间序列相干系数对所述时间序列干涉影像进行粗滤波，并根据粗滤波后的影像结合所述时间序列极相干系数分离出所述粗滤波后的影像中的稳定点，具体包括：
令所述时间序列平均幅度均值为时间序列幅度均值阈值，选择时间序列内平均幅度小于所述时间序列幅度阈值的像元作为一次粗滤波干涉影像；
设置时间序列幅度离差阈值，并基于所述一次粗滤波干涉影像，选择所述时间序列内幅度离差小于所述时间序列幅度离差阈值的像元作为二次粗滤波干涉影像；
设置时间序列相干系数阈值，并基于二次粗滤波干涉影像选择所述时间序列内平均相干系数大于所述时间序列相干系数阈值的像元作为最终粗滤波干涉影像，所述最终粗滤波干涉影像即为所述粗滤波后的影像。
设置极相干系数阈值；
基于所述粗滤波后的影像，提取所述极相干系数大于所述极相干系数阈值的永久散射点作为稳定点。
4.根据权利要求1所述的方法，其特征在于，所述根据所述稳定点相位数据集应用最小二乘算法确定所述稳定点相位数据集的高斯函数半宽度信息，具体包括：
设定所述稳定点相位数据集为(xf，f(xf))，xi为第f个稳定点的相位；f(xf)为第f个稳定点的相位的概率；f＝1,2,3……，M0；f(xf)为高斯函数，
Figure FDA0003602435960000021
a为高斯函数峰值，μg为高斯函数估计均值，S为高斯函数半宽度信息；
对所述高斯函数f(xf)取对数得到
Figure FDA0003602435960000022
令zf＝ln f(xf),
Figure FDA0003602435960000023
得出所述稳定点相位数据集的矩阵形式
Figure FDA0003602435960000024
将所述稳定点相位数据集的矩阵形式简化为Z＝XB，并基于最小二乘算法，得到矩阵B的广义最小二乘解；
基于所述矩阵B的广义最小二乘解得到所述稳定点相位数据集的高斯函数半宽度信息。
5.根据权利要求1或4所述的方法，其特征在于，所述根据所述高斯函数半宽度信息确定标准差最优补偿系数和标准差限差，具体包括：
基于所述高斯函数半宽度信息和基于极大似然估计的稳定点干涉相位样本标准差确定标准差最优补偿系数；
根据所述标准差最优补偿系数和所述高斯函数半宽度信息确定标准差限差。
6.根据权利要求5所述的方法，其特征在于，所述噪声识别规则包括：
(1)当滑动窗口中的永久散射体点的个数大于第一预设值时，且所述滑动窗口内的所述永久散射体点的相位标准差小于所述标准差限差时，采用W1公式进行噪声点估计；
W1＝(xps-μ)＞g(s,α)·λ1σs，其中，
Figure FDA0003602435960000031
λ1为标准差最优补偿系数；σ为基于极大似然估计的稳定点干涉相位样本标准差；xps为永久散射体点相位；μ为滑动窗内永久散射体点相位均值；g为Grubbs临界值，s为永久散射体点数目；α为显著性水平；σs为滑动窗口内PS点相位标准差；
(2)当滑动窗口中的永久散射体点的个数大于第一预设值时，且所述滑动窗口内的所述永久散射体点的相位标准差大于所述标准差限差时，采用W2公式进行噪声点估计；
Figure FDA0003602435960000032
其中，σg为标准差限差，
Figure FDA0003602435960000033
(3)当滑动窗口中的永久散射体点的个数小于第一预设值时，采用W3公式进行噪声点估计；
W3＝(xps-μ)＞g(s,α)·λ1σg。
7.根据权利要求1或6所述的方法，其特征在于，所述基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点，具体包括：
依据八邻域判别规则，将当前所述滑动窗口的邻域按形变区域分布的空间位置分割为八种类型，得到八种邻域分割区域；
依据所述噪声识别规则分别对当前所述滑动窗口内的永久散射体点和每一所述邻域分割区域内的永久散射体点进行噪声识别，得到当前所述滑动窗口内的初始噪声点集和每一所述邻域分割区域内的初始噪声点集；
对所述滑动窗口内的初始噪声点集和每一所述邻域分割区域内的初始噪声点集取交集，得到当前所述滑动窗口内的真实噪声点集。
8.根据权利要求1所述的方法，其特征在于，所述基于反距离和相位偏差双加权算法对当前所述滑动窗口中的所述真实噪声点进行相位重构，具体包括：
基于相位偏差最小补偿原则，对当前所述滑动窗口中的所述真实噪声点选取预设个数的参考点；
计算对每一所述真实噪声点与每一所述参考点之间的距离参数和相位偏差；
对每一所述真实噪声点，基于对应的归一化后的所述距离参数和归一化后的所述相位偏差计算权重值；
将所述权重值带入到噪声补偿公式中得到每一所述真实噪声点的估计相位。
9.根据权利要求8所述的方法，其特征在于，所述权重值的计算公式为：
Figure FDA0003602435960000041
其中，pe为第e个参考点的权值；le表示一真实噪声点与第e个参考点的归一化后的距离参数；de表示一真实噪声点与第e个参考点的归一化后的相位偏差；E表示参考点个数；
所述噪声补偿公式的表达式为：
Figure FDA0003602435960000042
其中，φnoise为噪声点相位的估计值，φe为第e个参考点的相位值。
10.根据权利要求1至9任一项所述的方法的系统，其特征在于，包括：
数据获取模块，用于获取GB-SAR监测的露天矿边坡的时间序列干涉影像；
粗滤波及形变分离模块，用于对所述时间序列干涉影像进行粗滤波，并对粗滤波后的影像进行形变分离，得到所述粗滤波后的影像中的稳定点；
第一计算模块，用于对分离出的各所述稳定点建立稳定点相位数据集，根据所述稳定点相位数据集应用最小二乘算法确定所述稳定点相位数据集的高斯函数半宽度信息；
第二计算模块，用于根据所述高斯函数半宽度信息确定标准差最优补偿系数和标准差限差；
噪声识别规则构建模块，用于基于所述标准差最优补偿系数、所述标准差限差、永久散射体点的相位均值、永久散射体点的相位标准差以及格拉布斯准则构建噪声识别规则；
噪声识别模块，用于确定所述时间序列干涉影像的滑动窗口，基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点；
噪声重构模块，用于基于反距离和相位偏差双加权算法对当前所述滑动窗口中的所述真实噪声点进行相位重构，实现GB-SAR干涉影像滤波；
滑动窗口移动模块，用于移动所述滑动窗口至下一所述滑动窗口，令下一所述滑动窗口为当前所述滑动窗口，并返回执行噪声识别模块中的“基于当前所述滑动窗口结合八邻域规则和所述噪声识别规则对当前所述滑动窗口中的永久散射体点进行离群值判别，得到当前所述滑动窗口中的真实噪声点”，直至遍历所述时间序列干涉影像的全部影像数据。
CN202210407680.3A 2022-04-19 2022-04-19 一种gb-sar干涉影像滤波方法及系统 Pending CN114779249A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210407680.3A CN114779249A (zh) 2022-04-19 2022-04-19 一种gb-sar干涉影像滤波方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210407680.3A CN114779249A (zh) 2022-04-19 2022-04-19 一种gb-sar干涉影像滤波方法及系统

Publications (1)

Publication Number Publication Date
CN114779249A true CN114779249A (zh) 2022-07-22

Family

ID=82430790

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210407680.3A Pending CN114779249A (zh) 2022-04-19 2022-04-19 一种gb-sar干涉影像滤波方法及系统

Country Status (1)

Country Link
CN (1) CN114779249A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116385314A (zh) * 2023-05-30 2023-07-04 武汉大学 面阵成像系统噪声去除方法及系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116385314A (zh) * 2023-05-30 2023-07-04 武汉大学 面阵成像系统噪声去除方法及系统
CN116385314B (zh) * 2023-05-30 2023-08-15 武汉大学 面阵成像系统噪声去除方法及系统

Similar Documents

Publication Publication Date Title
EP3229038B1 (en) Wavelet domain insar interferometric phase filtering method in combination with local frequency estimation
CN108415077B (zh) 边缘检测低序级断层识别方法
CN108663017A (zh) 一种监测城市地铁沿线地表沉降的方法
CN104020492B (zh) 一种三维地震资料的保边滤波方法
CN113960595A (zh) 一种地表形变监测方法及系统
Wang et al. A new approach to selecting coherent pixels for ground-based SAR deformation monitoring
CN109613531B (zh) 一种微变感知预警雷达的多阈值优化形变反演方法及系统
CN105118065A (zh) 小波域极化距离变换的极化sar图像变化检测方法
Jiang et al. Application of multitemporal InSAR covariance and information fusion to robust road extraction
CN112324422B (zh) 一种电成像测井缝洞识别方法、系统及孔隙结构表征方法
CN103645476A (zh) 一种合成孔径雷达差分干涉图序列的时空同质滤波方法
CN113091596B (zh) 一种基于多极化时序sar数据的地表形变监测方法
CN114779249A (zh) 一种gb-sar干涉影像滤波方法及系统
CN113281749A (zh) 一种顾及同质性的时序InSAR高相干点选取方法
CN114325706A (zh) 一种分布式散射体滤波方法
CN113204023B (zh) 联合ps目标与ds目标的双极化相位优化地表形变监测方法
Warner et al. Pine Island Glacier (Antarctica) velocities from Landsat7 images between 2001 and 2011: FFT-based image correlation for images with data gaps
CN115079172A (zh) 一种MTInSAR滑坡监测方法、设备及存储介质
CN111142165A (zh) 一种利用探地雷达获取含水层的水位信息的方法
US9563602B2 (en) Method of analyzing 3D geological structure using structure index
Shahzad et al. Nonlinear analysis of drainage systems to examine surface deformation: an example from Potwar Plateau (Northern Pakistan)
Xiang et al. PS Selection Method for and Application to GB‐SAR Monitoring of Dam Deformation
CN102314675A (zh) 基于小波高频的贝叶斯去噪方法
CN110874833A (zh) 基于超图匹配的sar图像变化检测方法
Bignami et al. Comparing and combining the capability of detecting earthquake damages in urban areas using SAR and optical data

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