CN106950556A - 基于分布式散射体时序干涉sar技术的遗产地形变监测方法 - Google Patents

Abstract

本发明公开了基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，首先确定遗产地的监测区域并收集相关资料；再对收集的SAR影像进行互配准、辐射校正和影像裁剪预处理，提取永久散射体目标和分布式散射体目标；去除分布式散射体目标的地形与平地相位，采用最大似然估计和LBFGS技术优化相位估计；再融合初始相位的永久散射体目标和经过优化相位的分布式目标，设置短时间/空间基线阈值生成小基线集干涉对，对所述小基线集干涉对进行相位解缠并去除地形、大气延迟造成的误差；最后求解平均年形变速率及形变历史变化值。本发明监测遗产地形变，克服郊区可提取点稀少问题，适用于从单一古建筑到整个遗产地多尺度稳定监测，利于遗产地影响因素的评估。

Description

基于分布式散射体时序干涉SAR技术的遗产地形变监测方法

技术领域

本发明涉及遗产地监测领域，特别是涉及基于分布式散射体时序干涉SAR技术的遗产地形变监测方法。

背景技术

遗产地，如文物古迹、自然奇观，是人类历史上具有重要意义和巨大价值的无以替代的财富。5000余年的历史文明以及辽阔的地域使得中国成为遗产大国，截止2015年7月列入《世界遗产名录》的有48项，数量位居世界第二，仅次于意大利。不幸的是，遗产地在自然灾害和人类活动的侵蚀下正在变得日益不稳定。自然灾害包括滑坡、地震、洪水、恶劣的天气以及气候突变。人类活动包括战争、资源的过度开采、城镇化建设和失控的旅游业。这些外界因素使得很多遗产地岌岌可危，监测并实施有效的保护势在必行。

遗产地不稳定的程度可以用形变来衡量，大致可以分为两类，一是长期的缓慢形变，尺度在毫米级或厘米级，主要由自然的侵蚀、地壳运动、人类活动引起的地表沉陷等引起；二是突发的剧烈形变，尺度在厘米级以上，主要由地震、滑坡等自然灾害以及人类破坏等引起。传统的形变测量方法如水准、GPS或物探技术均存在着以下不足：1)效率低、成本高、需要耗费大量的人力物力，受野外作业环境的严重限制；2)在空间维度上无法获得大面积形变场，时间维度上难以获取长时间序列形变，时空变化规律难以展现。

现有的时间序列差分干涉SAR技术(Multi-Temporal differentialInterferometric Synthetic Aperture Radar,MTInSAR)联合多基线重轨SAR影像测量地表微小形变，克服了传统方法的局限性，使得快速、大面积、时序监测成为了可能，而且理论上精度可以达到毫米级，因此更加适合于遗产地的监测。常用的MTInSAR技术包括：永久散射体技术(Persistent Scatterers,PSInSAR),小基线集技术(Small BAselines Subset,SBAS)，Stamps以及SqueeSAR技术。

但是现有的技术在遗产地监测中仍存在一些亟待解决的问题：1)对SAR影像数量的要求较高(通常要求为25景以上，如PSInSAR和SqueeSAR)；2)郊区干涉对的时间失相干现象严重，永久散射体(Persistent Scatterers,PS)稀少，可用于参数反演的点数量不足，影响相位解缠和形变估计的精度。SBAS方法虽然利用短时间/空间基线的干涉对来减少失相干的影响，以挖掘短时间内能够保持相干的像元信息来增加监测点的密度，但是多视操作降低了空间分辨率，相干性估计也容易出现偏差。

MTInSAR技术目前在郊区的研究主要集中在分布式散射体(DistributedScatterers，DS)上，它指的是占据一定面积均匀地物，在SAR影像上对应的像元呈现出相似的后向散射特征，裸地和低矮的植被都属于DS点。引入DS点能够显著提升遗产地可监测目标的数量，改进形变反演的结果。但是分布式散射体相干性差、信噪比低，需要进行复杂特殊的滤波处理。

由此可见，现有的MTInSAR技术在遗产地形变监测中仍然存在着缺陷，如何能创设一种更加适用的基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，成为业界追求的目标。

发明内容

本发明要解决的技术问题是提供一种基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，使其实现大范围、长时间、高精度的监测遗产地不稳定现象，从而克服现有的遗产地形变监测方法的不足。

为解决上述技术问题，本发明提供1、基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，包括如下步骤：

(一)确定遗产地的监测区域及资料收集

(11)确定遗产地的监测范围，包括遗产地和所述遗产地周边的缓冲带；

(12)调查所述遗产地的背景知识，并收集所述遗产地的时间序列SAR单视复数影像；

(二)SAR数据处理与参数反演

(21)对所述步骤(12)收集的时间序列SAR单视复数影像进行互配准、辐射校正和影像裁剪预处理；

(22)提取永久散射体目标和分布式散射体目标；

(23)去除分布式散射体目标的地形与平地相位，采用最大似然估计和LBFGS技术进行优化相位估计；

(24)融合初始相位的永久散射体目标和经过优化相位的分布式散射体目标，设置短时间/空间基线阈值生成小基线集干涉对，对所述小基线集干涉对进行相位解缠并去除地形、大气延迟造成的误差；

(25)求解平均年形变速率及形变历史变化值。

作为本发明的一种改进，所述方法还包括结果分析与解译步骤，具体为：利用野外实地调查测量结果对步骤(25)所得结果进行精度评定；并结合野外实地调查测量结果与步骤(25)所得结果共同分析所述遗产地的形变特征及引发不稳定的因素。

进一步改进，所述步骤(12)中所述遗产地的背景知识包括所处的环境、遭受的威胁、陆地景观、地表特征、空间尺度、遗产单体的几何形状和材料属性。

进一步改进，所述步骤(22)中提取永久散射体目标的具体方法为：根据单幅影像的光谱特征与多幅影像的振幅离差阈值法提取。

进一步改进，所述步骤(22)中提取分布式散射体目标的具体方法包括：

A、对所述时间序列SAR单视复数影像中像元的时序幅度值进行排序，转换得到像元幅度值累积分布函数，其像元幅度值累积分布函数的无偏估计表示为:

其中，F_N(X)为像元幅度值累积分布函数的无偏估计，X代表像元幅度值，N为配准后的SAR影像数量，x₁为时序幅度值排序中的第1个元素的幅度值，x_k为时序幅度值排序中的第k个元素的幅度值；x_k+1为时序幅度值排序中的第k+1个元素的幅度值；k取值范围为0…N-1。

B、再用KS检验或AD检验方法对相邻像元的无偏估计值进行判定，满足一定阈值的即为同质像元；

C、对于每个像元可获得其估计窗口内连通的所有同质像元，同质数量大于设定阈值的像元即为分布式散射体目标。

进一步改进，所述AD检验方法定义为：

其中，p,q代表两个像元，为像元p的累积分布函数的无偏估计值，为像元q的累积分布函数的无偏估计值，是为两点的经验累积分布函数的无偏估计值，N为配准后的SAR影像数量,x为像元的时序幅度值，A_N,N满足一定阈值，p,q即属于同质像元。

进一步改进，所述步骤(24)中对所述小基线集干涉对的相位解缠方法采用时空3D相位解缠方法。

进一步改进，针对第一代SAR数据进行轨道、地形相关误差去除的模型如下：

其中，x和y为距离向和方位向的坐标，h为高程；ε为随机相位误差；a_i为待估计参数。

进一步改进，所述步骤(25)中平均年形变速率是针对差分干涉对的相位模型采用最小二乘法求得；所述形变历史变化值是利用SVD算法求得，具体算法为：

假设第i幅干涉图，主辅影像获取时间分别为t_A和t_B，且t_B＞t_A，像元x处的差分干涉相位表示为：

式中φ(t_A,x)和φ(t_A,x)为在t_A和t_B时刻影像上的形变相位值，d(t_A,x)和d(t_B,x)是相对于参考时间t₀的雷达视线向形变量，参考时间形变量d(t₀,x)≡0，λ为雷达波长，用φ表示N个时刻SAR图像上的形变相位矩阵，Δφ表示M个差分干涉图上的相位矩阵，其矩阵形式表示为：

Δφ＝Aφ

其中，系数矩阵A[M×N]的每行对应一幅干涉图，每列对应一个时间上的SAR影像，用SVD算法求解得每个时间SAR影像上形变对应的相位，经转换可得每个时间上的形变值。

采用这样的设计后，本发明至少具有以下优点：

本发明基于分布式散射体的时序干涉SAR技术用于遗产地形变监测，使其能利用较少的SAR影像数据实现更大区域、更多地物的形变测量，从而推进遗产地的稳定监测与保护。

本发明属于一种适用于大范围、长时间、高精度监测遗产地不稳定现象的方法，特别适合于建筑物稀少的遗郊区。

附图说明

上述仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，以下结合附图与具体实施方式对本发明作进一步的详细说明。

图1是本发明基于分布式散射体时序干涉SAR技术的遗产地形变监测方法示意图。

具体实施方式

基于遗产地的监测需求和地表类型，本发明融合了SBAS和SqueeSAR算法，识别、提取DS点和PS点并进行联合分析，采用短时间/空间基线生成干涉对集合，最后去除地形、大气等误差生成形变结果，具体方法如下。

参照附图1所示，本发明基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，包括如下步骤：

一、确定遗产地的监测区域及资料收集

1.1、确定监测范围，应包括遗产地以及遗产地周边一定范围的缓冲带，该缓冲带有助于了解遗产地所处的环境，解释遗产地发生不稳定的原因和分析其发展规律。

1.2、调查确定遗产地的背景先验知识，包括其所处的环境、遭受的威胁(如滑坡等自然灾害或地面沉降等)，陆地景观、地表特征、空间尺度、遗产单体的几何形状和材料属性等数据。这些数据决定了要采取的监测方法和数据处理流程。

1.3、收集与该遗产地相关的时间序列SAR单视复数影像；

二、SAR数据处理与参数反演

2.1、对收集的时间序列SAR单视复数影像进行数据预处理，包括辐射校正、互配准、影像裁剪和主影像选择等。

2.2、提取永久散射体目标和分布式散射体目标；

根据SAR影像的后向散射特征和研究区域的地物类型提取永久散射体目标(PS点)和分布式散射体目标(DS点)。

由于PS点在SAR影像上表现为十分稳定的强散射，能够保持高相干性，对应的地物为遗产区的房屋、桥梁、道路以及裸露的岩石等。PS点具体的提取方法为：根据单幅影像的光谱特征与多幅影像的振幅离差阈值法提取。

DS点对应着分辨单元内所有散射体的后向散射系数大致相同的地物，多为郊区普遍存在的裸地和稀疏植被等，它们仅在部分干涉图保持一定的相干特性。

DS点的提取则需要对SAR影像上像元的后向散射特征进行统计分析，具体方法为：

1)对像元的时序幅度值进行排序，转换得到其累积分布函数；

2)用Kolmogorov-Smirnov(KS)检验或Anderson-Darling(AD)检验方法对相邻像元的无偏估计值进行判定，满足一定阈值的即为同质像元；

3)对于每个像元可获得其估计窗口内连通的所有同质像元，同质数量大于设定阈值的像元即为DS点。

其中，像元幅度值累积分布函数的无偏估计可表示为:

其中，F_N(X)为像元幅度值累积分布函数的无偏估计，X代表像元幅度值，N为配准后的SAR影像数量，x₁为时序幅度值排序中的第1个元素的幅度值，x_k为时序幅度值排序中的第k个元素的幅度值；x_k+1为时序幅度值排序中的第k+1个元素的幅度值；k取值范围为0…N-1。

无参数双尾AD检验方法可定义为：

p,q代表两个像元，为像元p的累积分布函数的无偏估计值，为像元q的累积分布函数的无偏估计值，是为两点的经验累积分布函数的无偏估计值，N为配准后的SAR影像数量,x为像元的时序幅度值，A_N,N满足一定阈值，p,q即属于同质像元。

2.3、去除分布式像元的地形与平地相位，然后进行时空三维的优化滤波；

因为DS点相干性差、信噪比低，需要进行时空滤波。首先像元p基于同质像元集为Ω的相干系数矩阵Γ(P)可改进为：

相干系数矩阵第k行、j列的元素γ_kj表示为第k景和第j景SAR影像组成干涉对的相干系数。主对角线为自身影像干涉且值为1。

基于相干系数矩阵和最大似然估计准则，假设第一幅SAR影像相位为0，N个相位值为最大似然估计值为：

其中，Λ＝exp(iφ)为N维向量，为复相干性矩阵估计值，ο为Hadamard乘积。采用LBFGS(Limited memory Broyden-Fletcher-Goldfarb-Shanno)算法求解，然后用解算后的相位取代初始影像的相位，即实现了优化估计。

2.4、融合永久散射体目标和分布式目标；

对PS点初始相位和DS点的优化相位进行差分干涉对重新组合。

2.5、设置短时间/空间基线阈值生成小基线集干涉对；

好的相干性是差分干涉技术测量形变的前提，而引起干涉对失相干最严重的因素是时间基线和空间基线。以单一主影像形成的干涉对集合因为时间/空间基线大、而且对影像数据量要求多，难以满足遗产地监测需求。以短时间/空间基线为标准形成的多基线干涉对能够克服失相干严重的问题，且能生成大量的干涉对降低对数据量的需求。

假设有N景影像，则小基线干涉对数量M满足：

差分干涉对的相位模型可表示为：

其中，a为常数项，Ρ_ξ、Ρ_η是大气效应和基线引起的干涉相位分别沿方位向ξ和距离向η的线性变化系数，Δh为高程误差，Δt为辅影像到主影像的时间基线，v表示线性形变速率，ε为残差，包含了大气残余、非线性形变等。Β_⊥为干涉对的垂直基线，λ为雷达波长，R为卫星传感器到目标的斜距，θ为主图像的入射角。其中ξ、η、Β_⊥、λ、R、θ、Δt为已知量，a、Ρ_ξ、Ρ_η、Δh、ν、ε为未知量。

2.6、相位解缠并去除地形、大气延迟等误差；

缠绕干涉图的相位值范围是(-π,π)，解决2π整数倍问题的过程称为相位解缠。这个过程是至关重要的，点密度过稀、地形起伏和形变较大等都能引起解缠错误。基于上述差分干涉对的相位模型，对于大量的小基线集干涉对可采用时空3D相位解缠方法，以对干涉对进行相位解缠并去除地形、大气延迟等误差。

第一代卫星影像的形变相位容易与轨道误差、大气延迟以及残余地形相位等混淆；长波段卫星的轨道误差能够达到30cm，在影像上呈现为二阶模型的变化趋势，同时轨道误差能引起地形估计的偏差，大气的湍流效应也与地形线性相关。

针对第一代SAR数据进行轨道、地形相关误差去除的模型如下：

其中，x和y为距离向和方位向的坐标，h为高程；ε为随机相位误差；a_i为待估计参数。

2.7、求解平均年形变速率及形变历史变化值；

对差分干涉对的相位模型用最小二乘法求解可得到平均年形变速率。采用奇异值去分解(Singular Value Decomposition SVD)的方法解算去过各种误差的解缠相位可获得形变的历史变化值。假设第i幅干涉图，主辅影像获取时间分别为t_A和t_B(t_B＞t_A)，像元x处的差分干涉相位可表示为：

式中φ(t_A,x)和φ(t_A,x)为在t_A和t_B时刻影像上的形变相位值，d(t_A,x)和d(t_B,x)是相对于参考时间t₀的雷达视线向(LOS)形变量，参考时间形变量d(t₀,x)≡0，λ为雷达波长。用φ表示N个时刻SAR图像上的形变相位矩阵，Δφ表示M个差分干涉图上的相位矩阵，则矩阵形式可表示为：

Δφ＝Aφ

其中，系数矩阵A[M×N]的每行对应一幅干涉图，每列对应一个时间上的SAR影像。用SVD算法求解可得每个时间SAR影像上形变对应的相位，然后可转换为每个时间上的形变值。

上述SAR数据处理融合了SBAS和SqueeSAR算法的优势，一方面联合提取与分析DS点和PS点，显著提高监测点的数量，扩大监测的范围，改善相位解缠与大气延迟去除的效果；另一方面采用小基线干涉对集合进行形变估计，在保留原始分辨率的同时降低对SAR影像数量的需求。

三、对干涉结果进行精度评定、解译以及分析遗产地的不稳定机制

3.1、野外调查测量与干涉结果精度评定；

野外调查测量是必不可少的工作，一方面是要确定监测点对应的地物，如一些树林中的废墟、遗址等，在普通光学影像上无法看见但能够被高分辨率SAR监测到的地物。另一方面是对形变结果进行实地验证，如发现古建筑物出现裂缝，地表明显塌陷等现象。条件允许，可获取实地测量数据(如水准，GPS等)对干涉形变结果进行绝对精度评定。

3.2、分析遗产地的形变特征及引发不稳定的因素与机制。

分析该遗产地的形变特征及不稳定的因素，需要结合遗产地的背景环境资料，如地质灾害、地下水开发或城市建设等进行，最终得出该遗产地的形变特征和不稳定因素。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰，均落在本发明的保护范围内。

Claims (9)

1.基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，其特征在于，包括如下步骤：

(一)确定遗产地的监测区域及资料收集

(11)确定遗产地的监测范围，包括遗产地和所述遗产地周边的缓冲带；

(12)调查所述遗产地的背景知识，并收集所述遗产地的时间序列SAR单视复数影像；

(二)SAR数据处理与参数反演

(21)对所述步骤(12)收集的时间序列SAR单视复数影像进行互配准、辐射校正和影像裁剪预处理；

(22)提取永久散射体目标和分布式散射体目标；

(23)去除分布式散射体目标的地形与平地相位，采用最大似然估计和LBFGS技术进行优化相位估计；

(24)融合初始相位的永久散射体目标和经过优化相位的分布式散射体目标，设置短时间/空间基线阈值生成小基线集干涉对，对所述小基线集干涉对进行相位解缠并去除地形、大气延迟造成的误差；

(25)求解平均年形变速率及形变历史变化值。

2.根据权利要求1所述的基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，其特征在于，所述方法还包括结果分析与解译步骤，具体为：利用野外实地调查测量结果对步骤(25)所得结果进行精度评定；并结合野外实地调查测量结果与步骤(25)所得结果共同分析所述遗产地的形变特征及引发不稳定的因素。

3.根据权利要求1所述的基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，其特征在于，所述步骤(12)中所述遗产地的背景知识包括所处的环境、遭受的威胁、陆地景观、地表特征、空间尺度、遗产单体的几何形状和材料属性。

4.根据权利要求1所述的基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，其特征在于，所述步骤(22)中提取永久散射体目标的具体方法为：根据单幅影像的光谱特征与多幅影像的振幅离差阈值法提取。

5.根据权利要求1所述的基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，其特征在于，所述步骤(22)中提取分布式散射体目标的具体方法包括：

A、对所述时间序列SAR单视复数影像中像元的时序幅度值进行排序，转换得到像元幅度值累积分布函数，其像元幅度值累积分布函数的无偏估计表示为:

$F_N\left(X\right)=\left\{\begin{array}{ccc}0,& if& X<x_1\\ \frac{k}{N},& if& x_k\&le;X<x_{k+1}\\ 1,& if& X\&GreaterEqual;x_N\end{array}\right.$

其中，F_N(X)为像元幅度值累积分布函数的无偏估计，X代表像元幅度值，N为配准后的SAR影像数量，x₁为时序幅度值排序中的第1个元素的幅度值，x_k为时序幅度值排序中的第k个元素的幅度值；x_k+1为时序幅度值排序中的第k+1个元素的幅度值；k取值范围为0…N-1。

B、再用KS检验或AD检验方法对相邻像元的无偏估计值进行判定，满足一定阈值的即为同质像元；

C、对于每个像元可获得其估计窗口内连通的所有同质像元，同质数量大于设定阈值的像元即为分布式散射体目标。

6.根据权利要求5所述的基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，其特征在于，所述AD检验方法定义为：

$A_{N,N}^2=\frac{N}{2}\underset{x\&Element;\{x_{p,i},x_{q,i}\}}{\&Sigma;}\frac{{\left({\hat{F}}_p\right(x)-{\hat{F}}_q(x\left)\right)}^2}{{\hat{F}}_{pq}\left(x\right)(1-{\hat{F}}_{pq}(x\left)\right)}$

其中，p,q代表两个像元，为像元p的累积分布函数的无偏估计值，为像元q的累积分布函数的无偏估计值，是为两点的经验累积分布函数的无偏估计值，N为配准后的SAR影像数量,x为像元的时序幅度值，A_N,N满足一定阈值，p,q即属于同质像元。

7.根据权利要求1所述的基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，其特征在于，所述步骤(24)中对所述小基线集干涉对的相位解缠方法采用时空3D相位解缠方法。

8.根据权利要求7所述的基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，其特征在于，针对第一代SAR数据进行轨道、地形相关误差去除的模型如下：

其中，x和y为距离向和方位向的坐标，h为高程；ε为随机相位误差；a_i为待估计参数。

9.根据权利要求1所述的基于分布式散射体时序干涉SAR技术的遗产地形变监测方法，其特征在于，所述步骤(25)中平均年形变速率是针对差分干涉对的相位模型采用最小二乘法求得；所述形变历史变化值是利用SVD算法求得，具体算法为：

假设第i幅干涉图，主辅影像获取时间分别为t_A和t_B，且t_B＞t_A，像元x处的差分干涉相位表示为：

${\mathrm{\&Delta;\&phi;}}_i\left(x\right)=\mathrm{\&phi;}(t_B,x)-\mathrm{\&phi;}(t_A,x)\&ap;\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&lsqb;d(t_B,x)-d(t_A,x)\&rsqb;$

式中φ(t_A,x)和φ(t_A,x)为在t_A和t_B时刻影像上的形变相位值，d(t_A,x)和d(t_B,x)是相对于参考时间t₀的雷达视线向形变量，参考时间形变量d(t₀,x)≡0，λ为雷达波长，用φ表示N个时刻SAR图像上的形变相位矩阵，Δφ表示M个差分干涉图上的相位矩阵，其矩阵形式表示为：

Δφ＝Aφ

其中，系数矩阵A[M×N]的每行对应一幅干涉图，每列对应一个时间上的SAR影像，用SVD算法求解得每个时间SAR影像上形变对应的相位，经转换可得每个时间上的形变值。