光學觀測中幾種太陽夾角計算方法及精度分析

虞炳文，婁廣國，蔡紅維，周小杰，楊 宏

(西昌衛星發射中心，四川 西昌 615000)

0 引言

在航天發射、低軌衛星長期在軌運行管理，以及天文觀測中，光學觀測是一個重要且不可或缺的的觀測方式，光學觀測具備別的測量控制設備所不具備的優勢，比如直觀，在光學影像中，可以很直觀的看見飛行器的外部狀態，另外，光學觀測的測角精度較高，可以作為飛行器在空間中位置的一個重要參考量。尤其在航天發射的起飛階段中，光學觀測設備是重要的輔助決策判決的手段，在發射場中得到了大量的運用。但是制約光學設備使用并發揮作用的限制條件也較為明顯，除了包括山體、云等遮蔽物之外，直面部分強光源的物體的直射，可能會對光學設備的感光模塊造成不可逆轉的破壞，其中較為明顯的，就是太陽以及月球的光線直射，因此，如何通過準確計算出飛行器與太陽的夾角，以達到有效規避強光直射的風險，避免太陽對光學設備的損壞，是文中將要論述的重點。

在計算飛行器與太陽夾角的過程中，有很多的算子可供選擇，可區分為簡單算子和復雜算子，簡單算子的優勢是只需要少量的參數即可計算出結果，但是精度相對較低，復雜算子的優勢是計算精度高，但是其需要裝填的參數較多，甚至有些參數需要專業設備才能測得，在日常使用中并不便捷。因此在日常使用中，通常會選用精度較高的簡單算子來替代復雜算子，本文將通過計算及精度分析，論證簡單算子中與復雜算子計算結果最為接近的算子組合，以達到在簡單便捷使用的同時，又能保證一定的精度。

1 計算太陽赤緯及太陽時角

首先需要介紹一下什么是太陽赤緯及太陽時角。

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;

}

2 太陽高度角及方位角的計算方法

太陽高度角和方位角，即在觀察者位置觀察太陽時的方位俯仰角，與當地的水平面相關，太陽高度角從地表橫切面起算，取值范圍為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;

}

3 飛行器與太陽夾角的計算方法

在計算出太陽高度角和方位角以后，同時已知飛行器的方位俯仰角，此時便能計算出飛行器與太陽之間的夾角。

計算夾角有兩種方法：一種是向量點乘的方法，另一種是三角函數推導的方法。

飛行器與太陽及觀測點之間的相互關系如圖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 飛行器與太陽夾角的精度分析

為比較各算法的計算結果，進行精度分析。

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算子的赤緯角一致性最好。

5 結束語

本次實驗實現了從計算太陽赤緯角和太陽時角，到計算太陽高度角和方位角，最后計算得到飛行器與太陽夾角。在C++環境中，實現了9種太陽赤緯角，6種太陽時角，2種太陽夾角以及相應的蒙氣差修正等過程的計算，并分析了9種太陽赤緯角算法，以及6種太陽時角算法對太陽夾角的影響，結果表明，所有簡單算法中，VSOP87的時角算子與Wang赤緯角算子的計算結果最為真實可靠。在日常使用的簡單計算中，推薦使用Wang赤緯角算子及VSOP87時角算子計算太陽夾角。