CN111239734B - 一种适用于深层黄土稳定面散射体的提取方法 - Google Patents

Abstract

本发明公开了一种适用于深层黄土稳定面散射体的提取方法，有效解决了深层黄土区域SAR影像难以获得稳定点信息的问题；通过联合非参数检验与区间估计方法，既避免了传统面散射体提取中非参数检验方法耗时过久配置要求过高的缺点，也避免了区间估计方法需要先验信息假设整体分布的缺点。本发明可保证在小数量集的条件下，仍能获取足够高密度的面散射体像素，为后续的滤波和形变解算等步骤提供支持，能够扩展时序InSAR技术在黄土区域的应用范围，提高合成孔径雷达在黄土区域的形变提取精度。

Description

一种适用于深层黄土稳定面散射体的提取方法

技术领域

本发明属于合成孔径雷达干涉测量技术领域，具体涉及一种适用于深层黄土稳定面散射体的提取方法。

背景技术

黄土是干旱与半干旱气候条件下形成的一种特殊的第四纪陆相沉淀物，由未经固结的粉土颗粒组成。黄土一般有灰黄色，质地均匀，不成层，空隙度高等特性，遇水浸润后易崩解并发生沉陷。由于黄土结构的特殊性，其力学强度会急剧降低，因此黄土区域容易产生地质灾害和工程灾害。我国的黄土区域仅黄土高原面积就有64万平方千米，约占我国领土面积的6.64％，所以监测黄土区域的地表形变无论是对于工程建设还是地质灾害研究都具有十分重要的意义。

面散射体是在一定时间内具有相对稳定回波的面状地物，例如裸地、低矮灌木等,面散射体广泛的存在于非城市区域。虽然单个像素的面散射体信噪比不如永久散射体(Persistent Scatterer,PS)，不足以支撑起PS-InSAR毫米级的监测精度。但是由于同一的面状地物具有类似的反射特性，若将整个面散射体视为一个整体，那么这个整体在长时序的后向散射信息的统计特性将会相当稳定，将具有可以媲美PS的信噪比。因此，基于分布式散射体(Distributed Scatterer,DS，也被称为面散射体)的时序InSAR形变监测方法不仅在技术上能够拓宽传统PS-InSAR方法的应用领域，在实践上也能有效的服务于黄土区域长时序形变监测，满足对我国大面积黄土区域灾前预警或灾后监测的需求，具有重要的研究价值。

但已有面散射体提取方法一般是基于统计学上的非参数检验中的同质性检验方法，例如双样本的K-S检验与A-D检验。非参数检验的优势在于，不用考虑两个样本所服从的是哪种总体分布，不需要对样本具有先验知识，仅仅从样本本身的数据来判断二者是否服从同一分布。但因此同质性检验也存在一定的缺陷：对于时序SAR影像来说，样本量太大，完成检验所需要的时间太久，对计算机处理性能也要求较高。而参数估计方法虽然能够快速提取面散射体，但由于需要对样本总体分布类型进行假设，若真实地物总体分布与假设不符，则参数估计所提取的面散射体将不对应真实地物的面散射体。

发明内容

针对现有技术中的上述不足，本发明提供的适用于深层黄土稳定面散射体的提取方法解决了深层黄土区域SAR影像难以获得稳定点信息的问题。

为了达到上述发明目的，本发明采用的技术方案为：一种适用于深层黄土稳定面散射体的提取方法，包括以下步骤：

S1、以深层黄土区域的SAR影像强度图序列中任意类型的像素A为中心构造一个估计窗口，并通过K-S检验对SAR影像中的待定点和像素A进行统计同质性像素检验；

其中，待定点为构造的估计窗口中除像素A以外的像素；

S2、根据统计同质性像素检验结果，构造像素A的统计同质性像素集合S_A；

S3、对像素A和统计同质性像素集合S_A中像素的强度数据进行区间估计，构造置信度为1-α的像素A的置信区间C_A；

其中，α为显著性水平；

S4、将SAR影像中所有落入置信区间C_A的像素加入集合U_A；

S5、移动构造的估计窗口，并根据集合U_A对置信区间C_A进行参数更新；

S6、遍历SAR影像强度图序列，根据步骤S1～S5，确定每种类型的像素对应的置信区间；

S7、设置数量阈值和相关系数阈值，根据双阈值法，将每种类型的像素对应的置信区间中的统计同质性像素个数高于数量阈值且其相关系数高于相关系数阈值的像素作为面散射体，完成面散射体的提取。

进一步地，所述步骤S1中通过K-S检验对SAR影像强度图序列中的待定点和像素A进行统计同质性像素检验的方法为：

A1、依次计算构造的估计窗口中的待定点与像素A进行K-S检验时的检验统计量D_N；

A2、判断检验统计量D_N是否服从Kolmogorov分布；

若是，则该待定点与像素A来自同一总体分布，通过统计同质性像素检验；

若否，未通过统计同质性像素检验。

进一步地，所述步骤A1中的检验统计量D_N的计算公式为：

式中，sup为由待定点组成的集合R中的上确界；

N为SAR影像强度图序列中的影像景数；

为第P_i个像素的累计经验分布函数，且第P_i个像素为任意待定点，i＝1,2,3,...,M，M为集合SA中的待定点总数；

所述步骤A2中的Kolmogorov分布的表达式为：

式中，P(·)为Kolmogorov累计分布函数；

t为Kolmogorov累计分布函数中的自变量。

进一步地，所述步骤S2中的统计同质性像素集合S_A中的像素包括通过统计同质性像素检验并与像素A相邻的像素和通过完成统计同质性像素检验的像素与像素A连通的像素。

进一步地，所述步骤S3中构造置信度为1-α的像素A的置信区间C_A的过程表达式为：

式中，P{·}为为符合{·}中条件的概率；

E为集合S_A中像素的强度均值；

σ为标准差；

μ(p)为像素p的强度值；

z_1-α/2为标准正态分布位于置信度1-α处的概率。

进一步地，所述步骤S5中移动估计窗口后，对置信区间C_A进行参数更新的方法具体为：

B1、判断移动后的估计窗口中的中心像素B是否属于集合U_A；

若是，则进入步骤B2；

若否，则进入步骤B3；

B2、依次判断移动后的估计窗口中的属于集合U_A的像素是否与中心像素B相邻或通过集合U_A的像素连通；

若是，则进入步骤B4；

若否，则进入步骤B5；

B3、继续移动估计窗口，并返回步骤B1；

B4、将对应的像素加入到集合S_B中，进入步骤B6；

B5、将对应的像素从集合U_A中剔除；；

B6、将S_A与S_B的像素重复步骤三和四，更新C_A的参数。

进一步地，所述步骤S7中的数量阈值和相关系数阈值的设置方法为：

根据待提取面散射体的SAR影像强度图序列集数据情况人工设定或将这个SAR影像强度图序列集中的统计同质性像素个数中值及其像素相关系数中值分别作为数量阈值和相关系数阈值。

本发明的有益效果为：

本发明提供了一种适用于深层黄土稳定面散射体的提取方法，有效解决了深层黄土区域SAR影像难以获得稳定点信息的问题；通过联合非参数检验与区间估计方法，既避免了传统面散射体提取中非参数检验方法耗时过久配置要求过高的缺点，也避免了区间估计方法需要先验信息假设整体分布的缺点。本发明可保证在小数量集的条件下，仍能获取足够高密度的面散射体像素，为后续的滤波和形变解算等步骤提供支持，能够扩展时序InSAR技术在黄土区域的应用范围，提高合成孔径雷达在黄土区域的形变提取精度。

附图说明

图1为本发明提供的适用于深层黄土稳定面散射体的提取方法流程图。

图2为本发明提供的根据置信区间参数方法流程图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，一种适用于深层黄土稳定面散射体的提取方法，包括以下步骤：

S1、以深层黄土区域的SAR影像强度图序列中任意类型的像素A为中心构造一个估计窗口，并通过K-S检验对SAR影像强度图序列中的待定点和像素A进行统计同质性像素(Statistically Homogeneous Pixel,SHP)检验；

其中，待定点为构造的估计窗口中除像素A以外的像素；

S2、根据统计同质性像素检验结果，构造像素A的统计同质性像素集合S_A；

S3、对像素A和统计同质性像素集合S_A中像素的强度数据进行区间估计，构造置信度为1-α的像素A的置信区间C_A(Confidence interval of A，A像素的置信区间)；

其中，α为显著性水平；

S4、将SAR影像强度图序列中所有落入置信区间C_A的像素加入集合U_A(Universe ofS_A，A像素SHP的全集)；

S5、移动构造的估计窗口，并根据集合U_A对置信区间C_A进行参数更新；

S6、遍历SAR影像强度图序列，根据步骤S1～S5，确定每种类型的像素对应的置信区间；

S7、设置数量阈值和相关系数阈值，根据双阈值法，将每种类型的像素对应的置信区间中的统计同质性像素个数高于数量阈值且其相关系数高于相关系数阈值的像素作为面散射体，完成面散射体的提取。

上述步骤S1中的K-S检验为双样本检验，用于判断任意待定点与像素A是否互为SHP，其检验的方法为：

A1、依次计算构造的估计窗口中的待定点与像素A进行K-S检验时的检验统计量D_N；

A2、判断检验统计量D_N是否服从Kolmogorov分布；

若是，则该待定点与像素A来自同一总体分布，通过统计同质性像素检验；

若否，未通过统计同质性像素检验。

在本发明中的面散射体探测中，每组数据就是每个像素的强度序列所对应的N维列向量中的N个元素(N景影像的情况下)，将这N哥元素从小到大排列。可以形成一个能代表这个像素统计特性的累计分布函数S_N(x)为：

式中，x_i表示将强度值从小到大排列后，序列中第i个值。所以，每一个像素都能够依此计算出一个累积经验分布函数，第P_i个像素的累计经验分布函数记为

K-S检验是通过计算两个累计经验分布函数差值的最大值D_N作为检验统计量从而来判断两个分布是否相似，该D_N为：

式中，sup为由待定点组成的集合R中的上确界；

N₁与N₂为两个样本集的样本个数；

在时序InSAR影像中，每个像素都是N维向量，样本个数都为N个，即N₁＝N₂＝N。所以上述步骤A1中，针对时序InSAR像素的K-S检验的检验统计量可以改写为：

式中，sup为由待定点组成的集合R中的上确界；

N为SAR影像强度图序列中的影像景数；

为第P_i个像素的累计经验分布函数，且第P_i个像素为任意待定点，i＝1,2,3,...,M，M为集合SA中的待定点总数；

如果D_N服从Kolmogorov分布的话，则零假设成立，即两个样本来自于同一总体分布，上述步骤A2中的Kolmogorov分布的表达式为：

式中，P(·)为Kolmogorov累计分布函数；

t为Kolmogorov累计分布函数中的自变量。

上述步骤S2中在构建集合S_A时，显著性水平α＝1-P(D_N≤t)，一般可设为0.05，即当α≤0.05时，接受原假设，即原假设成立，P₁和P₂两个像素的振幅值服从同一分布，从而认为P₁和P₂是统计同质性像素，即通过检验，能够作为面散射体的候选点；另外根据面散射体在真实地物上总是空间相邻的特征，反应在影像上对应像素也需要是连通的，所以要剔除与中心像素不连通的孤立像素；因此，步骤S2中的统计同质性像素集合S_A中的像素包括通过统计同质性像素检验并与像素A相邻的像素和通过完成统计同质性像素检验的像素与像素A连通的像素。

上述步骤S3中，对于在步骤S2中探测出的S_A来说，作为A的SHP的集合，理论上都属于地表同一类地物，则通过S_A中元素的值可以去估计总体可能属于的分布类型，一般来说则服从高斯分布；因此，上述步骤S3中构造置信度为1-α的像素A的置信区间C_A的过程表达式为：

式中，P{·}为为符合{·}中条件的概率；

E为集合S_A中像素的强度均值；

σ为标准差；

μ(p)为像素p的强度值；

z_1-α/2为标准正态分布位于置信度1-α处的概率。

其中，区间估计中显著性水平α通常也设为0.05，至此我们构造出了以像素A为特征点的95％置信度的区间C_A，即落入C_A范围内的像素，至少有95％可能性是像素A的SHP。

上述步骤S4中，当我们构造出C_A后，就不再需要利用双样本的K-S检验来判断某像素是否和像素A为SHP了，只需判断该像素是否落入C_A即可。而对整幅影像的判断，最终就能够提取出在整幅影像中与A为同一地物类型的像素集合U_A。

如步骤S2所言，U_A中与A连通的像素就可以认为是像素A的最终SHP，在步骤S5中，为了结果的精确性，随着被标记为A的最终SHP的像素的变化，可以对C_A的参数进行更新，以保证对A像素所属地物的各种参数的高精度估计，就能够确保最后提取的这类面散射体所对应的像素的精确性。我们对像素A的置信区间进行参数更新时，移动估计窗口后，对置信区间C_A进行参数更新的方法具体为：

B1、判断移动后的估计窗口中的中心像素B是否属于集合U_A；

若是，则进入步骤B2；

若否，则进入步骤B3；

B2、依次判断移动后的估计窗口中的属于集合U_A的像素是否与中心像素B相邻或通过集合U_A的像素连通；

若是，则进入步骤B4；

若否，则进入步骤B5；

B3、继续移动估计窗口，并返回步骤B1；

步骤B3中，若该像素不属于U_A，则说明它与像素A不是同一种地物，在和面就需要对该新的地物像素的SHP进行探测；

B4、将对应的像素加入到集合S_B中，进入步骤B6；

B5、将对应的像素从集合U_A中剔除；；

B6、将S_A与S_B的像素重复步骤三和四，更新C_A的参数。

在步骤B6中对C_A的参数进行更新时，因为像素B属于U_A，即B与A是同质性像素，所以S_A和S_B将会是高度相似的集合，可能只有少数像素增减的变化，对C_A的更新就是将S_A与S_B并集中的像素作为样本去求置信区间。

上述步骤S7中的数量阈值和相关系数阈值的设置方法为：

根据待提取面散射体的SAR影像强度图序列集数据情况人工设定或将这个SAR影像强度图序列集中的统计同质性像素个数中值及其像素相关系数中值分别作为数量阈值和相关系数阈值。另外也可以根据自己的需求进行增减阈值类型与数量，自由调节。

在本发明的一个实施例中，利用20景Sentinel 1-a影像进行了实验，验证了本发明提出的面散射体提取算法能够较好的提取研究区域的稳定目标，在较严格的SHP数目阈值与相干系数阈值下，点目标数量仍旧能够达到影像总像素个数的21.41％，远高于ADI阈值为0.5时提取的PS点比例的0.31％，说明本方法能够有效弥补黄土覆盖区域PS点密度不足的缺陷，从而完成时序形变解算任务。

Claims (7)

1.一种适用于深层黄土稳定面散射体的提取方法，其特征在于，包括以下步骤：

S1、以深层黄土区域的SAR影像强度图序列中任意类型的像素A为中心构造一个估计窗口，并通过K-S检验对SAR影像中的待定点和像素A进行统计同质性像素检验；

其中，待定点为构造的估计窗口中除像素A以外的像素；

S2、根据统计同质性像素检验结果，构造像素A的统计同质性像素集合S_A；

S3、对像素A和统计同质性像素集合S_A中像素的强度数据进行区间估计，构造置信度为1-α的像素A的置信区间C_A；

其中，α为显著性水平；

S4、将SAR影像中所有落入置信区间C_A的像素加入集合U_A；

S5、移动构造的估计窗口，并根据集合U_A对置信区间C_A进行参数更新；

S6、遍历SAR影像强度图序列，根据步骤S1～S5，确定每种类型的像素对应的置信区间；

S7、设置数量阈值和相关系数阈值，根据双阈值法，将每种类型的像素对应的置信区间中的统计同质性像素个数高于数量阈值且其相关系数高于相关系数阈值的像素作为面散射体，完成面散射体的提取。

2.根据权利要求1所述的适用于深层黄土稳定面散射体的提取方法，其特征在于，所述步骤S1中通过K-S检验对SAR影像强度图序列中的待定点和像素A进行统计同质性像素检验的方法为：

A1、依次计算构造的估计窗口中的待定点与像素A进行K-S检验时的检验统计量D_N；

A2、判断检验统计量D_N是否服从Kolmogorov分布；

若是，则该待定点与像素A来自同一总体分布，通过统计同质性像素检验；

若否，未通过统计同质性像素检验。

3.根据权利要求2所述的适用于深层黄土稳定面散射体的提取方法，其特征在于，所述步骤A1中的检验统计量D_N的计算公式为：

式中，sup为由待定点组成的集合R中的上确界；

N为SAR影像强度图序列中的影像景数；

为第P_i个像素的累计经验分布函数，且第P_i个像素为任意待定点，i＝1,2,3,...,M，M为集合S_A中的待定点总数；

所述步骤A2中的Kolmogorov分布的表达式为：

式中，P(·)为Kolmogorov累计分布函数；

t为Kolmogorov累计分布函数中的自变量。

4.根据权利要求1所述的适用于深层黄土稳定面散射体的提取方法，其特征在于，所述步骤S2中的统计同质性像素集合S_A中的像素包括通过统计同质性像素检验并与像素A相邻的像素和通过完成统计同质性像素检验的像素与像素A连通的像素。

5.根据权利要求1所述的适用于深层黄土稳定面散射体的提取方法，其特征在于，所述步骤S3中构造置信度为1-α的像素A的置信区间C_A的过程表达式为：

式中，P{·}为为符合{·}中条件的概率；

E为集合S_A中像素的强度均值；

σ为标准差；

μ(p)为像素p的强度值；

z_1-α/2为标准正态分布位于置信度1-α处的概率；

N为SAR影像强度图序列中的影像景数。

6.根据权利要求5所述的适用于深层黄土稳定面散射体的提取方法，其特征在于，所述步骤S5中移动估计窗口后，对置信区间C_A进行参数更新的方法具体为：

B1、判断移动后的估计窗口中的中心像素B是否属于集合U_A；

若是，则进入步骤B2；

若否，则进入步骤B3；

B2、依次判断移动后的估计窗口中的属于集合U_A的像素是否与中心像素B相邻或通过集合U_A的像素连通；

若是，则进入步骤B4；

若否，则进入步骤B5；

B3、继续移动估计窗口，并返回步骤B1；

B4、将对应的像素加入到集合S_B中，进入步骤B6；

B5、将对应的像素从集合U_A中剔除；

B6、将S_A与S_B的像素重复步骤S3和S4，更新C_A的参数。

7.根据权利要求1所述的适用于深层黄土稳定面散射体的提取方法，其特征在于，所述步骤S7中的数量阈值和相关系数阈值的设置方法为：

根据待提取面散射体的SAR影像强度图序列集数据情况人工设定或将这个SAR影像强度图序列集中的统计同质性像素个数中值及其像素相关系数中值分别作为数量阈值和相关系数阈值。

Images