CN104714243A - 近地轨道微小卫星所在位置地磁场强度的确定方法 - Google Patents

Abstract

近地轨道微小卫星所在位置地磁场强度的确定方法，涉及航天器所在位置地磁场强度的确定方法。为了解决目前还没有一种能够根据近地轨道微小卫星在J2000坐标系下的位置、速度与UTC解算其所在位置的地磁场强度与方向的方法的问题。本发明已知近地微小卫星在J2000坐标系下的速度、位置以及UTC，通过岁差转换矩阵P和地球自转转换矩阵R求得坐标转换矩阵W，求得微小卫星所在位置的地心经度λ、地心纬度与地心距R_d；利用地磁场球谐系数，由地磁场位函数理论，求得地磁场磁位并对磁位求偏导，将地磁场强度B在北东地坐标系下表达出来；解算出北东地坐标系与轨道坐标系间的夹角α，求得微小卫星所在位置的地磁场强度在轨道坐标系下的表达。本发明适用于航天器姿态控制技术领域。

Description

近地轨道微小卫星所在位置地磁场强度的确定方法

技术领域

本发明涉及航天器所在位置地磁场强度的确定方法，属于航天器姿态控制技术领域。

背景技术

解算地磁场强度对于近地轨道微小卫星姿态控制系统具有重要意义。在姿态确定方面，解算出特定坐标系下的地磁场方向矢量，可配合其他矢量(如太阳相对于微小卫星的位置矢量等)进行双矢量定姿；在姿态机动方面，可以通过磁力矩器将地磁力矩作为环境力矩加以利用，用以提供姿态控制力矩，或为飞轮卸载等。很多微小卫星的姿控系统是由星载GPS提供自身在J2000坐标系下的速度、位置以及时间信息(UTC)，进行姿态确定与机动。然而目前还没有一种能够根据近地轨道微小卫星在J2000坐标系下的位置、速度以及UTC信息解算其所在位置的地磁场强度与方向的方法。

发明内容

本发明为了解决目前还没有一种能够根据近地轨道微小卫星在J2000坐标系下的位置、速度以及UTC信息解算其所在位置的地磁场强度与方向的方法的问题。

近地轨道微小卫星所在位置地磁场强度的确定方法，包括以下步骤：

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

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

步骤三、利用国际地磁参考场(International Geomagnetic Reference Field，简称IGRF)2010年更新的地磁场球谐系数，根据地磁场位函数理论，将地磁场的磁位表达出来；

步骤四、对地磁场的磁位求偏导，得到地磁场的北向分量、东向分量与垂直分量，将地磁场强度B在北东地坐标系下表达出来；

步骤五、利用J2000坐标系下微小卫星的位置和速度，解算出北东地坐标系与轨道坐标系之间的夹角α，可得北东地坐标系与轨道坐标系之间的坐标转换矩阵，进而得到微小卫星所在位置的地磁场强度在轨道坐标系下的表达式。

本发明具有以下有益效果：本发明提出了一种根据近地轨道微小卫星在J2000坐标系下的位置、速度以及UTC信息，解算其所在位置的地磁场强度与方向，并将其在轨道坐标系中表示出来的方法。用本发明提供的方法解算出的轨道坐标系下微小卫星所在位置的地磁场强度(N＝6时)，与利用STK生成的相同轨道参数下的地磁场强度，有很高的吻合度，地磁场强度的指向精度平均误差不超过0.5°，最高不超过2.3°。

附图说明

图1近地轨道微小卫星所在位置地磁场强度的确定方法流程图；

图2地磁场B_{orbit_uni_x}与B_{orbit_STK_x}对比效果图；

图3地磁场B_{orbit_uni_y}与B_{orbit_STK_y}对比效果图；

图4地磁场B_{orbit_uni_z}与B_{orbit_STK_z}对比效果图；

图5地磁场强度单位矢量方向的指向误差效果图。

具体实施方式

具体实施方式一：近地轨道微小卫星所在位置地磁场强度的确定方法，于包括以下步骤：

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

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

步骤三、利用国际地磁参考场(International Geomagnetic Reference Field，简称IGRF)2010年更新的地磁场球谐系数，根据地磁场位函数理论，将地磁场的磁位表达出来；

步骤四、对地磁场的磁位求偏导，得到地磁场的北向分量、东向分量与垂直分量，将地磁场强度B在北东地坐标系下表达出来；

步骤五、利用J2000坐标系下微小卫星的位置和速度，解算出北东地坐标系与轨道坐标系之间的夹角α，可得北东地坐标系与轨道坐标系之间的坐标转换矩阵，进而得到微小卫星所在位置的地磁场强度在轨道坐标系下的表达式。

具体实施方式二：本实施方式的步骤一的具体操作步骤如下：

首先给出J2000坐标系与ECEF的定义：

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

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

由J2000坐标系到ECEF需进行四次坐标转换，依次为岁差转换，章动转换，地球自转转换，地球极移转换；本方法只考虑在由J2000坐标系到ECEF的坐标转换中起决定性作用的岁差转换矩阵(记为P)与地球自转转换矩阵(记为R)；

岁差转换矩阵P的求取如下：

首先利用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{hour}}{24}\end{array}---\left(1\right)$

其中：year，month，day分别为格林尼治年、月、日，hour，minute，second分别为UTC小时，分钟，秒；int(·)表示对括号里的式子取整；

利用儒略日JD求得从历元J2000.0起算的儒略世纪数T：

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

三个岁差角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{\θ}}_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(4\right);$

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

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

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

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

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

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

在忽略章动和地球极移影响的情况下，J2000坐标系经过岁差转换之后，再绕其z'轴逆时针转动θ'_GST，即与ECEF重合；J2000坐标系经过岁差转换后的坐标系三轴表示分别为x',y',z'；

地球自转转换矩阵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(8\right)$

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

W＝R·P    (9)

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

$\begin{array}{c}{r}_{J2000}={\left[\begin{array}{ccc}{r}_{J2000\_x}& r_{J2000\_y}& r_{J2000\_z}\end{array}\right]}^T\\ v_{J2000}={\left[\begin{array}{ccc}{v}_{J2000\_x}& v_{J2000\_y}& v_{J2000\_z}\end{array}\right]}^T\end{array}---\left(10\right)$

其中，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下，微小卫星相对于地心的位置表示为r_ECEF，微小卫星相对于地心的速度表示为v_ECEF：

$\begin{array}{c}{r}_{\mathrm{ECEF}}={\left[\begin{array}{ccc}{r}_{\mathrm{ECEF}\_x}& r_{\mathrm{ECEF}\_y}& r_{\mathrm{ECEF}\_z}\end{array}\right]}^T\\ r_{\mathrm{ECEF}}={\left[\begin{array}{ccc}{v}_{\mathrm{ECEF}\_x}& v_{\mathrm{ECEF}\_y}& v_{\mathrm{ECEF}\_z}\end{array}\right]}^T\end{array}---\left(11\right)$

其中，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得到：

                                              (12)。

r_ECEF＝W·r_J2000

v_ECEF＝W·v_J2000

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

具体实施方式三：本实施方式的步骤二的具体操作步骤如下：

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

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

$R_d=\sqrt{{r_{\mathrm{ECEF}\_x}}^2+{r_{\mathrm{ECEF}\_y}}^2+{r_{\mathrm{ECEF}\_z}}^2}$

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

对于地心经度，由于反三角函数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(14\right)$

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

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

具体实施方式四：本实施方式的步骤三的具体操作步骤如下：

由高斯球谐波分析法得到的地磁场磁位表达式：

其中，R_e为地球半径，R_d为地心距，为高斯公式系数，N为地磁场模型的阶数，为n阶m次施密特Legendre多项式。

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

具体实施方式五：本实施方式的

P_n,m(·)其求解过程如下：

$P_{n,m}\left(x\right)={(1-x^2)}^{m/2}{\stackrel{\&OverBar;}{P}}_{n,m}\left(x\right)$

${\stackrel{\&OverBar;}{P}}_{n,m}\left(x\right)={\mathrm{\&epsiv;}}_{n,m}{P}_n^m\left(x\right)$

$P_n^m\left(x\right)=\frac{d^m{P}_n\left(x\right)}{d{x}^m}---\left(16\right)$

${\mathrm{\&epsiv;}}_{n,m}=\sqrt{K_m\frac{(n-m)!}{(n+m)!}},K_m=\left\{\begin{array}{c}1(m=0)\\ 2(m\&GreaterEqual;1)\end{array}\right.$

其中，x为函数自变量，P_n(x)为Legendre多项式，其函数表达式为：

$P_n\left(x\right)=\frac{1}{2^{n}n!}\frac{d^n{(x^2-1)}^n}{d{x}^n}---\left(17\right).$

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

具体实施方式六：本实施方式的地磁场模型的阶数N根据实际需要选择阶数。

其他步骤和参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式的地磁场模型的阶数N取6。

其他步骤和参数与具体实施方式一至六之一相同。

具体实施方式八：本实施方式的步骤三中的高斯公式系数的确定方法如下：

高斯公式系数由IGRF官方网站每五年更新一次系数及其导数，

本模型(N＝6)采用2010年更新的高斯系数：

$\begin{array}{c}G=\left[\begin{array}{ccccccc}-29496.5& -1585.9& 0& 0& 0& 0& 0\\ -2396.6& 3026& 1668.6& 0& 0& 0& 0\\ 1339.7& -2326.3& 1231.7& 634.2& 0& 0& 0\\ 912.6& 809& 166.6& -357.1& 89.7& 0& 0\\ -231.1& 357.2& 200.3& -141.2& -163.1& -7.7& 0\\ 72.8& 68.6& 76& -141.4& -22.9& 13.1& -77.9\end{array}\right]\\ H=\left[\begin{array}{ccccccc}0& 4945.1& 0& 0& 0& 0& 0\\ 0& -2707.7& -575.4& 0& 0& 0& 0\\ 0& -160.5& 251.7& -536.8& 0& 0& 0\\ 0& 286.4& -211.2& 164.4& -309.2& 0& 0\\ 0& 44.7& 188.9& -118.1& 0.1& 100.9& 0\\ 0& -20.8& 44.2& 61.5& -66.3& 3.1& 54.9\end{array}\right]\end{array}---\left(18\right)$

这里G(i,j)表示G中的第i行，第j列元素，同理于H(i,j)。

其他步骤和参数与具体实施方式一至七之一相同。

具体实施方式九：本实施方式的步骤四的具体操作步骤如下：

地磁场强度B可以表示成标量磁位的负梯度：

其中为梯度算子符号；

首先定义北东地坐标系：坐标原点O_NED位于目标点处(本模型对应微小卫星质心所在位置)，x_NED即北轴，沿着天球经圈指向正北方向；y_NED即东轴，沿着天球纬圈指向正东方向；z_NED即地轴，垂直于天球表面并指向地心方向；此处天球是以地心为圆心，以地心到目标点(微小卫星所在位置)的距离为半径的假想球；

北东地坐标系与地心距离R_d、地心经度λ与地心纬度存在如下关系：

dz_NED＝-dR_d

其中，dx_NED，dy_NED，dz_NED分别为北东地坐标系三个轴单位向量的微元形式；

地磁场强度B在北东地坐标系下的分量表达式为：

$B_{z_{\mathrm{NED}}}=\frac{\&PartialD;V}{\&PartialD;R_d}$

其中，V即为是由公式(15)求得的磁位；B_xNED，B_yNED，B_zNED分别为地磁场强度B在北东地坐标系下的分量；设地磁场强度在北东地坐标系下的表达为B_NED，则：

$B_{\mathrm{NED}}={\left[\begin{array}{ccc}{B}_{x_{\mathrm{NED}}}& B_{y_{\mathrm{NED}}}& B_{z_{\mathrm{NED}}}\end{array}\right]}^T---\left(22\right).$

其他步骤和参数与具体实施方式一至八之一相同。

具体实施方式十：本实施方式的步骤五的具体操作步骤如下：

首先定义轨道坐标系：坐标原点O_orbit位于微小卫星质心，z_orbit轴指向地心；x_orbit轴位于微小卫星轨道平面内，与z_orbit轴垂直且指向微小卫星速度方向；y_orbit轴与其他两轴形成右手系；

北东地z_NED轴与轨道坐标系z_orbit轴重合，因而可通过一次绕北东地坐标系z_NED轴的旋转使北东地坐标系与轨道坐标系重合；设绕z_NED轴旋转的角度为α，α即为北东地坐标系与轨道坐标系之间的夹角，求出α即可得到由北东地坐标系到轨道坐标系的坐标转换矩阵；

利用微小卫星所在的子午面与微小卫星轨道面之间的几何关系，可求得北东地坐标系与轨道坐标系之间的夹角α：

$\mathrm{\&phi;}=\frac{(z_{J2000}\&times;r_{J2000})\&times;(-r_{J2000})}{\mathrm{norm}((z_{J2000}\&times;r_{J2000})\&times;(-r_{J2000}))}---\left(23\right)$

$\mathrm{\&alpha;}=a\cos ({{\stackrel{\&OverBar;}{v}}_{J2000}}^T\&CenterDot;\mathrm{\&phi;})$

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

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

$\mathrm{\&beta;}={(\mathrm{\&phi;}\&times;{\stackrel{\&OverBar;}{v}}_{J2000})}^T\&CenterDot;(-r_{J2000})---\left(24\right)$

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

由北东地坐标系到轨道坐标系的坐标转换矩阵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(25\right)$

坐标转换矩阵Q点乘北东地坐标系下的地磁场B_NED，即得到地磁场强度在轨道坐标系下的表达B_orbit：

B_orbit＝Q·B_NED    (26)。

在实际应用过程中，根据需要对其进行单位化或求模长处理。

其他步骤和参数与具体实施方式一至九之一相同。

实施例

设某近地轨道微小卫星于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)，以此作为输入，利用专利中的方法，设定N＝6，解算出对应时刻地磁场强度在轨道坐标系下的表达式，此表达式经单位化后记为B_{orbit_uni}＝[B_{orbit_uni_x} B_{orbit_uni_y} B_{orbit_uni_z}]^T，其中B_{orbit_uni_x}，B_{orbit_uni_y}，B_{orbit_uni_z}分别为B_{orbit_uni}在轨道坐标系下的三个分量。验证时间为2015年5月1日00：00：01至2015年5月1日13：00：01。将B_{orbit_uni}与利用STK生成的轨道坐标系下的地磁场强度单位矢量(生成路径：STK/Vectors(VVLH)/MagField(IGRF)/x/Magnitude y/Magnitudez/Magnitude)，记为B_{orbit_STK}＝[B_{orbit_STK_x} B_{orbit_STK_y} B_{orbit_STK_z}]^T，其中B_{orbit_STK_x},B_{orbit_STK_y},B_{orbit_STK_z}分别为B_{orbit_STK}在轨道坐标系下的三个分量，进行对比，三轴的对比情况分别如图2-图4所示；

以上图2-4中，黑色线条分别表示B_{orbit_uni_x}，B_{orbit_uni_y}，B_{orbit_uni_z}。灰色线条分别表示B_{orbit_STK_x}B_{orbit_STK_y}，B_{orbit_STK_z}。为了辨别，MATLAB仿真作图时黑色线条比灰色线条宽。由图可知，由专利中算法解算得到的地磁场强度，与STK在相同轨道条件下生成的地磁场强度具有很高的吻合度。

为了进一步验证该专利提供的算法解算出的轨道坐标系下地磁场强度的准确性，求取矢量B_{orbit_uni}，与矢量B_{orbit_STK}之间的夹角，此夹角体现由本算法解算出的地磁场强度矢量方向的指向误差，指向误差如下图5所示；

由图5可知，由本专利中提出的算法得到的地磁场强度方向的指向误差较小，平均指向误差角不超过0.5°，最高误差不超过2.3°。

Claims (10)

1.近地轨道微小卫星所在位置地磁场强度的确定方法，其特征在于包括以下步骤：

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

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

步骤三、利用国际地磁参考场2010年更新的地磁场球谐系数，根据地磁场位函数理论，将地磁场的磁位表达出来；

步骤四、对地磁场的磁位求偏导，得到地磁场的北向分量、东向分量与垂直分量，将地磁场强度B在北东地坐标系下表达出来；

步骤五、利用J2000坐标系下微小卫星的位置和速度，解算出北东地坐标系与轨道坐标系之间的夹角α，可得北东地坐标系与轨道坐标系之间的坐标转换矩阵，进而得到微小卫星所在位置的地磁场强度在轨道坐标系下的表达式。

2.根据权利要求1所述的近地轨道微小卫星所在位置地磁场强度的确定方法，其特征在于：步骤一的具体操作步骤如下：

岁差转换矩阵P的求取如下：

首先利用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{hour}}{24}\end{array}---\left(1\right)$

其中：year，month，day分别为格林尼治年、月、日，hour，minute，second分别为UTC小时，分钟，秒；int(·)表示对括号里的式子取整；

利用儒略日JD求得从历元J2000.0起算的儒略世纪数T：

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

三个岁差角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 \left({-z}_A\right)& 0\\ -\sin \left({-z}_A\right)& \cos \left({-z}_A\right)& 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 \left({-\mathrm{\&zeta;}}_A\right)& \sin \left({-\mathrm{\&zeta;}}_A\right)& 0\\ -\sin \left({-\mathrm{\&zeta;}}_A\right)& \cos \left({-\mathrm{\&zeta;}}_A\right)& 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重合；J2000坐标系经过岁差转换后的坐标系三轴表示分别为x′,y′,z′；

地球自转转换矩阵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得到：

W＝R·P                           (9)

在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下，微小卫星相对于地心的位置表示为r_ECEF，微小卫星相对于地心的速度表示为v_ECEF：

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

                                  (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得到：

r_ECEF＝W·r_J2000

                                    (12)。

v_ECEF＝W·v_J2000

3.根据权利要求1或2所述的近地轨道微小卫星所在位置地磁场强度的确定方法，其特征在于：步骤二的具体操作步骤如下：

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

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

$R_d=\sqrt{{r_{\mathrm{ECEF}\_x}}^2+{r_{\mathrm{ECEF}\_y}}^2+{r_{\mathrm{ECEF}\_z}}^2}$

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

对于地心经度，判断和转化过程如下：

$\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所述的近地轨道微小卫星所在位置地磁场强度的确定方法，其特征在于：步骤三的具体操作步骤如下：

由高斯球谐波分析法得到的地磁场磁位表达式：

其中，R_e为地球半径，R_d为地心距，为高斯公式系数，N为地磁场模型的阶数，为n阶m次施密特Legendre多项式。

5.根据权利要求4所述的近地轨道微小卫星所在位置地磁场强度的确定方法，其特征在于：P_n,m(·)其求解过程如下：

$P_{n,m}\left(x\right)={(1-x^2)}^{m/2}{\stackrel{\&OverBar;}{P}}_{n,m}\left(x\right)$

${\stackrel{\&OverBar;}{P}}_{n,m}\left(x\right)={\mathrm{\&epsiv;}}_{n,m}{P}_n^m\left(x\right)$

$P_n^m\left(x\right)=\frac{d^m{P}_n\left(x\right)}{{\mathrm{dx}}^m}---\left(16\right)$

${\mathrm{\&epsiv;}}_{n,m}=\sqrt{K_m\frac{(n-m)!}{(n+m)!}},K_m=\left\{\begin{array}{c}1(m=0)\\ 2(m\&GreaterEqual;1)\end{array}\right.$

其中，x为函数自变量，P_n(x)为Legendre多项式，其函数表达式为：

$P_n\left(x\right)=\frac{1}{2^{n}n!}\frac{d^n{(x^2-1)}^n}{{\mathrm{dx}}^n}---\left(17\right).$

6.根据权利要求5所述的近地轨道微小卫星所在位置地磁场强度的确定方法，其特征在于：地磁场模型的阶数N根据实际需要选择阶数。

7.根据权利要求6所述的近地轨道微小卫星所在位置地磁场强度的确定方法，其特征在于：地磁场模型的阶数N取6。

8.根据权利要求7所述的近地轨道微小卫星所在位置地磁场强度的确定方法，其特征在于：步骤三中的高斯公式系数的确定方法如下：

高斯公式系数由IGRF官方网站每五年更新一次系数及其导数，

本模型采用2010年更新的高斯系数：

$G=\left[\begin{array}{ccccccc}-29496.5& -1585.9& 0& 0& 0& 0& 0\\ -2396.6& 3026& 1668.6& 0& 0& 0& 0\\ 1339.7& -2326.3& 1231.7& 634.2& 0& 0& 0\\ 912.6& 809& 166.6& -357.1& 89.7& 0& 0\\ -231.1& 357.2& 200.3& -141.2& -163.1& -7.7& 0\\ 72.8& 68.6& 76& -141.4& -22.9& 13.1& -77.9\end{array}\right]$

                                                          (18)

$H=\left[\begin{array}{ccccccc}0& 4945.1& 0& 0& 0& 0& 0\\ 0& -2707.7& -575.4& 0& 0& 0& 0\\ 0& -160.5& 251.7& -536.8& 0& 0& 0\\ 0& 286.4& -211.2& 164.4& -309.2& 0& 0\\ 0& 44.7& 188.9& -118.1& 0.1& 100.9& 0\\ 0& -20.8& 44.2& 61.5& -66.3& 3.1& 54.9\end{array}\right]$

这里G(i,j)表示G中的第i行，第j列元素，同理于H(i,j)。

9.根据权利要求8所述的近地轨道微小卫星所在位置地磁场强度的确定方法，其特征在于：步骤四的具体操作步骤如下：

地磁场强度B可以表示成标量磁位的负梯度：

其中为梯度算子符号；

北东地坐标系与地心距离R_d、地心经度λ与地心纬度存在如下关系：

dz_NED＝-dR_d

其中，dx_NED，dy_NED，dz_NED分别为北东地坐标系三个轴单位向量的微元形式；

地磁场强度B在北东地坐标系下的分量表达式为：

$B_{z_{\mathrm{NED}}}=\frac{\&PartialD;V}{{\&PartialD;R}_d}$

其中，V即为是由公式(15)求得的磁位；分别为地磁场强度B在北东地坐标系下的分量；设地磁场强度在北东地坐标系下的表达为B_NED，则：

$B_{\mathrm{NED}}{\left[\begin{array}{ccc}{B}_{x_{\mathrm{NED}}}& B_{y_{\mathrm{NED}}}& B_{z_{\mathrm{NED}}}\end{array}\right]}^T---\left(22\right),$

10.根据权利要求9所述的近地轨道微小卫星所在位置地磁场强度的确定方法，其特征在于：步骤五的具体操作步骤如下：

北东地z_NED轴与轨道坐标系z_orbit轴重合，因而可通过一次绕北东地坐标系z_NED轴的旋转使北东地坐标系与轨道坐标系重合；设绕z_NED轴旋转的角度为α，α即为北东地坐标系与轨道坐标系之间的夹角；

利用微小卫星所在的子午面与微小卫星轨道面之间的几何关系，可求得北东地坐标系与轨道坐标系之间的夹角α：

$\mathrm{\&phi;}=\frac{(z_{j2000}\&times;r_{J2000})\&times;\left({-r}_{J2000}\right)}{\mathrm{norm}((z_{J2000}\&times;r_{J2000})\&times;\left({-r}_{J2000}\right))}---\left(23\right)$

$\mathrm{\&alpha;}=a\cos ({{\stackrel{\&OverBar;}{v}}_{J2000}}^T\&CenterDot;\mathrm{\&phi;})$

其中，φ为求解过程中的中间变量，norm(·)表示对括号内的矢量求模长，为单位化后的v_J2000，即 ${\stackrel{\&OverBar;}{v}}_{J2000}=\frac{{\left[\begin{array}{ccc}{v}_{J2000\_x}& v_{J2000\_y}& v_{J2000\_z}\end{array}\right]}^T}{\mathrm{norm}\left(v_{J2000}\right)},$ z_J2000为J2000坐标系下的z_J2000轴单位向量z_J2000＝[0 0 1]^T；

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

$\mathrm{\&beta;}={(\mathrm{\&phi;}\&times;{\stackrel{\&OverBar;}{v}}_{J2000})}^T\&CenterDot;(-r_{J2000})---\left(24\right)$

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

由北东地坐标系到轨道坐标系的坐标转换矩阵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(25\right)$

坐标转换矩阵Q点乘北东地坐标系下的地磁场B_NED，即得到地磁场强度在轨道坐标系下的表达B_orbit：

B_orbit＝Q·B_NED      (26)。