CN101881823B - InSAR区域网平差干涉参数定标与控制点加密方法 - Google Patents

Abstract

本发明公开了一种InSAR区域网平差干涉参数定标与控制点加密方法，它包括以下步骤：①获取各个干涉图的解缠结果；②根据获取的干涉图解缠结果，利用图像重叠区内的同名连接点和图像覆盖范围内的地面控制点，采用区域网平差方法解算轨道参数、基线参数和干涉相位偏置，并加密干涉处理用的控制点；③重建DEM和制作正射影像。本发明采用了区域网平差方法进行大面积、多套InSAR数据的干涉参数定标，可有效减少大面积多套InSAR数据干涉参数定标所需的地面控制点数量，提高干涉参数定标性能，降低不同干涉数据接边处反演高程值的差异，提高获取DEM及正射影像的整体一致性和连续性，并且能够加密出干涉处理用的大量地面控制点。

Description

InSAR区域网平差干涉参数定标与控制点加密方法

技术领域

本发明涉及InSAR干涉参数定标技术，特别涉及一种InSAR区域网平差干涉参数定标与控制点加密方法。

背景技术

InSAR(Interferometric Synthetic Aperture Radar，合成孔径雷达干涉测量)技术是随着信息技术、摄影测量技术、数字信号处理技术等相关技术的发展而迅速发展起来的一种高精度对地观测新技术。它在地形测绘、地表形变监测、冰川运动研究等方面都表现出快速、高精度、全天时、全天候等突出优势。InSAR技术已成为目前发展迅速、极具潜力的对地观测及测绘新技术之一。

在InSAR研究和应用方面，以往国内对InSAR技术的研究主要是针对单独干涉像对的干涉处理进行的，影响了我国InSAR技术的实际测图应用。

若采用处理单独干涉像对的方法来处理覆盖大面积的多套InSAR数据，将存在以下问题：(1)对于每对干涉数据，都需要足够数量、分布合理的地面控制点来进行轨道参数、基线参数和干涉相位偏置等参数的解算。因此，在整个大区域内完成所有干涉像对的干涉处理和DEM(Digital Elevation Model，数字高程模型)及正射影像的获取工作就需要获取大量、分布合理的地面控制点，从而需要消耗大量的人力和财力。在这种情况下，InSAR技术很难在境内无图区、稀疏控制区及境外地区的地形测绘中发挥作用。(2)对每个干涉像对单独进行干涉处理，由于没有考虑和利用数据之间的连接条件及相关关系，不同干涉数据获取的DEM及正射影像之间存在较大的相对误差，特别是在DEM接边处会出现明显的高程差异，在正射影像的接边处会出现明显的影像错位现象。

因此，为了有效减少大面积多套InSAR数据干涉处理所需的地面控制点数量，提高干涉参数定标性能，降低不同干涉数据接边处反演高程值的差异，提高获取DEM及正射影像的整体一致性和连续性，促进我国InSAR技术的发展和实际测图应用，本发明提出了InSAR区域网平差干涉参数定标与控制点加密方法。

发明内容

本发明的目的是为了提供一种InSAR区域网平差干涉参数定标与控制点加密方法，它与现有的单独像对InSAR数据干涉参数定标方法相比，具有所需地面控制点数量少、干涉参数定标性能文件、DEM接边处高程差异较小的突出优势。

为实现上述目的，本发明采用以下技术方案，它包括以下步骤：

①获取各个干涉图的解缠结果；

②根据获取的干涉图解缠结果，利用图像重叠区内的同名连接点和图像覆盖范围内的地面控制点，采用区域网平差方法解算轨道参数、基线参数和干涉相位偏置，并加密干涉处理用的控制点；

③重建DEM和制作正射影像。

在所述的步骤②中，采用整体解法或两类未知数交替趋近法进行区域网平差，并加密控制点。

在所述的步骤①中，获取各个干涉图的解缠结果包括如下步骤：

i)、对各个InSAR干涉像对进行高精度匹配；

ii)、干涉像对进行高精度匹配之后，对得到的数据进行重采样处理，并逐像元计算干涉相位，生成干涉图；

iii)、对得到的干涉图进行滤波处理；

iv)、计算干涉质量图；

v)、对干涉图进行相位解缠处理。

在所述的步骤i)中，采用影像金字塔实现InSAR干涉像对的粗匹配，然后采用相关系数法实现像元级匹配、最后采用最小二乘法实现子像元级匹配。

采用上述技术方案的本发明，采用了区域网平差方法进行大面积、多套InSAR数据的干涉参数定标，可有效减少大面积多套InSAR数据干涉参数定标所需的地面控制点数量，提高干涉参数定标性能，降低不同干涉数据接边处反演高程值的差异，提高获取DEM及正射影像的整体一致性和连续性，并且能够加密出干涉处理用的大量地面控制点。

附图说明

图1为本发明中InSAR区域网平差干涉参数定标与控制点加密流程图；

图2为本发明中InSAR区域网平差干涉参数定标的控制点分布示意图；

图3为本发明中以两套InSAR数据区域网平差干涉参数定标为例说明的控制点分布图；

图4为两套InSAR数据分别进行单独干涉参数定标时所需的控制点分布图；

图5为实验数据0001_04的控制点分布图；

图6为实验数据0001_03的控制点分布图；

图7为实验数据1001_04的控制点分布图；

图8为实验数据1001_03的控制点分布图；

图9为实验数据0001_04的连接点分布图；

图10为实验数据0001_03的连接点分布图；

图11为实验数据1001_04的连接点分布图；

图12为实验数据1001_03的连接点分布图；

图13为各套干涉数据分别进行单独干涉参数定标后连接点处的高差分布图；

图14为各套干涉数据进行区域网平差干涉参数定标后连接点处的高差分布图；

图15为对各干涉数据进行区域网平差干涉参数定标后获取的正射影像镶嵌结果图；

图16为对各干涉数据进行区域网平差干涉参数定标后获取的DEM镶嵌结果图。

具体实施方式

如图1所示，本发明的InSAR区域网平差干涉参数定标与控制点加密方法包括以下处理步骤：

①获取各个干涉图的解缠结果；

②根据获取的干涉图解缠结果，进行干涉参数定标与控制点加密。即利用各干涉数据重叠范围内的同名连接点和图像覆盖范围内地面控制点，采用区域网平差方法解算轨道参数、基线参数、干涉相位偏置等参数，并加密出干涉处理用的地面控制点；需要说明的是，上述的地面控制点包括平面高程控制点、高程控制点、平面控制点。

③重建DEM和制作正射影像。根据SAR影像的定向参数、相应地面点的高程值，采用直接法或间接法两种方法重建DEM并制作正射影像。本发明能够采用少量地面控制点对大面积多套InSAR数据进行干涉参数定标和控制点加密，其干涉参数定标性能稳健，获取的DEM及正射影像在接边处差异较小。

需要说明的是，上述获取各个干涉图解缠结果的方法有很多，本实施例提供以下一种，即包括：

i)、对各个InSAR干涉像对进行高精度匹配。该匹配的过程可人工进行匹配，也可以采用机器匹配。在利用机器进行匹配时，充分利用轨道参数、POS(Position and orientation system)数据等相关参数，采用影像金字塔匹配、相关系数匹配、最小二乘匹配等匹配策略实现各InSAR像对的自动、高精度匹配及不同像对之间的匹配，为干涉处理及连接点选取奠定基础。需要说明的是，上述的影像金字塔匹配、相关系数匹配、最小二乘匹配等匹配策略均为本领域普通技术人员所熟知的技术。

ii)、干涉图生成。将各干涉像对分别进行精确匹配之后，对包括强度和相位的相应数据进行重采样处理，并逐像元计算干涉相位，生成干涉图。

iii)、干涉图滤波。为了降低相位解缠的难度，提高相位解缠的可解性，采用空间域、频率域等滤波方法对干涉图进行滤波处理。上述的空间域、频率域等滤波方法也为本领域普通技术人员所熟知的技术。

iv)、干涉质量图计算。为了评价干涉质量，指导干涉图滤波及相位解缠，分别计算相干图、伪相干图等干涉质量图。上述计算干涉质量图也为本领域普通技术人员所熟知的现有技术。

v)、相位解缠。采用质量图区域生长、直接相位积分、最小二乘相位解缠、移动曲面相位拟合等策略对各干涉图进行相位解缠处理。上述的相位解缠过程也为本领域普通技术人员所熟知的技术。

在本发明中，步骤②InSAR区域网平差干涉参数定标和控制点加密为关键技术，下面对步骤②加以详细说明：

在步骤②的流程中，与传统单独像对干涉参数定标的区别主要体现在：充分利用不同像对之间的同名连接点，根据各像对中的平面控制点、高程控制点、平面高程控制点，采用区域网平差方法，同时整体解算出各干涉像对的基线参数、轨道参数、干涉相位偏置等参数，并加密出各连接点(像点)处对应的地面点坐标。

为了说明InSAR区域网平差中各类控制点及连接点的应用方案，图2给出了一种平面控制点、高程控制点、平面高程控制点和连接点均参与区域网平差解算的点位分布示意图。图2中，数据1、数据2、......数据6，均表示相应的InSAR干涉像对，平面控制点、高程控制点、平面高程控制点和连接点(加密点)都参与平差运算，可同时整体解算出相应各像对的干涉参数(包括基线参数、干涉相位偏置)及与所有连接点对应的地面点的地面坐标。

下面介绍基于区域网平差的干涉参数(包括基线参数和干涉相位偏置)定标和连接点处高程值加密的整体技术方案。

令R表示主天线相位中心到相应地面点的斜距，ΔR表示InSAR成像时两天线相位中心到同一地面点的斜距差，B表示两天线相位中心之间的基线长度，α表示基线与水平方向的夹角，H为主天线相位中心的高程，h为地面点的高程。则由InSAR的基本原理可得到如下关系式：

式(1)中，

对于机载双天线InSAR系统而言：

$\mathrm{\ΔR}=\frac{{\mathrm{\φ}}_0+\mathrm{\Δ\φ}}{2\mathrm{\π}}\mathrm{\λ};$

对于星载重复轨道InSAR系统而言：

$\mathrm{\ΔR}=\frac{{\mathrm{\φ}}_0+\mathrm{\Δ\φ}}{4\mathrm{\π}}\mathrm{\λ};$

其中，φ₀表示干涉相位偏置，Δφ表示解缠干涉相位，λ表示雷达波波长；

对式(1)进行变换，可得：

$F=B\sin (\mathrm{\θ}-\mathrm{\α})+\mathrm{\ΔR}-\frac{B^2}{2R}+\frac{{\mathrm{\ΔR}}^2}{2R}=0---\left(2\right)$

式(2)可简记为：

F(B，α，φ₀，h)＝0                         (3)

考虑到误差因素的影响，式(3)对基线长度、基线角、干涉相位偏置三个参数及高程值h线性化后对应的误差方程可表示为：

v＝F(B，α，φ₀，h)＝b₀ΔB+b₁Δα+b₂Δφ₀+b₃Δh-l      (4)

其系数分别为：

$\left\{\begin{array}{c}{b}_0=\frac{\∂F}{\∂B}=\sin (\arccos \frac{H-h}{R}-\mathrm{\α})-\frac{B}{R}\\ b_1=\frac{\∂F}{\∂\mathrm{\α}}=-B\cos (\arccos \frac{H-h}{R}-\mathrm{\α})\\ b_2=\frac{\∂F}{{\∂\mathrm{\φ}}_0}=\frac{\mathrm{\λ}}{2\mathrm{\π}}+\frac{\mathrm{\ΔR}}{R}\frac{\mathrm{\λ}}{2\mathrm{\π}}\\ b_3=\frac{\∂F}{\∂h}=\frac{1}{\sqrt{R^2-{(H-h)}^2}}B\cos (\arccos \frac{H-h}{R}-\mathrm{\α})\end{array}\right.$

常数项为：

$l=-F_0(B,\mathrm{\α},{\mathrm{\φ}}_0,h)=-B_0\sin (\arccos \frac{H-h}{R}-{\mathrm{\α}}_0)-\mathrm{\ΔR}+\frac{{B_0}^2}{2R}-\frac{{\mathrm{\ΔR}}^2}{2R}$

式(4)写成矩阵形式为：

$V=\left[\begin{array}{cc}A& B\end{array}\right]\left[\begin{array}{c}{\mathrm{\Δ}}_1\\ {\mathrm{\Δ}}_2\end{array}\right]-L---\left(5\right)$

其中：

V＝[v]；

A＝[b₀ b₁ b₂]；

B＝[b₃]；

Δ₁＝[ΔB Δα Δφ₀]^T；

Δ₂＝[Δh]；

L＝[l]；

为了简要说明基于区域网平差的干涉参数定标和控制点加密思想，我们假定采用的数据是具有一定影像重叠的两个干涉像对，所采用的控制点为平面高程控制点，需要解算的是两个干涉像对的干涉参数。如图3所示，若采用区域网平差干涉参数定标方法，则在这两个像对中，控制点、连接点都可根据式(4)列出相应的误差方程式。

以图3中的控制点1为例，所列出的误差方程形式为：

v_GCP1＝F(B₁，α₁，φ₁₀，h_GCP1)                     (6)

而对于图3中的连接点2，在不同的干涉像对中所列出的误差方程分别为：

v_1TP2＝F(B₁，α₁，φ₁₀，h_TP2)                      (7)

v_2TP2＝F(B₂，α₂，φ₂₀，h_TP2)                      (8)

依此类推，对于所有的控制点和连接点，都可列出相应的误差方程式；之后给定未知数的初值，根据误差方程组成法方程式进行各未知数改正量的答解；再根据解算出的未知数改正量对初值进行修正，并对上述计算过程进行迭代，直至满足给定的收敛条件。最后可得到各干涉参数的定标结果和所有连接点处的高程值加密结果。此时，需要定标的干涉参数共计2×3＝6个，重叠范围内的1个控制点可列1×2＝2个误差方程式，其余2个控制点可列2×1＝2个误差方程式，2个连接点可列2×2＝4个误差方程式并同时引入2个未知数，所以此时刚好能列2+2+4＝8个误差方程式来解算6+2＝8个未知数。此时，在整个影像覆盖范围内至少需要3个分布合理的地面控制点。并且当干涉像对数量增多时，也至少需要3个地面控制点。假定需要进行干涉参数定标的干涉像对数为100，且各数据间仅存在两度重叠，控制点均不位于影像重叠范围内，则需要定标的干涉参数数量为：n×3＝100×3＝300，由控制点可列出的误差方程式数量为：3×1＝3，由连接点列出的误差方程式数量为：m×2＝2m，由连接点引入的高程未知数数量为：m×1＝m，此时：300+m＝3+2m，m＝297。所以此时至少需要3个控制点和297个连接点就可以完成100个干涉像对的干涉参数定标，有效减少了所需的地面控制点数量。

但是如果采用各像对单独进行干涉参数定标方法，如图4所示，对于两套干涉数据，则要求每套数据中都至少布设分布合理的3个地面控制点，假定影像重叠区域内存在1个地面控制点，则总共需要至少5个地面控制点；航线数或者干涉像对数越大，所需的地面控制点数量就越大。当航线数或者干涉像对数为100时，假定各影像重叠范围内均存在且只存在1个地面控制点，则至少需要n×3-(n-1)×2+(n-1)＝2n+1＝2×100+1＝201个地面控制点；

综上所述，采用基于区域网平差的干涉参数定标和控制点加密方法可以有效减少对地面控制点数量的需求。

对于区域网平差干涉参数定标和控制点加密，可采取两种解算方案：

方案A：整体解法

对于每一个像点(分别对应于控制点和连接点)可以根据式(5)列出一组误差方程式，误差方程式中含有两类未知数Δ₁和Δ₂。其中，Δ₁对应于所有干涉像对的干涉参数(每个像对有3个干涉参数)，Δ₂对应于所有连接点的高程值。

相应的法方程式为：

$\left[\begin{array}{cc}{A}^{T}A& A^{T}B\\ B^{T}A& B^{T}B\end{array}\right]\left[\begin{array}{c}{\mathrm{\Δ}}_1\\ {\mathrm{\Δ}}_2\end{array}\right]=\left[\begin{array}{c}{A}^{T}L\\ B^{T}L\end{array}\right]---\left(9\right)$

对于大区域而言，当涉及的干涉像对数和连接点数很多时，误差方程式的总数十分可观。在解算过程中可先消去一类未知数而只求另一类未知数。考虑到一般选取的连接点数较多，其未知数Δ₂的个数远大于干涉参数Δ₁的个数，因此消去Δ₂以后，可得Δ₁的解为：

Δ₁＝[A^TA-A^TB(B^TB)^-1B^TA]^-1[A^TL-A^TB(B^TB)^-1B^TL]            (10)

方案B：两类未知数交替趋近法

首先把连接点处的高程近似值作为已知值，求出每个干涉像对的干涉参数，再利用干涉参数的新值计算各连接点的高程值，如此反复趋近，直至个干涉像对的干涉参数的改正值和连接点的高程改正值均小于某个限值时，迭代结束。

为验证本发明的效果，特作如下实验：

采用机载双天线InSAR系统获取的山东某地区多航带干涉数据进行了实验。在实验区内，包含了平地和山区等典型的地形要素。实验所采用的InSAR系统的部分相关参数如表1所示；

参数类型                  参数值

波(m)                     0.0312

波段                      X

方位向分辨率(m)           1.1

距离向分辨率(m)           1.25

绝对航高(m)               6190.0

多普勒中心频率(Hz)        0

极化方式                  HH

                表1

图5～图8分别为实验区两条航带的InSAR数据强度图。其中，图5和图6为航线0001中相邻两块数据的控制点分布情况；图7和图8为航线1001中相邻两块数据的控制点分布情况。表2为根据SAR强度影像中的明显特征点利用差分GPS到实地测量的高程值。

                   表2

图9～图12为各数据块中选取的连接点分布情况。表3为采用各航带单独定标方法对干涉相位偏置、基线长度和基线水平角进行定标的结果；表4为采用区域网平差方法对干涉参数进行定标的结果；从实验结果可以看出：对于干涉数据100103，由于覆盖范围内只有两个地面控制点，当采用干涉参数单独定标方法时，由于控制点数量不足，不能完成干涉参数定标工作；而采用区域网平差干涉参数定标与控制点加密方法，则可根据与其它数据区域的连接点条件，进行有效的干涉参数定标处理，并能进行高程控制点加密，进而有效地进行高程信息反演。

                       表3

                       表4

表5分别列出了采用单独定标和区域网平差定标的干涉参数反演连接点处高程值的差异情况。图13和图14分别表示单独定标和区域网平差定标后连接点处高差分布情况；从表5和图13、14可以看出，采用区域网平差干涉参数定标方法，可以有效减少接边处的高程差异。

                    表5

图15和图16为根据区域网平差定标结果得到的实验区的正射影像镶嵌结果和DEM镶嵌结果。

加密出的高程控制点坐标如表6所示。

点名    高程值(m)

1       52.302

2       50.460

3       55.817

4       47.016

5       56.465

6       47.312

7       47.994

8       49.311

9       54.844

10      45.100

11      49.946

12      67.35

13      59.239

14      54.533

15      56.983

16      56.501

17      56.041

18    45.415

19    59.606

20    52.042

21    49.541

22    51.295

23    49.580

24    55.928

25    59.064

26    62.578

27    59.840

28    60.737

30    63.912

31    63.026

33    56.933

34    63.877

35    50.700

36    66.838

37    95.817

38    98.639

68    56.317

    表6

Claims (4)

1.一种InSAR区域网平差干涉参数定标与控制点加密方法，其特征在于，它包括以下步骤：

①获取各个干涉图的解缠结果；

②根据获取的干涉图解缠结果，利用图像重叠区内的同名连接点和图像覆盖范围内的地面控制点，采用区域网平差方法解算轨道参数、基线参数和干涉相位偏置，并加密干涉处理用的地面控制点；

③重建DEM和制作正射影像。

2.根据权利要求1所述的InSAR区域网平差干涉参数定标与控制点加密方法，其特征在于：在所述的步骤②中，采用整体解法或两类未知数交替趋近法进行区域网平差，并加密控制点；

所述的整体解法是指；对应于每个控制点和连接点，列出一组误差方程式，误差方程式中含有两类未知数；对于大区域而言，当涉及的干涉像对数和连接点数很多时，在解算过程中先消去一类未知数而只求另一类未知数；

所述的两类未知数交替趋近法是指：首先把连接点处的高程近似值作为已知值，求出每个干涉像对的干涉参数，再利用干涉参数的新值计算各连接点的高程值，如此反复趋近，直至每个干涉像对的干涉参数的改正值和连接点的高程改正值均小于某个限值时，迭代结束。

3.根据权利要求1所述的InSAR区域网平差干涉参数定标与控制点加密方法，其特征在于：在所述的步骤①中，获取各个干涉图的解缠结果包括如下步骤：

i)、对各个InSAR干涉像对进行高精度匹配；

ii)、干涉像对进行高精度匹配之后，对得到的数据进行重采样处理，并逐像元计算干涉相位，生成干涉图；

iii)、对得到的干涉图进行滤波处理；

iv)、计算干涉质量图；

v)、对干涉图进行相位解缠处理。

4.根据权利要求3所述的InSAR区域网平差干涉参数定标与控制点加密方法，其特征在于：在所述的步骤i)中，采用影像金字塔实现InSAR干涉像对的粗匹配，然后采用相关系数法实现像元级匹配、最后采用最小二乘法实现子像元级匹配。

Images