CN109946682B

CN109946682B - 基于ICESat/GLAS的GF3数据基线估计方法

Info

Publication number: CN109946682B
Application number: CN201910267939.7A
Authority: CN
Inventors: 贾蓓蕾; 王志斌; 李真芳; 鲁焕; 张庆君
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2019-04-03
Filing date: 2019-04-03
Publication date: 2022-12-02
Anticipated expiration: 2039-04-03
Also published as: CN109946682A

Abstract

本发明公开了一种基于ICESat/GLAS的GF3数据基线估计方法，主要解决现有技术在实时轨道快速生成DEM中存在精度差的问题。其实现方案是：用ICESat/GLAS本身参数及外部DEM库高程信息对ICESat/GLAS高程数据进行筛选，并获取斜平面上TanDEM‑DEM库的高程值，将该高程值与SAR图像对干涉得到的斜平面高程作差，得到斜平面上的高程误差；用模糊高度、高程误差及基线误差间的关系，获得场景中基线误差；用基线误差对实时轨道误差进行校正。本发明克服了地形因素对基线估计的限制，提高了基线估计的精度，减小了实时轨道进行DEM重建时精度损失，可用于对InSAR信号的处理。

Description

基于ICESat/GLAS的GF3数据基线估计方法

技术领域

本发明属于信号处理技术领域，特别涉及高分三号GF3数据的基线估计方法，可用于对干涉合成孔径雷达InSAR图像信号的处理。

背景技术

数字高程模型DEM是干涉合成孔径雷达技术的主要产物，是干涉相位经过高程转换和地理编码得到的均匀分布在格网上的高程点，能够展示实际地形的起伏变化状态，在军事、经济方面都具有广阔的应用前景。基线是InSAR中的一个重要参数。它决定了有效干涉像对的选择，并对InSAR数据处理中的平地效应、多视平均等有很大影响，其精度还影响到最终提取的地表高程或形变的精度。研究高精度的干涉SAR基线估计方法对获取场景目标高质量的DEM至关重要，也是InSAR信号处理技术的一个研究热点。

基线，被定义为照射同一区域的两部天线相位中心之间的距离。双天线单航过或单天线双航过工作模式下的雷达平台扫过地面某一场景区域时，会对这一场景区域形成两幅SAR图像。由于星载平台的限制，单个平台难以同时安装两个天线。通常，单天线双航过模式用于获得星载干涉雷达的干涉图像。然而在这种模式下难以保证准确的基线参数，因为地面条件和散射特征可能由于两个图像之间的时间间隔而改变。要想获得高精度DEM，必须对基线进行准确计算和校正。为了满足InSAR生成DEM的要求，一般情况下基线误差要在毫米量级。

现有的基线估计方法主要分为三大类：卫星轨道状态矢量法；快速傅里叶变换FFT法；外部控制点法。卫星轨道状态矢量法估计基线精度较差，受轨道影响较大；FFT法对地形有一定的限制，不适用于地震等自然灾害的特殊场景；外部控制点法精度高，计算量大，且不能保证所有的数据都能获取控制点信息。

对InSAR时序信号的处理而言，基线估计通常是利用估计的基线误差量，对基线和轨道进行校正。常用的方法是卫星轨道状态矢量法和FFT法。其中基于卫星轨道状态矢量的基线估计方法中，是直接利用轨道的位置和速度矢量信息进行基线估计，该基线估计方法存在如下缺点：当两幅干涉SAR图像的轨道是实时轨道时，基线估计精度很差。在基于FFT的基线估计方法中，利用干涉图中蕴含的相位信息如干涉条纹、干涉相位差等来估算基线分量。该方法是一种盲估计方法，仅适用于平坦地区,且利用了轨道参数进行计算,因此估计出的基线精度也不高。

发明内容

本发明的目的在于针对上述现有技术的不足，提出一种基于ICESat/GLAS的GF3数据基线估计方法，以克服地形限制，减小轨道影响，进一步提高基线估计精度。

为实现上述目的，本发明的技术方案包括如下：

1.一种基于ICESat/GLAS的GF3数据基线估计方法，其特征在于，包括如下步骤：

(1)对ICESat卫星星载地学激光测高系统ICESat/GLAS的激光测高数据产品GLA14进行筛选，提取其在设定经纬度范围内的高程值；

(2)利用选取的高程数据，校正TanDEM-X卫星数字高程模型DEM库TanDEM-DEM的高程数据；

(3)通过卫星轨道以及参数文件中提供的参数信息，利用反定位的方法获得校正后的TanDEM-DEM在斜平面上高程值h₁；

(4)对两幅SAR图像进行干涉处理得到解缠相位信息，并将其转换为干涉图像在斜平面的高程值h₂；

(5)将TanDEM-DEM在斜平面上高程值h₁与两幅SAR图像干涉处理得到的斜平面高程值h₂作差，得到斜平面上的高程误差Δh；

(6)根据参数文件中提供的波长λ和斜距R参数，求出SAR图像观测场景中的模糊高度h_amb；

(7)利用上述获得的高程误差Δh、模糊高度h_amb，根据基线误差之间的关系获得场景中随着方位向线性变化的平行基线分量误差

(8)根据平行基线分量误差ΔB_||，对场景方位向的每一个位置都建立坐标系进行实时轨道误差校正，得到基线误差校正后的轨道参数。

本发明与现有的技术相比具有以下优点：

第一，本发明通过计算干涉SAR图像对的斜平面高程值与TanDEM-DEM高程值之间的高程误差以及模糊高度，获得基线误差，实现轨道校正，由于没有直接采用卫星轨道信息，因此在非精密轨道的条件下，估计的基线参数受轨道影响较小，更加准确。

第二，由于本发明估计的基线误差是一个线性变化量，相较于传统方法估计的基线误差是个常数，提高了基线估计的精度。

第三，本发明中由于利用精度较高的ICESat/GLAS和TanDEM-DEM高程库信息计算高程误差，因而在没有外部控制点的情况下，也能估计出较高的基线参数。

附图说明

图1为本发明的实现总流程图；

图2为本发明中进行ICESat/GLAS高程数据筛选的子流程图；

图3为用本发明方法与现有的轨道状态矢量法、FFT法的基线误差对比图；

图4为未进行基线误差校正直接生成的高程值与TanDEM-DEM高程值的对比图；

图5为利用本发明法进行基线误差校正后生成的高程值与TanDEM-DEM高程值的对比图。

具体实施方式

下面结合附图对本发明实施例及效果做进一步详细说明：

参照图1，本实例的实现步骤如下：

步骤1，对ICESat卫星星载地学激光测高系统ICESat/GLAS的激光测高数据产品GLA14进行筛选，提取其在设定经纬度范围内的高程值。

参照图2，本步骤的具体实现如下：

(1a)初步筛选：

从设定经纬度范围内的GLA14标准数据中筛选高程数据，这些高程数据同时满足地表高程可用性标志为“0”、饱和度校正标志为“0”或“1”、姿态质量指标参数为“00”、反射率参数“<1”、云量描述参数为“15”、大气散射增益参数“<100”这六个条件；

(1b)第二次筛选：

将(1a)中从ICESat/GLAS中得到的高程数据与在相同经纬度上TanDEM-X卫星数字高程模型DEM库TanDEM-DEM的高程数据作差，剔除误差小于10m的高程数据，最终得到筛选后的ICESat/GLAS高程数据；

(1c)坐标转换：

将(1b)筛选出的ICESat/GLAS高程数据在海洋环流测量卫星TOPEX/Poseidon椭球坐标系下的坐标值，转换为在世界大地测量系统1984WGS84椭球坐标系下的坐标值，转换公式如下：

H_G＝H+dh

其中H为坐标转换之前的高程值，H_G为转换后的高程值，dh＝-cos²Bda-sin²Bdb为高程变换修正值，在dh中，B代表纬度，da、db分别表示WGS84椭球和TOPEX/Poseidon椭球长轴和短轴之间的差异值，转换误差在mm量级。

步骤2，利用选取的高程数据，校正TanDEM-X卫星数字高程模型DEM库TanDEM-DEM的高程数据。

(2a)绘制相同经纬度范围内ICESat/GLAS和TanDEM-DEM的高程数据，统计二者的均方根误差和标准差；

(2b)根据统计出的均方根误差和标准差，对所选经纬度范围内的TanDEM-DEM高程数据补加标准差，实现校正，经过校正之后，可以获取精度更高的TanDEM-DEM高程库信息，这是对基线误差校正建模的基础，同时，受轨道误差影响较小，使得估计的基线更加准确。

步骤3，通过卫星轨道以及参数文件中提供的参数信息，利用反定位的方法获得校正后的TanDEM-DEM在斜平面上高程值h₁。

(3a)将TanDEM-DEM高程数据在地平面纬度、地平面经度和高程三维坐标系中的坐标值，按如下关系转换为在WGS84坐标系中的坐标值：

其中X,Y,Z分别代表WGS84坐标系中的坐标，B,L,H分别代表地平面纬度、地平面经度和高程三维坐标系的坐标，N＝a(1-e²sin²B)^-1/2是曲率半径，a是N中的半长轴，e²表示第一偏心率；

(3b)通过聚焦多普勒中心得到TanDEM-DEM高程数据在斜平面的方位向位置，通过高分三号卫星数据的参数文件中所提供的距离向采样间隔和斜距信息，得到TanDEM-DEM高程数据在斜平面上的距离向位置，根据已确定的方位向位置和距离向位置，获得TanDEM-DEM高程数据在斜平面上的高程矩阵；

(3c)对高程矩阵中不准确的点进行插值，获得完整的高程值h₁。

步骤4，对两幅SAR图像进行干涉处理得到解缠相位信息，并将其转换为干涉图像在斜平面的高程值h₂。

(4a)通过图像配准、相位生成、相位滤波以及相位解缠处理，得到解缠相位

(4b)通过如下转换关系式，将解缠相位

转换为主、辅星天线的斜距差ΔR：

其中λ为波长；

(4c)构建定位方程：

根据主、辅星天线到地面的斜距关系以及主天线的多普勒方程，获得以下方程组：

其中，R₁为主星天线到地面的斜距，R₂为辅星天线到地面的斜距，且R₂＝R₁+ΔR，f_dc1为主星成像多普勒中心频率，

为主星天线相位中心矢量，

为辅星天线相位中心矢量，

表示地面目标矢量，λ为波长，

为主星速度矢量，

为地面目标点的速度矢量；

将方程组中

这五个矢量以坐标形式表示，且考虑到地面目标点为静止目标点，即

将上述三个方程表述为如下定位方程形式：

其中，(x₁,y₁,z₁)为主星天线相位中心的坐标，(x₂,y₂,z₂)为辅星天线相位中心的坐标，(x,y,z)为地面目标的坐标，(V_x,V_y,V_z)为主星天线速度；

(4d)利用牛顿迭代法对(4c)得到的定位方程进行求解：

(4d1)设非线性方程组

令x^(k)表示非线性方程组的一个近似根；

对非线性方程组在x^(k)处，按照多元函数的泰勒公式展开，并取线性相关得到如下牛顿迭代公式：

其中：

(4d2)将(4c)得到的定位方程代入上述牛顿迭代公式，解方程组，得到地面目标的坐标(x,y,z)，其中，z为目标点(x,y)的斜平面高程h₂。

步骤5，将TanDEM-DEM在斜平面上高程值h₁与两幅SAR图像干涉处理得到的斜平面高程值h₂作差，得到斜平面上的高程误差：Δh＝h₁-h₂。

步骤6，根据参数文件中提供的波长λ和斜距R参数，求出SAR图像观测场景中的模糊高度：

其中：λ是雷达波长，R是雷达到目标的斜距，θ是下视角，B_⊥是垂直有效基线。

步骤7，利用上述获得的高程误差Δh、模糊高度h_amb，根据基线误差之间的关系获得场景中随着方位向线性变化的平行基线分量误差：

至此，获得场景中的基线误差，由于模糊高度h_amb是沿距离向线性变化的，因此由上述公式得到的基线误差ΔB_||也是线性变化的，而传统估计方法中基线误差是个常数，用线性变化的基线误差校正基线比用常数误差校正基线更加精确。

步骤8，根据平行基线分量误差ΔB_||，对场景方位向的每一个位置都建立坐标系进行实时轨道误差校正，得到基线误差校正后的轨道参数。

(8a)基线误差的根源在于轨道误差，对基线误差沿着轨道状态矢量的三个坐标方向进行分解，得到每个方向的误差分量；

(8b)将原始轨道在每个方向的轨道速度和位置信息加上该方向的轨道误差分量，获得校正后的轨道速度和位置信息。

本发明的效果可通过以下实验进一步说明：

一.实验参数：

本实验中使用的两幅干涉SAR图像的基本参数信息如下表I：

表I中国北京地区的GF-3数据基本参数信息

二.实验步骤

第一步，对两幅干涉SAR图像进行干涉处理，得到干涉图像的斜平面高程信息，再通过地理编码，生成在地平面上数字高程模型DEM产品：

(1.1)对两幅SAR图像进行干涉处理得到解缠相位信息，并将其转换为干涉图像在斜平面的高程矩阵；

(1.2)将斜平面每个像点的高程值经地理编码处理投影到地平面像点上，得到最终的DEM产品；

(1.2.1)设置误差阈值T_X和T_Y，其中T_X和T_Y分别为地平面坐标X和Y的误差阈值；

(1.2.2)设定地平面坐标初值(X₀,Y₀)，并将几何构象模型表达为关于地面点平面坐标的线性化方程组形式，即

C·Δ_G-L＝0

式中，C为系数矩阵，Δ_G为改正量向量，L为几何构像模型中函数F₁(x,y,X,Y,Z)和F₂(x,y,X,Y,Z)的初值向量，其中，(x,y)为斜平面像点坐标，(X,Y)为地平面像点坐标，Z为每个像点的高程值；

(1.2.3)设误差向量

权阵为P，得到误差方程式：

V＝C·Δ_G-L·P，

(1.2.4)取权阵P为单位阵，得到对应的法方程式：

C^TCΔ_G-C^TL＝0

其中C^T表示系数矩阵C的转置矩阵；

根据上述法方程式得到地平面坐标的改正量向量Δ_G为：

Δ_G＝(C^TC)^-1C^TL，

(1.2.5)将改正量向量Δ_G与地平面坐标初值(X₀,Y₀)相加，得到修正后的地面点平面坐标(X,Y)：

其中ΔX为Δ_G在X方向的分量，ΔY为Δ_G在Y方向的分量。

(1.2.6)将修正后的平面坐标(X,Y)作为新的初值(X₀,Y₀)，返回(1.2.2)，直至ΔX＜T_X、ΔY＜T_Y，得到最终的(X,Y)坐标值，生成地平面上数字高程模型DEM产品，如图4(b)所示。

第二步，利用本发明提出的方法计算平行基线分量误差，校正基线和轨道。

第三步，利用校正之后的轨道信息进行DEM重建。

(3.1)对完成基线和轨道校正后的两幅SAR图像再次进行干涉处理得到解缠相位信息，并将其转换为干涉图像在斜平面的高程矩阵；

(3.2)将斜平面每个像点的高程值经地理编码处理投影到地平面像点上，得到最终的DEM产品，如图5(b)所示。

第四步，结果对比。

(4.1)将本发明方法的基线误差与现有的轨道状态矢量法、FFT法的基线误差进行对比，对比结果如图3所示；

(4.2)将干涉图像未进行基线误差校正直接生成的DEM与TanDEM-DEM库进行比对，对比结果如图4所示，其中图4(a)为TanDEM-DEM库的DEM，4(b)为干涉图像未进行基线误差校正直接生成的DEM；

(4.3)将DEM重建结果与TanDEM-DEM库进行对比，对比图结果如图5所示，其中图5(a)为TanDEM-DEM库的DEM，5(b)为DEM重建结果。

基线的精度将直接影响DEM重建的精度，因此可以通过DEM重建后的结果来验证基线估计的结果，下面结合实验结果对本发明的技术效果再做说明：

从图3可见，轨道状态矢量法的精度最差，校正结果存在线性误差，因此不能完全校正基线误差，并且FFT法比轨道状态矢量法精度高，但误差校正后，依然存在少许的基线误差；本发明方法进行基线误差校正之后，得到的基线误差曲线在零附近，意味着本发明提出的方法可以较准确地校正基线误差；

从图4可见，未进行基线误差校正生成的DEM与TanDEM-DEM库的DEM相比，在经度和纬度两个方向上均存在较大误差；

图5可见，使用本发明提出的方法进行基线校正后生成的DEM与TanDEM-DEM库的DEM相比，在经度和纬度两个方向上都不存在调制误差和斜坡误差，进一步体现了本发明的提高基线估计精度和DEM产品精度的优越性。

Claims

2.根据权利要求1所述的方法，其特征在于：(1)中对ICESat卫星星载地学激光测高系统ICESat/GLAS的激光测高数据产品GLA14进行筛选，实现步骤如下：

(1a)初步筛选：

(1b)第二次筛选：

(1c)坐标转换：

将(1b)筛选出的ICESat/GLAS高程数据在海洋环流测量卫星TOPEX/Poseidon椭球坐标系下的坐标值，转换为在世界大地测量系统1984 WGS84椭球坐标系下的坐标值，转换公式如下：

H_G＝H+dh

其中H为坐标转换之前的高程值，H_G为转换后的高程值，dh＝-cos² Bda-sin² Bdb为高程变换修正值，在dh中，B代表纬度，da、db分别表示WGS84椭球和TOPEX/Poseidon椭球长轴和短轴之间的差异值。

3.根据权利要求1所述方法，其特征在于：(2)中利用选取的高程数据，校正TanDEM-X卫星数字高程模型DEM库TanDEM-DEM的高程数据，实现步骤如下：

(2b)根据统计出的均方根误差和标准差，对所选经纬度范围内的TanDEM-DEM高程数据补加标准差，实现校正。

4.根据权利要求1所述的方法，其特征在于：(3)中利用反定位的方法获得校正后的TanDEM-DEM在斜平面上高程值h₁，实现步骤如下：

其中X,Y,Z分别代表WGS84坐标系中的坐标，B,L,H分别代表地平面纬度、地平面经度和高程三维坐标系的坐标，N＝a(1-e² sin² B)^-1/2是曲率半径，a是N中的半长轴，e²表示第一偏心率；

5.根据权利要求1所述的方法，其特征在于：(4)中对两幅SAR图像进行干涉处理得到解缠相位信息，并将其转换为干涉图像在斜平面的高程值h₂，是通过图像配准、相位生成、相位滤波以及相位解缠处理，得到解缠相位，再进行相位与高程的转换，获得干涉图像在斜平面的高程值h₂。

6.根据权利要求1所述的方法，其特征在于：(6)中根据参数文件中提供的波长λ和斜距R参数，求出SAR图像观测场景中的模糊高度h_amb，通过以下公式得到：

7.根据权利要求1所述的方法，其特征在于：(8)中根据平行基线分量误差ΔB_||，对场景方位向的每一个位置都建立坐标系进行实时轨道误差校正，其实现如下：

(8a)对基线误差沿着轨道状态矢量的三个坐标方向进行分解，得到每个方向的误差分量；