CN106023157A

CN106023157A - 一种基于sar图像的山区地表微形变信息提取方法

Info

Publication number: CN106023157A
Application number: CN201610304277.2A
Authority: CN
Inventors: 陈彦; 陈浩; 童玲
Original assignee: University of Electronic Science and Technology of China
Current assignee: University of Electronic Science and Technology of China
Priority date: 2016-05-10
Filing date: 2016-05-10
Publication date: 2016-10-12
Anticipated expiration: 2036-05-10
Also published as: CN106023157B

Abstract

本发明公开了一种基于SAR图像的山区地表微形变信息提取方法，通过对SAR影像组进行精配准并提取振幅与相位信息，然后根据像元点的振幅离差指数阈值法和归一化同极化相位差阈值法选取永久散射点；同时利用SAR影像组和高程数据进行差分干涉处理得到差分相位图；再利用二次滤波算法选取分布式散射点，并对分布式散射点的差分相位值进行优化；再将永久散射点和分布式散射点合并为稳定散射点，构造稳定散射点的Delaunay三角网；最后利用网络平差法对稳定散射点的差分相位进行求解，最终得到地表微形变信息。

Description

一种基于SAR图像的山区地表微形变信息提取方法

技术领域

本发明属于图像处理技术领域，更为具体地讲，涉及一种基于SAR图像的山区地表微形变信息提取方法。

背景技术

永久散射体差分干涉测量技术(PS-InSAR)是一种利用时间序列SAR影像组提取地表形变信息的技术。PS-InSAR技术能够全天时全天候对地表进行监测而不受天气影响，因而为大范围的地表形变监测提供了快速、有效且经济的手段。中国地域辽阔，自然环境复杂，PS-InSAR技术能够在观测条件有限的情况下实现地表信息提取。相对于传统的位移计、水准仪等实地测量，该技术能够极大的提高覆盖范围、降低成本、提高空间分辨率。PS-InSAR提取地表形变信息的步骤包括：配准、差分干涉处理、永久散射点的提取、相位分解等。由于处理过程环环相扣，任意一个环节的误差都将直接被带入下一步处理中，因此，需要每一步都达到较高的精度，才能保证最后结果的准确性。

PS-InSAR技术是一种对永久散射体进行差分干涉处理，从而提取出整幅影像地表形变信息的技术。永久散射体指的是研究区域中具有高强度稳定散射特性的地物目标，可能是人工建筑、裸露岩石等。在SAR图像上，永久散射体所对应的像元具有较高的亮度和信噪比，在时间轴上呈现出很好的稳定性。可以通过对时间序列SAR影像组进行分析从而选取出永久散射点，永久散射点选取算法的准确度和精度将直接影响PS-InSAR的技术的精度。

目前关于永久散射点选取算法虽然较多，但存在着许多不足之处。传统的永久散射点选取算法受到地表类型限制，人工建筑、裸露岩石砂砾等地面目标具有较强的相关性因而能够被选取为永久散射点，然而植被覆盖较多的山区往往提取不到足够的永久散射点。同时算法对时间序列影像的数量要求高，传统PS-InSAR算法往往需要十几甚至二十几幅图像才能保证永久散射点的可靠性，因此，必须对现有的技术进行改进，以实现提高山区稳定散射点的数量，同时降低对SAR图像数量的要求。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于SAR图像的山区地表微形变信息提取方法，通过提取稳定散射点，以实现山区的地表微形变信息提取。

为实现上述发明目的，本发明为一种基于SAR图像的山区地表微形变信息提取方法，其特征在于，包括以下步骤：

(1)、对时间序列SAR影像组进行预处理

在预先采集的N幅时间序列SAR影像组中，根据综合相关函数选取一幅SAR图像作为公共主影像，其余的为从影像；再利用sarscape软件对每幅从影像按主影像进行配准，并提取每幅SAR图像的振幅与相位信息，再为每幅SAR图像各生成一幅振幅图与相位图；

(2)、差分干涉

将时间序列SAR影像组进行干涉处理并去平，得到干涉相位图；同时将外部高程数据转化为高程相位分量；然后将干涉相位图的相位与高程相位分量进行差分处理，得到差分干涉相位图；

(3)、利用振幅离差阈值法与归一化同极化相位差阈值法选取永久散射点

(3.1)、根据振幅图，计算每个像元点的振幅离差指数D_A(i,j)：

D_{A} (i, j) = \frac{σ_{A} (i, j)}{m_{A} (i, j)}

其中，(i,j)表示像元点的坐标，m_A(i,j)、σ_A(i,j)分别为振幅的均值与方差：

m_{A} (i, j) = \frac{1}{N} Σ_{k = 1}^{N} A_{k} (i, j), σ_{A} = \sqrt{\frac{1}{N} Σ_{k = 1}^{N} {(A_{k} (i, j) - m_{A} (i, j))}^{2}}

其中，A_k(i,j)表示第k幅振幅图的振幅值，N表示振幅图或相位图的总个数；

(3.2)、根据相位图，计算每个像元点的归一化同极化相位差χ(i,j)：

其中分别为第k相位图的HH极化与VV极化相位值；

(3.3)、根据振幅离差指数D_A(i,j)和归一化同极化相位差χ(i,j)选取永久散射点；

(4)、利用双重自适应性滤波算法选取分布式散射点

(4.1)、选取中心像元点

根据振幅图，计算每个像元点的同极化振幅幅比ρ(i,j)：

ρ (i, j) = \frac{1}{N} Σ_{k = 1}^{N} \frac{A_{k}^{H H} (i, j)}{A_{k}^{V V} (i, j)}

其中，表示第k幅振幅图在HH极化振幅值，表示第k幅振幅图在VV极化振幅值；

设定振幅离差指数阈值阈值D'_T、同极化振幅幅比阈值将振幅离差指数满足D_A(i,j)＜D'_T且同极化振幅幅比值满足的像元点设置为中心像元点；

(4.2)、利用统计学中的KS检验以中心像元点为中心、滤波窗口内，判定中心像元点的同质点；

(4.3)、确定同质点族：在滤波窗口内，舍去没有与中心像元点直接或间接相连的同质点，将中心像元点与其直接或间接相连的同质点作为一个同质点族；

(4.4)、判断同质点族大小是否超过阈值T，如果同质点族大小超过阈值T，则将该同质点归入一次候选分布式散射点集CDS1中；

(4.5)、计算同质点族的振幅均值，并作为CDS1的振幅值；

(4.6)、对CDS1重复步骤(4.2)～(4.4)，再次判定同质点族，当同质点族大小超过阈值T时，归入二次候选分布式散射点集CDS2中；

(4.7)、根据差分干涉相位图计算出CDS2同质点族的差分相位均值，并作为CDS2的差分相位值；同时计算CDS2中各像元点的振幅离差指数，将振幅离差指数D_A(i,j)小于预设阈值D″_T的像元点设置为分布式散射点；

(5)、形变信息提取

将永久散射点与分布式散射点合并为稳定散射点集，建立Delaunay三角网，利用网络平差法对稳定散射点的差分相位进行相位分解，从而提取出山区地表微形变信息。

本发明的发明目的是这样实现的：

本发明一种基于SAR图像的山区地表微形变信息提取方法，通过对SAR影像组进行精配准并提取振幅与相位信息，然后根据像元点的振幅离差指数阈值法和归一化同极化相位差阈值法选取永久散射点；同时利用SAR影像组和高程数据进行差分干涉处理得到差分相位图；再利用二次滤波算法选取分布式散射点，并对分布式散射点的差分相位值进行优化；再将永久散射点和分布式散射点合并为稳定散射点，构造稳定散射点的Delaunay三角网；最后利用网络平差法对稳定散射点的差分相位进行求解，最终得到地表微形变信息。

同时，本发明基于SAR图像的山区地表微形变信息提取方法还具有以下有益效果：

本发明采用的分布式散射点选取算法和传统算法相比，在降低计算量的同时，提高了分布式散射点的置信度；本发明采用的永久散射点算法在很大程度上降低了对影像数量的要求。将这两种算法运用于PS-InSAR技术中，将PS-InSAR技术的运用领域从城市地区扩展到了地形复杂、植被覆盖的山区。

附图说明

图1是本发明基于SAR图像的山区地表微形变信息提取方法流程图。

图2是主影像HH极化的振幅与相位图。

图3是提取到的PS点。

图4是两次滤波处理的示意图。

图5是提取到的DS点。

图6是PS点与DS点共同构成的点集合。

图7是根据网络平差法求解的线性形变速率分布图。

图8是根据网络平差法求解的高程修正量分布图。

图9是经过地理校正后的线性形变速率图。

图10是研究区域2013年四月至九月地表形变累积量分布图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

为了方便描述，先对具体实施方式中出现的相关专业术语进行说明：

SAR(Synthetic Aperture Radar)：合成孔径雷达；

KS(Kolmogorov-Smirnov)：柯尔莫哥洛夫-斯摩洛夫检验。

在本实施例中，如图1所示，本发明一种基于SAR图像的山区地表微形变信息提取方法，其特征在于，包括以下步骤：

S1、对时间序列SAR影像组进行预处理

在预先采集的N幅时间序列SAR影像组中，根据综合相关函数选取一幅SAR图像作为公共主影像，其余的为从影像；再利用sarscape软件对每幅从影像按主影像进行配准，并提取每幅SAR图像的振幅与相位信息，再为幅SAR图像各生成一幅振幅图与相位图；

在本实施例中，主影像生成的振幅图与相位图如图2所示，振幅图中白色的亮点代表高振幅值，图中白色亮条纹是由于山体叠掩效应造成的。

S2、差分干涉

S3、利用振幅离差阈值法与归一化同极化相位差阈值法选取永久散射点

S3.1、根据振幅图，计算每个像元点的振幅离差指数D_A(i,j)：

D_{A} (i, j) = \frac{σ_{A} (i, j)}{m_{A} (i, j)}

m_{A} (i, j) = \frac{1}{N} Σ_{k = 1}^{N} A_{k} (i, j), σ_{A} = \sqrt{\frac{1}{N} Σ_{k = 1}^{N} {(A_{k} (i, j) - m_{A} (i, j))}^{2}}

S3.2、根据相位图，计算每个像元点的归一化同极化相位差χ(i,j)：

其中分别为第k相位图的HH极化与VV极化相位值；

S3.3、根据振幅离差指数D_A(i,j)和归一化同极化相位差χ(i,j)选取永久散射点；

在本实施例中，设置阈值σ₁＝0.2、σ₂＝0.8，ε₁＝0.2、ε₂＝0.25；

将归一化同极化相位差χ(i,j)＜σ₁，D_A(i,j)＜ε₁的像元点设置为永久散射点，以及将χ(i,j)＞σ₂，D_A(i,j)＜ε₂的像元点也设置为永久散射点；其中，χ(i,j)还满足χ(i,j)∈[0,1]；经过上述二次选取后，得到最终的永久散射点如图3所示，其中，永久散射点分布稀疏不均匀，数量为242个，约占整个影像的0.069％。

S4、利用双重自适应性滤波算法选取分布式散射点

S4.1、选取中心像元点

根据振幅图，计算每个像元点的同极化振幅幅比ρ(i,j)：

ρ (i, j) = \frac{1}{N} Σ_{k = 1}^{N} \frac{A_{k}^{H H} (i, j)}{A_{k}^{V V} (i, j)}

S4.2、利用统计学中的KS检验以中心像元点为中心、7×7滤波窗口内，判定中心像元点的同质点；

S4.3、确定同质点族：在滤波窗口内，舍去没有与中心像元点直接或间接相连的同质点，将中心像元点与其直接或间接相连的同质点作为一个同质点族,如图4(a)所示；

S4.4、判断同质点族大小是否超过阈值T＝20，如果同质点族大小超过阈值T＝20，则将该同质点归入一次候选分布式散射点集CDS1中；

S4.5、计算同质点族的振幅均值，并作为CDS1的振幅值，这样减小分布式散射点散射特性不稳定带来的误差；

S4.6、对CDS1重复步骤S4.2～S4.4，再次判定同质点族，当同质点族大小超过阈值T＝20时，归入二次候选分布式散射点集CDS2中；

在本实施例中，如图4所示，对比图4(a)与图4(b)，同一个像元点的同质点族可能不同，这是由于图4(b)中的像元点得到了优化，因而更接近真实情况；

S4.7、根据差分干涉相位图计算出CDS2同质点族的差分相位均值，并作为CDS2的差分相位值，这样可以减小系统噪声带来的误差；同时计算CDS2中各像元点的振幅离差指数，将振幅离差指数D_A(i,j)小于预设阈值D″_T的像元点设置为分布式散射点，如图5所示，选取的分布式散射点有650个，约占整幅影像的0.18％，为永久散射点的2.5倍；

S5、形变信息提取

将永久散射点与分布式散射点合并为稳定散射点集，建立Delaunay三角网。稳定散射点集构成的Delaunay三角网如图6所示，从图中可以看出稳定散射点的分布较为均匀，数量和密度上都满足了PS-InSAR技术的要求。

利用网络平差法对稳定散射点的差分相位进行相位分解，从而提取出山区地表微形变信息。图7与图8分别为影像区域的线性形变速率分布图和高程修正量分布图。线性形变速率代表着影像区域地形的变化趋势，而高程修正量可以用来修正高程数据，提高高程数据的精度。最后需要对线性形变速率图进行地理校正，如图9所示，可以根据地理位置的经纬度信息直接求出对应的线性形变速率。图10为研究区域2013年四月至九月地表形变累积量分布图，从图中可以直观的看出这五个月来影像区域的地表形量。

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

Claims

1.一种基于SAR图像的山区地表微形变信息提取方法，其特征在于，包括以下步骤：

(1)、对时间序列SAR影像组进行预处理

(2)、差分干涉

(3.1)、根据振幅图，计算每个像元点的振幅离差指数D_A(i,j)：

D_{A} (i, j) = \frac{σ_{A} (i, j)}{m_{A} (i, j)}

m_{A} (i, j) = \frac{1}{N} Σ_{k = 1}^{N} A_{k} (i, j), σ_{A} = \sqrt{\frac{1}{N} Σ_{k = 1}^{N} {(A_{k} (i, j) - m_{A} (i, j))}^{2}}

(3.2)、根据振幅图，计算每个像元点的归一化同极化相位差χ(i,j)：

其中分别为第k相位图的HH极化与VV极化相位值；

(4)、利用双重自适应性滤波算法选取分布式散射点

(4.1)、选取中心像元点

根据振幅图，计算每个像元点的同极化振幅比ρ(i,j)：

ρ (i, j) = \frac{1}{N} Σ_{k = 1}^{N} \frac{A_{k}^{H H} (i, j)}{A_{k}^{V V} (i, j)}

设定振幅离差指数阈值阈值D'_T、同极化振幅比阈值将振幅离差指数满足D_A(i,j)＜D'_T且同极化振幅比值满足的像元点设置为中心像元点；

(4.5)、计算同质点族的振幅均值，并作为CDS1的振幅值；

(4.7)、根据差分干涉相位图计算出CDS2同质点族的差分相位均值，并作为CDS2的差分相位值；同时计算CDS2中各像元点的振幅离差指数，将振幅离差指数D_A(i,j)小于预设阈值D'_T'的像元点设置为分布式散射点；

(5)、形变信息提取

2.根据权利要求1所述的基于SAR图像的山区地表微形变信息提取方法，其特征在于，所述步骤(3.3)中选取永久散射点的方法为：

设置阈值σ₁、σ₂，ε₁、ε₂；

将归一化同极化相位差χ(i,j)＜σ₁，D_A(i,j)＜ε₁的像元点设置为永久散射点，以及将χ(i,j)＞σ₂，D_A(i,j)＜ε₂的像元点也设置为永久散射点；其中，χ(i,j)可满足χ(i,j)∈[0,1]。

3.根据权利要求1所述的基于SAR图像的山区地表微形变信息提取方法，其特征在于，所述步骤(4.2)中滤波窗口的取值为7×7。