CN104729457A - 太阳相对近地轨道微小卫星位置的确定方法 - Google Patents

Abstract

太阳相对近地轨道微小卫星位置的确定方法，涉及太阳相对于航天器的位置的确定方法。为了解决目前还没有一种能够根据近地轨道微小卫星在J2000坐标系下的位置、速度与UTC信息解算太阳相对于微小卫星的位置的方法的问题。本发明已知近地微小卫星在J2000坐标系下的速度、位置以及UTC，通过岁差转换矩阵P和地球自转转换矩阵R求得坐标转换矩阵W，求得微小卫星所在位置的地心经度λ、地心纬度与地心距R_d；计算出太阳高度角h和太阳方位角A，将太阳相对于微小卫星的位置在北东地坐标系下表示出来；解算出北东地坐标系与轨道坐标系间的坐标转换矩阵Q，将太阳相对于微小卫星的位置在轨道坐标系下表示出来。本发明适用于航天器姿态控制技术领域。

Description

太阳相对近地轨道微小卫星位置的确定方法

技术领域

本方法涉及太阳相对于航天器位置的确定方法，属于航天器姿态控制技术领域。

背景技术

近地轨道微小卫星姿态控制系统常利用星载GPS输出的J2000坐标系下卫星的速度、位置以及时间信息进行姿态确定，利用GPS的输出信息可以解算出太阳相对于微小卫星的位置，此相对位置可以作为矢量之一，配合另一矢量(如地磁场方向矢量)进行双矢量定姿，或配合太阳敏感器共同为微小卫星的姿态确定与姿态机动提供参考。

发明内容

本发明为了解决目前还没有一种能够根据近地轨道微小卫星在J2000坐标系下的位置、速度以及协调世界时UTC(Coordinate Universal Time)信息解算太阳相对于近地轨道微小卫星的位置的方法的问题。

太阳相对于近地轨道微小卫星的位置的确定方法，包括以下步骤：

步骤一、在已知由星载GPS输出的微小卫星在J2000坐标系下的速度、位置，以及格林尼治年、月、日信息与协调世界时UTC(Coordinate Universal Time)的前提下，通过岁差转换矩阵P和地球自转转换矩阵R求得由J2000坐标系到地心地固坐标系(Earth-Centered,Earth-Fixed Coordinate System，简称ECEF)的坐标转换矩阵W，进而将微小卫星的速度、位置矢量在ECEF下表示出来；

步骤二、根据微小卫星在ECEF下的位置，利用反三角函数以及经纬度象限条件求得微小卫星所在位置的地心经度λ、地心纬度

步骤三、利用步骤二的结果，计算出积日N，日角θ、太阳赤纬角δ、太阳时角τ、地方时S_d，进而计算出太阳高度角h、太阳方位角A；

步骤四、利用步骤三求得的太阳高度角h、太阳方位角A，在已知太阳相对于微小卫星的距离的前提下，将太阳相对于微小卫星的位置在北东地坐标系下表示出来；

步骤五、利用步骤一中J2000坐标系下微小卫星的位置和速度信息，解算出北东地坐标系与轨道坐标系之间的夹角，利用夹角得到北东地坐标系到轨道坐标系之间的坐标转换矩阵，将在北东地坐标系下表示的太阳相对于微小卫星的位置在轨道坐标系中表示出来。

本发明具有以下有益效果：本发明提出了一种根据近地轨道微小卫星在J2000坐标系下的位置、速度以及UTC信息，解算太阳相对于微小卫星的位置，并将此相对位置在轨道坐标系中表示出来的方法。用本发明提供的方法解算出的轨道坐标系下太阳相对于微小卫星的位置，与利用STK生成的相同轨道参数下的太阳相对于微小卫星的位置，有较高的吻合度，太阳相对于微小卫星的位置矢量指向精度评价误差不超过0.2°，最高不超过0.4°。

附图说明

图1太阳相对于近地轨道微小卫星的位置矢量解算过程；

图2太阳相对于近地轨道微小卫星位置r_{sat_sun_uni_x}与r_{sat_sun_STK_x}对比效果图；

图3太阳相对于近地轨道微小卫星位置r_{sat_sun_uni_y}与r_{sat_sun_STK_y}对比效果图；

图4太阳相对于近地轨道微小卫星位置r_{sat_sun_uni_z}与r_{sat_sun_STK_z}对标效果图；

图5太阳相对于近地轨道微小卫星位置矢量指向误差效果图。

具体实施方式

具体实施方式一：本实施方式的本方法提出的太阳相对于微小卫星的位置矢量计算过程如图1所示。

步骤一、在已知由星载GPS输出的微小卫星在J2000坐标系下的速度、位置，以及格林尼治年、月、日信息与协调世界时UTC(Coordinate Universal Time)的前提下，通过岁差转换矩阵P和地球自转转换矩阵R求得由J2000坐标系到地心地固坐标系(Earth-Centered,Earth-Fixed Coordinate System，简称ECEF)的坐标转换矩阵W，进而将微小卫星的速度、位置矢量在ECEF下表示出来；

步骤二、根据微小卫星在ECEF下的位置，利用反三角函数以及经纬度象限条件求得微小卫星所在位置的地心经度λ、地心纬度

步骤三、利用步骤二的结果，计算出积日N，日角θ、太阳赤纬角δ、太阳时角τ、地方时S_d，进而计算出太阳高度角h、太阳方位角A；

步骤四、利用步骤三求得的太阳高度角h、太阳方位角A，在已知太阳相对于微小卫星的距离的前提下，将太阳相对于微小卫星的位置在北东地坐标系下表示出来；

步骤五、利用步骤一中J2000坐标系下微小卫星的位置和速度信息，解算出北东地坐标系与轨道坐标系之间的夹角，利用夹角得到北东地坐标系到轨道坐标系之间的坐标转换矩阵，将在北东地坐标系下表示的太阳相对于微小卫星的位置在轨道坐标系中表示出来。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中：由J2000坐标系到ECEF坐标系需进行四次坐标转换，依次为岁差转换，章动转换，地球自转转换，地球极移转换。本方法只考虑在由J2000坐标系到ECEF坐标系的坐标转换中起决定性作用的岁差转换矩阵(记为P)，与地球自转转换矩阵(记为R)；

(一)求得儒略世纪数T：

首先利用格林尼治年、月、日与UTC求取儒略日JD：

$\begin{array}{c}\mathrm{JD}=367*\mathrm{year}-\mathrm{int}\left(\frac{7*(\mathrm{year}+\mathrm{int}\left(\frac{\mathrm{month}+9}{12}\right))}{4}\right)+\mathrm{int}\left(\frac{275*\mathrm{month}}{9}\right)\\ +\mathrm{day}+1721013.5+\frac{\frac{(\frac{\sec \mathrm{ond}}{60}+\min \mathrm{ute})}{60}+\mathrm{hout}}{24}\end{array}---\left(1\right)$

其中year，month，hour，minute，second分别为格林尼治年，月，日，小时，分钟，秒，实际应用时由GPS实时提供给微小卫星姿态控制系统；int()表示对括号内的数值取整；利用儒略日求得儒略世纪数：

$T=\frac{\mathrm{JD}-2451545}{36525}---\left(1\right)$

(二)求取岁差转换矩阵P与地球自转转换矩阵R：

岁差转换矩阵的求取：

三个岁差角z_A，θ_A，ζ_A的计算公式为(单位为弧度)：

z_A＝0.011180860T+5.308*10^-6T²+8.9*10^-8T³

θ_A＝0.009717173T-2.068*10^-6T²-2.02*10^-7T³   (2)

ζ_A＝0.011180860T+1.464*10^-6T²+8.7*10^-8T³

根据坐标转换矩阵的原理，岁差矩阵可以表示为：

$P=\left[\begin{array}{ccc}\cos (-z_A)& \sin (-z_A)& 0\\ -\sin (-z_A)& \cos (-z_A)& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}\cos \left({\mathrm{\θ}}_A\right)& 0& -\sin \left({\mathrm{\θ}}_A\right)\\ 0& 1& 0\\ \sin \left({\mathrm{\θ}}_A\right)& 0& \cos \left({\mathrm{\θ}}_A\right)\end{array}\right]\left[\begin{array}{ccc}\cos (-{\mathrm{\ζ}}_A)& \sin (-{\mathrm{\ζ}}_A)& 0\\ -\sin (-{\mathrm{\ζ}}_A)& \cos (-{\mathrm{\ζ}}_A)& 0\\ 0& 0& 1\end{array}\right]---\left(3\right)$

地球自转旋转矩阵的求取：

将J2000坐标系转换到ECEF，关键在于求取地球自转旋转矩阵，求解地球自转旋转矩阵关键在于求取格林尼治恒星时(Greenwich Sidereal Time)θ_GST，θ_GST可由公式(2)求得的儒略世纪数T表达出来：

θ_GST＝67310.54841+3164400184.812866T+0.0093104T²-6.2*10^-6T³   (4)

此公式求得的θ_GST是以秒(second)为单位，根据时间与角度、角度与弧度之间的换算关系：

将θ_GST转换成以弧度为单位：

${\mathrm{\θ}}_{\mathrm{GST}}^{\′}={\mathrm{\θ}}_{\mathrm{GST}}*\frac{15}{3600}*\frac{\mathrm{\π}}{180}---\left(6\right)$

在忽略章动和地球极移影响的情况下，J2000坐标系经过岁差转换之后，再绕其Z'轴(J2000坐标系经过岁差转换后的坐标系三轴表示分别为X',Y',Z')逆时针转动θ'_GST，即与ECEF重合。地球自转转换矩阵R可表示为：

$R=\left[\begin{array}{ccc}\cos \left({\mathrm{\θ}}_{\mathrm{GST}}^{\′}\right)& \sin \left({\mathrm{\θ}}_{\mathrm{GST}}^{\′}\right)& 0\\ -\sin \left({\mathrm{\θ}}_{\mathrm{GST}}^{\′}\right)& \cos \left({\mathrm{\θ}}_{\mathrm{GST}}^{\′}\right)& 0\\ 0& 0& 1\end{array}\right]---\left(7\right)$

由J2000坐标系到ECEF的坐标转换矩阵W由地球自转转换矩阵R点乘岁差转换矩阵P得到：

在J2000坐标系下，微小卫星相对于地心的位置表示为r_J2000，微小卫星相对于地心的速度表示为v_J2000：

r_J2000＝[r_{J2000_x} r_{J2000_y} r_{J2000_z}]^T

                                (9)

v_J2000＝[v_{J2000_x} v_{J2000_y} v_{J2000_z}]^T

其中r_{J2000_x}，r_{J2000_y}，r_{J2000_z}分别为r_J2000在J2000坐标系下的三个分量，v_{J2000_x}，v_{J2000_y}，v_{J2000_z}分别为v_J2000在J2000坐标系下的三个分量。

在ECEF下，微小卫星相对于地心的位置在ECEF中表示为r_ECEF，微小卫星相对于地心的速度表示为v_ECEF：

r_ECEF＝[r_{ECEF_x} r_{ECEF_y} r_{ECEF_z}]^T

                                (10)

v_ECEF＝[v_{ECEF_x} v_{ECEF_y} v_{ECEF_z}]^T

其中r_{ECEF_x}，r_{ECEF_y}，r_{ECEF_z}分别为r_ECEF在ECEF下的分量；v_{ECEF_x}，v_{ECEF_y}，v_{ECEF_z}分别为v_ECEF在ECEF下的分量。r_ECEF与v_ECEF可分别由J2000坐标系到ECEF的坐标转换矩阵W点乘r_J2000、v_J2000得到：

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：根据微小卫星相对地球的位置r_ECEF在ECEF下分量r_{ECEF_x}，r_{ECEF_y}，r_{ECEF_z}，利用反三角函数以及经纬象限条件求解微小卫星所在位置的地心经度λ、地心纬度

在本方法中，设定地心经度λ的取值范围为[-π,π]，东经为正，西经为负；地心纬度的取值范围为北纬为正，南纬为负。

对于地心经度，由于反三角函数atan求解后的结果与期望得到的地心经度不是一一对应的关系，需要根据象限条件判断后进行转换，判断和转换过程如下：

$\left\{\begin{array}{c}{r}_{\mathrm{ECEF}\_x}>0\left\{\begin{array}{c}{r}_{\mathrm{ECEF}\_y}\&GreaterEqual;0\&RightArrow;\mathrm{\&lambda;}=a\tan \left(\frac{r_{\mathrm{ECEF}\_y}}{r_{\mathrm{ECEF}\_x}}\right)\\ r_{\mathrm{ECEF}\_y}<0\&RightArrow;\mathrm{\&lambda;}=a\tan \left(\frac{r_{\mathrm{ECEF}\_y}}{r_{\mathrm{ECEF}\_x}}\right)+2\mathrm{\&pi;}\end{array}\right.\\ r_{\mathrm{ECEF}\_x}<0\&RightArrow;\mathrm{\&lambda;}=a\tan \left(\frac{r_{\mathrm{ECEF}\_y}}{r_{\mathrm{ECEF}\_x}}\right)+\mathrm{\&pi;}\\ r_{\mathrm{ECEF}\_x}=0\left\{\begin{array}{c}{r}_{\mathrm{ECEF}\_y}>0\&RightArrow;\mathrm{\&lambda;}=\frac{\mathrm{\&pi;}}{2}\\ r_{\mathrm{ECEF}\_y}\&le;0\&RightArrow;\mathrm{\&lambda;}=\frac{3\mathrm{\&pi;}}{2}\end{array}\right.\end{array}\right.---\left(13\right)$

由于规定了东经为正，西经为负，经过上式判断后，若λ＞π，则减去2π，得到的结果即为实际的地心经度(此时地心经度为负值，位于西经某处)。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：(一)日角θ：

$\mathrm{\&theta;}=\frac{2\mathrm{\&pi;}(N-N_0)}{365.2422}---\left(14\right)$

其中N为积日，即为从当年1月1日起开始计算的天数(例如：每年的1月1日积日为第1日N＝1，2月1日为第32日，N＝32，以此类推)。N₀的计算公式为：

$N_0=79.6764+0.2422*(\mathrm{year}-1985)-\mathrm{int}\left(\frac{\mathrm{year}-1985}{4}\right)---\left(15\right)$

(二)太阳赤纬角δ利用日角θ计算得到(单位为弧度)：

$\mathrm{\&delta;}=\left(\begin{array}{c}0.3723+23.2567\sin \mathrm{\&theta;}+0.1149\sin 2\mathrm{\&theta;}-0.1712\sin 3\mathrm{\&theta;}\\ -0.758\cos \mathrm{\&theta;}+0.3656\cos 2\mathrm{\&theta;}+0.0201\cos 3\mathrm{\&theta;}\end{array}\right)*\left(\frac{\mathrm{\&pi;}}{180}\right)---\left(16\right)$

(三)太阳时角τ：

真太阳时S由地方时S_d与时差E_t计算得到：

$S=S_d+\frac{E_t}{60}---\left(17\right)$

地方时S_d可由GPS输出的UTC时间：小时(hour)，分钟(minute)秒(second)与步骤二中解算出的地心经度λ计算得到：

$S_d=\mathrm{hour}+\frac{\min \mathrm{ute}}{60}+\frac{\sec \mathrm{ond}}{3600}+\mathrm{\&lambda;}\frac{180}{\mathrm{\&pi;}*15}---\left(18\right)$

由于本模型中微小卫星位于西经位置时λ为负，根据每天24小时的实际情况，对S_d做如下处理：

$\left\{\begin{array}{c}{S}_d<0\&DoubleRightArrow;S_d=S_d+24\\ \mathrm{else}\&DoubleRightArrow;S_d=S_d\end{array}\right.---\left(19\right)$

本模型将UTC近似作为平太阳时，时差E_t为平太阳时与真太阳时之间的差值，通过日角θ解算出来：

E_t＝0.0028-1.9857sinθ+9.9059sin2θ-7.0924cosθ-0.6882cos2θ   (20)

太阳时角τ的计算公式为(计算结果单位为弧度)：

$\mathrm{\&tau;}=15*(S-12)*\left(\frac{\mathrm{\&pi;}}{180}\right)---\left(21\right)$

(四)太阳高度角h与太阳方位角A：

太阳高度角h是指太阳直射光线与地平面间的夹角，太阳方位角A是指太阳直射光线在地平面上的投影线与地平面正南向所夹的角，规定以正南向为0°，向西为正值，向东为负值，太阳高度角h与太阳方位角A同太阳赤纬角δ、太阳时角τ，以及卫星的地心纬度存在如下关系：

由于反三角函数的值域与期望角度值不是一一对应的关系，因此在求解方位角A时，在程序中需进行象限判断：

$\left\{\begin{array}{c}\\ \sin A>0\&DoubleRightArrow;A=a\cos \left(\cos A\right)\\ \sin A<0\&DoubleRightArrow;A=2\mathrm{\&pi;}-a\cos \left(\cos A\right)\\ \sin A=0\left\{\begin{array}{c}\cos A>0\&DoubleRightArrow;A=0\\ \cos A\&le;0\&DoubleRightArrow;A=\mathrm{\&pi;}\end{array}\right.\end{array}\right.---\left(23\right)$

太阳高度角h求解公式为：

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：在北东地坐标系下，太阳相对于微小卫星的位置矢量r_{sat_sun}为：

r_{sat_sun}＝[r_{sat_sun_x} r_{sat_sun_y} r_{sat_sun_z}]^T   (25)

其中r_{sat_sun_x}，r_{sat_sun_y}，r_{sat_sun_z}分别为太阳相对于微小卫星位置矢量的北向分量、东向分量与向地心方向的分量。北向、东向与向地心方向的分量与太阳与微小卫星的距离R_{sun_sat}，太阳方位角A与太阳高度角h之间的关系为：

r_{sat_sun_x}＝-R_{sun_sat}*cosA*cosh

r_{sat_sun_y}＝-R_{sun_sat}*sinA*cosh   (26)

r_{sat_sun_z}＝-R_{sun_sat}*sinh

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：

本模型中定义轨道坐标系为z_orbit轴指向地心，x_orbit轴位于轨道平面内，与z_orbit轴垂直且指向速度方向，y_orbit轴与其他两轴形成右手系，由此可知北东地与轨道坐标系z_orbit轴重合，因而可通过一次绕z_orbit轴的旋转得到两者的坐标转换矩阵，利用卫星所在的子午面与轨道面之间的几何关系，可求得北东地坐标系与轨道坐标系之间的夹角α：

其中χ为求解过程中的中间变量，norm()表示对括号内的矢量求模长；为单位化后的v_J2000，即z_J2000为J2000坐标系下的z_J2000轴单位向量z_J2000＝[0 0 1]^T。

北东地坐标系向轨道坐标系转换的过程中，北东地坐标系绕自身z_NED轴旋转的方向(正向/负向)需进一步判定，将判定量记为β：

式中χ由公式(28)求得，若β＞0，则α即为公式(29)所得结果，若β≤0，则α为公式(29)所得结果的相反数；

由北东地坐标系到轨道坐标系的坐标转换矩阵Q为：

$Q=\left[\begin{array}{ccc}\cos \mathrm{\&alpha;}& \sin \mathrm{\&alpha;}& 0\\ -\sin \mathrm{\&alpha;}& \cos \mathrm{\&alpha;}& 0\\ 0& 0& 1\end{array}\right]---\left(29\right)$

坐标转换矩阵Q点乘北东地坐标系下的太阳相对于微小卫星的位置矢量r_{sat_sun}，即得r_{sat_sun}在轨道坐标系中的表达r_{sat_sun_orbit}：

在实际应用中，可根据需要对r_{sat_sun_orbit}进行单位化处理。

实施例

设某近地轨道微小卫星于2015年5月1日00时00分01秒入轨，入轨参数为：

$\left[\begin{array}{c}{r}_{J2000-x}\\ r_{J2000\_y}\\ r_{J2000-z}\end{array}\right]=\left[\begin{array}{c}180877.444\\ 6875758.257\\ 571.469\end{array}\right]\left(m\right),\left[\begin{array}{c}{v}_{J2000\_x}\\ v_{J2000\_y}\\ v_{J2000\_z}\end{array}\right]=\left[\begin{array}{c}992.206\\ 25.474\\ 7547.628\end{array}\right](m/s)$

在此轨道初始条件下，利用STK(Satellite Tool Kit)生成每秒更新的J2000坐标系下微小卫星的位置(生成路径：STK/Cartesian Position/J2000/x y z)和速度(生成路径：STK/Cartesian Velocity/J2000/x y z)，以此作为输入，利用专利中的方法，解算出对应时刻太阳相对于微小卫星的位置在轨道坐标系下的表达式，此表达式经单位化后记为。r_{sat_sun_orbit_unb}＝[r_{sat_sun_orbit_uni_x} r_{sat_sun_orbit_uni_y} r_{sat_sun_orbit_uni_z}]^T，其中r_{sat_sun_orbit_uni_x}，r_{sat_sun_orbit_uni_y}，r_{sat_sun_orbit_uni_z}分别为r_{sat_sun_orbit_uni}在轨道坐标系下的三个分量。验证时间为2015年5月1日00：00：01至2015年5月1日13：00：01。将r_{sat_sun_orbit_uni}与利用STK生成的轨道坐标系下的太阳相对于微小卫星(单位化后的)位置矢量(生成路径：STK/Vectors(VVLH)/Sun/x/Magnitude y/Magnitude z/Magnitude)，记为r_{sat_sun_STK}＝[r_{sat_sun_STK_x} r_{sat_sun_STK_y} r_{sat_sun_STK_z}]^T，其中r_{sat_sun_STK_x}，r_{sat_sun_STK_y}，r_{sat_sun_STK_z}分别为r_{sat_sun_STK}在轨道坐标系下的三个分量，进行对比，三轴的对比情况分别如图2-图4所示；

以上图2-4中，黑色线条分别表示r_{sat_sun_orbit_uni_x}，r_{sat_sun_orbit_uni_y}，r_{sat_sun_orbit_uni_z}。灰色线条分别表示r_{sat_sun_STK_x}，r_{sat_sun_STK_y}，r_{sat_sun_STK_z}。为了辨别，MATLAB仿真作图时黑色线条比灰色线条宽。由图可知，由专利中算法解算得到的太阳相对于微小卫星的位置，与STK在相同轨道条件下生成的太阳相对于微小卫星的位置具有很高的吻合度。

为了进一步验证该专利提供的算法解算出的太阳相对于微小卫星位置的准确性，求取矢量r_{sat_sun_orbit_uni}与矢量r_{sat_sun_STK}之间的夹角，此夹角体现由本算法解算出的太阳相对于微小卫星位置矢量方向的指向误差，指向误差效果图如下图5所示；

由图5可知，由本专利中提出的算法得到的太阳相对于微小卫星的位置矢量的指向误差较小，平均指向误差不超过0.2°，最高误差不超过0.4°。

Claims (6)

1.太阳相对近地轨道微小卫星位置的确定方法，其特征在于它按以下步骤实现：

步骤一、已知由星载GPS输出的微小卫星在J2000坐标系下的速度、位置，以及格林尼治年、月、日信息与协调世界时UTC，通过岁差转换矩阵P和地球自转转换矩阵R求得由J2000坐标系到ECEF的坐标转换矩阵W，将微小卫星的速度、位置矢量在ECEF下表示；

其中，所述J2000坐标系：坐标原点O_J2000位于地心，X_J2000轴指向J2000.0平春分点，Z_J2000轴向北，指向J2000.0平赤道的极点，Y_J2000轴与X_J2000、Y_J2000轴构成右手系；J2000.0指的是2000年1月1日12：00：00TDB；

ECEF坐标系：坐标原点O_ECEF位于地心，x_ECEF轴在赤道面上并指向格林尼治子午线方向，z_ECEF轴垂直于赤道面指向北极，y_ECEF轴与x_ECEF轴、z_ECEF轴形成右手系；

步骤二、根据微小卫星在ECEF下的位置，利用反三角函数以及经纬度象限条件求得微小卫星所在位置的地心经度λ、地心纬度

步骤三、利用步骤二的结果，计算出积日N，日角θ、太阳赤纬角δ、太阳时角τ、地方时S_d，进而计算出太阳高度角h、太阳方位角A；

步骤四、利用步骤三求得的太阳高度角h、太阳方位角A，在已知太阳相对于微小卫星的距离的前提下，将太阳相对于微小卫星的位置在北东地坐标系下表示；

步骤五、利用步骤一中J2000坐标系下微小卫星的位置和速度信息，解算出北东地坐标系与轨道坐标系之间的夹角，利用夹角得到北东地坐标系到轨道坐标系之间的坐标转换矩阵，将在北东地坐标系下表示的太阳相对于微小卫星的位置在轨道坐标系中表示。

2.根据权利要求1所述的太阳相对近地轨道微小卫星位置的确定方法，其特征在于步骤一具体为：

(一)求得儒略世纪数T：

首先利用格林尼治年、月、日与UTC求取儒略日JD：

$\begin{array}{c}\mathrm{JD}=360*\mathrm{year}-\mathrm{int}\left(\frac{7*(\mathrm{year}+\mathrm{int}\left(\frac{\mathrm{nonth}+9}{12}\right))}{4}\right)+\mathrm{int}\left(\frac{275*\mathrm{month}}{9}\right)\\ +\mathrm{day}+1721013.5+\frac{\frac{(\frac{\sec \mathrm{ond}}{60}+\min \mathrm{ute})}{60}+\mathrm{hour}}{24}\end{array}---\left(1\right)$

其中，year，month，hour，minute，second分别为格林尼治年，月，日，小时，分钟，秒；利用儒略日求得儒略世纪数：

$T=\frac{\mathrm{JD}-2451545}{36525}---\left(2\right)$

(二)求取岁差转换矩阵P与地球自转转换矩阵R：

岁差转换矩阵的求取：

三个岁差角z_A，θ_A，ζ_A的计算公式为：

z_A＝0.011180860T+5.308*10^-6T²+8.9*10^-8T³

θ_A＝0.009717173T-2.068*10^-6T²-2.02*10^-7T³         (3)

ζ_A＝0.011180860T+1.464*10^-6T²+8.7*10^-8T³

根据坐标转换矩阵P的原理，岁差转换矩阵可以表示为：

$P=\left[\begin{array}{ccc}\cos (-z_A)& \sin (-z_A)& 0\\ -\sin (-z_A)& \cos (-z_A)& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}\cos \left({\mathrm{\&theta;}}_A\right)& 0& -\sin \left({\mathrm{\&theta;}}_A\right)\\ 0& 1& 0\\ \sin \left({\mathrm{\&theta;}}_A\right)& 0& \cos \left({\mathrm{\&theta;}}_A\right)\end{array}\right]\left[\begin{array}{ccc}\cos (-{\mathrm{\&zeta;}}_A)& \sin (-{\mathrm{\&zeta;}}_A)& 0\\ -\sin \left({\mathrm{\&zeta;}}_A\right)& \cos (-{\mathrm{\&zeta;}}_A)& 0\\ 0& 0& 1\end{array}\right]---\left(4\right)$

地球自转转换矩阵R的求取：

格林尼治恒星θ_GST由公式(2)求得的儒略世纪数T表达出来：

θ_GST＝67310.54841+3164400184.812866T+0.0093104T²-6.2*10^-6T³    (5)

求得的θ_GST是以秒为单位，根据时间与角度、角度与弧度之间的换算关系：

将θ_GST转换成以弧度为单位：

${\mathrm{\&theta;}}_{\mathrm{GST}}^{\&prime;}={\mathrm{\&theta;}}_{\mathrm{GST}}*\frac{15}{3600}*\frac{\mathrm{\&pi;}}{180}---\left(7\right)$

在忽略章动和地球极移影响的情况下，J2000坐标系经过岁差转换之后，再绕其Z'轴逆时针转动θ'_GST，即与ECEF重合；地球自转转换矩阵R表示为：

$R=\left[\begin{array}{ccc}\cos \left({\mathrm{\&theta;}}_{\mathrm{GST}}^{\&prime;}\right)& \sin \left({\mathrm{\&theta;}}_{\mathrm{GST}}^{\&prime;}\right)& 0\\ -\sin \left({\mathrm{\&theta;}}_{\mathrm{GST}}^{\&prime;}\right)& \cos \left({\mathrm{\&theta;}}_{\mathrm{GST}}^{\&prime;}\right)& 0\\ 0& 0& 1\end{array}\right]---\left(8\right)$

由J2000坐标系到ECEF的坐标转换矩阵W由地球自转转换矩阵R点乘岁差转换矩阵P得到：

在J2000坐标系下，微小卫星相对于地心的位置表示为r_J2000，微小卫星相对于地心的速度表示为v_J2000：

r_J2000＝[r_{J2000_x} r_{J2000_y} r_{J2000_z}]^T

                                     (10)

v_J2000＝[v_{J2000_x} v_{J2000_y} v_{J2000_z}]^T

其中r_{J2000_x}，r_{J2000_y}，r_{J2000_z}分别为r_J2000在J2000坐标系下的三个分量，v_{J2000_x}，v_{J2000_y}，v_{J2000_z}分别为v_J2000在J2000坐标系下的三个分量；

在ECEF下，微小卫星相对于地心的位置在ECEF中表示为r_ECEF，微小卫星相对于地心的速度表示为v_ECEF：

r_ECEF＝[r_{ECEF_x} r_{ECEF_y} r_{ECEF_z}]^T

A                         (11)

v_ECEF＝[v_{ECEF_x} v_{ECEF_y} v_{ECEF_z}]^T

其中r_{ECEF_x}，r_{ECEF_y}，r_{ECEF_z}分别为r_ECEF在ECEF下的分量；v_{ECEF_x}，v_{ECEF_y}，v_{ECEF_z}分别为v_ECEF在ECEF下的分量；r_ECEF与v_ECEF可分别由J2000坐标系到ECEF的坐标转换矩阵W点乘r_J2000、v_J2000得到：

3.根据权利要求2所述的太阳相对近地轨道微小卫星位置的确定方法，其特征在于步骤二具体为：

根据微小卫星相对地球的位置r_ECEF在ECEF下分量r_{ECEF_x}，r_{ECEF_y}，r_{ECEF_z}，利用反三角函数以及经纬象限条件求解微小卫星所在位置的地心经度λ、地心纬度

$\mathrm{\&lambda;}=a\tan \left(\frac{r_{\mathrm{ECEF}\_y}}{r_{\mathrm{ECEF}\_x}}\right)$

                       (13)

设定地心经度λ的取值范围为[-π,π]，东经为正，西经为负；地心纬度的取值范围为北纬为正，南纬为负；

根据象限条件判断后进行转换，判断和转换过程如下：

$\left\{\begin{array}{c}{r}_{\mathrm{ECEF}\_x}>0\left\{\begin{array}{c}{r}_{\mathrm{ECEF}\_y}\&GreaterEqual;0\&RightArrow;\mathrm{\&lambda;}=a\tan \left(\frac{r_{\mathrm{ECEF}\_y}}{r_{\mathrm{ECEF}\_x}}\right)\\ r_{\mathrm{ECEF}\_y}<0\&RightArrow;\mathrm{\&lambda;}=a\tan \left(\frac{r_{\mathrm{ECEF}\_y}}{r_{\mathrm{ECEF}\_x}}\right)+2\mathrm{\&pi;}\end{array}\right.\\ r_{\mathrm{ECEF}\_x}<0\&RightArrow;\mathrm{\&lambda;}=a\tan \left(\frac{r_{\mathrm{ECEF}\_y}}{r_{\mathrm{ECEF}\_x}}\right)+\mathrm{\&pi;}\\ r_{\mathrm{ECEF}\_x}=0\left\{\begin{array}{c}{r}_{\mathrm{ECEF}\_y}>0\&RightArrow;\mathrm{\&lambda;}=\frac{\mathrm{\&pi;}}{2}\\ r_{\mathrm{ECEF}\_y}\&le;0\&RightArrow;\mathrm{\&lambda;}=\frac{3\mathrm{\&pi;}}{2}\end{array}\right.\end{array}\right.---\left(14\right)$

由于规定了东经为正，西经为负，经过上式判断后，若λ>π，则减去2π，得到的结果即为实际的地心经度。

4.根据权利要求3所述的太阳相对近地轨道微小卫星位置的确定方法，其特征在于步骤三具体为：

(一)日角θ：

$\mathrm{\&theta;}=\frac{2\mathrm{\&pi;}(N-N_0)}{365.2422}---\left(15\right)$

其中N为积日，即为从当年1月1日起开始计算的天数，N₀为计算日角θ的中间变量，其计算公式为：

$N_0=79.6764+0.2422*(\mathrm{year}-1985)-\mathrm{int}\left(\frac{\mathrm{year}-1985}{4}\right)---\left(16\right)$

(二)太阳赤纬角δ利用日角θ计算得到：

$\mathrm{\&delta;}=\left(\begin{array}{c}0.3723+23.2567\sin \mathrm{\&theta;}+0.1149\sin 2\mathrm{\&theta;}-0.1712\sin 3\mathrm{\&theta;}\\ -0.758\cos \mathrm{\&theta;}+0.3656\cos 2\mathrm{\&theta;}+0.0201\cos 3\mathrm{\&theta;}\end{array}\right)*\left(\frac{\mathrm{\&pi;}}{180}\right)---\left(17\right)$

(三)太阳时角τ：

真太阳时S由地方时S_d与时差E_t计算得到：

$S=S_d+\frac{E_t}{60}---\left(18\right)$

地方时S_d可由GPS输出的UTC时间：小时、分钟秒与步骤二中解算出的地心经度λ计算得到：

$S_d=\mathrm{hour}+\frac{\min \mathrm{ute}}{60}+\frac{\sec \mathrm{ond}}{3600}+\mathrm{\&lambda;}\frac{180}{\mathrm{\&pi;}*15}---\left(19\right)$

微小卫星位于西经位置时λ为负，根据每天24小时的实际情况，对S_d做如下处理：

$\left\{\begin{array}{c}{S}_d<0\&DoubleRightArrow;S_d=S_d+24\\ \mathrm{else}\&DoubleRightArrow;S_d=S_d\end{array}\right.---\left(20\right)$

将UTC近似作为平太阳时，时差E_t为平太阳时与真太阳时之间的差值，通过日角θ解算出来：

E_t＝0.0028-1.9857sinθ+9.9059sin2θ-7.0924cosθ-0.6882cos2θ     (21)

太阳时角τ的计算公式为：

$\mathrm{\&tau;}=15*(S-12)*\left(\frac{\mathrm{\&pi;}}{180}\right)---\left(22\right)$

(四)太阳高度角h与太阳方位角A的计算：

太阳高度角h是指太阳直射光线与地平面间的夹角，太阳方位角A是指太阳直射光线在地平面上的投影线与地平面正南向所夹的角，规定以正南向为0°，向西为正值，向东为负值，太阳高度角h与太阳方位角A同太阳赤纬角δ、太阳时角τ，以及卫星的地心纬度存在如下关系：

$\sin A=\frac{\cos \mathrm{\&delta;}*\sin \mathrm{\&tau;}}{\cosh }---\left(23\right)$

由于反三角函数的值域与期望角度值不是一一对应的关系，因此在求解太阳方位角A时，在程序中需进行象限判断：

$\left\{\begin{array}{c}\sin A>0\&DoubleRightArrow;A=a\cos \left(\cos A\right)\\ \sin A<0\&DoubleRightArrow;A=2\mathrm{\&pi;}-a\cos \left(\cos A\right)\\ \sin A=0\left\{\begin{array}{c}\cos A>0\&DoubleRightArrow;A=0\\ \cos A\&le;0\&DoubleRightArrow;A=\mathrm{\&pi;}\end{array}\right.\end{array}\right.---\left(24\right)$

太阳高度角h求解公式为：

5.根据权利要求4所述的太阳相对近地轨道微小卫星位置的确定方法，其特征在于步骤四具体为：

在北东地坐标系下，太阳相对于微小卫星的位置矢量r_{sat_sun}为：

r_{sat_sun}＝[r_{sat_sun_x} r_{sat_sun_y} r_{sat_sun_z}]^T           (26)

其中r_{sat_sun_x}，r_{sat_sun_y}，r_{sat_sun_z}分别为太阳相对于微小卫星位置矢量的北向分量、东向分量与向地心方向的分量，北向、东向向地心方向的分量与太阳与微小卫星的距离R_{sun_sat}，太阳方位角A与太阳高度角h之间的关系为：

r_{sat_sun_x}＝-R_{sun_sat}*cos A*cos h

r_{sat_sun_y}＝-R_{sun_sat}*sin A*cos h            (27)。

r_{sat_sun_z}＝-R_{sun_sat}*sin h

6.根据权利要求7所述的太阳相对近地轨道微小卫星位置的确定方法，其特征在于所述步骤五中具体为：

定义轨道坐标系为z_orbit轴指向地心，x_orbit轴位于轨道平面内，与z_orbit轴垂直且指向速度方向，y_orbit轴与其他两轴形成右手系，北东地坐标系与轨道坐标系z_orbit轴重合，通过一次绕z_orbit轴的旋转得到两者的坐标转换矩阵，利用卫星所在的子午面与轨道面之间的几何关系，可求得北东地坐标系与轨道坐标系之间的夹角α：

$\mathrm{\&chi;}=\frac{(Z_{J2000}\&times;r_{J2000})\&times;(-r_{J2000})}{\mathrm{norm}((Z_{J2000}\&times;r_{J2000})\&times;(-r_{J2000}))}---\left(28\right)$

其中χ为求解过程中的中间变量，norm()表示对括号内的矢量求模长；为单位化后的v_J2000，即z_J2000为J2000坐标系下的z_J2000轴单位向量z_J2000＝[0 0 1]^T；

北东地坐标系向轨道坐标系转换的过程中，北东地坐标系绕自身z_NED轴旋转的方向需进一步判定，将判定量记为β：

式中χ由公式(28)求得，若β>0，则α即为公式(29)所得结果，若β≤0，则α为公式(29)所得结果的相反数；

由北东地坐标系到轨道坐标系的坐标转换矩阵Q为：

$Q=\left[\begin{array}{ccc}\cos \mathrm{\&alpha;}& \sin \mathrm{\&alpha;}& 0\\ -\sin \mathrm{\&alpha;}& \cos \mathrm{\&alpha;}& 0\\ 0& 0& 1\end{array}\right]---\left(30\right)$

坐标转换矩阵Q点乘北东地坐标系下的太阳相对于微小卫星的位置矢量r_{sat_sun}，即得r_{sat_sun}在轨道坐标系中的表达r_{sat_sun_orbit}：