CN102129066A - 一种宽幅星载sar快速地理编码方法 - Google Patents

Abstract

本发明将严格符合星地空间几何关系的斜距多普勒定位模型和基于GCP发展的遥感图像多项式校正模型进行了整合，得到一套快速星载SAR图像地理编码方法框架，不但简化了传统星载SAR地理编码的流程，提高宽幅星载SAR地理编码的效率，而且有效消除了SAR成像算法中不精准多普勒中心频率带来的几何变形。本发明在单像素地理编码中，给出了斜距多普勒非线性方程组的代数解析解，解决了传统斜距-多普勒(RD)定位模型使用数值迭代方法面临的初始值选取的问题。本发明是一种新的快速星载SAR地理编码的方法，能广泛应用于星载SAR图像地理编码。

Description

一种宽幅星载SAR快速地理编码方法

技术领域

本发明属于遥感图像处理技术领域，涉及遥感图像地理编码方法，特别涉及一种宽幅星载SAR快速地理编码方法。

背景技术

随着我国综合国力和空间技术的发展，星载合成孔径雷达作为一种工作在微波波段的主动式遥感器，能够不受日照和天气条件的限制完成对地观测，因此在农、林、地质、环境、水文、海洋、灾害、测绘与军事等领域的应用受到越来越多的重视。

作为遥感应用的重要支柱之一，高精度的星载SAR图像地理编码产品作为地面处理系统最重要的产品，直接面向用户。面对星载SAR观测带更宽，观测数据庞大以及应用更广泛的需求，具有高精度高效率的宽幅星载SAR自动地理编码方法成为工程研究的热点。

目前常用的星载SAR地理编码方法可以分为两类：

一种是从SAR图像和对应的数字地图中抽取地面控制点(GCP)进行自动匹配校正来完成图像地理编码的方法。该类方法图像定位精度高，但实现自动匹配的技术鲁棒性不高，尤其在海面、沙漠等缺乏有效地面控制点的地区难以实现。该类方法往往是从事光学遥感的研究人员，借助摄影测量的共线方程或者多项式近似方程来实现的。

另一种是利用卫星平台及SAR等系统参数直接进行SAR图像地理编码。该类方法不需要GCP，易于实现自动地理编码。这些方法包括基于星地坐标转换的地理编码方法和基于斜距多普勒定位模型的地理编码方法。

1、基于星地坐标转换的方法

基于星地坐标转换的方法，利用星地坐标之间的6个坐标系关系，借助观测视角的中间变量，不需要迭代运算完成星载SAR图像地理编码。但该方法步骤繁琐，而且没有考虑成像算法带来的图像变形，对宽幅的星载SAR图像地理编码使用性差。

2、基于斜距多普勒定位模型的方法

基于斜距多普勒定位模型的方法不但是一种严格符合星地空间几何关系的精确的SAR图像地理编码方法，而且基于多普勒方程的约束条件，有效消除成像算法中不精准多普勒中心频率带来的图像几何变形。但该方法由三个较复杂的非线性方程构成。传统方法均使用数值迭代，而当初始值选取不合适时，结果可能不收敛，同时迭代处理需要消耗过多的时间，严重影响SAR地理编码的效率。而且基于斜距多普勒模型建立图像映射关系过程中会出现病态方程，导致星载SAR地理编码自动流程失败。

发明内容

本发明的目的在于提供一种基于斜距多普勒定位模型和多项式校正模式的宽幅星载SAR快速地理编码方法，利用地固坐标系下目标速度为零的特性，简化斜距多普勒模型构成的非线性方程组，转化成一元四次方程，利用费拉里的求根方法，给出斜距多普勒非线性方程组的解析解，并结合多项式校正模型，进行宽幅星载SAR地理编码。

为实现上述目的，本发明的主要步骤如下：

步骤S21：提取星载SAR实图像产品中辅助文件的相关参数；

步骤S22：拟合卫星轨道曲线，获得卫星位置速度和图像方位时刻的关系；

步骤S23：根据地图图像产品用途选择地图投影模型参数、地图投影模型变形系数和地图图像像素分辨率；

步骤S24：单像素地理编码；

步骤S25：在斜距图像上，沿距离向均匀选择M组，沿方位向均匀选择若干N组，共M×N组像素点，重复步骤S24的单像素地理编码；

步骤S26：根据多项式校正模型，用步骤S25得到的若干组斜距图像像素和地图图像像素，计算多项式的参数；

步骤S27：根据步骤S26得到的多项式参数，依次计算地图图像每个像素对应的图像像素，然后使用双线性插值法对图像进行重采样。

步骤S28：格式化输出。

所述的宽幅星载SAR快速地理编码方法中，所述步骤S25中的M组多于4组，N组多于3组，具体可以根据星载SAR斜距图像的幅宽进行调整，若幅宽小于或等于50×50公里，M组为6组，N组为4组，共6×4＝24组像素即可。

所述的宽幅星载SAR快速地理编码方法中，所述步骤S26中的多项式校正模型为三次多项式校正模型，具体如下：

i＝a₀+a₁E+a₂N+a₃EN+a₄E²+a₅N²+a₆E²N+a₇EN²+a₈E³+a₉N³

j＝b₀+b₁E+b₂N+b₃EN+b₄E²+b₅N²+b₆E²N+b₇EN²+b₈E³+b₉N³

在S24单像素地理编码中，给出斜距多普勒非线性方程组的解析解，具体步骤如下：

步骤S11：获取建立斜距多普勒非线性方程组的相关参数

步骤S12：利用地固坐标系下目标速度为零的特性，简化斜距多普勒模型构成的非线性方程组，通过分别将斜距方程和多普勒中心频率方程联立及多普勒中心频率方程和地球模型方程联立，转化成目标矢量y轴(R_ty)的一元四次方程。

步骤S13：利用费拉里求根方法，求解步骤S12构建的一元四次方程，并根据视角方向选择R_ty的有效解。

步骤S14：将步骤S13的解代入表达式，并根据视角方向选择目标矢量x轴(R_tx)的有效解。

步骤S15：将步骤S13的解和步骤S14的解代入多普勒中心频率方程，得到目标矢量z轴(R_tz)的解。

本发明的有益效果：

本发明中推导的斜距多普勒非线性方程组的解析解，解决了传统迭代方法需要设置初始值的问题，避免了迭代运算和不合适的初始值导致无法收敛的缺陷，减少星载SAR地理编码模块的输入接口，提高可靠性。

本发明采用的斜距多普勒模型，避免了繁琐的六坐标转换，星地几何关系更加准确；处理精度不但不会受到卫星姿态参数的影响，而且充分利用成像处理中的多普勒参数，有效校正成像处理使用有误差的多普勒参数带来的成像几何变形。

本发明采用的三次多项式校正模型，既保证了斜距多普勒模型的定位精度，又简化了斜距多普勒模型逆运算的处理流程，提高了地理编码效率。

本发明结合了斜距多普勒模型和三次多项式校正模型，不但充分考虑成像处理带来的几何畸变，而且处理流程模块化强，效率高，适合数据量大、自动处理的需要，有利于计算机或其它专用设备实现。

附图说明

图1是本发明基于斜距多普勒定位模型和多项式校正模型的宽幅星载SAR快速地理编码方法的完整流程图；

图2是本发明星载SAR图像地理编码示意图；

图3是本发明斜距多普勒(RD)模型组成的非线性方程组的解析解推导流程图；

图4是本发明求解斜距多普勒(RD)非线性方程组的流程图；

图5是本发明单像素地理编码流程图；

图6是本发明基于斜距多普勒定位模型和多项式校正模型的宽幅星载SAR快速地理编码方法的基本步骤图；

图7、图8和图9是本发明实施例的效果示意图，其中：

图7a是SAR的斜距图像(图像大小为距离向5088×方位向18432)，图7b是SAR地理编码后的地图图像(图像大小为指东9539×指北7454)；

图9a和图9b是本发明的SAR地图图像和google earth上的光学地图在相同地图投影坐标系下比较。

图中符号说明如下：

λ：SAR的工作波长；f_d：多普勒中心频率；R：斜距；R_sx：卫星x方向位置；R_sy：卫星y方向位置；R_sz：卫星z方向位置；V_sx：卫星x方向速度；V_sy：卫星y方向速度；V_sz：卫星z方向速度；R_tx：目标x方向位置；R_ty：目标y方向位置；R_tz：目标z方向位置。

具体实施方式

星载SAR地理编码是将SAR原始回波信号经过成像处理后得到的斜距图像，按照某种通用的地图投影方式，投影到以地图坐标系(指东指北)的图像过程。在该过程中，需要进行图像像素的精确定位，并对图像几何畸变进行校正，从而得到符合制图标准的具有地理信息的图像，以便于人们对SAR图像进行理解和判读。

本发明提供的基于斜距多普勒定位模型和多项式校正模式的宽幅星载SAR快速地理编码方法，利用地固坐标系下目标速度为零的特性，简化斜距多普勒模型构成的非线性方程组，转化成一元四次方程，利用费拉里的求根方法，给出斜距多普勒非线性方程组的解析解，并结合多项式校正模型，进行宽幅星载SAR地理编码。

首先给出权利要求5中斜距多普勒非线性方程组的解析解的具体步骤为：

步骤S11：获取参数。获取建立斜距多普勒非线性方程组的参数，包括：图像产品某像素(i，j)对应的斜距值R，多普勒中心频率fd，卫星位置矢量 ${\stackrel{\→}{R}}_s={[R_{\mathrm{sx}},R_{\mathrm{sy}},R_{\mathrm{sz}}]}^T,$ 卫星速度矢量 ${\stackrel{\→}{V}}_s={[V_{\mathrm{sx}},V_{\mathrm{sy}},V_{\mathrm{sz}}]}^T,$ 以及图像产品对应的波长λ，拟采用的地球模型[R_e，R_p]。

步骤S12：计算一元四次方程系数变量的数值大小。将S11中的斜距多普勒非线性方程组的参数代入如下表达式，计算得到一元四次方程 $\mathrm{AA}\·R_{\mathrm{ty}}^4+\mathrm{BB}\·R_{\mathrm{ty}}^3+\mathrm{CC}\·R_{\mathrm{ty}}^2+\mathrm{DD}\·R_{\mathrm{ty}}+\mathrm{EE}=0$ 的系数变量AA，BB，CC，DD，EE。

K₁＝R²

$K_2=-\frac{\mathrm{\λ}\·R\·\mathrm{fd}}{2\·V_{\mathrm{sz}}}$ $T_1=R_e^2$

$K_3=\frac{V_{\mathrm{sx}}}{V_{\mathrm{sz}}}$ $T_2=\frac{R_e}{R_p}(R_{\mathrm{sz}}+\frac{\mathrm{\λ}\·R\·\mathrm{fd}}{2\·V_{\mathrm{sz}}}+\frac{V_{\mathrm{sx}}}{V_{\mathrm{sz}}}{R}_{\mathrm{sx}}+\frac{V_{\mathrm{sy}}}{V_{\mathrm{sz}}}{R}_{\mathrm{sy}})$

$K_4=\frac{V_{\mathrm{sy}}}{V_{\mathrm{sz}}}$ $T_3=\frac{V_{\mathrm{sx}}}{V_{\mathrm{sz}}}\·\frac{R_e}{R_p}$

A＝R_sx-R_tx

B＝R_sy-R_ty， $T_4=\frac{V_{\mathrm{sy}}}{V_{\mathrm{sz}}}\·\frac{R_e}{R_p},$

$M=(K_3^2+1)T_3{T}_4-(T_3^2+1)K_3{K}_4$

$N=R_{\mathrm{sx}}(T_3^2+1)(K_3^2+1)-T_2{T}_3(K_3^2+1)+K_3{K}_4{R}_{\mathrm{sy}}(T_3^2+1)-K_2{K}_3(T_3^2+1)$

$P_1=-(T_3^2+T_4^2+1)$

Q₁＝2T₂T₄

$R_1=T_1{T}_3^2-T_2^2+T_1$

$P_2=-(K_3^2+K_4^2+1)$

$Q_2=2R_{\mathrm{sy}}(K_3^2+K_4^2+1)-2K_2{K}_4$

$R_2=-(K_3^2+K_4^2+1)R_{\mathrm{sy}}^2+2K_2{K}_4{R}_{\mathrm{sy}}+K_1{K}_3^2-K_2^2+K_1,$

$F_1=M^2-{(K_3^2+1)}^2{P}_1-{(T_3^2+1)}^2{P}_2$

$F_2=2\mathrm{MN}-{(K_3^2+1)}^2{Q}_1-{(T_3^2+1)}^2{Q}_2$

$F3=N^2-{(K_3^2+1)}^2{R}_1-{(T_3^2+1)}^2{R}_2$

$F_4=2(K_3^2+1)(T_3^2+1)$

$\mathrm{AA}=F_4^2{P}_1{P}_2-F_1^2$

$\mathrm{BB}=F_4^2(P_1{Q}_2+P_2{Q}_1)-2F_1{F}_2$

$\mathrm{CC}=F_4^2(R_1{R}_2+Q_1{Q}_2+P_2{R}_1)-2F_1{F}_3-F_2^2$

$\mathrm{DD}=F_4^2(Q_1{R}_2+R_1{Q}_2)-2F_2{F}_3$

$\mathrm{EE}=F_4^2{R}_1{R}_2-F_3^2.$

步骤S13：利用费拉里求根方法，求解步骤S12构建的一元四次方程 $\mathrm{AA}\·R_{\mathrm{ty}}^4+\mathrm{BB}\·R_{\mathrm{ty}}^3+\mathrm{CC}\·R_{\mathrm{ty}}^2+\mathrm{DD}\·R_{\mathrm{ty}}+\mathrm{EE}=0,$ 根据视角方向选择R_ty的有效解。

步骤S14：将步骤S13的解代入如下表达式，并根据视角方向选择R_tx的有效解，

$R_{\mathrm{tx}}=\frac{T_2{T}_3-T_3{T}_4{R}_{\mathrm{ty}}}{T_3^2+1}\±\frac{\sqrt{-(T_3^2+T_4^2+1)R_{\mathrm{ty}}^2+2T_2{T}_4{R}_{\mathrm{ty}}+T_1{T}_3^2-T_2^2+T_1}}{T_3^2+1}$

步骤S15：将步骤S13的解和步骤S14的解代入如下表达式，得到R_tz的解，

$R_{\mathrm{tz}}=R_{\mathrm{sz}}+\frac{\mathrm{\λ}\·R\·\mathrm{fd}}{2\·V_{\mathrm{sz}}}+\frac{V_{\mathrm{sx}}}{V_{\mathrm{sz}}}(R_{\mathrm{sx}}-R_{\mathrm{tx}})+\frac{V_{\mathrm{sy}}}{V_{\mathrm{sz}}}(R_{\mathrm{sy}}-R_{\mathrm{ty}})$

所述根据斜距多普勒模型，并结合多项式校正模型，进行宽幅星载SAR地理编码的具体步骤为：

步骤S21：提取星载SAR实图像产品中辅助文件的相关参数，包括：图像产品的初始斜距R_min，采样频率f_s，多普勒中心频率沿距离像素变化的参数：常数项f_d0，一次项f_d1，二次项f_d2，图像产品初始时刻t_min，脉冲重复频率prf，波长λ，图像产品对应时刻内若干组卫星位置和速度矢量t_n：[R_sxn，R_syn，R_szn，V_sxn，V_syn，V_szn]，以及地球模型[R_e，R_p]。

步骤S22：拟合卫星轨道曲线，获得卫星位置速度和图像方位时刻的关系。

步骤S23：根据地图图像产品用途选择地图投影模型参数，包括：投影中心点经度，投影中心点纬度，投影中心点指东值，投影中心点指北值，以及地图投影模型变形系数。同时还需要根据地图图像产品用途选择地图图像像素分辨率。

步骤S24：单像素地理编码。根据权力要求5的方法求解斜距多普勒非线性方程组，实现对图像某像素(i，j)的地理编码。

步骤S25：在斜距图像上，根据权力要求2，3的原则选择若干组像素重复步骤S24的单像素地理编码。

步骤S26：根据权利要求4设置如下多项式校正模型，用步骤S25得到的若干组斜距图像像素(i，j)和地图图像像素(E，N)，计算多项式的参数a_i，b_i。

i＝a₀+a₁E+a₂N+a₃EN+a₄E²+a₅N²+a₆E²N+a₇EN²+a₈E³+a₉N³

j＝b₀+b₁E+b₂N+b₃EN+b₄E²+b₅N²+b₆E²N+b₇EN²+b₈E³+b₉N³

步骤S27：根据步骤S26得到的多项式参数，依次计算地图图像每个像素(E_n，N_n)对应的(i_n，j_n)，然后对图像进行重采样。

步骤S28：格式化输出。将步骤S23得到的地图图像和地图投影模型以及地图图像像素分辨率等信息，按照国际geotiff格式标准打包，形成geotiff图。

所述步骤S24中的单像素地理编码具体步骤为：

步骤S241：根据权利要求5的方法求解斜距多普勒模型组成的非线性方程组，进行像素定位，得到斜距图像某像素(i，j)对应的地固坐标系下的位置矢量(x_t，y_t，z_t)；

步骤S242：计算经纬度，将地固坐标系的位置矢量(x_t，y_t，z_t)转换为地理经纬度(B，L)；

步骤S243：根据步骤S23设定的地图投影模型，将地理经纬度(B，L)转换为地图坐标系下的指东指北坐标(E′，N′)；

步骤S244：根据步骤S23设定的地图图像像素分辨率，将地图投影后的指东指北坐标(E′，N′)转换为地图图像坐标(E，N)。

所述步骤S25中选择像素的原则为：距离向多于4组，方位向多于3组，并且保证均匀性，可以根据星载SAR斜距图像的幅宽进行调整。若幅宽小于等于50公里×50公里，沿距离向均匀选择6组，沿方位向均匀选择4组，共6×4＝24组像素即可。

实施例

本发明以2006年10月12日山东地区的ALOS的1.1级产品为例，结合附图作进一步的说明。本发明的具体实施步骤如图6所示：

步骤S21：提取星载SAR实图像产品中辅助文件的相关参数，包括以下参数以及若干组卫星位置速度矢量。

本步骤的实施条件需要根据ALOS的产品格式说明，提取相关的卫星参数。

步骤S22：拟合卫星轨道曲线，获得卫星位置速度和图像方位时刻的关系。

步骤S23：根据地图图像产品用途选择地图投影模型参数，包括：投影中心点经度，投影中心点纬度，投影中心点指东值，投影中心点指北值，以及地图投影模型变形系数。同时还需要根据地图图像产品用途选择地图图像像素分辨率ΔE，ΔN。

选择UTM投影方式：投影中心点经度(37度)，投影中心点纬度(118.5度)，投影中心点指东值(4.0e5)，投影中心点指北值(1.0e5)，变形系数(0.99960)。

选择地图图像像素分辨率ΔE＝9，ΔN＝9。

步骤S24：单像素地理编码。如图5所示，具体步骤由步骤S241，S242，S243，S244组成：

步骤S241：求解斜距多普勒模型组成的非线性方程组，进行像素定位，计算斜距图像某像素(i，j)对应的地固坐标系下的位置矢量(x_t，y_t，z_t)。如图3所示，具体步骤由步骤S11，S12，S13，S14，S15组成：

步骤S11：获取像素(i，j)对应的斜距多普勒非线性方程组的参数。

f_d＝f_d0+f_d1·i+f_d2·i²

$t_j=t_{\min }+\frac{j}{\mathrm{prf}}$

$R_{\mathrm{sx}}=x_0+x_1\·t_j+x_2\·t_j^2+x_3\·t_j^3$

$R_{\mathrm{sy}}=y_0+y_1\·t_j+y_2\·t_j^2+y_3\·t_j^3$

$R_{\mathrm{sz}}=z_0+z_1\·t_j+z_2\·t_j^2+z_3\·t_j^3$

$V_{\mathrm{sx}}=x_1+{2x\·}_2\·t_j+3\·x_3\·t_j^2$

$V_{\mathrm{sy}}=y_1+{2y\·}_2\·t_j+3\·y_3\·t_j^2$

$V_{\mathrm{sz}}=z_1+{2z\·}_2\·t_j+3\·z_3\·t_j^2$

步骤S12：计算中间变量及一元四次方程系数

K₁＝R²

$K_2=-\frac{\mathrm{\λ}\·R\·\mathrm{fd}}{2\·V_{\mathrm{sz}}}$ $T_1=R_e^2$

$K_3=\frac{V_{\mathrm{sx}}}{V_{\mathrm{sz}}}$ $T_2=\frac{R_e}{R_p}(R_{\mathrm{sz}}+\frac{\mathrm{\λ}\·R\·\mathrm{fd}}{2\·V_{\mathrm{sz}}}+\frac{V_{\mathrm{sx}}}{V_{\mathrm{sz}}}{R}_{\mathrm{sx}}+\frac{V_{\mathrm{sy}}}{V_{\mathrm{sz}}}{R}_{\mathrm{sy}}$

$K_4=\frac{V_{\mathrm{sy}}}{V_{\mathrm{sz}}}$ $T_3=\frac{V_{\mathrm{sx}}}{V_{\mathrm{sz}}}\·\frac{R_e}{R_p}$

A＝R_sx-R_tx

B＝R_sy-R_ty $T_4=\frac{V_{\mathrm{sy}}}{V_{\mathrm{sz}}}\·\frac{R_e}{R_p},$

$M=(K_3^2+1)T_3{T}_4-(T_3^2+1)K_3{K}_4$

$N=R_{\mathrm{sx}}(T_3^2+1)(K_3^2+1)-T_2{T}_3(K_3^2+1)+K_3{K}_4{R}_{\mathrm{sy}}(T_3^2+1)-K_2{K}_3(T_3^2+1)$

$P_1=-(T_3^2+T_4^2+1)$

Q₁＝2T₂T₄

$R_1=T_1{T}_3^2-T_2^2+T_1$

$P_2=-(K_3^2+K_4^2+1)$

$Q_2=2R_{\mathrm{sy}}(K_3^2+K_4^2+1)-2K_2{K}_4$

$R_2=-(K_3^2+K_4^2+1)R_{\mathrm{sy}}^2+2K_2{K}_4{R}_{\mathrm{sy}}+K_1{K}_3^2-K_2^2+K_1,$

$F_1=M^2-{(K_3^2+1)}^2{P}_1-{(T_3^2+1)}^2{P}_2$

$F_2=2\mathrm{MN}-{(K_3^2+1)}^2{Q}_1-{(T_3^2+1)}^2{Q}_2$

$F3=N^2-{(K_3^2+1)}^2{R}_1-{(T_3^2+1)}^2{R}_2,$

$F_4=2(K_3^2+1)(T_3^2+1)$

$\mathrm{AA}=F_4^2{P}_1{P}_2-F_1^2$

$\mathrm{BB}=F_4^2(P_1{Q}_2+P_2{Q}_1)-2F_1{F}_2$

$\mathrm{CC}=F_4^2(R_1{R}_2+Q_1{Q}_2+P_2{R}_1)-2F_1{F}_3-F_2^2$

$\mathrm{DD}=F_4^2(Q_1{R}_2+R_1{Q}_2)-2F_2{F}_3$

$\mathrm{EE}=F_4^2{R}_1{R}_2-F_3^2.$

步骤S13：利用费拉里求根方法，求解步骤S12构建的一元四次方程 $\mathrm{AA}\·R_{\mathrm{ty}}^4+\mathrm{BB}\·R_{\mathrm{ty}}^3+\mathrm{CC}\·R_{\mathrm{ty}}^2+\mathrm{DD}\·R_{\mathrm{ty}}+\mathrm{EE}=0,$ 根据视角方向选择R_ty的有效解。

步骤S14：将步骤S13的解代入如下表达式，并根据视角方向选择R_tx的有效解，

$R_{\mathrm{tx}}=\frac{T_2{T}_3-T_3{T}_4{R}_{\mathrm{ty}}}{T_3^2+1}\±\frac{\sqrt{-(T_3^2+T_4^2+1)R_{\mathrm{ty}}^2+2T_2{T}_4{R}_{\mathrm{ty}}+T_1{T}_3^2-T_2^2+T_1}}{T_3^2+1}$

步骤S15：将步骤S13的解和步骤S14的解代入如下表达式，得到R_tz的解，

$R_{\mathrm{tz}}=R_{\mathrm{sz}}+\frac{\mathrm{\λ}\·R\·\mathrm{fd}}{2\·V_{\mathrm{sz}}}+\frac{V_{\mathrm{sx}}}{V_{\mathrm{sz}}}(R_{\mathrm{sx}}-R_{\mathrm{tx}})+\frac{V_{\mathrm{sy}}}{V_{\mathrm{sz}}}(R_{\mathrm{sy}}-R_{\mathrm{ty}})$

步骤S242：计算经纬度，将地固坐标系的位置矢量(x_t，y_t，z_t)转换为地理经纬度(B，L)；

步骤S243：根据步骤S23设定的地图投影模型，将地理经纬度(B，L)转换为地图坐标系下的指东指北坐标(E′，N′)；

步骤S244：根据步骤S23设定的地图图像像素分辨率，将地图投影后的指东指北坐标(E′，N′)转换为地图图像坐标(E，N)。

步骤S25：在斜距图像上，沿距离向均匀选择6组，沿方位向均匀选择4组，共6×4＝24组像素点，重复步骤S24的单像素地理编码。

步骤S26：根据多项式校正模型，用步骤S25得到的若干组斜距图像像素(i，j)和地图图像像素(E，N)，计算多项式的参数。

i＝a₀+a₁E+a₂N+a₃EN+a₄E²+a₅N²+a₆E²N+a₇EN²+a₈E³+a₉N³

j＝b₀+b₁E+b₂N+b₃EN+b₄E²+b₅N²+b₆E²N+b₇EN²+b₈E³+b₉N³

步骤S27：根据步骤S26得到的多项式参数，依次计算地图图像每个像素(E_n，N_n)对应的(i_n，j_n)，然后使用双线性插值法对图像进行重采样。

步骤S28：格式化输出。将步骤S27得到的地图图像和步骤S23设置的地图投影模型以及地图图像像素分辨率等信息，按照国际geotiff格式标准打包，形成geotiff图。

图7是地理编码后前后对比图。图7a是SAR的斜距图像(图像大小为距离向5088×方位向18432)，图7b是SAR地理编码后的地图图像(图像大小为指东9539×指北7454)。

下面说明该发明具体实施后的定位精度和运行效率。

在地图图像上选取30个特征点(图8中的十字)，与google earth上的光学地图的实际地理经纬度进行比较。

表1：地图图像定位精度表

从表1可以看出，利用本发明进行的星载SAR地理编码图像，定位精度最小11.75米，最大46.18米，平均达到26.67米。

在奔腾IV3.0GHz，内存2G，操作系统XP的台式机上，运行该星载SAR地理编码程序，所用28.5秒时间，具有执行效率高的特点。

图9给出本发明的SAR地图图像(图9a)和google earth上的光学地图(图9b)在相同地图投影坐标系下比较。从红线可以看出，两幅图像的地理位置匹配准确。

Claims (5)

1.一种宽幅星载SAR快速地理编码方法，主要步骤如下：

步骤S21：提取星载SAR实图像产品中辅助文件的相关参数；

步骤S22：拟合卫星轨道曲线，获得卫星位置速度和图像方位时刻的关系；

步骤S23：根据地图图像产品用途选择地图投影模型参数、地图投影模型变形系数和地图图像像素分辨率；

步骤S24：单像素地理编码；

步骤S25：在斜距图像上，沿距离向均匀选择M组，沿方位向均匀选择若干N组，共M×N组像素点，重复步骤S24的单像素地理编码；

步骤S26：根据多项式校正模型，用步骤S25得到的若干组斜距图像像素和地图图像像素，计算多项式的参数；

步骤S27：根据步骤S26得到的多项式参数，依次计算地图图像每个像素对应的图像像素，然后使用双线性插值法对图像进行重采样；

步骤S28：格式化输出。

2.根据权利要求1所述的宽幅星载SAR快速地理编码方法，其中，所述步骤S25中的M组多于4组，N组多于3组。

3.根据权利要求1或2所述的宽幅星载SAR快速地理编码方法，其中，所述步骤S25中的M组为6组，N组为4组。

4.根据权利要求1所述的宽幅星载SAR快速地理编码方法，其中，所述步骤S26的多项式校正模型为三次多项式：

i＝a₀+a₁E+a₂N+a₃EN+a₄E²+a₅N²+a₆E²N+a₇EN²+a₈E³+a₉N³

j＝b₀+b₁E+b₂N+b₃EN+b₄E²+b₅N²+b₆E²N+b₇EN²+b₈E³+b₉N³

5.在权利要求1所述的步骤S24单像素地理编码中，给出斜距多普勒非线性方程组的解析解，具体步骤如下：

步骤S11：获取建立斜距多普勒非线性方程组的相关参数；

步骤S12：利用地固坐标系下目标速度为零的特性，简化斜距多普勒模型构成的非线性方程组，通过分别将斜距方程和多普勒中心频率方程联立及多普勒中心频率方程和地球模型方程联立，转化成目标矢量y轴的一元四次方程；

步骤S13：利用费拉里求根方法，求解步骤S12构建的一元四次方程，并根据视角方向选择目标矢量y轴的有效解；

步骤S14：将步骤S13的解代入表达式，并根据视角方向选择目标矢量x轴的有效解；

步骤S15：将步骤S13的解和步骤S14的解代入多普勒中心频率方程，得到目标矢量z轴的解。

Images