CN101414003B - 一种基于星地坐标转换的星载sar图像地理编码方法 - Google Patents

Abstract

本发明提出了一种基于星地坐标转换的星载SAR图像地理编码方法，它是利用星地坐标之间的6个坐标系关系，巧妙的借用观测视角这一中间变量，在不需要迭代运算的前提下完成对星载SAR图像的高精度快速地理编码。传统的星载SAR图像地理编码方法大多需要进行迭代运算，处理效率低，在计算中是采用近似处理，影响定位精度，特别是卫星速度误差，它是影响传统星载SAR图像地理编码结果精度的最大因素，而本方法却不受卫星速度误差的影响，在星历参数非理想时，本方法的定位精度远远优于传统定位的方法，因而更适合于实际工程应用。本方法是一种全新的星载SAR图像定位方法，对实现星载SAR图像快速高精度地理编码具有重要意义。

Description

一种基于星地坐标转换的星载SAR图像地理编码方法

技术领域

本发明涉及一种图像地理编码方法，特别涉及一种基于星地坐标转换的星载SAR图像地理编码方法。属于图像处理技术领域。

背景技术

合成孔径雷达(SAR)卫星近些年来发展迅速，由于SAR卫星不受天气、地理、时间等因素的限制，能够对地进行全天时的观测，且具有一定的穿透力，因而被广泛的应用于军事侦察、资源探测、海洋观测、生态监测、快速救援等方面。

随着应用需求的不断增加，星载SAR图像产品的等级也不断提高，其中对星载SAR图像产品进行地理编码是一个重要的星载SAR图像处理内容，具有广泛的应用前景，比如：快速营救、精确打击、资源探测、图像拼接等。对于星载SAR图像，目前常用的星载SAR图像地理编码基本可以分为两类：

1、基于距离-多普勒(R-D)定位模型的方法

基于距离-多普勒(R-D)定位模型的方法是一种精确的SAR图像地理编码方法，其核心是利用了3个方程求解SAR图像中的像素在转动地心坐标系下的坐标，进而获取该像素位置所对应的经纬度，3个方程分别为地球椭球模型方程、斜距方程及多普勒方程。由于这3个方程表示形式较为复杂，故最后的问题归结到求解非线性方程组的问题。一般情况下，非线性方程组不具有显示解，需要利用数值分析的方法求解近似值，而当初值选取不合适时，结果可能不收敛，同时迭代处理需要消耗过多的时间，严重影响了SAR图像的地理编码效率。

2、基于近似模型的方法

正由于R-D定位模型需要迭代运算，处理效率低，在实际应用中往往对处理模型进行修改，在不影响应用的前提下对模型进行各种近似，进而获取某种条件下的显示解，这种方法虽然提高了处理效率，但处理精度受到了很大的影响，难以实现星载SAR图像的高精度地理编码。

本方法解决了上述的问题，通过建立星地坐标之间的6组转换关系，利用雷达天线视角这一中间变量，将问题转换为求解一元四次方程的问题，在不损失精度的前提下显著地提高了处理的效率，是一种全新的星载SAR图像地理编码方法，对实现星载SAR图像高精度实时地理编码具有重要意义。

发明内容

本发明的目的是提供一种基于星地坐标转换的星载SAR图像地理编码方法，它是通过星地坐标转换矩阵，将雷达天线视角作为未知量带入星地坐标转换矩阵，将其转换到转动地心坐标系中，并带入地球椭球模型方程进行求解，通过求解一元四次方程，求解出雷达天线视角大小，再将雷达天线视角带入星地坐标转换矩阵即可求解SAR图像中某像素所对应的经纬度，依次对星载SAR图像中的像素点重复上述操作，即可快速、精度的完成星载SAR图像的地理编码。

为了更好对本方法进行介绍，首先需要介绍星地6个坐标系之间的转换关系。

首先说明星地转换6个坐标系的定义

1.不转动的地心坐标系E_o

坐标原点：地球球心

Z轴：沿地球的自转轴指向正北极

X轴：在赤道平面内，指向春分点

Y轴：在赤道平面内，使该坐标系构成右手直角坐标系

不转动的地心坐标系是在不转动的地心赤道参考系(惯性参考系)中建立的直角坐标系。

2.转动的地心坐标系E_g

坐标原点：地球球心

Z轴：沿地球的自转轴指向正北极

X轴：在赤道平面内，通过格林威治子午线上半分支

Y轴：在赤道平面内，使该坐标系构成右手坐标系

3.卫星轨道平面坐标系E_v

坐标原点：卫星椭圆轨道的一个焦点(即地球球心)

Z轴：垂直于卫星轨道平面，正向指向卫星的角动量矢量方向

Y轴：在卫星平面内，正向指向近心点

X轴：在卫星轨道平面内，使该坐标系构成右手直角坐标系

4.卫星平台坐标系E_r

坐标原点：卫星质心

Z轴：垂直于卫星轨道平面，正向指向卫星的角动量矢量方向

X轴：在卫星轨道平面内，陀螺平台住纵轴方向(卫星的设计分性方向)

Y轴：在卫星轨道平面内，使该坐标系构成右手直角坐标系

5.卫星星体坐标系E_e

坐标原点：卫星质心；

X轴：沿卫星星体纵轴方向(卫星的真实飞行方向)

Y轴、Z轴：沿卫星星体的另外两个惯性主轴方向

6.天线坐标系E_a

坐标原点：天线相位中心点；

X轴：正向指向卫星的真实飞行方向；

Y轴：沿天线瞄准线，指向地球方向为正向；

Z轴：右手准则给出，使该坐标系构成右手直角坐标系。

其次说明星地六个坐标系之间的转换关系

坐标系转换示意图如上图所示，A_og为A_go的逆矩阵(A_go·A_og＝I)，以下类同。

1.转动的地心坐标系E_g/不转动的地心坐标系E_o

E_g＝A_go·E_o

不转动的地心坐标系E_o绕轴逆时针转过一个春分点的格林威治时角H_G就得到转动的地心坐标系E_g。

$A_{\mathrm{go}}=\left(\begin{array}{ccc}\cos H_G& \sin H_G& 0\\ -\sin H_G& \cos H_G& 0\\ 0& 0& 1\end{array}\right)$

2.不转动的地心坐标系E_o/轨道平面坐标系E_v

E_o＝A_ov·E_v

不转动的地心坐标系E_o经三次旋转得到轨到平面坐标系E_v。第一次，将不转动的地心坐标系绕Z轴逆时针旋转一个角Ω；第二次，将得到的坐标系再绕X轴逆时针旋转一个角度i；第三次，再将所得坐标系绕Z轴逆时针旋转一个角度ω，最后得到轨道平面坐标系E_v。

$A_{\mathrm{ov}}=\left(\begin{array}{ccc}\cos \mathrm{\Ω}& -\sin \mathrm{\Ω}& 0\\ \sin \mathrm{\Ω}& \cos \mathrm{\Ω}& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos i& -\sin i\\ 0& \sin i& \cos i\end{array}\right)\left(\begin{array}{ccc}\cos \mathrm{\ω}& -\sin \mathrm{\ω}& 0\\ \sin \mathrm{\ω}& \cos \mathrm{\ω}& 0\\ 0& 0& 1\end{array}\right)$

3.轨道平面坐标系E_v/卫星平台坐标系E_r

E_v＝A_vr·E_r

卫星平台坐标系E_r绕Z轴逆时针旋转一个角度90°+θ-γ得到卫星轨道平面坐标系E_v。其中，θ是卫星的真近心角，γ是卫星的航迹角。

真近心角θ由开普勒方程及真近心角与偏心角E的对应关系求得

$E-e\sin E=\sqrt{\frac{\mathrm{\μ}}{a^3}}(t-\mathrm{\τ})$

$\tan \frac{\mathrm{\θ}}{2}=\sqrt{\frac{1+e}{1-e}}\tan \frac{E}{2}$

航迹角γ

$\tan \mathrm{\γ}=\frac{e\sin \mathrm{\θ}}{1+e\cos \mathrm{\θ}}$ |γ|≤90°

$A_{\mathrm{vr}}=\left(\begin{array}{ccc}-\sin (\mathrm{\θ}-\mathrm{\γ})& -\cos (\mathrm{\θ}-\mathrm{\γ})& 0\\ \cos (\mathrm{\θ}-\mathrm{\γ})& -\sin (\mathrm{\θ}-\mathrm{\γ})& 0\\ 0& 0& 1\end{array}\right)$

其中，μ为引力场常数，取3.986013e14，t为星上某一时刻。

4.卫星平台坐标系E_r/卫星形体坐标系E_e

E_r＝A_re·E_e

卫星星体坐标系E_e经三次旋转得到卫星平台坐标系E_r。第一次，将卫星星体坐标系E_e绕X轴顺时针旋转一个横滚角度θ_r；第二次，将得到的坐标系绕Z轴顺时针旋转一个俯仰角度θ_p；第三次，再将所得到的坐标系绕Y轴逆时针旋转一个角度θ_y，最后得到卫星平台坐标系E_r。

$A_{\mathrm{re}}=\left(\begin{array}{ccc}\cos {\mathrm{\θ}}_y& 0& -\sin {\mathrm{\θ}}_y\\ 0& 1& 0\\ \sin {\mathrm{\θ}}_y& 0& \cos {\mathrm{\θ}}_y\end{array}\right)\left(\begin{array}{ccc}\cos {\mathrm{\θ}}_p& -\sin {\mathrm{\θ}}_p& 0\\ \sin {\mathrm{\θ}}_p& \cos {\mathrm{\θ}}_p& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos {\mathrm{\θ}}_r& -\sin {\mathrm{\θ}}_r\\ 0& \sin {\mathrm{\θ}}_r& \cos {\mathrm{\θ}}_r\end{array}\right)$

5.卫星星体坐标系E_e/天线坐标系E_a

E_e＝E_ea·E_a

天线坐标系E_a绕X轴逆时针旋转一个角度θ_L得到卫星星体坐标系E_e。

$A_{\mathrm{ea}}=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos {\mathrm{\θ}}_L& \sin {\mathrm{\θ}}_L\\ 0& -\sin {\mathrm{\θ}}_L& \cos {\mathrm{\θ}}_L\end{array}\right)$

注意：这些变换矩阵的转置矩阵就是逆矩阵。

综上所述，本发明一种基于星地坐标转换的星载SAR图像地理编码方法，其具体操作步骤如下：

步骤一：依照星载SAR图像数据产品说明书，对相关信息副产品进行读取，获取所需的星历参数，包括：偏心率e、半长轴a、升交点赤径Ω、近地点俯角ω、轨道倾角i、波前斜距R_{_min}、信号采样率f_s、脉冲重复频率prf、过近地点时刻τ、光速c、格林威治时角H_G；

本步骤的实施条件在于需要提供SAR图像相关信息的副产品，卫星星历参数都将包含在该副产品内，目前国际上标准SAR图像数据都提供相关信息的副产品。

步骤二：依照星载SAR图像数据产品说明书，读入星载SAR图像数据，计算图像每一列所对应的斜距大小，相关计算公式如下，其中，R_{_j}代表第j个距离门对应的斜距，R_{_min}代表第一个距离门对应的斜距，f_s为采样率，c为光速，j为距离门编号：

$R_{\_j}=R_{\_\min }+j\·\frac{c}{2\·f_s}$

本步骤的实施条件在于获取相关SAR图像信息，尤其是波前斜距R_{_min}、采样率f_s和SAR图像大小。

步骤三：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，依照星载SAR图像数据产品说明书，对相关信息副产品进行读取，获取其对应的星上时间，计算该时刻下卫星位置矢量在转动地心坐标系下的坐标，相关计算公式如下，其中，(x_gs，y_gs，z_gs)代表卫星位置矢量在转动地心坐标系下的坐标，A_go·A_ov分别代表了不转动地心坐标系到转动地心坐标系的转换矩阵和轨道坐标系到不转动地心坐标系的转换矩阵，θ代表真近心角，M代表平均真近心角，r代表极矢径，a代表半长轴，μ代表地球引力常数，τ代表过近心点时刻，t代表星上时间，e代表偏心率：

$\left(\begin{array}{c}{x}_{\mathrm{gs}}\\ y_{\mathrm{gs}}\\ z_{\mathrm{gs}}\end{array}\right)=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·\left(\begin{array}{c}r\cos \mathrm{\θ}\\ r\sin \mathrm{\θ}\\ 0\end{array}\right)$

$\mathrm{\θ}=M+e\·(2-\frac{1}{4}{e}^2+\frac{5}{96}{e}^4)\·\sin M+e^2\·(\frac{5}{4}-\frac{11}{24}{e}^2)\·\sin 2M$

$+e^3\·(\frac{13}{12}-\frac{43}{64}{e}^2)\·\sin 3M+\frac{103}{96}\·e^4\·\sin 4M+\frac{1097}{960}\·e^5\·\sin 5M$

$M=(t-\mathrm{\τ})\·\sqrt{(\mathrm{\μ}/a^3)}$

r＝a·(1-e²)/(1+e·cosθ)

本步骤的实施条件在于获取相关SAR图像信息，尤其是偏心率e、半长轴a、升交点赤径Ω、近地点俯角ω、轨道倾角i、脉冲重复频率prf、过近地点时刻τ、光速c、格林威治时角H_G。

步骤四：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，依照星载SAR图像数据产品说明书，对相关信息副产品进行读取，获取其对应的星上时间、三轴姿态控制角(θ_y，θ_p，θ_r)、偏心率e、半长轴a、升交点赤径Ω、近地点俯角ω、轨道倾角i、脉冲重复频率prf、过近地点时刻τ、光速c、格林威治时角H_G，计算该时刻下卫星天线相位中心位置矢量(x_e，y_e，z_c)在转动地心坐标系下的坐标(x_ge，y_ge，z_ge)，相关计算公式如下，其中，(x_e，y_e，z_e)代表卫星天线相位中心在卫星星体坐标系下的坐标，(x_ge，y_ge，z_ge)代表卫星天线相位中心在转动地心坐标系下的坐标，A_go、A_ov、A_vr、A_re分别代表不转动地心坐标系到转动地心坐标系的转换矩阵、轨道坐标系到不转动地心坐标系的转换矩阵、卫星平台坐标系到轨道坐标系的转换矩阵、卫星星体坐标系到卫星平台坐标系的转换矩阵：

$\left(\begin{array}{c}{x}_{\mathrm{ge}}\\ y_{\mathrm{ge}}\\ z_{\mathrm{ge}}\end{array}\right)=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}\·\left(\begin{array}{c}{x}_e\\ y_e\\ z_e\end{array}\right)$

本步骤的实施条件在于获取相关SAR图像信息，尤其是偏心率e、半长轴a、升交点赤径Ω、近地点俯角ω、轨道倾角i、脉冲重复频率prf、过近地点时刻τ、光速c、格林威治时角H_G，同时还需要获取卫星天线相位中心坐标。

步骤五：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，依照星载SAR图像数据产品说明书，对相关信息副产品进行读取，获取其对应的星上时间、三轴姿态控制角(θ_y，θ_p，θ_r)、偏心率e、半长轴a、升交点赤径Ω、近地点俯角ω、轨道倾角i、脉冲重复频率prf、过近地点时刻τ、光速c、格林威治时角H_G，计算该时刻下卫星星体坐标系到转动地心坐标系的转换矩阵A_ge，相关计算公式如下，其中，A_ge代表由卫星星体坐标系到转动地心坐标系的转换矩阵，A_go、A_ov、A_vr、A_re分别代表不转动地心坐标系到转动地心坐标系的转换矩阵、轨道坐标系到不转动地心坐标系的转换矩阵、卫星平台坐标系到轨道坐标系的转换矩阵、卫星星体坐标系到卫星平台坐标系的转换矩阵：

$A_{\mathrm{ge}}=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right)$

本步骤的实施条件在于获取相关SAR图像信息，尤其是偏心率e、半长轴a、升交点赤径Ω、近地点俯角ω、轨道倾角i、过近地点时刻τ、光速c、格林威治时角H_G。

步骤六：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，计算其对应的斜距，求解一元四次方程组，解出此时的雷达天线视角，相关计算公式如下，其中，R_{_j}代表第j个距离门对应的斜距，R_{_min}代表第一个距离门对应的斜距，f_s为采样率，c为光速，j为距离门编号，A代表地球长半轴，B代表地球短半轴，θ_L代表雷达天线视角：

$R_{\_j}=R_{\_\min }+j\·\frac{c}{2\·f_s}$

$({m_1}^2+{m_4}^2)x^4+2(m_1{m}_2+m_4{m}_5)x^3+({2m}_1{m}_3+{m_2}^2-{m_4}^2+{m_5}^2)x^2$

$+2(m_2{m}_3-m_4{m}_5)x^3+({m_3}^2-{m_5}^2)=0$

m₁＝R_{_j} ²·(a₁₂ ²+a₂₂ ²+k·a₃₂ ²-a₁₃ ²-a₂₃ ²-k·a₃₃ ²)

m₂＝2·R_{_j}·(a₁₂·b₁+a₂₂·b₂+k·a₃₂·b₃)

m₃＝R_{_j} ²·(a₁₃ ²+a₂₃ ²+k·a₃₃ ²)+(b₁ ²+b₂ ²+b₃ ²-A²)

m₄＝2·R_{_j} ²·(a₁₂a₁₃+a₂₂a₂₃+k·a₃₂a₃₃)

m₅＝2·R_{_j}·(a₁₃b₁+a₂₃b₂+k·a₃₃b₃)

$k={\left(\frac{A}{B}\right)}^2$

x＝cosθ_L

本步骤的实施条件在于需要考虑边界条件，避免求解值超出(0，π/2)范围。

步骤七：将求解出的雷达天线视角带入星地坐标转换矩阵，求解星载SAR图像某一像素点(i，j)在转动地心坐标系下的坐标，进而求出该像素点所对应的经纬度(Λ，Φ)，相关计算公式如下，其中，A_go、A_ov、A_vr、A_re、A_ea分别代表不转动地心坐标系到转动地心坐标系的转换矩阵、轨道坐标系到不转动地心坐标系的转换矩阵、卫星平台坐标系到轨道坐标系的转换矩阵、卫星星体坐标系到卫星平台坐标系的转换矩阵、天线坐标系到卫星星体坐标系的转换矩阵，R_{_j}代表第j个距离门对应的斜距，(x_gs，y_gs，z_gs)代表卫星位置矢量在转动地心坐标系下的坐标，(x_go，y_go，z_go)代表天线指向点在转动地心坐标系下的坐标：

$\left(\begin{array}{c}{x}_{\mathrm{go}}\\ y_{\mathrm{go}}\\ z_{\mathrm{go}}\end{array}\right)=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}\·A_{\mathrm{ea}}\left(\begin{array}{c}0\\ R_{\_j}\\ 0\end{array}\right)+\left(\begin{array}{c}{x}_{\mathrm{gs}}\\ y_{\mathrm{gs}}\\ z_{\mathrm{gs}}\end{array}\right)+A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}\left(\begin{array}{c}{x}_e\\ y_e\\ z_e\end{array}\right)$

$\tan \mathrm{\Λ}=\frac{y_{\mathrm{go}}}{z_{\mathrm{go}}}$

$\tan \mathrm{\Φ}=\frac{z_{\mathrm{go}}}{\sqrt{{x_{\mathrm{go}}}^2+{y_{\mathrm{go}}}^2+{z_{\mathrm{go}}}^2}}$

本步骤的实施条件在于获取相关SAR图像信息，尤其是偏心率e、半长轴a、升交点赤径Ω、近地点俯角ω、轨道倾角i、过近地点时刻τ、光速c、格林威治时角H_G。

步骤八：每完成一个像素点的地理编码，计算i-M和j-N的值，其中M和N分别代表SAR图像的行数和列数，判断此时j-N是否小于0，如果小于0，则从步骤六进行重复操作；如果大等于0，则判断此时i-M是否小于0，如果小于0，则转入图像下一行处理，从步骤2进行重复操作，否则星载SAR图像编码完成。本步骤的实施条件在于硬件设施所能开辟的内存大小不能小于SAR图像大小。

本发明一种基于星地坐标转换的星载SAR图像地理编码方法的优点在于：

(1)本方法能够利用星地转换关系，通过求解雷达天线视角进一步求解图像像素所对应的经纬度，避免了直接求解所需的迭代运算，显著提高了处理效率。

(2)本方法能够利用星地转换关系，通过求解一元四次方程来得到精确的地理编码结果，显著提高了处理结果的精度。

(3)本方法可避免迭代运算，因而避免了传统方法中由于采用迭代运算而导致结果不收敛的缺陷。

(4)本方法计算结果精度不受卫星速度影响。传统地理编码方法中，卫星速度是影响定位精度的最大因素，而采用本方法将不受卫星速度误差的影响，即当星历参数不理想时，采用本方法所得到的地理编码误差将远小于传统方法的地理编码误差。

(5)本方法计算效率不受初始迭代值影响。传统地理编码方法中，其编码效率受初始迭代值影响，当初始迭代值选取不合适时，其编码效率将很低，甚至出现结果不收敛、效率为0的状况，而本方法可避免这种状况，在任何初始条件下都具有更高的处理效率。

(6)本方法对星载SAR图像采用逐行处理方式进行图像地理编码，由于星载SAR图像行与行之间的处理过程完全独立，故可对星载SAR图像进行分块处理，进一步提高处理效率。

附图说明

图1基于星地坐标转换的星载SAR图像地理编码方法的流程图

图2SAR图像说明示意图

图3计算卫星位置矢量在转动地心坐标系下坐标的流程图

图4计算卫星天线相位中心位置矢量在转动地心坐标系的流程图

图5计算转动地心坐标系到卫星星体坐标系的转换矩阵的流程图

图6计算求解雷达天线视角的流程图

图7仿真场景摆放示意图

图中符号说明如下：

A_go 不转动地心坐标系到转动地心坐标系的转换矩阵

A_ov 轨道坐标系到不转动地心坐标系的转换矩阵

A_vr 卫星平台坐标系到轨道坐标系的转换矩阵

A_re 卫星星体坐标系到卫星平台坐标系的转换矩阵

A_ea 天线坐标系到卫星星体坐标系的转换矩阵

R_{_min} 波前距离，即第一个距离门所对应的斜距

R_{_j}   第j个距离门所对应的斜距

t₁  第1行SAR图像数据对应的星上时间

t_i  第i行SAR图像数据对应的星上时间

(x_go，y_go，z_go)　目标像素点在转动地心坐标系下的坐标

具体实施方式

见图1、图2、图3、图4、图5、图6、图7所示，

本发明一种基于星地坐标转换的星载SAR图像地理编码方法，以一幅大小为2048x2048的SAR图像为例，该方法具体步骤如下：

步骤一：读入卫星星历参数，包括：偏心率e＝0.003、半长轴a＝7071140.0m、升交点赤径Ω＝69.517572°、近地点俯角ω＝0.0°、轨道倾角i＝99.038314°、波前斜距R_{_min}＝787796.306m、信号采样率f_s＝75.0MHz、脉冲重复频率prf＝2300.0Hz、过近地点时刻τ＝-9.6465s、光速c＝3.0e8m/s、格林威治时角H_G＝0.0s。

本步骤的实施条件在于需要提供SAR图像相关信息的副产品，卫星星历参数都将包含在该副产品内，目前国际上标准SAR图像数据都提供相关信息的副产品。

步骤二：读入星载SAR图像，计算图像每一列所对应的斜距大小。

星载SAR成像机制在距离向是一种距离分辨机制，其每一距离门(即星载SAR图像每一列，如图2所示)中的像素对应相同的斜距，即像素所代表的目标距天线相位中心具有相同的距离。星载SAR图像第j列所代表的斜距可按下式进行计算，其中，R_{_min}代表第一个距离门所对应的斜距(即图像的第一列)，其中，f_s＝75.0MHz，R_{_min}＝787796.306m，

$R_{\_j}=R_{\_\min }+j\·\frac{c}{2\·f_s}---\left(1\right)$

本步骤的实施条件在于获取相关SAR图像信息，尤其是波前斜距R_{_min}＝787796.306m、采样率f_s＝75.0MHz和SAR图像大小。

步骤三：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，计算该时刻下卫星位置矢量在转动地心坐标系下的坐标，流程图见图3，以i＝1024，j＝1024为例，

转动地心坐标系下的坐标(x_gs，y_gs，z_gs)求取方法如下式：

$\left(\begin{array}{c}{x}_{\mathrm{gs}}\\ y_{\mathrm{gs}}\\ z_{\mathrm{gs}}\end{array}\right)=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·\left(\begin{array}{c}r\cos \mathrm{\θ}\\ r\sin \mathrm{\θ}\\ 0\end{array}\right)---\left(2\right)$

其中，r为极矢径，θ为真近心角，

$\mathrm{\θ}=M+e\·(2-\frac{1}{4}{e}^2+\frac{5}{96}{e}^4)\·\sin M+e^2\·(\frac{5}{4}-\frac{11}{24}{e}^2)\·\sin 2M+$

                                          (3)

其中，M代表真近心角， $M=(t-\mathrm{\τ})\·\sqrt{(\mathrm{\μ}/a^3)}$ ，μ为引力常数，取3.986013e14。r＝a·(1-e²)/(1+e·cosθ)＝7049927.6994m(4)

本步骤的实施条件在于获取相关SAR图像信息，尤其是偏心率e＝0.003、半长轴a＝7071140.0m、升交点赤径Ω＝69.517572°、近地点俯角ω＝0.0°、轨道倾角i＝99.038314°、波前斜距R_{_min}＝787796.306m、信号采样率f_s＝75.0MHz、脉冲重复频率prf＝2300.0Hz、过近地点时刻τ＝-9.6465s、光速c＝3.0e8m/s、格林威治时角H_G＝0.0s。

步骤四：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，计算该时刻下卫星天线相位中心在转动地心坐标系下的坐标，流程图见图4。

首先读取卫星天线相位中心在卫星星体坐标系下的坐标(x_e，y_e，z_e)及三轴姿态控制角(θ_y，θ_p，θ_r)，一般情况下，该坐标取(0，0，0)，并利用坐标转换矩阵A_go、A_ov、A_vr、A_re将其转换到转动地心坐标系下的坐标(x_ge，y_ge，z_ge)，以i＝1024，j＝1024为例，则(x_e，y_e，z_e)＝(0.0，0.0，0.0)，(θ_y，θ_p，θ_r)＝(3.79°，0.0°，0.0°)

$\left(\begin{array}{c}{x}_{\mathrm{ge}}\\ y_{\mathrm{ge}}\\ z_{\mathrm{ge}}\end{array}\right)=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{ve}}\·A_{\mathrm{re}}\·\left(\begin{array}{c}{x}_e\\ y_e\\ z_e\end{array}\right)---\left(5\right)$

本步骤的实施条件在于获取相关SAR图像信息，尤其是偏心率e＝0.003、半长轴a＝7071140.0m、升交点赤径Ω＝69.517572°、近地点俯角ω＝0.0°、轨道倾角i＝99.038314°、波前斜距R_{_min}＝787796.306m、信号采样率f_s＝75.0MHz、脉冲重复频率prf＝2300.0Hz、过近地点时刻τ＝-9.6465s、光速c＝3.0e8m/s、格林威治时角H_G＝0.0s，同时还需要获取卫星天线相位中心坐标(x_e，y_e，z_e)＝(0.0，0.0，0.0)及三轴姿态控制角(θ_y，θ_p，θ_r)＝(3.79°，0.0°，0.0°)。

步骤五：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，计算该时刻下卫星星体坐标系到转动地心坐标系的转换矩阵，流程图见图5，以i＝1024，j＝1024为例，

为了便于表示，将卫星星体坐标系到转动地心坐标系的转换矩阵用A_ge表示，则

$A_{\mathrm{ge}}=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right)---\left(6\right)$

本步骤的实施条件在于获取相关SAR图像信息，尤其是偏心率e＝0.003、半长轴a＝7071140.0m、升交点赤径Ω＝69.517572°、近地点俯角ω＝0.0°、轨道倾角i＝99.038314°、波前斜距R_{_min}＝787796.306m、信号采样率f_s＝75.0MHz、脉冲重复频率prf＝2300.0Hz、过近地点时刻τ＝-9.6465s、光速c＝3.0e8m/s、格林威治时角H_G＝0.0s。

步骤六：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，根据公式1计算其对应的斜距，求解一元四次方程组，解出此时的雷达天线视角，流程图见图6，以i＝1024，j＝1024为例，

在转动地球坐标系下，卫星位置矢量、天线相位中心位置矢量、目标到天线相位中心的距离矢量、目标点矢量形成了闭合的四边形，故其矢量和为0，故有，

$\left(\begin{array}{c}{x}_{\mathrm{go}}\\ y_{\mathrm{go}}\\ z_{\mathrm{go}}\end{array}\right)=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}\·A_{\mathrm{ea}}\left(\begin{array}{c}0\\ R_{\_j}\\ 0\end{array}\right)+\left(\begin{array}{c}{x}_{\mathrm{gs}}\\ y_{\mathrm{gs}}\\ z_{\mathrm{gs}}\end{array}\right)+A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}\left(\begin{array}{c}{x}_e\\ y_e\\ z_e\end{array}\right)---\left(7\right)$

其中，对于SAR图像每一行而言，后两项为固定值，为了便于分析，作如下进一步化简，

$\left(\begin{array}{c}{x}_{\mathrm{go}}\\ y_{\mathrm{go}}\\ z_{\mathrm{go}}\end{array}\right)=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right)\·\left(\begin{array}{ccc}1& 0& 0\\ 0& {\cos \mathrm{\θ}}_L& {\sin \mathrm{\θ}}_L\\ 0& {-\sin \mathrm{\θ}}_L& \cos {\mathrm{\θ}}_L\end{array}\right)\left(\begin{array}{c}0\\ R_{\_j}\\ 0\end{array}\right)+\left(\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right)---\left(8\right)$

其中，θ_L即为雷达天线视角，(x_go，y_go，z_go)为目标点在转动地心坐标系下的坐标。

由于(x_go，y_go，z_go)为转动地心坐标系下的坐标，故其满足地球椭球方程，

$\frac{{x_{\mathrm{go}}}^2+{y_{\mathrm{go}}}^2}{A^2}+\frac{{z_{\mathrm{go}}}^2}{B^2}=1---\left(9\right)$

A、B分别代表了地球的长、短半轴，分别取A＝6378140.0mB＝6356755.0m，将(x_go，y_go，z_go)用θ_L表示，并代入表达式(8)，则可化简为关于cosθ_L的一元四次方程，

(m₁ ²+m₄ ²)x⁴+2(m₁m₂+m₄m₅)x³+(2m₁m₃+m₂ ²-m₄ ²+m₅ ²)x²

                                            (10)

+2(m₂m₃-m₄m₅)x³+(m₃ ²-m₅ ²)＝0

其中，

m₁＝R_{_j} ²·(a₁₂ ²+a₂₂ ²+k·a₃₂ ²-a₁₃ ²-a₂₃ ²-k·a₃₃ ²)

m₂＝2·R_{_j}·(a₁₂·b₁+a₂₂·b₂+k·a₃₂·b₃)

m₃＝R_{_j} ²·(a₁₃ ²+a₂₃ ²+k·a₃₃ ²)+(b₁ ²+b₂ ²+b₃ ²-A²)

m₄＝2·R_{_j} ²·(a₁₂a₁₃+a₂₂a₂₃+k·a₃₂a₃₃)

m₅＝2·R_{_j}·(a₁₃b₁+a₂₃b₂+k·a₃₃b₃)

$k={\left(\frac{A}{B}\right)}^2$

x＝cosθ_L

则问题最终转换为求解一元四次方程，其求解过程可参考各种数学手册即可，最终解出θ_L＝30.0°。

本步骤的实施条件在于需要考虑边界条件，避免求解值超出(0，π/2)范围。

步骤七：将求解出的雷达天线视角带入星地坐标转换矩阵，求解星载SAR图像某一像素点(i，j)在转动地心坐标系下的坐标。

此时由于cosθ_L和sinθ_L已经求出，将其代入公式(7)即可求解出(x_go，y_go，z_go)，则其对应的经度为：

纬度为：

本步骤的实施条件在于获取相关SAR图像信息，尤其是偏心率e＝0.003、半长轴a＝7071140.0m、升交点赤径Ω＝69.517572°、近地点俯角ω＝0.0°、轨道倾角i＝99.038314°、波前斜距R_{_min}＝787796.306m、信号采样率f_s＝75.0MHz、脉冲重复频率prf＝2300.0Hz、过近地点时刻τ＝-9.6465s、光速c＝3.0e8m/s、格林威治时角H_G＝0.0s。

步骤八：判断此时j是否小于图像列数，如果小于图像列数，则从步骤六进行重复操作；如果大于图像列数，则判断此时i是否小于图像行数，如果小于图像行数，则转入图像下一行处理，从步骤二进行重复操作，否则星载SAR图像编码完成。

本步骤的实施条件在于硬件设施所能开辟的内存大小不能小于SAR图像大小。

为了说明该方法的有效性，进行如下仿真试验，仿真参数如表1。点阵摆放如图7所示，在地球表面摆放3x3的点阵共9个点，每一个点所代表的经纬度如表2：

  表1 仿真参数

 

仿真参数	所选参数值
		地球长半轴(m)	6378140.0
地球短半轴(m)	6356755.0
		地球平均半径(m)	6371140.0
半长轴(m)	7071140.0
		近地点幅角(度)	0.0
轨道倾角(度)	99.038314
		升交点赤经(度)	69.517572
偏心率	0.003

      表2  点阵摆放经纬度(单位：度)

 

(0.2027425，73.0686908)	(0.2027425，73.0715556)	(0.2027425，73.0744204)
			(0.2056073，73.0686908)	(0.2056073，73.0715556)	(0.2056073，73.0744204)
(0.2084721，73.0686908)	(0.2084721，73.0715556)	(0.2084721，73.0744204)

(注：第一个参数代表纬度，第二个参数为经度，后面表格也遵循此标准)

 表3  对3x3点阵进行地理编码后的经纬度(单位：度)

 

(0.2027394，73.0686905)	(0.2027444，73.0715576)	(0.2027406，73.0744167)
			(0.2056027，73.0686900)	(0.2056076，73.0715568)	(0.2056055，73.0744235)
(0.2084660，73.0686895)	(0.2084708，73.0715559)	(0.2084686，73.0744222)

          表4  地理编码误差大小(单位：m)

  表5  两种地理编码方法对比结果

 

	基于星地坐标转换的地理编码方法	基于R-D模型的地理编码方法
			所用时间(s)	67	405
地理编码结果误差的均值(m)	(-0.248，0.028)	(-0.275，-0.102)
			地理编码结果误差的标准差(m)	(0.255，0.214)	(-0.251，0.220)

从表3中可以看出，用基于星地坐标转换的方法对目标进行地理编码，其结果可以精确到小数点第6位，具有非常高的精度；从表4中可以看出，每一个点的地理编码误差不超过0.4m，也验证本方法的精确性；从表5中可以看出，对一个2048x2048的SAR图像进行地理编码，达到同样数量极的地理编码精度，基于R-D模型的地理编码方法消耗的时间约为基于星地坐标转换的地理编码方法的6～7倍，且SAR图像越大，两种方法消耗的时间之比也将越来越大，这也验证了本方法的快速性。

Claims (1)

1.一种基于星地坐标转换的星载SAR图像地理编码方法，其特征在于，它包括如下步骤：

步骤一：依照星载SAR图像数据产品说明书，对相关信息副产品进行读取，获取所需的星历参数，包括：偏心率e、半长轴a、升交点赤径Ω、近地点俯角ω、轨道倾角i、第一个距离门对应的斜距R_{_min}、信号采样率f_s、脉冲重复频率prf、过近地点时刻τ、光速c、格林威治时角H_G；

步骤二：依照星载SAR图像数据产品说明书，读入星载SAR图像数据，计算图像每一列所对应的斜距大小，相关计算公式如下，其中，R_{_j}代表第j个距离门对应的斜距，R_{_min}代表第一个距离门对应的斜距，f_s为采样率，c为光速，j为距离门编号：

$R_{\_j}=R_{\_\min }+j\·\frac{c}{2\·f_s}$

步骤三：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，依照星载SAR图像数据产品说明书，对相关信息副产品进行读取，获取其对应的星上时间，计算该时刻下卫星位置矢量在转动地心坐标系下的坐标，相关计算公式如下，其中，(x_gs，y_gs，z_gs)代表卫星位置矢量在转动地心坐标系下的坐标，A_go·A_ov分别代表了不转动地心坐标系到转动地心坐标系的转换矩阵和轨道坐标系到不转动地心坐标系的转换矩阵，θ代表真近心角，M代表平均真近心角，r代表极矢径，a代表半长轴，μ代表地球引力常数，τ代表过近地点时刻，t代表星上时间，e代表偏心率：

$\left(\begin{array}{c}{x}_{\mathrm{gs}}\\ y_{\mathrm{gs}}\\ z_{\mathrm{gs}}\end{array}\right)=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·\left(\begin{array}{c}r\cos \mathrm{\θ}\\ r\sin \mathrm{\θ}\\ 0\end{array}\right)$

$\mathrm{\θ}=M+e\·(2-\frac{1}{4}{e}^2+\frac{5}{96}{e}^4)\·\sin M+e^2\·(\frac{5}{4}-\frac{11}{24}{e}^2)\·\sin 2M$

$+e^3\·(\frac{13}{12}-\frac{43}{64}{e}^2)\·\sin 3M+\frac{103}{96}\·e^4\·\sin 4M+\frac{1097}{960}\·e^5\·\sin 5M$

$M=(t-\mathrm{\τ})\·\sqrt{(\mathrm{\μ}/a^3)}$

r＝a·(1-e²)/(1+e·cosθ)

步骤四：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，依照星载SAR图像数据产品说明书，对相关信息副产品进行读取，获取其对应的星上时间、三轴姿态控制角(θ_y，θ_p，θ_r)、偏心率e、半长轴a、升交点赤径Ω、近地点俯角ω、轨道倾角i、脉冲重复频率prf、过近地点时刻τ、光速c、格林威治时角H_G，计算该时刻下卫星天线相位中心位置矢量(x_e，y_e，z_e)在转动地心坐标系下的坐标(x_ge，y_ge，z_ge)，相关计算公式如下，其中，(x_e，y_e，z_e)代表卫星天线相位中心在卫星星体坐标系下的坐标，(x_ge，y_ge，z_ge)代表卫星天线相位中心在转动地心坐标系下的坐标，A_go、A_ov、A_vr、A_re分别代表不转动地心坐标系到转动地心坐标系的转换矩阵、轨道坐标系到不转动地心坐标系的转换矩阵、卫星平台坐标系到轨道坐标系的转换矩阵、卫星星体坐标系到卫星平台坐标系的转换矩阵：

$\left(\begin{array}{c}{x}_{\mathrm{ge}}\\ y_{\mathrm{ge}}\\ z_{\mathrm{ge}}\end{array}\right)=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}\·\left(\begin{array}{c}{x}_e\\ y_e\\ z_e\end{array}\right)$

步骤四的实施条件在于获取相关SAR图像信息以及卫星天线相位中心坐标(x_e，y_e，z_e)＝(0.0，0.0，0.0)及三轴姿态控制角(θ_y，θ_p，θ_r)＝(3.79°，0.0°，0.0°)；

步骤五：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，依照星载SAR图像数据产品说明书，对相关信息副产品进行读取，获取其对应的星上时间、三轴姿态控制角(θ_y，θ_p，θ_r)、偏心率e、半长轴a、升交点赤径Ω、近地点俯角ω、轨道倾角i、脉冲重复频率prf、过近地点时刻τ、光速c、格林威治时角H_G，计算该时刻下卫星星体坐标系到转动地心坐标系的转换矩阵A_ge，相关计算公式如下，其中，A_ge代表由卫星星体坐标系到转动地心坐标系的转换矩阵，A_go、A_ov、A_vr、A_re分别代表不转动地心坐标系到转动地心坐标系的转换矩阵、轨道坐标系到不转动地心坐标系的转换矩阵、卫星平台坐标系到轨道坐标系的转换矩阵、卫星星体坐标系到卫星平台坐标系的转换矩阵：

$A_{\mathrm{ge}}=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right)$

步骤六：对星载SAR图像中某一像素点(i，j)，其中i为行数、j为列数，计算其对应的斜距，求解一元四次方程组，解出此时的雷达天线视角，相关计算公式如下，其中，R_{_j}代表第j个距离门对应的斜距，R_{_min}代表第一个距离门对应的斜距，f_s为采样率，c为光速，j为距离门编号，A代表地球长半轴，B代表地球短半轴，θ_L代表雷达天线视角：

$R_{\_j}=R_{\_\min }+j\·\frac{c}{2\·f_s}$

(m₁ ²+m₄ ²)x⁴+2(m₁m₂+m₄m₅)x³+(2m₁m₃+m₂ ²-m₄ ²+m₅ ²)x²

+2(m₂m₃-m₄m₅)x³+(m₃ ²-m₅ ²)＝0

m₁＝R_{_j} ²·(a₁₂ ²+a₂₂ ²+k·a₃₂ ²-a₁₃ ²-a₂₃ ²-k·a₃₃ ²)

m₂＝2·R_{_j}·(a₁₂·b₁+a₂₂·b₂+k·a₃₂·b₃)

m₃＝R_{_j} ²·(a₁₃ ²+a₂₃ ²+k·a₃₃ ²)+(b₁ ²+b₂ ²+b₃ ²-A²)

m₄＝2·R_{_j} ²·(a₁₂a₁₃+a₂₂a₂₃+k·a₃₂a₃₃)

m₅＝2·R_{_j}·(a₁₃b₁+a₂₃b₂+k·a₃₃b₃)

$k={\left(\frac{A}{B}\right)}^2$

x＝cosθ_L

参数b₁、b₂、b₃由下列公式推导：

$\left(\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right)=\left(\begin{array}{c}{x}_{\mathrm{gs}}\\ y_{\mathrm{gs}}\\ z_{\mathrm{gs}}\end{array}\right)+A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}\left(\begin{array}{c}{x}_e\\ y_e\\ z_e\end{array}\right)$

步骤六的实施条件是考虑边界条件，避免cosθ_L的求解值超出(0，π/2)范围；

步骤七：将求解出的雷达天线视角带入星地坐标转换矩阵，求解星载SAR图像某一像素点(i，j)在转动地心坐标系下的坐标，像素点(i，j)在转动地心坐标系下的坐标即天线指向点在转动地心坐标系下的坐标(x_go，y_go，z_go)，进而求出该像素点所对应的经纬度(Λ，Φ)，相关计算公式如下，其中，A_go、A_ov、A_vr、A_re、A_ea分别代表不转动地心坐标系到转动地心坐标系的转换矩阵、轨道坐标系到不转动地心坐标系的转换矩阵、卫星平台坐标系到轨道坐标系的转换矩阵、卫星星体坐标系到卫星平台坐标系的转换矩阵、天线坐标系到卫星星体坐标系的转换矩阵，R_{_j}代表第j个距离门对应的斜距，(x_gs，y_gs，z_gs)代表卫星位置矢量在转动地心坐标系下的坐标，(x_go，y_go，z_go)代表天线指向点在转动地心坐标系下的坐标：

$\left(\begin{array}{c}{x}_{\mathrm{go}}\\ y_{\mathrm{go}}\\ z_{\mathrm{go}}\end{array}\right)=A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}\·A_{\mathrm{ea}}\left(\begin{array}{c}0\\ R_{\_j}\\ 0\end{array}\right)+\left(\begin{array}{c}{x}_{\mathrm{gs}}\\ y_{\mathrm{gs}}\\ z_{\mathrm{gs}}\end{array}\right)+A_{\mathrm{go}}\·A_{\mathrm{ov}}\·A_{\mathrm{vr}}\·A_{\mathrm{re}}\left(\begin{array}{c}{x}_e\\ y_e\\ z_e\end{array}\right)$

$\tan \mathrm{\Λ}=\frac{y_{\mathrm{go}}}{z_{\mathrm{go}}}$

$\tan \mathrm{\Φ}=\frac{z_{\mathrm{go}}}{\sqrt{{x_{\mathrm{go}}}^2+{y_{\mathrm{go}}}^2+{z_{\mathrm{go}}}^2}}$

步骤八：每完成一个像素点的地理编码，计算i-M和j-N的值，其中M和N分别代表SAR图像的行数和列数，判断此时j-N是否小于0，如果小于0，则从步骤六进行重复操作；如果大等于0，则判断此时i-M是否小于0，如果小于0，则转入图像下一行处理，从步骤2进行重复操作，否则星载SAR图像编码完成。

Images