CN109116351A - 一种星载InSAR定位解析算法 - Google Patents

Abstract

本发明公开了一种星载InSAR定位解析算法，包括以下步骤：a、根据主雷达天线相位中心在地心坐标系的位置和速度，在主雷达方位向上建立以主雷达天线相位中心为原点的局部坐标系；b、根据主雷达天线相位中心、辅雷达天线相位中心的位置和速度，计算获得待定位点P在局部坐标系下的坐标；c、根据主雷达天线相位中心在地心坐标系下的坐标和待定位点P在局部坐标系下的坐标，获得待定位点P在地心坐标系下的坐标。本发明的算法利用更加严密的算法模型，而且几何意义明显，可以同时保证算法编码精度，并使编码效率明显提高。星载InSAR定位得到DSM时，最终使得算法精度与Mora模型法相同，但效率提高5倍。

Description

一种星载InSAR定位解析算法

技术领域

本发明涉及定位算法，尤其涉及一种星载InSAR定位解析算法。

背景技术

合成孔径雷达干涉(InSAR)，是一种应用于测绘和遥感的雷达技术。它是 利用合成孔径雷达对同一地区观测的两幅复数值影像数据进行相干处理，以获 取地表高程信息的技术。InSAR定位解析得到高精度的数字高程模型(DEM)， 是InSAR数据处理的一个重要的环节，其目的是利用基线矢量、轨道、速度和 绝对干涉相位等信息求解地面点三维坐标。已有的SAR定位算法大体可分为两 类，第一类是高程平面分开答解方法，即先求出高程，再求解出平面坐标，最 终得到地面点的三维坐标，典型的有Ambiguity方法、Rodriguez法和Schwabish法等。第二类方法是将地面点的三维坐标作为未知数，利用一定的 模型进行求解，典型的有Mora模型法、距离向闭合定位模型法。

无论是Ambiguity、Rodriguez还是Schwabish方法，在进行平面高程分开 求解时均进行了一定程度的近似。这种近似使得定位模型不严密，从而导致定 位结果存在较大误差。Mora模型法采用严密定位模型进行编码，其优点是定位 过程不损失精度，缺点是方程无法直接答解，但只能通过迭代方式进行，编码 效率不高。距离向闭合定位模型法采用直接法进行求解，模型严密，但其缺点 是模型没有明显的几何意义，同时逐点求解方向矢量导致求解效率较低。

发明内容

为了解决解决现有技术条件下，星载InSAR定位得到DSM时，或定位误差较 大、或计算效率较低的问题，本发明提供了一种星载InSAR定位解析算法。

为了解决以上技术问题，本发明采用的技术方案是：一种星载InSAR定位 解析算法，包括以下步骤：

a、根据主雷达天线相位中心在地心坐标系的位置和速度，在主雷达方位 向上建立以主雷达天线相位中心为原点的局部坐标系；

b、根据主雷达天线相位中心、辅雷达天线相位中心的位置和速度，计算 获得待定位点P在局部坐标系下的坐标；

c、根据主雷达天线相位中心在地心坐标系下的坐标和待定位点P在局部坐 标系下的坐标，获得待定位点P在地心坐标系下的坐标。

进一步地，所述步骤a中，设主雷达天线相位中心在地心坐标系的位置为速度为定义以主雷达天线相位中心为原点的局 部坐标系O-X′Y′Z′，Y′轴为X′轴为Z′轴为X′×Y′，遵循右手定 律。

进一步地，所述步骤b中，待定位点P在局部坐标系O-X′Y′Z′下的坐标的计 算方法具体为:

根据Mora模型法获得以下三个方程：

其中，X₁，Y₁，Z₁表示主雷达天线相位中心的坐标，X₂，Y₂，Z₂表示辅雷达天 线相位中心的坐标，X_t，Y_t，Z_t表示待定位点P坐标。R₁表示主雷达天线相位中 心到待定位点P距离，R₂表示辅雷达天线相位中心到待定位点P距离，λ表示雷 达发射信号波长，φ表示干涉相位，ρ表示模糊高度，f_d表示多普勒频率；

将上述方程应用于局部坐标系O-X′Y′Z′，设是成像时，主雷达天线相位 中心与辅雷达天线相位中心的连线,B_X、B_Y、B_Z是基线矢量分别在O-X′Y′Z′ 坐标系下三轴的投影；此时，X₁＝Y₁＝Z₁＝0，X₂＝B_X，Y₂＝B_Y，Z₂＝B_Z，则方程 (1)，(2)变为：

(B_X-X_t)²+(B_Y-Y_t)²+(B_Z-Z_t)²＝(R₁+λφ/2ρπ)² (5)

方程(3)表示以主星APC速度为轴的一个圆锥，圆锥的张角表示为 γ,

将方程(3)带入方程(1)，方程(1)变成一个周长固定的圆锥的底部，圆锥的 边长为R₁，圆的半径为R₀＝R₁sinγ。点P必然在这个圆上，其Y′坐标为R₁cosγ；

将Y′＝R₁cosγ带入方程(4)和方程(5)，可得

设则方程(7)变成

方程(6)、(8)表示在X′OZ′坐标系下两个圆相交，其交点分别为A₁和A₂，其 中A₁为虚解，A₂为正解；

易得，X_t＝R₀cos(α-β)，Z_t＝R₀sin(α-β)，其中α＝atan(B_Z/B_X)，

由此可得待定位点P在局部坐标系O-X′Y′Z′下的坐标 (X_t,Y_t,Z_t)＝[R₁cos(α-β),R₁cosγ,R₁sin(α-β)]。

进一步地，所述步骤c中，待定位点P在地心坐标系下的坐标(X_P,Y_P,Z_P)的计 算方法具体为：

其中，(x_a,y_a,z_a)，(x_b,y_b,z_b)和(x_c,y_c,z_c)分别为局部坐标系O-X′Y′Z′的X′、Y′ 和Z′三轴在地心坐标系下的单位矢量；(X₁,Y₁,Z₁)是主雷达APC在地心坐标系下 的坐标。

进一步地，所述星载InSAR定位解析算法还包括：步骤d、按照步骤a-c获 得目标地面每一点在地心坐标系下的坐标，并按照其构建高精度的数字高程模 型。

本发明在首先在方位向上建立局部坐标系并推导出地面待测点在局部坐标 下的三维坐标，然后建立局部坐标系与地心坐标系之间的转换关系，并求解出 该点在地心坐标系下的坐标。本发明的算法利用更加严密的算法模型，而且几 何意义明显，可以同时保证算法编码精度，并使编码效率明显提高。星载 InSAR定位得到DSM时，最终使得算法精度与Mora模型法相同，但效率提高5 倍。

附图说明

图1为局部坐标系定义示意图。

图2为多普勒方程与距离方程联合示意图。

图3为X′OZ′坐标系下两个圆相交的示意图。

图4为试验地区雷达强度图像。

图5为试验地区DEM数据图。

图6为试验过程中控制点分布情况示意图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细的说明。

卫星在飞行过程中，其速度和位置会不断发生变化，但是对于同一距离 向，其基线矢量、主雷达天线相位中心(APC)速度、位置以及多普勒中心频率 相同。本发明思路就是在雷达图像每一距离向上建立局部坐标系，求解距离向 上每一点在局部坐标系中的坐标，然后将其转换到地心坐标系中，从而求得地 面点的三维坐标。

a、局部坐标系定义：

设主雷达APC在地心坐标系的位置为速度为定义Y′轴为定义X′轴为Z′轴为X′×Y′,符号 ×表示矢量叉乘，遵循右手定律，如图1所示。

R₂是辅雷达APC到地面点P的距离，是主辅雷达成像时APC的连线,B_X、 B_Y、B_Z是基线矢量分别在O-X′Y′Z′坐标系下三轴的投影。

b、求解地面点P在局部坐标系下的坐标：

Mora法可用以下三个方程表示：

其中，X₁，Y₁，Z₁表示雷达主星坐标，X₂，Y₂，Z₂表示雷达辅星坐标， X_t，Y_t，Z_t表示地面点坐标。R₁表示主星到地面点距离，R₂表示辅星到地面点距 离，λ表示雷达发射信号波长，φ表示干涉相位，ρ表示模糊高度，f_d表示多 普勒频率。

本发明算法的原理为：在局部坐标系下，X₁＝Y₁＝Z₁＝0，X₂＝B_X，Y₂＝B_Y， Z₂＝B_Z，则方程(1)，(2)变为：

(B_X-X_t)²+(B_Y-Y_t)²+(B_Z-Z_t)²＝(R₁+λφ/2ρπ)² (5)

方程(3)表示以主星APC速度为轴的一个圆锥，圆锥的张角表示为 γ,如图2所示。

将方程(3)带入方程(1)，方程(1)变成一个周长固定的圆锥的底部，圆锥的 边长为R₁，圆的半径为R₀＝R₁sinγ。点P必然在这个圆上，其Y′坐标为R₁cosγ。

将Y′＝R₁cosγ带入方程(4)和方程(5)，可得

设则方程(7)变成

方程(6)、(8)表示在X′OZ′坐标系下两个圆相交，如图3所示，其交点分别为 A₁和A₂，其中A₁为虚解，A₂为正解。

易得，X_t＝R₀cos(α-β)，Z_t＝R₀sin(α-β)。

其中α＝atan(B_Z/B_X)，

由此可得地面点P在局部坐标系下的坐标 (X_t,Y_t,Z_t)＝[R₁cos(α-β),R₁cosγ,R₁sin(α-β)]。

从图2和图3可以看出，在局部坐标下，当γ＝0时，此时的原理示意图与经 典InSAR模型 相同，也即经典定位模型是本模型的一个特例，当多普勒中心频率f_d＝0时，本模型退化为经 典InSAR二维模型。

c、局部坐标系向全局坐标系的转换：

设局部坐标系O-X′Y′Z′三轴在地心坐标系下的单位矢量分别为(x_a,y_a,z_a)， (x_b,y_b,z_b)和(x_c,y_c,z_c)，地面点P在局部坐标系下的坐标为(X_t,Y_t,Z_t)，则其在地心 坐标系下的坐标(X_P,Y_P,Z_P)为：

其中，是主雷达APC在地心坐标系下的坐标。

d、获得目标地面每一点在地心坐标系下的坐标，并按照其构建高精度的 数字高程模型。

一、试验与结果分析：

为了验证算法的正确性与效率，采用了基于回波仿真模式的仿真数据进行 了试验。仿真地区为河北尚义地区，场景大小10km×10km、X波段、轨道高度 515km、波束下视角41°，波长为0.032m、影像采样间隔为3.2×1.6m(方位向×距 离向)，像素个数为4272×6207(方位向×距离向)，基线长度约为740m，其强 度图像和DEM数据如图4和图5所示。

二、定位精度试验：

场地中共布设有100个控制点，其分布如图6所示，其中横轴是距离向，纵 轴是方位向。从图中可以看出，控制点分布基本是均匀分布。

定位精度试验主要是验证三种不同地理编码方法的定位精度，分别用三种 方法计算控制点的坐标，并与实际坐标进行比较，统计三轴的中误差和总定位 误差如表1。从表1可以看出，三种方法定位精度基本一致，差别在毫米级，误 差基本可以忽略不计。

表1三种方法定位精度统计 单位(m)

三、编码效率试验：

以整个场景数据为单位，逐点进行地理编码，然后统计三种方法编码所需 时间。测试环境为：PⅢXEON 2.8G处理器、内存2GB、硬盘1TB，其所消耗 时间如表2所示。

表2编码所花费时间 单位(s)

可以看出，相同场景情况下，Mora模型法消耗时间最长，直接地理编码法 消耗时间最短，这是因为后者只需要在整个距离向上定义一次局部坐标系即 可，其余就是局部坐标系向全局坐标系转换所耗费时间。距离向闭合定位模型 法消耗时间在两者之间，这是因为该方法需要对每个点计算一次斜距矢量因此耗费了时间。

上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本 技术领域的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或 替换，也均属于本发明的保护范围。

Claims (5)

1.一种星载InSAR定位解析算法，其特征在于：包括以下步骤：

a、根据主雷达天线相位中心在地心坐标系的位置和速度，在主雷达方位向上建立以主雷达天线相位中心为原点的局部坐标系；

b、根据主雷达天线相位中心、辅雷达天线相位中心的位置和速度，计算获得待定位点P在局部坐标系下的坐标；

c、根据主雷达天线相位中心在地心坐标系下的坐标和待定位点P在局部坐标系下的坐标，获得待定位点P在地心坐标系下的坐标。

2.根据权利要求1所述的星载InSAR定位解析算法，其特征在于：所述步骤a中，设主雷达天线相位中心在地心坐标系的位置为速度为定义以主雷达天线相位中心为原点的局部坐标系O-X′Y′Z′，Y′轴为X′轴为Z′轴为X′×Y′，遵循右手定律。

3.根据权利要求2所述的星载InSAR定位解析算法，其特征在于：所述步骤b中，待定位点P在局部坐标系O-X′Y′Z′下的坐标的计算方法具体为：

根据Mora模型法获得以下三个方程：

其中，X₁，Y₁，Z₁表示主雷达天线相位中心的坐标，X₂，Y₂，Z₂表示辅雷达天线相位中心的坐标，X_t，Y_t，Z_t表示待定位点P坐标。R₁表示主雷达天线相位中心到待定位点P距离，R₂表示辅雷达天线相位中心到待定位点P距离，λ表示雷达发射信号波长，φ表示干涉相位，ρ表示模糊高度，f_d表示多普勒频率；

将上述方程应用于局部坐标系O-X′Y′Z′，设是成像时，主雷达天线相位中心与辅雷达天线相位中心的连线,B_X、B_Y、B_Z是基线矢量分别在O-X′Y′Z′坐标系下三轴的投影；此时，X₁＝Y₁＝Z₁＝0，X₂＝B_X，Y₂＝B_Y，Z₂＝B_Z，则方程(1)，(2)变为：

(B_X-X_t)²+(B_Y-Y_t)²+(B_Z-Z_t)²＝(R₁+λφ/2ρπ)² (5)

方程(3)表示以主星APC速度为轴的一个圆锥，圆锥的张角表示为γ，

将方程(3)带入方程(1)，方程(1)变成一个周长固定的圆锥的底部，圆锥的边长为R₁，圆的半径为R₀＝R₁sinγ。点P必然在这个圆上，其Y′坐标为R₁cosγ；

将Y′＝R₁cosγ带入方程(4)和方程(5)，可得

设则方程(7)变成

方程(6)、(8)表示在X′OZ′坐标系下两个圆相交，其交点分别为A₁和A₂，其中A₁为虚解，A₂为正解；

易得，X_t＝R₀cos(α-β)，Z_t＝R₀sin(α-β)，其中α＝atan(B_Z/B_X)，

由此可得待定位点P在局部坐标系O-X′Y′Z′下的坐标

(X_t,Y_t,Z_t)＝[R₁cos(α-β),R₁cosγ,R₁sin(α-β)]。

4.根据权利要求3所述的星载InSAR定位解析算法，其特征在于：所述步骤c中，待定位点P在地心坐标系下的坐标(X_P,Y_P,Z_P)的计算方法具体为：

其中，(x_a,y_a,z_a)，(x_b,y_b,z_b)和(x_c,y_c,z_c)分别为局部坐标系O-X′Y′Z′的X′、Y′和Z′三轴在地心坐标系下的单位矢量；(X₁,Y₁,Z₁)是主雷达APC在地心坐标系下的坐标。

5.根据权利要求4所述的星载InSAR定位解析算法，其特征在于：所述星载InSAR定位解析算法还包括：步骤d、按照步骤a-c获得目标地面每一点在地心坐标系下的坐标，并按照其构建高精度的数字高程模型。