CN103868502B - 一种大椭圆卫星地球背景成像区域昼夜分布的确定方法 - Google Patents

Abstract

本发明公开了一种大椭圆遥感卫星地球背景成像区域昼夜分布的定量确定方法，包括：建立地心固联坐标系、传感器坐标系和图像坐标系，构建图像坐标系与地心固联坐标系之间的转换关系；确定当前时刻t，获取当前时刻下卫星、太阳在地心固联坐标系中的坐标，以及传感器参数；判断当前时刻下，成像面上点映射到物面上的点是否为地球上的点，并计算该点在上地球上的点，计算地球上的点在当前时刻下的太阳高度角，并根据太阳高度角判断这个点的昼夜情况；对成像面上的所有点都做上述处理，判断成像面上各点在地球上对应点的昼夜情况。本发明建立在坐标系变换和基本公式的基础上，建立了数学模型，可以定量确定大椭圆遥感卫星图像地球背景昼夜分布情况。

Description

一种大椭圆卫星地球背景成像区域昼夜分布的确定方法

技术领域

本发明属于计算机视觉与航空航天技术交叉的科学技术领域，更具体地，涉及一种大椭圆卫星地球背景成像区域昼夜分布的确定方法。

背景技术

遥感技术在灾害监测、矿物探测、环境监测等领域具有明显的优势和广泛的应用前景。但是，遥感探测特别是航天遥感成本高，传感器的制造周期长，阻碍目前遥感成像技术的不断发展。随着计算机技术的快速发展，特别是虚拟现实技术的发展，遥感成像仿真技术极为有效的解决这类问题。自然界中的一切物体，都会反射太阳辐射或发射自身热辐射。遥感成像技术是研究如何探测的反射和发射辐射，并将辐射能量通过电子系统转换为图像信号的科学技术。遥感成像仿真是基于遥感成像机理，在某种成像条件（大气条件、传感器平台）下，模拟场景中的地物如何反射太阳辐射和发射自身热辐射，将辐射能量映射到图像信号，最终生成遥感图像。显然，太阳的照射对于遥感图像仿真起着重要的作用。为了更好地对地球进行遥感探测，我们需要研究不同地球背景的辐射统计特性随昼夜时相的变化关系。因此，研究遥感图像地球背景成像区域昼夜分布对调整优化遥感器系统参数和数据获取方案具有重要的现实意义。

大椭圆轨道，与地球静止轨道不同，远地点距离地面高度几万公里，近地点距离地面几百公里。地球静止轨道（GEO）卫星与地面相对静止，固定在赤道上空，其星下点轨迹是一个点，运动周期与地球自转周期相等，它的轨道离心率和轨道倾角均为零，距地面高度约为3.6万公里。地球静止轨道卫星携带的传感器始终指向地球上固定一片区域，所以其地球背景成像区域昼夜分布只需考虑地球自转和地球与太阳之间相对位置的变化。

大椭圆轨道（HEO）卫星大不相同，卫星与地球相对位置实时变化。大椭圆轨道卫星所携带的传感器与地球静止轨道卫星的传感器工作模式不同，大椭圆卫星要求传感器相机始终指向北极点，所以大椭圆卫星的传感器指向地球上的区域也是实时变化的。

综上所述，对于大椭圆轨道卫星而言，地球自转，地球、太阳和卫星三者相对位置的实时变化以及传感器始终指向北极点等因素都是造成地球背景成像区域昼夜分布实时变化的原因。因此，确定大椭圆卫星地球背景成像区域的昼夜分布是非常复杂的。国内外没有相关文献提到这种确定方法。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种大椭圆遥感卫星地球背景成像区域昼夜分布的定量确定方法，包括：

（1）建立地心固联坐标系、传感器坐标系和图像坐标系，获取当前时刻t的太阳位置、卫星位置以及传感器参数，确定当前时刻t下成像面到物面的转换关系；

（2）判断当前时刻t下，成像面上点（u，v）映射到物面上的点是否为地球上的点，如果是地球上的点，则根据成像面到物面的转换关系求像素点（u，v）在物面上的点（X，Y，Z），如果不是地球上的点，则处理成像面上下一个像素点；

（3）计算像素点（u，v）在物面上的点（X，Y，Z）在当前时刻t的太阳高度角，并根据该太阳高度角的正负来判断这个点的昼夜情况，即太阳高度角为正，则为白天，否则为晚上；

（4）遍历地球背景成像区域的成像面上的所有像素点，计算其在物面上的对应点，并判断物面上的每个对应点的昼夜情况，得到地球背景成像区域的昼夜分布图。

优先地，所述步骤（1）中各坐标系定义如下：

地心固联坐标系：原点O为地心，X轴指向本初子午线，Z轴指向正北，按照右手螺旋法则确定Y轴方向，地心固联坐标系下的点坐标为(X,Y,Z)；

传感器坐标系：原点O_c与卫星位置O_s重合，传感器指向Z_c的方向为卫星位置指向北极点的方向，Z_c方向与卫星与太阳连线方向叉乘得到Y_c轴方向，按照右手螺旋法则确定X_c轴；

图像坐标系：包括图像像素坐标系和图像物理坐标系，图像像素坐标系的原点O₁定义为卫星传感器光轴与图像平面的交点，其横轴u和纵轴v分别与传感器坐标系的X_c和Y_c平行，其像素坐标以(u,v)表示，(u₀,v₀)为原点O₁在图像像素坐标系下的坐标；图像物理坐标系与图像像素坐标系平行，其像素坐标以(x,y)表示，图像物理坐标系中每一个像素点在x轴与y轴方向上的物理尺寸为dx、dy；

成像面中的点(u,v)与其在物面中的对应点(X,Y,Z)满足如下关系：

$\left\{\begin{array}{c}{X}^2+Y^2+Z^2=R^2\\ \frac{X-x_s}{m}=\frac{Y-y_s}{n}=\frac{Z-z_s}{l}\end{array},\right.$

$\stackrel{\→}{l}=(m,n,l),$

$\left(\begin{array}{c}\stackrel{\→}{l}\\ 1\end{array}\right)=M_1^{-1}\left(\begin{array}{c}\stackrel{\→}{l_c}\\ 1\end{array}\right)=M_1^{-1}{M}^{\′}\left(\begin{array}{c}u\\ v\\ 1\end{array}\right),$

其中，M₁为地心固联坐标系到传感器坐标系的转换矩阵，

$M_1=T({-x}_s,-y_s,-z_s)\left[\begin{array}{cccc}{\mathrm{\μ}}_x& {\mathrm{\μ}}_y& {\mathrm{\μ}}_z& 0\\ {\mathrm{\υ}}_x& {\mathrm{\υ}}_y& {\mathrm{\υ}}_z& 0\\ {\mathrm{\η}}_x& {\mathrm{\η}}_y& {\mathrm{\η}}_z& 0\\ 0& 0& 0& 1\end{array}\right],$

$T({-x}_s,{-y}_s,{-z}_s)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ -x_s& -y_s& {-z}_s& 1\end{array}\right],$

$\left\{\begin{array}{c}{\mathrm{\η}}_x=-x_s\\ {\mathrm{\η}}_y=-y_s\\ {\mathrm{\η}}_z=R-z_s\end{array}\right.,$

$\left\{\begin{array}{c}{\mathrm{\υ}}_x=(z_s-z_0)y_s-(R-z_s)(y_0-y_s)\\ {\mathrm{\υ}}_y=(z_0-z_s)x_s+(R-z_s)(x_0-x_s)\\ {\mathrm{\υ}}_z=-x_s,y_0+y_s{x}_0\end{array}\right.,$

$\left\{\begin{array}{c}{\mathrm{\μ}}_x=v_y{\mathrm{\η}}_z-v_z{\mathrm{\η}}_y\\ {\mathrm{\μ}}_y=v_z{\mathrm{\η}}_x-v_x{\mathrm{\η}}_z\\ {\mathrm{\μ}}_z=v_x{\mathrm{\η}}_y-v_y{\mathrm{\η}}_x\end{array}\right.,$

(x_s,y_s,z_s)为传感器坐标系原点O_s在地心固联坐标系中的坐标，(x₀,y₀,z₀)为太阳在地心固联坐标系中的坐标，R为地球半径；

其中 $M^{\′}=\left[\begin{array}{ccc}\tan \left(\mathrm{dixAngleCol}\right)& 0& 0\\ 0& \tan \left(\mathrm{dixAngleRow}\right)& 0\\ 0& 0& 1\\ 0& 0& 1\end{array}\right],$

dixAngleRow,dixAngleCol为相机的行列视场像元角。

进一步地，所述步骤（2）中判断成像面上点（u，v）映射到物面上的点是否为地球上的点，以及计算像面上的像素点（u，v）在物面上地球上的点（X，Y，Z），具体为：

（2.1）根据式 $\left(\begin{array}{c}\stackrel{\→}{l}\\ 1\end{array}\right)=M_1^{-1}{M}^{\′}\left(\begin{array}{c}u\\ v\\ 1\end{array}\right)$ 求得向量 $\stackrel{\→}{l}=(m,n,l);$

（2.2）判断式 $\left\{\begin{array}{c}{X}^2+Y^2+Z^2=R^2\\ \frac{X-x_s}{m}=\frac{Y-y_s}{n}=\frac{Z-z_s}{l}\end{array}\right.$ 是否有解，令则X＝mt+x_s,Y＝nl+y_s,Z＝lt+z_s，带入方程X²+Y²+Z²＝R²得一元二次方程组：(m+n+l)t²+2×(mx_s+ny_s+lz_s)×t+x² _s+y² _s+z² _s-R²＝0；当Δ＜0时，方程组无解，即点（u，v）映射到物面上的点不是地球上的点，，当Δ≥0时，方程组有解，点（u，v）映射到物面上的点不是地球上的点；其中Δ＝(2×(mx_s+ny_s+lz_s))²-4×(m+n+l)×(x² _s+y² _s+z² _s-R²)

（2.3）求解（2.2）中的一元二次方程组，得到成像面上的像素点（u，v）在物面上地球上的点（X，Y，Z）。

进一步地，所述步骤（3）具体包括：

（3.1）根据该物面上的点（X，Y，Z）及当前时刻t，确定地理纬度，太阳赤纬角δ，太阳时角τ；

（3.2）计算

（3.3）根据该物面上的点的sinh判断该点的昼夜情况，当sinh＞0时，为白天，当sinh＜0时，为晚上。

本发明首先通过定义地心固联坐标系、传感器坐标系和图像坐标系来阐述大椭圆轨道卫星大视场复杂的成像过程，然后用矩阵变换来阐述这三个坐标系的转换关系。这样就建立了像面信息和物面信息的映射关系，接着就可以求得成像面上的每个像素点代表地球上的点，进而用基本公式法求得这个地球上的点的太阳高度角，进而判断昼夜情况。本发明建立在坐标系变换和基本公式的基础上，建立了数学模型，可以定量确定大椭圆遥感卫星图像地球背景昼夜分布情况。

附图说明

图1是本发明总体流程图；

图2是各坐标系的定义示意图；

图3是传感器坐标系与传感器图像坐标系的关系示意图；

图4是卫星传感器坐标系和图像坐标系的关系示意图；

图5是大椭圆卫星春分日2：30地球背景昼夜分布的仿真图；

图6是大椭圆卫星春分日6：00地球背景昼夜分布的仿真图；

图7是大椭圆卫星夏至日2：30地球背景昼夜分布的仿真图；

图8是大椭圆卫星夏至日6：00地球背景昼夜分布的仿真图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明总体流程图如图1所示。下面结合附图和实例对本发明作进一步详细的说明。

（1）定义大椭圆轨道卫星传感器成像相关的坐标系：地心固联坐标系、传感器坐标系和图像坐标系，推导出三者的互相转换关系，并求出从地心固联坐标系到图像坐标系的转换矩阵。如图2所示是各坐标系的定义。

1)坐标系定义

①地心固联坐标系：原点O为地心，X轴指向本初子午线，Z轴指向正北，按照右手螺旋法则确定Y轴方向。设该坐标系下的点坐标为(X,Y,Z)。

②传感器坐标系：考虑太阳光照对传感器成像的影响，调整卫星相机姿态，使太阳在传感器坐标系的X_cO_cZ_c平面内，所以传感器坐标系原点O_c与卫星位置O_s重合，传感器指向Z_c的方向为卫星位置指向北极点的方向；Z_c方向与卫星与太阳连线方向叉乘得到Y_c轴方向，按照右手螺旋法则确定X_c轴。

③图像坐标系：卫星相机所成的图像以像素为单位，每个像素的坐标是该像素所在的行数和列数，需要建立以物理单位(例如km)表示的图像坐标系。存在两个图像坐标系，一个以像素所在行列数表示的坐标系，其像素坐标以(u,v)表示，一个以物理单位表示的坐标系，其像素坐标以(x,y)表示。

2)坐标系之间相互关系

①地心固联坐标系与传感器坐标系之间的关系

已知卫星位置在地心固联坐标系的坐标为O_s(x_s,y_s,z_s)，北极点在地心固联坐标系下的坐标为N(0,0,R)，其中R为地球半径。设传感器坐标系各坐标轴X_c，Y_c，Z_c的方向向量为 $\stackrel{\→}{\mathrm{\μ}}({\mathrm{\μ}}_x,{\mathrm{\μ}}_y,{\mathrm{\μ}}_z),\stackrel{\→}{\mathrm{\υ}}({\mathrm{\υ}}_x,{\mathrm{\υ}}_y,{\mathrm{\υ}}_z),\stackrel{\→}{\mathrm{\η}}({\mathrm{\η}}_x,{\mathrm{\η}}_y,{\mathrm{\η}}_z),$ 则轴的方向向量为： $\stackrel{\→}{\mathrm{\η}}=\stackrel{\→}{O_c{Z}_c}=\stackrel{\→}{O_{s}N}=-x_s\stackrel{\→}{i}-y_s\stackrel{\→}{j}+(R-z_s)\stackrel{\→}{k}$

其中：

$\left\{\begin{array}{c}{\mathrm{\η}}_x=-x_s\\ {\mathrm{\η}}_y=-y_s\\ {\mathrm{\η}}_z=R-z_s\end{array}\right.---\left(1\right)$

已知太阳位置在地心固联坐标系的坐标为S(x₀,y₀,z₀)，设卫星与太阳的连线方向为则有：

$\stackrel{\→}{O_{s}S}=(x_0-x_s)\stackrel{\→}{i}+(y_0-y_s)\stackrel{\→}{j}+(z_0-z_s)\stackrel{\→}{k}$

轴的方向向量由方向与叉乘得到：

$\stackrel{\→}{\mathrm{\υ}}=\stackrel{\→}{O_c{Y}_c}=\stackrel{\→}{O_c{Z}_c}\×\stackrel{\→}{O_{s}S}=\left|\begin{array}{ccc}\stackrel{\→}{i}& \stackrel{\→}{j}& \stackrel{\→}{k}\\ -x_s& -y_s& R-z_s\\ x_0-x_s& y_0-y_s& z_0-z_s\end{array}\right|$

则

$\left\{\begin{array}{c}{\mathrm{\υ}}_x=(z_s-z_0)y_s-(R-z_s)(y_0-y_s)\\ {\mathrm{\υ}}_y=(z_0-z_s)x_s+(R-z_s)(x_0-x_s)\\ {\mathrm{\υ}}_z=-x_s,y_0+y_s{x}_0\end{array}\right.---\left(2\right)$

由右手准则得传感器坐标系的轴的方向矢量为：

$\stackrel{\→}{\mathrm{\μ}}=\stackrel{\→}{O_c{X}_c}=\stackrel{\→}{O_c{Y}_c}\×\stackrel{\→}{O_c{Z}_c}=\left|\begin{array}{ccc}\stackrel{\→}{i}& \stackrel{\→}{j}& \stackrel{\→}{k}\\ {\mathrm{\υ}}_x& {\mathrm{\υ}}_y& {\mathrm{\υ}}_z\\ {\mathrm{\η}}_x& {\mathrm{\η}}_y& {\mathrm{\η}}_z\end{array}\right|$

则有：

$\left\{\begin{array}{c}{\mathrm{\μ}}_x=v_y{\mathrm{\η}}_z-v_z{\mathrm{\η}}_y\\ {\mathrm{\μ}}_y=v_z{\mathrm{\η}}_x-v_x{\mathrm{\η}}_z\\ {\mathrm{\μ}}_z=v_x{\mathrm{\η}}_y-v_y{\mathrm{\η}}_x\end{array}\right.---\left(3\right)$

则由正交变换可得到地心固联坐标系与卫星坐标系之间的转换关系为：

首先进行平移变换，使卫星坐标系的原点O_s(x_s,y_s,z_s)落到地心固联坐标系的原点。该地心固联坐标系到传感器坐标系的平移变换为T(-x_s,-y_s,-z_s)。

$T({-x}_s,{-y}_s,{-z}_s)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ -x_s& -y_s& {-z}_s& 1\end{array}\right]---\left(4\right)$

再进行旋转变换得地心固联坐标系到传感器坐标系的转换矩阵M₁，根据式(1)～(3)有：

$M_1=T({-x}_s,-y_s,-z_s)\left[\begin{array}{cccc}{\mathrm{\μ}}_x& {\mathrm{\μ}}_y& {\mathrm{\μ}}_z& 0\\ {\mathrm{\υ}}_x& {\mathrm{\υ}}_y& {\mathrm{\υ}}_z& 0\\ {\mathrm{\η}}_x& {\mathrm{\η}}_y& {\mathrm{\η}}_z& 0\\ 0& 0& 0& 1\end{array}\right]---\left(5\right)$

北极点在地心固联坐标系下的坐标为N(0,0,R)，传感器坐标系轴方向为卫星位置指向北极点，可求得北极点在传感器坐标系下的坐标为 $(\mathrm{0,0},\sqrt{x_s^2+y_s^2+{(z_s-R)}^2}).$

②传感器坐标系与图像坐标系之间的关系

图像坐标系有两种表现形式：图像像素坐标系和图像物理坐标系，如图3所示为图像像素坐标系与图像物理坐标系之间的关系示意图。在图像上定义直角坐标系u，v，每一象素的坐标(u,v)分别是该象素在数组中的列数和行数，(u,v)是以象素为单位的图像坐标系的坐标，由于(u,v)只表示象素位于数组中的列数和行数，并没有用物理单位来表示出该象素在图像中的位置，需要再建立以物理单位表示的图像坐标系，该坐标以图像内某一点O₁为原点，x轴和y轴分别与u轴和v轴平行，在x,y坐标系中，原点O₁定义为卫星传感器光轴与图像平面的交点，该点一般位于图像中心，但由于制作的原因，也会有所偏离，若O₁在u，v坐标系下的坐标为(u₀,v₀)，每一个象素点在x轴与y轴方向上的物理尺寸为dx、dy，则图像中任意一个象素在两个坐标系下的坐标转换关系：

$\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{\mathrm{dx}}& 0& u_0\\ 0& \frac{1}{\mathrm{dy}}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]---\left(6\right)$

卫星传感器成像几何关系可由图4所示，其中O点为卫星传感器光心，X_c轴和Y_c轴与图像的x轴和y轴平行，Z_c轴为卫星传感器的光轴，与图像平面垂直；光轴与图像平面的交点即为图像坐标系的原点，OO₁为卫星传感器焦距f。

空间任意一点P在图像上的成像位置可用针孔模型近似表示，即任意点P在图像上的投影位置P在图像上的投影位置为光心O与P点的连线OP与图像平面的交点，这种关系也称为中心投影或透视投影。由变换原理知北极点为传感器成像面的中心，即为O₁点。

由比例关系有如下关系式：

$Z_c\left(\begin{array}{c}x\\ y\\ 1\end{array}\right)=\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left(\begin{array}{c}{X}_c\\ Y_c\\ Z_c\\ 1\end{array}\right)---\left(7\right)$

由式（6）、（7），则有：

$Z_c\left(\begin{array}{c}u\\ v\\ 1\end{array}\right)=M\left(\begin{array}{c}{X}_c\\ Y_c\\ Z_c\\ 1\end{array}\right)---\left(8\right)$

其中相机参数矩阵M为：

$M=\left[\begin{array}{ccc}\frac{1}{\mathrm{dx}}& 0& u_0\\ 0& \frac{1}{\mathrm{dy}}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]=\left[\begin{array}{cccc}\frac{f}{\mathrm{dx}}& 0& u_0& 0\\ 0& \frac{f}{\mathrm{dx}}& v_0& 0\\ 0& 0& 1& 0\end{array}\right]---\left(9\right)$

3）运动轨迹在传感器成像面上的映射

地心固联坐标系中任意一点(X,Y,Z)，联立式5)、(8)，得到该点在成像面上的坐标(u,v)，满足：

$Z_c\left(\begin{array}{c}u\\ v\\ 1\end{array}\right)\mathrm{MM}\begin{array}{}\end{array}\left(\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right)---\left(10\right)$

其中M是相机参数矩阵，由相机内部参数决定，M₁是转换矩阵，与太阳和卫星位置有关，而太阳和卫星位置是实时变化的，所以M₁也是实时变化的。

4）成像面到地球上的点的映射

已知背景上的点在成像面上的像素坐标为(u,v)，根据相机的行列视场像元角dixAngleRow、dixAngleCol，可计算出该点在传感器坐标系下指向相应地球上的点向量的m_c,n_c,l_c，值，则有

$\left(\begin{array}{c}{m}_c\\ n_c\\ l_c\\ 1\end{array}\right)=M^{\′}\left(\begin{array}{c}u\\ v\\ 1\end{array}\right)$

$M^{\′}=\left[\begin{array}{ccc}\tan \left(\mathrm{dixAngleCol}\right)& 0& 0\\ 0& \tan \left(\mathrm{dixAngleRow}\right)& 0\\ 0& 0& 1\\ 0& 0& 1\end{array}\right]---\left(12\right)$

根据传感器坐标系到地心固联坐标系的转换矩阵式，可计算出该向量在地心固联坐标系下为则有

$\left(\begin{array}{c}\stackrel{\→}{l}\\ 1\end{array}\right)=M_1^{-1}\left(\begin{array}{c}\stackrel{\→}{l_c}\\ 1\end{array}\right)=M_1^{-1}{M}^{\′}\left(\begin{array}{c}u\\ v\\ 1\end{array}\right)---\left(13\right)$

计算由方向和传感器位置Os(xs,ys,zs)确定的直线与地球的交点A(X_a,Y_a,Z_a)，满足方程：

$\left\{\begin{array}{c}{X}^2+Y^2+Z^2=R^2\\ \frac{X-x_s}{m}=\frac{Y-y_s}{n}=\frac{Z-z_s}{l}\end{array}\right.---\left(14\right)$

解上述方程组时，可以令则X_a＝mt+x_s,Y_a＝nl+y_s,Z_a＝lt+z_s，带入第一个方程得一元二次方程组：(m+n+l)t²+2×(mx_s+ny_s+lz_s)×t+x² _s+y² _s+z² _s-R²＝0。当Δ＜0时，方程组无解，即点A不是地球上的点，当Δ≥0时，方程组有解，即可求得由传感器成像面上(u,v)点发出的视线与地球的交点A(X_a,Y_a,Z_a)。

（2）使用STK软件仿真大椭圆轨道卫星，导出不同时刻的太阳和卫星位置信息。

根据大椭圆轨道卫星的轨道参数和传感器参数，使用STK软件可以仿真大椭圆轨道卫星，利用STK软件的Report Graph Manager模块可以导出不同时刻大椭圆轨道卫星和太阳的位置信息，即卫星和太阳在地心固联坐标系的坐标。

（3）由以上两个步骤中的坐标系转换矩阵，太阳和卫星位置及传感器相关参数，可以建立物面信息和像面信息的映射关系。遍历成像面上的点，计算地球背景区域的像素点（u，v）映射到地球上的点（X，Y，Z）。

（4）求地球上的点（X，Y，Z）在当前时刻的太阳高度角，进而判断这个点的昼夜情况。

当已知观测点的位置和具体时间，就可以求得此时此地的太阳高度角时，再计算过程中主要涉及到下列一个基本公式和四个基本参数。其中基本公式为：

1）地理纬度观测地的纬度，已知观测点的地心固联坐标系下的坐标就可以容易求得观测地的经度和纬度。南半球纬度范围：北半球纬度范围：

2）太阳赤纬角（δ）:地球中心和太阳中心的连线与地球赤道平面的夹角称为太阳赤纬角。它在春分和秋分时刻等于零，,而在夏至和冬至时刻有极值，分别为±23°26.5′；

3）太阳时角（τ）:通过太阳的时圈与观测者子午圈之间的两面角。地球自转1周为360°，相应的时间为24小时，即每小时相应的时角为15°,每4min的时角为1°。太阳正午12时太阳时角为0°，其他时刻太阳时角的数值等于距离正午的时间乘以15，上午时角为负值,下午时角为正值；

4）太阳高度角（h）:地球表面上某点和太阳的连线与地平线之间的夹角,正午时太阳高度角最大，日出日落时太阳高度角为0。当某点太阳高度角为正时，这个点为白天，否则处于晚上。

从上述公式中可以看出，只要得到太阳赤纬（δ），地理纬度和太阳时角（τ），就可以求出太阳高度角h，进而可以判断观测地的昼夜情况。在这三个未知量中，观测点确定之后是很容易求得地理纬度的，下面主要介绍太阳赤纬（δ）和太阳时角（τ）的计算，用到的是公式法，即利用传统的公式计算。

1)太阳赤纬角计算

太阳赤纬角在周年运动中任何时刻的具体值都是严格已知的，即δ是时间函数，可以用以下表达式表述，

δ＝0.3723+23.2567sinθ+0.1149sin2θ-0.1712sin3θ-0.758cosθ+0.3656cos2θ+0.0201cos3θ

式中δ的单位是°，θ是日角，θ＝2πt/365.2422

这里的t又由两部分组成，t＝N-N_o

式中N为积日，就是日期在年内的顺序号，例如，1月1日的积日为1，平年12月31日的积日为365，闰年12月31日的积日为366。

而N_o的计算式：N_o＝79.6764+0.2422*(年份-1985)-(int)[(年份-1985)/4]

2)太阳时角计算

真正的太阳在黄道上的运动不是匀速的，而是时快时慢，因此，真太阳日的长短也就各不相同但人们的实际生活需要一种均匀不变的时间单位，这就需要寻找一个假想的太阳，它以均匀的速度在运行这个假想的太阳就称为平太阳，其周日的持续时间称平太阳日，由此而来的小时称为平太阳时。

从我们已知的北京时间转换到太阳时角需要三个步骤

①把北京时S转换为地方时S_d，即S_d＝S+(JD-120°)×4/60，式中JD是观测点的经度，北京时的标准经度是120°，乘4是将角度转换为分钟，除60是将分钟转换为小时。

②通过时差订正，把地方时S_d转换为真太阳时So，即So＝S_d+E_t/60，式中E_t是时差，单位是分钟，所以需要除60，将分钟转换成小时，而E_t可以下面公式求得：E_t＝0.0028-1.9857sinθ+9.9059sin2θ-7.0924cosθ-0.6882cos2θ

③把真太阳时So转换为太阳时角τ，即τ＝15×(S_o-12)，式中So的单位为小时，所以需要乘15，将小时转换为角度°。

最后根据计算出来的太阳高度角的正负判断此时此刻的昼夜情况，即当sinh＞0时，为白天，当sinh＜0时，为晚上。

按照以上步骤实施得仿真结果图如图所示，图像大小为8000行*12000列，灰色区域是表示白天，黑色区域表示晚上，白色圈为北极圈。其中，图5为春分日2：30时刻的仿真结果图，图6为春分日6：00时刻的仿真结果图，从图5、6中可以看出，虽然春分日北极圈内都是昼夜平分，但是成像区域不一样。同理，图7、8是夏至日2：30和6：00时刻仿真结果图（极昼）。

本发明可以通过基于坐标系转换和计算太阳高度角的方法来定量确定大椭圆遥感卫星地球背景昼夜分布情况。本发明的特点是理论建模，算法简单。

本发明不仅局限于上述具体实施方式，本领域一般技术人员根据本发明公开的内容，可以采用其它多种具体实施方式实施本发明，因此，凡是采用本发明的设计结构和思路，做一些简单的变化或更改的设计，都落入本发明保护的范围。

Claims (4)

1.一种大椭圆遥感卫星地球背景成像区域昼夜分布的确定方法，其特征在于，所述方法包括：

（1）建立地心固联坐标系、传感器坐标系和图像坐标系，获取当前时刻t的太阳位置、卫星位置以及传感器参数，确定当前时刻t下成像面到物面的转换关系；

（2）判断当前时刻t下，成像面上点（u，v）映射到物面上的点是否为地球上的点，如果是地球上的点，则根据成像面到物面的转换关系求像素点（u，v）在物面上的点（X，Y，Z），如果不是地球上的点，则处理成像面上下一个像素点；

（3）计算像素点（u，v）在物面上的点（X，Y，Z）在当前时刻t的太阳高度角，并根据该太阳高度角的正负来判断这个点的昼夜情况，即太阳高度角为正，则为白天，否则为晚上；

（4）遍历地球背景成像区域的成像面上的所有像素点，计算其在物面上的对应点，并判断物面上的每个对应点的昼夜情况，得到地球背景成像区域的昼夜分布图。

2.根据权利要求1所述的方法，其特征在于，所述步骤（1）中各坐标系定义如下：

地心固联坐标系：原点O为地心，X轴指向本初子午线，Z轴指向正北，按照右手螺旋法则确定Y轴方向，地心固联坐标系下的点坐标为(X,Y,Z)；

传感器坐标系：原点O_c与卫星位置O_s重合，传感器指向Z_c的方向为卫星位置指向北极点的方向，Z_c方向与卫星与太阳连线方向叉乘得到Y_c轴方向，按照右手螺旋法则确定X_c轴；

图像坐标系：包括图像像素坐标系和图像物理坐标系，图像像素坐标系的原点O₁定义为卫星传感器光轴与图像平面的交点，其横轴u和纵轴v分别与传感器坐标系的X_c和Y_c平行，其像素坐标以(u,v)表示，(u₀,v₀)为原点O₁在图像像素坐标系下的坐标；图像物理坐标系与图像像素坐标系平行，其像素坐标以(x,y)表示，图像物理坐标系中每一个像素点在x轴与y轴方向上的物理尺寸为dx、dy；

成像面中的点(u,v)与其在物面中的对应点(X,Y,Z)满足如下关系：

$\left\{\begin{array}{c}{X}^2+Y^2+Z^2=R^2\\ \frac{X-x_s}{m}=\frac{Y-y_s}{n}=\frac{Z-z_s}{l}\end{array},\right.$

$\stackrel{\→}{l}=(m,n,l),$

$\left(\begin{array}{c}\stackrel{\→}{l}\\ 1\end{array}\right)=M_1^{-1}\left(\begin{array}{c}\stackrel{\→}{l_c}\\ 1\end{array}\right)=M_1^{-1}{M}^{\′}\left(\begin{array}{c}u\\ v\\ 1\end{array}\right),$

其中，M₁为地心固联坐标系到传感器坐标系的转换矩阵，

$M_1=T({-x}_s,-y_s,-z_s)\left[\begin{array}{cccc}{\mathrm{\μ}}_x& {\mathrm{\μ}}_y& {\mathrm{\μ}}_z& 0\\ {\mathrm{\υ}}_x& {\mathrm{\υ}}_y& {\mathrm{\υ}}_z& 0\\ {\mathrm{\η}}_x& {\mathrm{\η}}_y& {\mathrm{\η}}_z& 0\\ 0& 0& 0& 1\end{array}\right],$

$T({-x}_s,{-y}_s,{-z}_s)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ -x_s& -y_s& {-z}_s& 1\end{array}\right],$

$\left\{\begin{array}{c}{\mathrm{\η}}_x=-x_s\\ {\mathrm{\η}}_y=-y_s\\ {\mathrm{\η}}_z=R-z_s\end{array}\right.,$

$\left\{\begin{array}{c}{\mathrm{\υ}}_x=(z_s-z_0)y_s-(R-z_s)(y_0-y_s)\\ {\mathrm{\υ}}_y=(z_0-z_s)x_s+(R-z_s)(x_0-x_s)\\ {\mathrm{\υ}}_z=-x_s,y_0+y_s{x}_0\end{array}\right.,$

$\left\{\begin{array}{c}{\mathrm{\μ}}_x=v_y{\mathrm{\η}}_z-v_z{\mathrm{\η}}_y\\ {\mathrm{\μ}}_y=v_z{\mathrm{\η}}_x-v_x{\mathrm{\η}}_z\\ {\mathrm{\μ}}_z=v_x{\mathrm{\η}}_y-v_y{\mathrm{\η}}_x\end{array}\right.,$

(x_s,y_s,z_s)为传感器坐标系原点O_s在地心固联坐标系中的坐标，(x₀,y₀,z₀)为太阳在地心固联坐标系中的坐标，R为地球半径；

其中 $M^{\′}=\left[\begin{array}{ccc}\tan \left(\mathrm{dixAngleCol}\right)& 0& 0\\ 0& \tan \left(\mathrm{dixAngleRow}\right)& 0\\ 0& 0& 1\\ 0& 0& 1\end{array}\right],$

dixAngleRow,dixAngleCol为相机的行列视场像元角。

3.如权利要求2所述的方法，其特征在于，所述步骤（2）中判断成像面上点（u，v）映射到物面上的点是否为地球上的点，以及计算像面上的像素点（u，v）在物面上地球上的点（X，Y，Z），具体为：

（2.1）根据式 $\left(\begin{array}{c}\stackrel{\→}{l}\\ 1\end{array}\right)=M_1^{-1}{M}^{\′}\left(\begin{array}{c}u\\ v\\ 1\end{array}\right)$ 求得向量 $\stackrel{\→}{l}=(m,n,l);$

（2.2）判断式 $\left\{\begin{array}{c}{X}^2+Y^2+Z^2=R^2\\ \frac{X-x_s}{m}=\frac{Y-y_s}{n}=\frac{Z-z_s}{l}\end{array}\right.$ 是否有解，令则X＝mt+x_s,Y＝nl+y_s,Z＝lt+z_s，带入方程X²+Y²+Z²＝R²得一元二次方程组：(m+n+l)t²+2×(mx_s+ny_s+lz_s)×t+x² _s+y² _s+z² _s-R²＝0；当Δ＜0时，方程组无解，即点（u，v）映射到物面上的点不是地球上的点，，当Δ≥0时，方程组有解，点（u，v）映射到物面上的点不是地球上的点；其中Δ＝(2×(mx_s+ny_s+lz_s))²-4×(m+n+l)×(x² _s+y² _s+z² _s-R²)；

（2.3）求解（2.2）中的一元二次方程组，得到成像面上的像素点（u，v）在物面上地球上的点（X，Y，Z）。

4.如权利要求1至3任一项所述的方法，其特征在于，所述步骤（3）具体包括：

（3.1）根据该物面上的点（X，Y，Z）及当前时刻t，确定地理纬度，太阳赤纬角δ，太阳时角τ；

（3.2）计算

（3.3）根据该物面上的点的sinh判断该点的昼夜情况，当sinh＞0时，为白天，当sinh＜0时，为晚上。