CN116299459A - 一种基于地形和信号反射强度的InSAR掩膜方法 - Google Patents

Abstract

本发明涉及星载图像处理技术领域，具体涉及一种基于地形和信号反射强度的InSAR掩膜方法，包括：根据卫星轨道参数和DEM数据，生成考虑地形因素的第一掩膜；定义一个滑动窗口，根据窗口内所有像元的平均幅度信息，生成考虑信号反射强度的第二掩膜；通过第一掩膜将SAR原始影像中的阴影像元和叠掩像元筛除，并确定疑似阴影像元和疑似叠掩像元；通过第二掩膜将SAR原始影像中疑似叠掩像元区域、疑似阴影像元区域、非疑似叠掩像元区域和非疑似阴影像元区域的弱反射像元筛除。本发明可同时剔除阴影像元、叠掩像元和弱反射像元，能够提取到更高质量的点目标，提高结果精度。

Description

一种基于地形和信号反射强度的InSAR掩膜方法

技术领域

本发明涉及星载图像处理技术领域，更具体的说是涉及一种基于地形和信号反射强度的InSAR掩膜方法。

背景技术

星载合成孔径雷达干涉测量(Interferometry Synthetic Aperture Radar,InSAR)技术作为一种快速高效获取地形形变信息的技术手段，因为其具有全天候、高精度、低成本等优点，在过去二十多年间，该技术取得了迅速的发展。2001年，Ferretti提出永久散射体雷达干涉测量(Permanent Scatterer InSAR,PSInSAR)技术，将InSAR技术带入了时间序列研究阶段，该技术从时序SAR影像中挑选在时间纬度上幅度保持稳定的高相干点进行干涉测量，较好地克服时空失相干和大气延迟等影响，实现精确可靠反演地表信息。

但是因为其侧视成像方式，单从一个方向来获取地表信息，容易出现阴影和叠掩区域，若在阴影或者叠掩区域取点会对结果造成较大误差影响，所以在PSInSAR处理之前必须将阴影和叠掩区域进行去除。为提取阴影和叠掩像元，已有众多学者先后提出基于幅度和相干系数阈值检测、基于成像几何关系检测、恒虚警检测、局部频率检测等多种方法。如在2011年，张继贤学者提出的利用DEM数据与卫星轨道数据来剔除阴影和叠掩像元的方法得到广泛应用，但该方法会受到DEM的精度或者地形的影响，存在漏剔的现象，不仅如此，因为在地球表面的反射和散射机制中存在许多干扰源，所以该方法还会保留下较多的反射强度非常弱的像元。保留较多的弱反射像元还会降低最终PS点质量，从而影响InSAR结果的精度。

发明内容

有鉴于此，本发明提供了一种基于地形和信号反射强度的InSAR掩膜方法，可同时剔除阴影像元、叠掩像元和弱反射像元，能够提取到更高质量的点目标，提高结果精度。

为了实现上述目的，本发明采用如下技术方案：

一种基于地形和信号反射强度的InSAR掩膜方法，包括：

根据卫星轨道参数和DEM数据，生成考虑地形因素的第一掩膜；

定义一个滑动窗口，根据窗口内所有像元的平均幅度信息，生成考虑信号反射强度的第二掩膜；

通过所述第一掩膜将SAR原始影像中的阴影像元和叠掩像元筛除，并确定疑似阴影像元和疑似叠掩像元；

通过所述第二掩膜将SAR原始影像中疑似叠掩像元区域、疑似阴影像元区域、非疑似叠掩像元区域和非疑似阴影像元区域的弱反射像元筛除。

进一步的，所述根据卫星轨道参数和DEM数据，生成考虑地形因素的第一掩膜，包括：

计算DEM数据中每个像元的坡向角β_n，利用卫星轨道参数中的轨道方位角Ω_s将DEM数据中像元的坡向角转换为基于卫星飞行轨道的方位向坡向角β_s；

判断方位向坡向角β_s是否大于180°，将大于180°的像元分为迎坡面，将小于180°的像元分为背坡面。

计算雷达俯角β和像元的距离向坡度角θ_r，如果迎坡面上像元的距离向坡度角θ_r与雷达俯角β之和大于等于90°，则视为叠掩像元，否则，设定阈值γ_layover，且30°≤γ_layover＜90°，如果像元的距离向坡度角θ_r与雷达俯角β的和值处于γ_layover到90°之间，则视为疑似叠掩像元；

如果背坡面上像元的雷达俯角β减去距离向坡度角θ_r的差值小于等于0°，则视为阴影像元；否则，设定阈值γ_shadow，且0°＜γ_shadow≤60°，如果像元的雷达俯角β减去距离向坡度角θ_r的差值处于0与γ_shadow之间，则视为疑似阴影像元。

进一步的，在所述根据卫星轨道参数和DEM数据，生成考虑地形因素的第一掩膜之前，还包括：使用克里金方法对DEM数据插值，使DEM数据的分辨率与SAR影像分辨率一致，且DEM数据像元与SAR影像像元一一对应。

进一步的，迎坡面和背坡面的确定过程为：

读取卫星飞行轨道的方位角Ω_s，并计算DEM数据中每个像元的坡向角β_n；

根据卫星飞行轨道方位角Ω_s和像元坡向角β_n计算方位向坡向角β_s，计算方式如下：

当β_n＞Ω_s时，β_s＝β_n-Ω_s；

当β_n≤Ω_s时，β_s＝360-Ω_s+β_n；

如果像元的方位向坡向角β_s大于180°，则认为是背坡面的像元，否则，是迎坡面的像元。

进一步的，雷达俯角β的计算公式如下：

其中，R_s为在地心坐标系中卫星传感器的位置矢量，R_T为在地心坐标系中地物的位置矢量，R为在地心坐标系中传感器到地物的斜距；

像元的距离向坡度角θ_r的计算公式如下：

其中，θ为像元的坡度角，β_s为像元的方位向坡向角。

进一步的，像元的坡度角θ和坡向角β_n的确定过程为：

确定待计算像元及其临近像元的分布关系，分布关系如下：

e5	e2	e6
			e1	e	e3
e8	e4	e7

其中，e为待计算坡向角β_n和坡度角θ的像元，e1、e2…e8为待计算像元的临近像元的高程；

利用下式计算像元e的坡向角β_n：

利用下式计算像元e的坡度角θ：

进一步的，所述定义一个滑动窗口，根据滑动窗口内所有像元的平均幅度信息，生成考虑信号反射强度的第二掩膜，包括：

使用配准之后的多景SAR影像，计算各幅影像中相同位置上所有像元的平均幅度，生成平均幅度影像；

定义一个可设置大小的滑动窗口，在所述平均幅度影像中，计算窗口内所有像元幅度的最大值和平均值，并与窗口中心位置像元的幅度值进行比较；

当窗口中心位置位于疑似叠掩像元或疑似阴影像元时，比较窗口中心位置像元的幅度值是否大于等于窗口内所有像元幅度的最大值，若大于，则认为该像元为强反射像元，否则，认为该像元为弱反射像元；

当窗口中心位置位于非疑似叠掩像元或非疑似阴影像元时，比较窗口中心位置像元的幅度值是否大于等于窗口内所有像元幅度的平均值，若大于，则认为该像元为强反射像元，否则，认为该像元为弱反射像元。

进一步的，所述平均幅度影像的像元幅度A_mu计算公式为：

其中，A₁、A₂、A_n为各幅SAR影像中相同位置上的像元的幅度值，n为影像的个数；

像元的幅度A_n的计算公式为：

其中，I为像元的虚部，R为像元的实部。

经由上述的技术方案可知，与现有技术相比，本发明公开提供了一种基于地形和信号反射强度的InSAR掩膜方法，首先，根据DEM数据和卫星轨道参数，生成第一掩膜，通过第一掩膜剔除叠掩像元和阴影像元，同时，为了尽可能消除叠掩和阴影像元对结果影响，设定阈值，提取疑似叠掩像元和疑似阴影像元。然后定义滑动窗口大小，通过计算平均幅度，生成第二掩膜，通过第二掩膜在疑似叠掩和疑似阴影区域采用窗口中心位置幅度值大于等于窗口内幅度最大值作为筛选条件，在非叠掩和非阴影区域采用窗口中心位置幅度值大于等于窗口内幅度的平均值作为筛选条件，挑选强反射像元和弱反射像元，最后将所有弱反射像元、阴影像元、叠掩像元从SAR原始影像中掩去，保留下来的像元参与后续的PSInSAR处理。本发明不仅考虑了地形影响，还考虑了反射强度对最终结果精度的影响，并提出了在疑似叠掩和疑似阴影像元位置与非疑似叠掩和非疑似阴影像元位置使用不同判别尺度生成掩膜的方法，在保证了最大程度去除地形因素对结果影响的前提下，还尽可能多地保留像元个数，提高了数据的精度和可靠性，从而也保证了最终PS点的质量。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。

图1为本发明提供的基于地形和信号反射强度的InSAR掩膜方法的流程图；

图2为本发明提供的第一掩膜的生成流程图；

图3为本发明提供的第二掩膜的生成流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明实施例公开了一种基于地形和信号反射强度的InSAR掩膜方法，包括以下步骤：

根据卫星轨道参数和DEM数据，生成考虑地形因素的第一掩膜；

定义一个滑动窗口，根据窗口内所有像元的平均幅度信息，生成考虑信号反射强度的第二掩膜；

通过第一掩膜将SAR原始影像中的阴影像元和叠掩像元筛除，并确定疑似阴影像元和疑似叠掩像元；

通过第二掩膜将SAR原始影像中疑似叠掩像元区域、疑似阴影像元区域、非疑似叠掩像元区域和非疑似阴影像元区域的弱反射像元筛除。

将SAR原始影像中的阴影像元、叠掩像元和弱反射像元筛除之后，只利用保留下来的像元参与到PSInSAR的后续处理。

在一个具体实施例中，如图2所示，第一掩膜的生成过程包括：

计算DEM数据中每个像元的坡向角β_n，利用卫星轨道参数中的轨道方位角Ω_s将DEM数据中像元的坡向角转换为基于卫星飞行轨道的方位向坡向角β_s；

判断方位向坡向角β_s是否大于180°，将大于180°的像元分为迎坡面，将小于180°的像元分为背坡面。

计算雷达俯角β和像元的距离向坡度角θ_r，如果迎坡面上像元的距离向坡度角θ_r与雷达俯角β之和大于等于90°，则视为叠掩像元，否则，设定阈值γ_layover，且30°≤γ_layover＜90°，如果像元的距离向坡度角θ_r与雷达俯角β的和值处于γ_layover到90°之间，则视为疑似叠掩像元，否则，为非疑似叠掩像元；

如果背坡面上像元的雷达俯角β减去距离向坡度角θ_r的差值小于等于0°，则视为阴影像元；否则，设定阈值γ_shadow，且0°＜γ_shadow≤60°，如果像元的雷达俯角β减去距离向坡度角θ_r的差值处于0与γ_shadow之间，则视为疑似阴影像元，否则，为非疑似阴影像元。

具体来说，迎坡面和背坡面的确定过程为：

读取卫星飞行轨道的方位角Ω_s，并计算DEM数据中每个像元的坡向角β_n；

根据卫星飞行轨道方位角Ω_s和像元坡向角β_n计算方位向坡向角β_s，计算方式如下：

当β_n＞Ω_s时，β_s＝β_n-Ω_s；

当β_n≤Ω_s时，β_s＝360-Ω_s+β_n；

如果像元的方位向坡向角β_s大于180°，则认为是背坡面的像元，否则，是迎坡面的像元。

雷达俯角β的计算公式如下：

其中，R_s为在地心坐标系中卫星传感器的位置矢量，R_T为在地心坐标系中地物的位置矢量，R为在地心坐标系中传感器到地物的斜距；

像元的距离向坡度角θ_r的计算公式如下：

其中，θ为像元的坡度角，β_s为像元的方位向坡向角；像元的坡度角θ和坡向角β_n的确定过程为：

确定待计算像元及其临近像元的分布关系，分布关系如下：

e5	e2	e6
			e1	e	e3
e8	e4	e7

其中，e为待计算坡向角β_n和坡度角θ的像元，e1、e2…e8为待计算像元的临近像元的高程；

利用下式计算像元e的坡向角β_n：

利用下式计算像元e的坡度角θ：

更有利的，在生成考虑地形因素的第一掩膜之前，还包括：使用克里金方法对DEM数据插值，使DEM数据的分辨率与SAR影像分辨率一致，且DEM数据像元与SAR影像像元一一对应。

在另一个具体实施例中，如图3所示，第二掩膜的生成过程包括：

使用配准之后的多景SAR影像，计算各幅影像中相同位置上所有像元的平均幅度，生成平均幅度影像；其中，多景SAR影像采用已有方法进行配准，举例来说，其配准过程为：利用轨道参数并结合斜距方程、多普勒方程、椭球方程计算得到主影像在各个辅影像上的偏移量，然后利用偏移量通过最小二乘方法计算变换模型参数，最后通过重采样与插值得到配准后的多景SAR影像；

定义一个可设置大小的滑动窗口，如3×3或5×5窗口；在平均幅度影像中，计算窗口内所有像元幅度的最大值和平均值，并与窗口中心位置像元的幅度值进行比较；

当窗口中心位置位于疑似叠掩像元或疑似阴影像元时，比较窗口中心位置像元的幅度值是否大于等于窗口内所有像元幅度的最大值，若大于，则认为该像元为强反射像元，否则，认为该像元为弱反射像元；

当窗口中心位置位于非疑似叠掩像元或非疑似阴影像元时，比较窗口中心位置像元的幅度值是否大于等于窗口内所有像元幅度的平均值，若大于，则认为该像元为强反射像元，否则，认为该像元为弱反射像元。

其中，平均幅度影像的像元幅度A_mu计算公式为：

其中，A₁、A₂、A_n为各幅SAR影像中相同位置上的像元的幅度值，n为影像的个数；

像元的幅度A_n的计算公式为：

其中，I为像元的虚部，R为像元的实部。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (8)

1.一种基于地形和信号反射强度的InSAR掩膜方法，其特征在于，包括以下步骤：

根据卫星轨道参数和DEM数据，生成考虑地形因素的第一掩膜；

定义一个滑动窗口，根据窗口内所有像元的平均幅度信息，生成考虑信号反射强度的第二掩膜；

通过所述第一掩膜将SAR原始影像中的阴影像元和叠掩像元筛除，并确定疑似阴影像元和疑似叠掩像元；

通过所述第二掩膜将SAR原始影像中疑似叠掩像元区域、疑似阴影像元区域、非疑似叠掩像元区域和非疑似阴影像元区域的弱反射像元筛除。

2.根据权利要求1所述的基于地形和信号反射强度的InSAR掩膜方法，其特征在于，所述根据卫星轨道参数和DEM数据，生成考虑地形因素的第一掩膜，包括：

计算DEM数据中每个像元的坡向角β_n，利用卫星轨道参数中的轨道方位角Ω_s将DEM数据中像元的坡向角转换为基于卫星飞行轨道的方位向坡向角β_s；

判断方位向坡向角β_s是否大于180°，将大于180°的像元分为迎坡面，将小于180°的像元分为背坡面。

计算雷达俯角β和像元的距离向坡度角θ_r，如果迎坡面上像元的距离向坡度角θ_r与雷达俯角β之和大于等于90°，则视为叠掩像元，否则，设定阈值γ_layover，且30°≤γ_layover＜90°，如果像元的距离向坡度角θ_r与雷达俯角β的和值处于γ_layover到90°之间，则视为疑似叠掩像元；

如果背坡面上像元的雷达俯角β减去距离向坡度角θ_r的差值小于等于0°，则视为阴影像元；否则，设定阈值γ_shadow，且0°＜γ_shadow≤60°，如果像元的雷达俯角β减去距离向坡度角θ_r的差值处于0与γ_shadow之间，则视为疑似阴影像元。

3.根据权利要求2所述的基于地形和信号反射强度的InSAR掩膜方法，其特征在于，在所述根据卫星轨道参数和DEM数据，生成考虑地形因素的第一掩膜之前，还包括：使用克里金方法对DEM数据插值，使DEM数据的分辨率与SAR影像分辨率一致，且DEM数据像元与SAR影像像元一一对应。

4.根据权利要求2所述的基于地形和信号反射强度的InSAR掩膜方法，其特征在于，迎坡面和背坡面的确定过程为：

读取卫星飞行轨道的方位角Ω_s，并计算DEM数据中每个像元的坡向角β_n；

根据卫星飞行轨道方位角Ω_s和像元坡向角β_n计算方位向坡向角β_s，计算方式如下：

当β_n＞Ω_s时，β_s＝β_n-Ω_s；

当β_n≤Ω_s时，β_s＝360-Ω_s+β_n；

如果像元的方位向坡向角β_s大于180°，则认为是背坡面的像元，否则，是迎坡面的像元。

5.根据权利要求2所述的基于地形和信号反射强度的InSAR掩膜方法，其特征在于，雷达俯角β的计算公式如下：

其中，R_s为在地心坐标系中卫星传感器的位置矢量，R_T为在地心坐标系中地物的位置矢量，R为在地心坐标系中传感器到地物的斜距；

像元的距离向坡度角θ_r的计算公式如下：

其中，θ为像元的坡度角，β_s为像元的方位向坡向角。

6.根据权利要求5所述的基于地形和信号反射强度的InSAR掩膜方法，其特征在于，像元的坡度角θ和坡向角β_n的确定过程为：

确定待计算像元及其临近像元的分布关系，分布关系如下：

e5 e2 e6 e1 e e3 e8 e4 e7

其中，e为待计算坡向角β_n和坡度角θ的像元，e1、e2…e8为待计算像元的临近像元的高程；

利用下式计算像元e的坡向角β_n：

利用下式计算像元e的坡度角θ：

7.根据权利要求1所述的基于地形和信号反射强度的InSAR掩膜方法，其特征在于，所述定义一个滑动窗口，根据滑动窗口内所有像元的平均幅度信息，生成考虑信号反射强度的第二掩膜，包括：

使用配准之后的多景SAR影像，计算各幅影像中相同位置上所有像元的平均幅度，生成平均幅度影像；

定义一个可设置大小的滑动窗口，在所述平均幅度影像中，计算窗口内所有像元幅度的最大值和平均值，并与窗口中心位置像元的幅度值进行比较；

当窗口中心位置位于疑似叠掩像元或疑似阴影像元时，比较窗口中心位置像元的幅度值是否大于等于窗口内所有像元幅度的最大值，若大于，则认为该像元为强反射像元，否则，认为该像元为弱反射像元；

当窗口中心位置位于非疑似叠掩像元或非疑似阴影像元时，比较窗口中心位置像元的幅度值是否大于等于窗口内所有像元幅度的平均值，若大于，则认为该像元为强反射像元，否则，认为该像元为弱反射像元。

8.根据权利要求7所述的基于地形和信号反射强度的InSAR掩膜方法，其特征在于，所述平均幅度影像的像元幅度A_mu计算公式为：

其中，A₁、A₂、A_n为各幅SAR影像中相同位置上的像元的幅度值，n为影像的个数；

像元的幅度A_n的计算公式为：

其中，I为像元的虚部，R为像元的实部。

Images