光学观测中几种太阳夹角计算方法及精度分析

虞炳文，娄广国，蔡红维，周小杰，杨 宏

(西昌卫星发射中心，四川 西昌 615000)

在航天发射、低轨卫星长期在轨运行管理，以及天文观测中，光学观测是一个重要且不可或缺的的观测方式，光学观测具备别的测量控制设备所不具备的优势，比如直观，在光学影像中，可以很直观的看见飞行器的外部状态，另外，光学观测的测角精度较高，可以作为飞行器在空间中位置的一个重要参考量。尤其在航天发射的起飞阶段中，光学观测设备是重要的辅助决策判决的手段，在发射场中得到了大量的运用。但是制约光学设备使用并发挥作用的限制条件也较为明显，除了包括山体、云等遮蔽物之外，直面部分强光源的物体的直射，可能会对光学设备的感光模块造成不可逆转的破坏，其中较为明显的，就是太阳以及月球的光线直射，因此，如何通过准确计算出飞行器与太阳的夹角，以达到有效规避强光直射的风险，避免太阳对光学设备的损坏，是文中将要论述的重点。

在计算飞行器与太阳夹角的过程中，有很多的算子可供选择，可区分为简单算子和复杂算子，简单算子的优势是只需要少量的参数即可计算出结果，但是精度相对较低，复杂算子的优势是计算精度高，但是其需要装填的参数较多，甚至有些参数需要专业设备才能测得，在日常使用中并不便捷。因此在日常使用中，通常会选用精度较高的简单算子来替代复杂算子，本文将通过计算及精度分析，论证简单算子中与复杂算子计算结果最为接近的算子组合，以达到在简单便捷使用的同时，又能保证一定的精度。

首先需要介绍一下什么是太阳赤纬及太阳时角。

1.1 太阳赤纬和时角的基本概念

要梳理清楚什么是太阳赤纬和时角，需要首先建立天球这个天文学概念。

1.1.1 天球坐标系的概念

首先需要理解的是，天球是一个无限大的球体，距离在这里没有意义，在天球坐标系中，天球的中心可以是任何的星球，虽然一般情况下，会将地球放在天球的中心。换言之，地球在天球的中心位置，此时的天球可以理解成，地球居民在地球上的观察视角，在某个位置，如果足够高，四周没有遮挡，观察者所能看见的就是半个天穹，即以观察者所在位置为切点的，与地球球体相切的一个无限延展的平面，理论上，因为观察者所观察的是整个宇宙，距离无穷远，此时观察者所站立的地球，其本身的半径相比较无穷大而言，可以忽略不计，因此，此时整个宇宙的所有星球，都可以区分为两种位置，一种是平面以上，一种是平面以下。

因此，可以简单地认为，观察者在任何一个位置，都能观察到其所在切平面以上的，即半个宇宙的所有星球，因为很难去准确知晓某个星球距离地球的位置，所以在天球坐标系中需要将距离的概念淡化，而淡化距离概念的最简单直接的方法，就是将所有的星球与地球的距离都假设成相同，即将所有的星球都放在同一个球面上，当忽略了距离，此时想要标识某个星球的位置，就只需要用到两个参数值，即一个用来标识其纵向上的位置，一个用来标识其横向上的位置的数值。

在天球坐标系中，选用何种标准来标识一个方位和一个俯仰上的值，其最简单直接的想法便是，既然都是球体，日常用来标识地球球体的经纬度体系，可否用来标识天球的度量体系？原理上是可以的，事实上也正是如此的。因此在天球坐标系的度量体系中，实际上选用了一套和地球大地坐标系同样的度量方式，区别在于地球坐标系中使用的经纬度和高程，而在天球坐标系中因为淡化了距离的概念，将所有的星球都假设在了同一距离的球面上，所以并不需要高程这一参数值，因为需要高程的地方所需的一个基本特征是能够测量高程，而在宇宙尺度下，距离是很难被获取的。将星球的距离都假设在同一个平面上的过程，也就是常说的投影。天球坐标系仅仅参考了地球大地坐标系的经纬度值，在天球坐标系中，这两个参数值被称之为赤经和赤纬。

在天球坐标系中，天球的中心是地球的中心。因为观察者虽然并不在地球中心观察，但是因为观察者所在的地球表面，其距离地球中心的距离，相比较天球无限大的距离来说，可以忽略不记，因此可以简单认为，观察者所在的位置就是地球中心。

1.1.2 天球坐标系

同样一个天球，根据定义的经过球心的水平面不同，可以区分为多种天球坐标系，在此重点介绍常见的天球坐标系有三种。

首先是经过球心的平面为赤道面时，该处所指赤道便是地球赤道经过往外无限延伸，最终势必将会与天球相交，此时便形成了第二赤道坐标系，此时的赤道面称之为天赤道，经地心，作一条与天赤道垂直的无限延长的线，与天球相交于两点，即图 1中所示P点(天北极)与P′点(天南极)。

图1 天球示意图

然后是将地球与太阳放在同一个平面上，所形成的平面，称之为黄道面，以此为基准，便是黄道坐标系，仅靠地球与太阳两点不能确认一个平面，还需要一到两个点，这另外两个点的选择，便是春分点和秋分点，即太阳直射点穿过赤道时的两个点。经过地心，作一条与黄道面垂直的直线，与天球相交于两点，便是图 1中的Q(黄北极)和Q′(黄南极)点。黄道面与赤道面相交所形成的夹角，便是黄赤交角。

最后还需要介绍的便是以观察者所在位置，所作的与地球相交的平面，称之为地平面，此时的地平面，在宇宙维度内，可以认为是经过地球地心的一个平面，过地心作一条垂直于地平面的直线，与天球相交于两点，即观察者的头顶与脚底，见图1中Z(天顶)和Z′(天底)。

1.1.3 对天球上的常用点线面的理解

对图 1中所涉及的标注进行解释。

F点指的是飞行器，G点指的时太阳，飞行器在此时，也可以当成一个映射在天球球面上的物体，M点为地球上观察者的位置，过天北极P点与天南极P′的赤经圈称之为天子午圈，天子午圈可以有很多条，但是通过天顶Z的，只有一个圈，称之为本初子午圈，即图 1中所画弧线PZP′所在圆。

C点为春分点，需要重点介绍一下，在后续计算过程中，很多的运算的数值，都与此概念相关。春分点时黄道面和赤道面的相交点，即太阳直射点从南向北经过赤道的那个点，这个点并不是一个实际上地球地理坐标点，是无法用经纬度表示的，春分点是两个平面的交点，始终处在两个平面的相交处。天文学上，将春分点作为很多天文概念的起算点，比如恒星时，比如太阳赤经。

太阳赤纬角为地心与太阳的连线与天赤道面的夹角，即图 1中δ角。

太阳赤经角，从春分点C向东，一直到太阳在天赤道上的投影线之间的夹角，即春分点所在的赤经圈与太阳所在的赤经圈之间的夹角，称之为赤经角，即图 1中的β角。

太阳时角，是从本初子午圈与天赤道相交的点向西，即太阳所在的赤经圈与本初子午圈之间的夹角，即图 1中的γ角。

最终需要求解的是太阳夹角，即从地球上某个位置观察飞行器和太阳时，太阳和飞行器之间的夹角，即图 1中的α角。

在此需要特别注意的是，上述涉及的天文概念中，前缀往往都带着参考的星球的名称，比如太阳，之所以这么解释，是因为可以不是太阳，可以是宇宙中任何一个天体，比如银河系的中心作为参考的银河坐标系。在本文中所涉及的未特别声明的天文概念，均以太阳作为参考。

1.2 太阳赤纬的计算方法

要计算太阳夹角，首先需要解决的问题，便是如何准确计算天球坐标系下，太阳赤纬及太阳时角的值。首先计算太阳赤纬的值。

计算太阳赤纬值的方法有很多，常见的简单算法有Cooper算法、Spencer算法、Stine算法、Bourges算法、Wang算法、Yu算法，上述六种算法均基于天文观测数据，通过公式拟合提出的数学公式，另外还可以通过对天文观测数据的数值拟合算法，求出拟合度更高的计算公式，只是通过该方法提出的公式，其系数会因年而异，所以不具备普适性，另外，李文提出通过傅里叶展开法提出对拟合公式的改进算法。除数学拟合算法之外，还有一种基于VSOP87行星理论，演算推导算法的方法。

下面结合C++代码语言，对太阳赤纬各算法进行逐个分析计算。

1.2.1 Bourges算子

Bourges是由Bourges提出的，Bourges算子的公式如下所示。其中ω为弧度参考系数，t为等效计日序数。太阳赤纬角的单位为弧度。

δ=0.3723+23.2567sin(ωt)+0.1149sin(2ωt)-0.1712sin(3ωt)-0.7580cos(ωt)+0.3656cos(2ωt)+0.02010cos(3ωt)

(1)

弧度参考系数如下所示：

(2)

计日序数的计算如下所示。其中dn为计日序数，即本年度的第几天。

t=dn-1-n0

(3)

n0的的计算公式如下，其中n表示年份，如本年度为2022年，则n=2022。

n0=78.801+[0.2422(n-1969)]-int[0.25(n-1969)]

(4)

Bourges算法实现的代码如下:

double Bourges(cDateTime date)

{

int n=getDaysOfYear(date);

int temp=0.25*(date.year-1969);

double n0=78.801+0.2422*(date.year-1969)-temp;

double t=n-1-n0;//等效计日序数

double w=2*M_PI/365.2422;//参考弧度系数

double dellta=0.3723+23.2567*sin(w*t)+0.1149*sin(2*w*t)-0.1712*sin(3*w*t)-0.7580*cos(w*t)+0.3656*cos(2*w*t)+0.0201*cos(3*w*t);

return dellta/180.0*M_PI;//弧度

}

1.2.2 Cooper算子

Cooper算子是由Cooper提出的，Cooper算子的公式如下所示。公式中n为该年度中的天数累计值，如1月11日，则n=11。太阳赤纬角的返回值单位是弧度。

(5)

δ为太阳赤纬角。代码实现如下：

double Cooper(cDateTime date)

{

int n=getDaysOfYear(date);

double dellta;

dellta=23.45*sin(2*M_PI*(284.0+n)/365.0);

return dellta/180.0*M_PI;//弧度

}

1.2.3 Spencer算子

Spencer算子是由Spencer提出的，其中太阳赤纬角的单位为弧度，θ为日角。算法公式如下所示：

δ=0.006918-0.399912cos(θ)+0.070257sin(θ)-0.006758cos(2θ)+0.000907sin(2θ)-0.002697cos(3θ)+0.00148sin(3θ)

日角的计算公式如下所示：

(6)

Spencer算子的实现代码如下：

double Spencer(cDateTime date)

{

int n=getDaysOfYear(date);

double gamma=2*M_PI*(n-1)/365.0;//日角

double dellta=0.006918-0.399912*cos(gamma)+0.070257*sin(gamma)-0.006758*cos(2*gamma)+0.000907*sin(2*gamma)-0.002697*cos(3*gamma)+0.00148*sin(3*gamma);

return dellta;//弧度

}

1.2.4 Yu算子

Yu算子是由Yu提出来的，是在Spencer算子的基础上做的简化，计算得到的太阳赤纬角单位为弧度。θ为日角，公式如下:

δ=0.006918-0.399912cos(θ)+0.070257sin(θ)-0.006758cos(2θ)+0.000907sin(2θ)

(7)

其中的日角的计算公式如下:

(8)

实现Yu算法的代码如下:

double Yu(cDateTime date)

{

int n=getDaysOfYear(date);

double theT=2*M_PI*(n-1)/365;//日角

double dellta=0.006918-0.399912*cos(theT)+0.070257*sin(theT)-0.006758*cos(2*theT)+0.000907*sin(2*theT);

return dellta;//弧度

}

1.2.5 Stine算子

Stine算子是由Stine提出的，其中太阳赤纬角的单位为弧度，其中n为计日序数，即本年度的第几天。

Stine算法的实现代码如下：

double Stine(cDateTime date)

{

int n=getDaysOfYear(date);

double dellta=asin(0.39795*cos(2*M_PI*(n-173)/365.242));

return dellta;//弧度

}

1.2.6 Wang算子

Wang算子是由王炳忠在Bourges的基础上，针对计日系数，修改了部分系数，并将起算年份从1969改为了1985，使得其更加贴近真实值。计算所得的太阳赤纬角的单位为弧度。公式如下所示，与Bourges算子相同。太阳赤纬角的单位为弧度。

δ=0.3723+23.2567sin(ωt)+0.1149sin(2ωt)-0.1712sin(3ωt)-0.7580cos(ωt)+0.3656cos(2ωt)+0.02010cos(3ωt)

(9)

其中的弧度参考系数也与Bourges相同。如下所示:

(10)

以及等效计日序数的算法也是相同的，如下:

t=dn-1-n0

(11)

与Bourges不同的是对n0的计算，如下所示：

n0=79.6764+[0.2422(n-1985)]-int[0.25(n-1985)]

(12)

Wang算法的实现代码如下所示:

double Wang(cDateTime date)

{

int n=getDaysOfYear(date);

int year=date.year;

int temp=0.25*(year-1985);

double n0=79.6764+0.2422*(year-1985)-temp;

double t=n-1-n0;//等效计日序数

double w=2*M_PI/365.2422;//参考弧度系数

double dellta=0.3723+23.2567*sin(w*t)+0.1149*sin(2*w*t)-0.1712*sin(3*w*t)-0.7580*cos(w*t)+0.3656*cos(2*w*t)+0.0201*cos(3*w*t);

return dellta/180.0*M_PI;//弧度

}

1.2.7 数值拟合法

数值拟合法是由李文提出的，根据2015年-2018年的天文年历，提出使用积日因子β作为多项式的中间系数，进行拟合计算。得到适2015-2018年的年份的改进公式。

太阳赤纬角的单位为弧度，其中太阳赤纬角的计算公式如下，其中的β为积日因子，ak为系数。

(13)

积日因子的求解公式如下所示。其中dn为积日系数。

(14)

n0的计算公式如下所示：

n0=79.6764+0.2422(n-1985)-int(n-1985)

(15)

数值拟合算法的代码如下所示:

double NumericalFitting(cDateTime date)

{

int n=getDaysOfYear(date);

int year=date.year;

int temp=year-1985;

double n0=79.6764+0.2422*(year-1985)-temp;

double b=2*M_PI*(n-1-n0)/365.2422;

int L=(year-2015)%4;

double dellta=0;

if(L==0)

{

dellta=calDellta(b,0.38835,22.911,-0.49055,-3.1217,0.043485,-0.19959,0.12879,0.045704,-0.040516,0.010031,-1.0946e-3,4.5142e-5);

}else if(L==1){

dellta=calDellta(b,0.38879,22.909,-0.49009,-3.1203,0.041657,-0.19950,0.12971,0.045246,-0.040482,0.010056,-1.1011e-3,4.5598e-5);

}else if(L==2){

dellta=calDellta(b,0.38769,22.909,-0.49277,-3.1178,0.046562,-0.20471,0.12953,0.047321,-0.041560,0.010309,-1.1302e-3,4.6935e-5);

}else if(L==3){

dellta=calDellta(b,0.38702,22.910,-0.49060,-3.1193,0.044708,-0.20270,0.12943,0.046646,-0.041166,0.010207,-1.1175e-3,4.6298e-5);

}

return dellta/180.0*M_PI;//弧度

}

double calDellta(double b, double a0,double a1, double a2, double a3, double a4, double a5, double a6, double a7, double a8, double a9, double a10, double a11)

{

double dellta=a0*pow(b,0)+a1*pow(b,1)+a2*pow(b,2)+a3*pow(b,3)+a4*pow(b,4)+a5*pow(b,5)+a6*pow(b,6)+a7*pow(b,7)+a8*pow(b,8)+a9*pow(b,9)+a10*pow(b,10)+a11*pow(b,11);

return dellta;

}

1.2.8 傅里叶展开法

傅里叶展开法是由李文提出的，在数值拟合法的基础上，进行的改进算法，通过观察连续的太阳赤纬角的变化规律，发现太阳赤纬角的变化具备周期性，每4年一个循环，因此将算法的周期，从1年修改为4年，定义算法中的年份为每个周期中的第三年，并且将计日序数从之前的1年中的第几天，改成整个4年周期中的第几天，即总的计日序数。

得到的太阳赤纬角的值，单位为弧度，其中β为积日因子。计算公式如下所示:

(16)

积日因子的计算公式如下所示，其中dn为积日系数。

(17)

式中,n0表达式如下所示:

n0=79.6764+0.2422(n-1985)-int(n-1985)

(18)

傅里叶展开法的计算公式如下所示:

double FourierBaseNF(cDateTime date)

{

int n=getDaysOfYear(date);

int year=date.year;

int y=0;//周期内的第3个年份

int dn=0;//周期内的总计日序数

if(year>=2015)

{

y=((year-2015)/4)*4+2017;

}else{

y=2013-((2015-year)/4)*4;

}

int temp=(year-3)%4;//这是周期内的第几年

if(temp==0)

{

dn=n;

}else if(temp==1)

{

if(QDate::isLeapYear(year-1))

{

dn=dn+366;

}else {

dn=dn+365;

}

dn=dn+n;

}else if(temp==2)

{

if(QDate::isLeapYear(year-2))

{

dn=dn+366;

}else {

dn=dn+365;

}

if(QDate::isLeapYear(year-1))

{

dn=dn+366;

}else {

dn=dn+365;

}

dn=dn+n;

}else if(temp==3)

{

if(QDate::isLeapYear(year-3))

{

dn=dn+366;

}else {

dn=dn+365;

}

if(QDate::isLeapYear(year-2))

{

dn=dn+366;

}else {

dn=dn+365;

}

if(QDate::isLeapYear(year-1))

{

dn=dn+366;

}else {

dn=dn+365;

}

dn=dn+n;

}

double n0=79.6764+0.2422*(y-1985)-int(y-1985);

double b=2*M_PI*(dn-1-n0)/365.2422;

double dellta=calFourierBaseNF(b,0.3783,-0.5624,0.3654,0.0156,-0.007662,-0.0005366,23.25,0.1082,-0.1705,-0.002773,0.003393);

return dellta/180.0*M_PI;//弧度

}

double calFourierBaseNF(double b,double a0,double a1,double a2,double a3,double a4,double a5,double b1,double b2,double b3,double b4,double b5)

{

double dellta=a0+a1*cos(b)+a2*cos(2*b)+a3*cos(3*b)+a4*cos(4*b)+a5*cos(5*b)+b1*sin(b)+b2*sin(2*b)+b3*sin(3*b)+b4*sin(4*b)+b5*sin(5*b);

return dellta;

}

1.2.9 VSOP87行星理论

VSOP87理论是由Bretagnon和Francou在1987年提出的，是VSOP82理论的进一步完善，VSOP理论是用来描述太阳系范围内的行星的轨道，在很长时间内的周期性变化的一种半分析理论。

VSOP形形理论除了可以算出精确的行星轨道参数，还可以直接计算出行星的位置，在此所言的位置，包含各种坐标系在内，如黄道坐标系。

计算VSOP87计算太阳赤纬角的过程，大致可分为七步。

1)计算儒略日数，用J表示。

2)计算相对儒略实际数，用T表示，如下所示:

(19)

3)计算太阳几何平黄经，用L表示，如下所示:

L=280.466456+36000.76982779×T+0.003032028×

(20)

4)计算平黄赤交角，用MB表示，如下所示:

(21)

5)计算太阳平近点角，用MSs表示，如下所示:

MSs=357.52191+35999.0503×T-

(22)

6)计算太阳真黄经，用SYS表示，如下所示:

SYS=L+(1.9146-0.004817×T-0.000014×T2)

(23)

7)计算太阳地心赤纬角，用Dec表示，如下所示:

(24)

上式中涉及的radeg为弧度与度的转换量纲，取radeg值为57.295 779 513 07。

用代码实现上述公式如下:

double VSOP87(cDateTime date)

{

//儒略日数

double J=date2JD(date);

//计算相对儒略世纪数

double T=(J- 2451545.0)/36525;

//太阳几何平黄经

double L=280.466456+36000.76982779*T+0.003032028*pow(T,2)+pow(T,3)/49931-pow(T,5)/15299;

//平黄赤交角

double MB=23.4392911111-(46.815/3600)*T-(0.00059/3600)*pow(T,2)+(0.001813/3600)*pow(T,3);

//太阳平近点角

double MSs=357.52191+35999.0503*T-0.0001559*pow(T,2)-0.00000048*pow(T,3);

//太阳真黄经

double SYS=L+(1.9146-0.004817*T-0.000014*pow(T,2))*sin(MSs/radeg)+(0.019993-0.000101*T)*sin(2*MSs/radeg)+0.00029*sin(3*MSs/radeg);

//太阳地心赤纬

double Dec=asin(sin(MB/radeg)*sin(SYS/radeg));

return Dec;

}

1.3 太阳时角的计算方法

太阳时角是以太阳为参考，其值代表了观察者所在的子午圈与太阳所在的子午圈之间的夹角。太阳时角的单位通常可以是角度值，如度分秒，此时的取值范围为0至360度，也可以是时间值，如时分秒，此时的取值范围为0至24小时。当用角度表示时，其表示的含义可以理解为两个子午圈之间的夹角，当使用时分秒表示时，其含义可以理解为太阳在时角时间之前通过观察者正上空。如时角为3.5h，则代表在3.5小时之前，太阳在3.5小时之前从观察者所在的正上空通过。

太阳时角的计算采用下列公式。其中tS为真太阳时。

(25)

真太阳时的计算采用下面中的公式计算，其中t为平太阳时，单位为小时。在天文观测中，很多天文概念都会加一个“平”字，其含义往往是代表着未加误差修正的简单值，加上修正以后，通常称之为“真”。eot即对时间的修正值。太阳时角的计算，其重点就在于误差修正值eot的计算。

(26)

下面将介绍6种计算太阳时角的算子，分别为Lamm算子，Spencer算子，Whillier算子，Wloof算子，Yu算子，以及VSOP87行星理论。其区别主要便在误差修正值的计算上。

1.3.1 Lamm算子

Lamm算子是由Lamm提出的，其特点是考虑了平年和闰年的区别，其公式如下。以四年为一个周期，其中n为全周期中的计日序数。Ak和Bk为系数值，在此采用杜春旭提出的参数值进行计算。计算所得eot值单位为分钟值，需要转化为小时进行下一步计算。

(27)

实现代码如下所示:

double Lamm(cDateTime mDateTime, double lon)

{

int n=getDaysOfYear(mDateTime);

int n0=mDateTime.year%4;

if(n0==0)

{

n=n;

}else if(n0==1)

{

n=n+366;

}else if(n0==2)

{

n=n+366+365;

}else if(n==3)

{

n=n+366+365+365;

}

double t =getHours(mDateTime);

double thet =2*M_PI*n/365.25;

//eot为分钟的单位

double eot=0.00020870*cos(thet*0)+0.0092869*cos(thet*1)-0.052258*cos(thet*2)-0.0013077*cos(thet*3)-0.0021867*cos(thet*4)-0.000151*cos(thet*5)

-0.12229*sin(thet*1)-0.15698*sin(thet*2)-0.0051602*sin(thet*3)-0.0029823*sin(thet*4)-0.00023463*sin(thet*5);

double ts=0;

ts=t+eot/60+(lon-120)/15;

//转为弧度。1个小时对应15°。

double w=M_PI/12*(ts-12);

return w;

}

1.3.2 Spencer算子

Spencer算子除了可以计算太阳赤纬，也提出了计算太阳时角的公式。其计算公式如下。其中θ为日角，单位为弧度。

eot=229.18[0.000075+0.001868cos(θ)-0.032077sin(θ)-0.014615cos(2θ)-0.04089sin(2θ)]

(28)

日角的计算公式如下，周期为1年，n为计日序数。

(29)

公式计算所得的eot值单位为分钟，需要转换为小时再进行下一步计算。

实现代码如下所示:

double SpencerTimeAngle(cDateTime mDateTime, double lon)

{

int n=getDaysOfYear(mDateTime);

double t =getHours(mDateTime);

double thet =2*M_PI*(n-1)/365;

double eot=229.18*(0.000075+0.001868*cos(thet)-0.032077*sin(thet)-0.014615*cos(2*thet)-0.04089*sin(2*thet));

double ts=0;

ts=t+eot/60+(lon-120)/15;

double w=M_PI/12*(ts-12);

return w;

}

1.3.3 Whillier算子

Whillier算子由Whillier提出，其中θ为日角，单位为弧度。计算所得eot单位为分钟，需要转换成小时。

eot=9.87sin(2θ)-7.53cos(θ)-1.5sin(θ)

(30)

日角的计算公式如下，周期为1年，n为计日序数。

(31)

实现代码如下:

double Whillier(cDateTime mDateTime, double lon)

{

int n=getDaysOfYear(mDateTime);

double t =getHours(mDateTime);

double thet =2*M_PI*(n-81)/364;

double eot=9.87*sin(2*thet)-7.53*cos(thet)-1.5*sin(thet);

double ts=0;

ts=t+eot/60+(lon-120)/15;

double w=M_PI/12*(ts-12);

return w;

}

1.3.4 Wloof算子

Wloof算子由Wloof提出，其中θ为日角，单位为弧度。计算所得eot单位为分钟，需要转换成小时。误差修正值计算公式如下:

eot=0.258cos(θ)-7.416sin(θ)-3.648cos(2θ)-9.228sin(2θ)

(32)

其中的日角公式如下，周期为1年，n为计日序数。

(33)

代码实现如下:

double Wloof(cDateTime mDateTime, double lon)

{

int n=getDaysOfYear(mDateTime);

double t =getHours(mDateTime);

double thet =2*M_PI*(n-1)/365.242;

double eot=0.258*cos(thet)-7.416*sin(thet)-3.648*cos(2*thet)-9.228*sin(2*thet);

double ts=0;

ts=t+eot/60+(lon-120)/15;

double w=M_PI/12*(ts-12);

return w;

}

1.3.5 Yu算子

Yu算子由Yu提出，除了可计算赤纬，可提出了计算太阳时角的公式，其中θ为日角，单位为弧度。计算所得eot单位为分钟，需要转换成小时。误差修正值计算公式如下:

eot=0.0172+0.4281cos(θ)-7.351sin(θ)-3.3495cos(2θ)-9.3619sin(2θ)

其中的日角公式如下，周期为1年，n为计日序数。

代码实现如下:

double YuTimeAngle(cDateTime mDateTime, double lon)

{

int n=getDaysOfYear(mDateTime);

double t =getHours(mDateTime);

double thet =2*M_PI*n/365;

double eot=0.0172+0.4281*cos(thet)-7.351*sin(thet)-3.3495*cos(2*thet)-9.3619*sin(2*thet);

double ts=0;

ts=t+eot/60+(lon-120)/15;

double w=M_PI/12*(ts-12);

return w;

}

1.3.6 VSOP87行星理论

关于VSOP87行星理论上文已有介绍，利用该理论计算的过程如下。

1)计算相对儒略世纪数，用T表示:

(34)

2)计算黄道与月球平轨道升交点黄经，用mw表示:

mw=125.04452-1934.136261×T+

(35)

3)计算太阳几何平黄经，用L表示:

L=280.466456+36000.76982779×T+

(36)

4)计算月球几何平黄经，用L1表示:

L1=218.3164591+481267.88134236×T-0.0013268×T2+0.0000019×T3

(37)

5)计算黄经章动，用NHJ表示:

(38)

6)计算平黄赤交角，用MB表示:

(39)

7)计算任意时刻格林尼治恒星时，并加上黄经章动修正，得到视恒星时用Q0表示，需要将Q0值转换到0～360°范围内：

Q0=280.4606183+360.98564736629×(J-2451545)+

(40)

8)计算太阳平近点角，用MSs表示:

MSs=357.52191+35999.0503×T-

(41)

9)计算太阳真黄经，用SYS表示:

SYS=L+(1.9146-0.004817×T-0.000014×T2)

(42)

10)计算太阳赤经，用RM表示:

(43)

11)计算当地太阳时角，用ω，需要将ω值转换到0～360°范围内，计算结果为度值，通常需要转换成弧度值，方便后续计算:

(44)

ω=Q0+lon-120-RM

(45)

实现上述公式的代码如下:

double VSOP87TimeAngle(cDateTime date,double lon)

{

//儒略日数

double J=date2JD(date);

//计算相对儒略世纪数

double T=(J- 2451545.0)/36525;

//黄道与月球平轨道升交点黄经

double mw=125.04452-1934.136261*T+0.0020708*pow(T,2)+pow(T,4)/450000;

//太阳几何平黄经

double L=280.466456+36000.76982779*T+0.003032028*pow(T,2)+pow(T,3)/49931-pow(T,5)/15299;

//月球几何平黄经

double L1=218.3164591+481267.88134236*T-0.0013268*pow(T,2)+0.0000019*pow(T,3);

//黄经章动

double NHJ=(-17.2/3600)*sin(mw/radeg)-(1.32/3600)*sin(2*L/radeg)-(0.23/3600)*sin(2*L1/radeg)+(0.21/3600)*sin(2*mw/radeg);

//平黄赤交角

double MB=23.4392911111-(46.815/3600)*T-(0.00059/3600)*pow(T,2)+(0.001813/3600)*pow(T,3);

//任意时刻格林尼治视恒星时，加黄经章动修正,得到视恒星时

double Q0=280.4606183+360.98564736629*(J-2451545)+0.000387933*pow(T,2)-pow(T,3)/38710000+NHJ*MB;

int n=Q0/360;

Q0=Q0-n*360;

//太阳平近点角

double MSs=357.52191+35999.0503*T-0.0001559*pow(T,2)-0.00000048*pow(T,3);

//太阳真黄经

double SYS=L+(1.9146-0.004817*T-0.000014*pow(T,2))*sin(MSs/radeg)+(0.019993-0.000101*T)*sin(2*MSs/radeg)+0.00029*sin(3*MSs/radeg);

//太阳地心赤经

double RM=radeg*atan2(cos(MB/radeg)*sin(SYS/radeg),cos(SYS/radeg));

//当地时角

double w =Q0+lon-120-RM;

int temp=w/360;

w=w-temp*360;

return w/180*M_PI;

}

太阳高度角和方位角，即在观察者位置观察太阳时的方位俯仰角，与当地的水平面相关，太阳高度角从地表横切面起算，取值范围为0～90°，水平角从正北方向起算，取值范围为0～360°。

计算太阳高度角时，需要考虑蒙气差的影响。

2.1 太阳高度角的计算

太阳高度角的计算公式如下，其中ψ为观察者所处地点的纬度值，δ为太阳赤纬角，e为地球椭圆扁率：

(46)

其中地球扁率公式为：

(47)

太阳高度角和方位角计算的代码如下:

cAE SunAltitudeAngle::calSunAltitudeAngle(cLocation mLocation, double hourAngle, double Dec)

{

cAE mSunAE;

double e=(earth_a-earth_b)/earth_a;

double cos_lat=cos(mLocation.dLatitude/radeg);

double cos_w=cos(hourAngle);

double cos_Dec=cos(Dec);

double sin_lat=sin(mLocation.dLatitude/radeg);

double sin_Dec=sin(Dec);

double up=cos_lat*cos_w*cos_Dec+sin_lat*sin_Dec;

double down=sqrt(1+(1/pow(1-e,2)-1)*pow(cos_Dec,2)*pow(sin_lat,2)+(pow(1-e,2)-1)*pow(cos_lat,2)*pow(sin_Dec,2));

double cos_phi=up/down;

double phi=acos(cos_phi);

double E=M_PI_2-phi;

double R0=MongolianQiDiff(1,E,16,880);//计算蒙气差，单位为角秒

E=E+(R0/3600)/radeg;

mSunAE.E=E;

double cos_A=(sin(E)*sin_lat-sin_Dec)/(cos(E)*cos_lat);

mSunAE.A=acos(cos_A);

return mSunAE;

}

2.2 太阳方位角的计算

太阳方位角的计算在计算得出高度角的基础上进行，公式如下，其中，α为太阳高度角，ψ为观察者所处纬度值，δ为太阳赤纬角。

(48)

2.3 蒙气差计算公式

蒙气差修正，即对因大气折射、温度、气压带来的误差进行修正。蒙气差的单位为角秒。在此采用李文提出的计算公式，以公式适用范围为依据，将蒙气差计算按俯仰角进行划分。下列公式输入值α为俯仰角，单位为弧度，输出值为修正值，输出为角秒。

2.3.1 蒙气差对大气折射的修正

2.3.1.1 俯仰角0～14°范围的修正

在此范围内，使用傅里叶展开法进行计算修正，修正系数采用真实记录值进行计算。计算公式如下。

(49)

实现代码如下：

double SunAltitudeAngle::Fourier(double alpha,double a[6],double b[5],double w)

{

double z=M_PI_2-alpha;

double R0=a[0]+a[1]*cos(tan(z)*w)+a[2]*cos(2*tan(z)*w)+a[3]*cos(3*tan(z)*w)+a[4]*cos(4*tan(z)*w)+a[5]*cos(5*tan(z)*w)+b[0]*sin(tan(z)*w)+b[1]*sin(2*tan(z)*w)+b[2]*sin(3*tan(z)*w)+b[3]*sin(4*tan(z)*w)+b[4]*sin(5*tan(z)*w);

return R0;

}

2.3.1.2 俯仰角15～45°范围的修正

在此范围内，采用数值拟合法进行修正，修正系数采用真实记录值进行计算，修正公式如下:

(50)

实现代码如下:

double SunAltitudeAngle::NumericalFitting(double alpha, double a[8])

{

double z=M_PI_2-alpha;

double R0=a[0]*pow(tan(z),0)+a[1]*pow(tan(z),1)+a[2]*pow(tan(z),2)+a[3]*pow(tan(z),3)+a[4]*pow(tan(z),4)+a[5]*pow(tan(z),5)+a[6]*pow(tan(z),6)+a[7]*pow(tan(z),7);

return R0;

}

2.3.1.3 俯仰角46～90°范围的修正

高俯仰范围内的修正方法较多，常见的有5种算法。

1)算法一的计算公式如下:

(51)

2)算法二的计算公式如下:

(51)

3)算法三的计算公式如下:

(52)

4)算法四的计算公式如下:

(53)

5)算法五的计算公式如下:

(54)

实现上述公式的代码如下:

if(type==1)

{

R0=60.2*tan(z);

}else if(type==2)

{

R0=(1819.08+194.89*alpha+1.47*pow(alpha,2)-0.042*pow(alpha,3))/(1+0.41*alpha+0.067*pow(alpha,2)+0.000085*pow(alpha,3));

}else if(type==3)

{

R0=1.02/tan((alpha+10.3/(alpha+5.11))*M_PI/180);

}else if(type==4)

{

R0=60.097*tan(z)+0.0109*pow(tan(z),2)-0.073*pow(tan(z),3)+0.002*pow(tan(z),4);

}else if(type==5)

{

R0=60.103*tan(z)-0.066*pow(tan(z),3)-0.00015*pow(tan(z),5);

}

2.3.2 蒙气差对温度、气压的修正

蒙气差还会受温度和湿度的影响，因此需要对蒙气差进行附加值修正。输入值为蒙气差值，单位为角秒，输出值为蒙气差修正值，单位为角秒。蒙气差修正公式如下。其中在高仰角(通常指大于45°)的情况下，大气折射对温度也会有一个影响，因此需要对温度再进行一次修正，其中μ为对温度的附加修正。

R=R0(1+μA+B)

(55)

其中:μ的修正公式如下:

(56)

低仰角时，无需对温度进行修正。蒙气差修正公式如下:

R=R0(1+A+B)

(57)

对温度的修正公式如下:

(58)

对气压的修正公式如下:

(59)

上述公式实现代码如下:

double SunAltitudeAngle::MongolianQiDiffPlus(double alpha,double R,double T,double P)

{

double R0=0;

double A=-0.00383*T/(1+0.00363*T);

double B=P/101324.72-1;

if(alpha<(45/radeg))

{

double a[6]={1.00011366,-8.95854338e-4,1.77241695e-3,-9.29720341e-5,2.00679856e-6,-1.5325798e-8};

double z=M_PI_2-alpha;

double mu=a[0]*pow(tan(z),0)+a[1]*pow(tan(z),1)+a[2]*pow(tan(z),2)+a[3]*pow(tan(z),3)+a[4]*pow(tan(z),4)+a[5]*pow(tan(z),5);

R0=R*(1+mu*A+B);

}else{

R0=R*(1+A+B);

}

return R0;

}

在计算出太阳高度角和方位角以后，同时已知飞行器的方位俯仰角，此时便能计算出飞行器与太阳之间的夹角。

计算夹角有两种方法：一种是向量点乘的方法，另一种是三角函数推导的方法。

飞行器与太阳及观测点之间的相互关系如图2。F点为飞行器，S点为太阳，P点为观察点，太阳夹角即图中所示∠FPS。

图2 太阳夹角示意图

计算时，为方便计算，取模的长度为1，α1为飞行器点的方位角，β1为飞行器的俯仰角，α2为太阳的方位角，β2为太阳的俯仰角。则F点方向可以假定一个模长为1点F′，其坐标可以表示为F′(cosβ1sinα1,cosβ1cosα1)，同理，在S点方向可以假定一个模长为1的点S′，其坐标可以表示为S′(cosβ2sinα2,cosβ2cosα2)。正东方向取为X轴正方向，正北方向取为Y轴正方向。

3.1 向量点乘法

向量点乘的计算公式如下:

(60)

根据F′与S′的值，代码实现如下:

double SolarAngle::dotProduct(cAE SunAE, cAE CraftAE)

{

cXYZ Sun_XYZ;

cXYZ CraftAE_XYZ;

double alpha1=SunAE.A;

double beta1=SunAE.E;

double alpha2=CraftAE.A;

double beta2=CraftAE.E;

Sun_XYZ.X=cos(beta1)*cos(alpha1);

Sun_XYZ.Y=cos(beta1)*sin(alpha1);

Sun_XYZ.Z=sin(beta1);

CraftAE_XYZ.X=cos(beta2)*cos(alpha2);

CraftAE_XYZ.Y=cos(beta2)*sin(alpha2);

CraftAE_XYZ.Z=sin(beta2);

double thet=acos(Sun_XYZ.X*CraftAE_XYZ.X+Sun_XYZ.Y*CraftAE_XYZ.Y+Sun_XYZ.Z*CraftAE_XYZ.Z);

return thet;

}

3.2 三角函数推导

在PF与PS两条线段上，从P点出发，取一个长度为1的线段，记作F′和S′，将F′和S′相连，成一个等腰三角形，计算F′S′长度，从而求出∠F′PS′，即太阳夹角。

经推导，计算公式结果如下:

(61)

代码实现如下:

double SolarAngle::trigonometric(cAE SunAE, cAE CraftAE)

{

double alpha1=SunAE.A;

double beta1=SunAE.E;

double alpha2=CraftAE.A;

double beta2=CraftAE.E;

double thet=2*asin(sqrt(2-2*cos(beta1)*cos(beta2)*cos(alpha1-alpha2)-2*sin(beta1)*sin(beta2))/2);

return thet;

}

为比较各算法的计算结果，进行精度分析。

4.1 精度分析方法

用以进行数值精度鉴定的参考数据，通常有两种来源，一是通过实际观察获取的实测记录值，另一种来源是精度更高的计算值。在本次计算赤经赤纬数值的精度鉴定中，采用第二种方法，以精度更高的计算值来作为参考值的办法。采用复杂算法SPA算法的计算值作为数据的参考值，SPA算法是综合考虑了海拔、压强、温度、倾斜度、方位角旋转、大气折射、时区、地球自转时间与地球时间之差等各方面因素的算法，其计算精度可达±0.000 3°，可以用来作为参考值对简单算法进行精度分析。

因此在对太阳夹角的数值分析时，选用SPA算法作为参考值。误差分析选用均值、方差、均方根三项指标进行分析。

飞行器飞行轨迹选用样本，按照每秒一个点，每秒在方位、俯仰方向上各新增0.000 1°，全时长为300秒设计。飞行时间设定为2018年8月8日上午10点整。起飞地点设置为东经102.241 897 39，北纬27.902 341 42。

数据样本为全飞行轨迹中，每个点的太阳夹角与参考值的差值。

均值即将样本值求均值。公式如下:

(62)

方差为将样本值与样本的均值相减，分别求平方后求和，再将平方和求均值。

(63)

均方根，将所有样本值分别求平方，再求平方和以后求平方和的均值，最后开平方。

(64)

分析对象为：9种太阳赤纬角算法，6种太阳时间算法。

4.2 精度分析结果

4.2.1 均值分析结果

均值分析结果如表1所示。

表1 均值计算结果

图3为折线图形式。横坐标为赤纬算子，每一条折现代表一种时角算法。

图3 均值折线图

均值能反映统计样本的大致的平均水准，均值越小，说明样本值整体越小。从表1均值计算结果分析可得以下结论。

1)不同的赤纬角算子下，VSOP87的时角计算结果偏差最小。

2)不同的时角算子下，Wang算法的赤纬角计算结果偏差最小。

4.2.2 方差分析结果

方差分析结果如表2所示。

表2 方差计算结果

图4为折线图形式。横坐标为赤纬算子，每一条折现代表一种时角算法。

图4 方差折线图

方差分析反映的是自变量对数值的影响。即对于均值的偏离度的一个分析，即数值的稳定程度。数值越小，说明数值越稳定，起伏越小。

1)在不同的赤纬角算子，VSOP87的时角计算结果波动最小。

2)在不同的时角算子下，Wang算子的赤纬角计算结果波动最小。

4.2.3 均方根分析结果

均方根分析结果如表3所示。

表3 均方根计算结果

图5为折线图形式。横坐标为赤纬算子，每一条折现代表一种时角算法。

图5 均方根折线图

均方根和均值一起分析，反映了样本的一致性情况。均方根和均值越接近，说明样本数据一致性越好。

1)在不同的赤纬角算子下，VSOP87的时角一致性最好。

2)在不同的时角算子下，Wang算子的赤纬角一致性最好。

本次实验实现了从计算太阳赤纬角和太阳时角，到计算太阳高度角和方位角，最后计算得到飞行器与太阳夹角。在C++环境中，实现了9种太阳赤纬角，6种太阳时角，2种太阳夹角以及相应的蒙气差修正等过程的计算，并分析了9种太阳赤纬角算法，以及6种太阳时角算法对太阳夹角的影响，结果表明，所有简单算法中，VSOP87的时角算子与Wang赤纬角算子的计算结果最为真实可靠。在日常使用的简单计算中，推荐使用Wang赤纬角算子及VSOP87时角算子计算太阳夹角。