光學(xué)觀測(cè)中幾種太陽(yáng)夾角計(jì)算方法及精度分析

虞炳文，婁廣國(guó)，蔡紅維，周小杰，楊 宏

(西昌衛(wèi)星發(fā)射中心，四川 西昌 615000)

0 引言

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

在計(jì)算飛行器與太陽(yáng)夾角的過(guò)程中，有很多的算子可供選擇，可區(qū)分為簡(jiǎn)單算子和復(fù)雜算子，簡(jiǎn)單算子的優(yōu)勢(shì)是只需要少量的參數(shù)即可計(jì)算出結(jié)果，但是精度相對(duì)較低，復(fù)雜算子的優(yōu)勢(shì)是計(jì)算精度高，但是其需要裝填的參數(shù)較多，甚至有些參數(shù)需要專業(yè)設(shè)備才能測(cè)得，在日常使用中并不便捷。因此在日常使用中，通常會(huì)選用精度較高的簡(jiǎn)單算子來(lái)替代復(fù)雜算子，本文將通過(guò)計(jì)算及精度分析，論證簡(jiǎn)單算子中與復(fù)雜算子計(jì)算結(jié)果最為接近的算子組合，以達(dá)到在簡(jiǎn)單便捷使用的同時(shí)，又能保證一定的精度。

1 計(jì)算太陽(yáng)赤緯及太陽(yáng)時(shí)角

首先需要介紹一下什么是太陽(yáng)赤緯及太陽(yáng)時(shí)角。

1.1 太陽(yáng)赤緯和時(shí)角的基本概念

要梳理清楚什么是太陽(yáng)赤緯和時(shí)角，需要首先建立天球這個(gè)天文學(xué)概念。

1.1.1 天球坐標(biāo)系的概念

首先需要理解的是，天球是一個(gè)無(wú)限大的球體，距離在這里沒(méi)有意義，在天球坐標(biāo)系中，天球的中心可以是任何的星球，雖然一般情況下，會(huì)將地球放在天球的中心。換言之，地球在天球的中心位置，此時(shí)的天球可以理解成，地球居民在地球上的觀察視角，在某個(gè)位置，如果足夠高，四周沒(méi)有遮擋，觀察者所能看見(jiàn)的就是半個(gè)天穹，即以觀察者所在位置為切點(diǎn)的，與地球球體相切的一個(gè)無(wú)限延展的平面，理論上，因?yàn)橛^察者所觀察的是整個(gè)宇宙，距離無(wú)窮遠(yuǎn)，此時(shí)觀察者所站立的地球，其本身的半徑相比較無(wú)窮大而言，可以忽略不計(jì)，因此，此時(shí)整個(gè)宇宙的所有星球，都可以區(qū)分為兩種位置，一種是平面以上，一種是平面以下。

因此，可以簡(jiǎn)單地認(rèn)為，觀察者在任何一個(gè)位置，都能觀察到其所在切平面以上的，即半個(gè)宇宙的所有星球，因?yàn)楹茈y去準(zhǔn)確知曉某個(gè)星球距離地球的位置，所以在天球坐標(biāo)系中需要將距離的概念淡化，而淡化距離概念的最簡(jiǎn)單直接的方法，就是將所有的星球與地球的距離都假設(shè)成相同，即將所有的星球都放在同一個(gè)球面上，當(dāng)忽略了距離，此時(shí)想要標(biāo)識(shí)某個(gè)星球的位置，就只需要用到兩個(gè)參數(shù)值，即一個(gè)用來(lái)標(biāo)識(shí)其縱向上的位置，一個(gè)用來(lái)標(biāo)識(shí)其橫向上的位置的數(shù)值。

在天球坐標(biāo)系中，選用何種標(biāo)準(zhǔn)來(lái)標(biāo)識(shí)一個(gè)方位和一個(gè)俯仰上的值，其最簡(jiǎn)單直接的想法便是，既然都是球體，日常用來(lái)標(biāo)識(shí)地球球體的經(jīng)緯度體系，可否用來(lái)標(biāo)識(shí)天球的度量體系？原理上是可以的，事實(shí)上也正是如此的。因此在天球坐標(biāo)系的度量體系中，實(shí)際上選用了一套和地球大地坐標(biāo)系同樣的度量方式，區(qū)別在于地球坐標(biāo)系中使用的經(jīng)緯度和高程，而在天球坐標(biāo)系中因?yàn)榈司嚯x的概念，將所有的星球都假設(shè)在了同一距離的球面上，所以并不需要高程這一參數(shù)值，因?yàn)樾枰叱痰牡胤剿璧囊粋€(gè)基本特征是能夠測(cè)量高程，而在宇宙尺度下，距離是很難被獲取的。將星球的距離都假設(shè)在同一個(gè)平面上的過(guò)程，也就是常說(shuō)的投影。天球坐標(biāo)系僅僅參考了地球大地坐標(biāo)系的經(jīng)緯度值，在天球坐標(biāo)系中，這兩個(gè)參數(shù)值被稱之為赤經(jīng)和赤緯。

在天球坐標(biāo)系中，天球的中心是地球的中心。因?yàn)橛^察者雖然并不在地球中心觀察，但是因?yàn)橛^察者所在的地球表面，其距離地球中心的距離，相比較天球無(wú)限大的距離來(lái)說(shuō)，可以忽略不記，因此可以簡(jiǎn)單認(rèn)為，觀察者所在的位置就是地球中心。

1.1.2 天球坐標(biāo)系

同樣一個(gè)天球，根據(jù)定義的經(jīng)過(guò)球心的水平面不同，可以區(qū)分為多種天球坐標(biāo)系，在此重點(diǎn)介紹常見(jiàn)的天球坐標(biāo)系有三種。

首先是經(jīng)過(guò)球心的平面為赤道面時(shí)，該處所指赤道便是地球赤道經(jīng)過(guò)往外無(wú)限延伸，最終勢(shì)必將會(huì)與天球相交，此時(shí)便形成了第二赤道坐標(biāo)系，此時(shí)的赤道面稱之為天赤道，經(jīng)地心，作一條與天赤道垂直的無(wú)限延長(zhǎng)的線，與天球相交于兩點(diǎn)，即圖 1中所示P點(diǎn)(天北極)與P′點(diǎn)(天南極)。

圖1 天球示意圖

然后是將地球與太陽(yáng)放在同一個(gè)平面上，所形成的平面，稱之為黃道面，以此為基準(zhǔn)，便是黃道坐標(biāo)系，僅靠地球與太陽(yáng)兩點(diǎn)不能確認(rèn)一個(gè)平面，還需要一到兩個(gè)點(diǎn)，這另外兩個(gè)點(diǎn)的選擇，便是春分點(diǎn)和秋分點(diǎn)，即太陽(yáng)直射點(diǎn)穿過(guò)赤道時(shí)的兩個(gè)點(diǎn)。經(jīng)過(guò)地心，作一條與黃道面垂直的直線，與天球相交于兩點(diǎn)，便是圖 1中的Q(黃北極)和Q′(黃南極)點(diǎn)。黃道面與赤道面相交所形成的夾角，便是黃赤交角。

最后還需要介紹的便是以觀察者所在位置，所作的與地球相交的平面，稱之為地平面，此時(shí)的地平面，在宇宙維度內(nèi)，可以認(rèn)為是經(jīng)過(guò)地球地心的一個(gè)平面，過(guò)地心作一條垂直于地平面的直線，與天球相交于兩點(diǎn)，即觀察者的頭頂與腳底，見(jiàn)圖1中Z(天頂)和Z′(天底)。

1.1.3 對(duì)天球上的常用點(diǎn)線面的理解

對(duì)圖 1中所涉及的標(biāo)注進(jìn)行解釋。

F點(diǎn)指的是飛行器，G點(diǎn)指的時(shí)太陽(yáng)，飛行器在此時(shí)，也可以當(dāng)成一個(gè)映射在天球球面上的物體，M點(diǎn)為地球上觀察者的位置，過(guò)天北極P點(diǎn)與天南極P′的赤經(jīng)圈稱之為天子午圈，天子午圈可以有很多條，但是通過(guò)天頂Z的，只有一個(gè)圈，稱之為本初子午圈，即圖 1中所畫弧線PZP′所在圓。

C點(diǎn)為春分點(diǎn)，需要重點(diǎn)介紹一下，在后續(xù)計(jì)算過(guò)程中，很多的運(yùn)算的數(shù)值，都與此概念相關(guān)。春分點(diǎn)時(shí)黃道面和赤道面的相交點(diǎn)，即太陽(yáng)直射點(diǎn)從南向北經(jīng)過(guò)赤道的那個(gè)點(diǎn)，這個(gè)點(diǎn)并不是一個(gè)實(shí)際上地球地理坐標(biāo)點(diǎn)，是無(wú)法用經(jīng)緯度表示的，春分點(diǎn)是兩個(gè)平面的交點(diǎn)，始終處在兩個(gè)平面的相交處。天文學(xué)上，將春分點(diǎn)作為很多天文概念的起算點(diǎn)，比如恒星時(shí)，比如太陽(yáng)赤經(jīng)。

太陽(yáng)赤緯角為地心與太陽(yáng)的連線與天赤道面的夾角，即圖 1中δ角。

太陽(yáng)赤經(jīng)角，從春分點(diǎn)C向東，一直到太陽(yáng)在天赤道上的投影線之間的夾角，即春分點(diǎn)所在的赤經(jīng)圈與太陽(yáng)所在的赤經(jīng)圈之間的夾角，稱之為赤經(jīng)角，即圖 1中的β角。

太陽(yáng)時(shí)角，是從本初子午圈與天赤道相交的點(diǎn)向西，即太陽(yáng)所在的赤經(jīng)圈與本初子午圈之間的夾角，即圖 1中的γ角。

最終需要求解的是太陽(yáng)夾角，即從地球上某個(gè)位置觀察飛行器和太陽(yáng)時(shí)，太陽(yáng)和飛行器之間的夾角，即圖 1中的α角。

在此需要特別注意的是，上述涉及的天文概念中，前綴往往都帶著參考的星球的名稱，比如太陽(yáng)，之所以這么解釋，是因?yàn)榭梢圆皇翘?yáng)，可以是宇宙中任何一個(gè)天體，比如銀河系的中心作為參考的銀河坐標(biāo)系。在本文中所涉及的未特別聲明的天文概念，均以太陽(yáng)作為參考。

1.2 太陽(yáng)赤緯的計(jì)算方法

要計(jì)算太陽(yáng)夾角，首先需要解決的問(wèn)題，便是如何準(zhǔn)確計(jì)算天球坐標(biāo)系下，太陽(yáng)赤緯及太陽(yáng)時(shí)角的值。首先計(jì)算太陽(yáng)赤緯的值。

計(jì)算太陽(yáng)赤緯值的方法有很多，常見(jiàn)的簡(jiǎn)單算法有Cooper算法、Spencer算法、Stine算法、Bourges算法、Wang算法、Yu算法，上述六種算法均基于天文觀測(cè)數(shù)據(jù)，通過(guò)公式擬合提出的數(shù)學(xué)公式，另外還可以通過(guò)對(duì)天文觀測(cè)數(shù)據(jù)的數(shù)值擬合算法，求出擬合度更高的計(jì)算公式，只是通過(guò)該方法提出的公式，其系數(shù)會(huì)因年而異，所以不具備普適性，另外，李文提出通過(guò)傅里葉展開(kāi)法提出對(duì)擬合公式的改進(jìn)算法。除數(shù)學(xué)擬合算法之外，還有一種基于VSOP87行星理論，演算推導(dǎo)算法的方法。

下面結(jié)合C++代碼語(yǔ)言，對(duì)太陽(yáng)赤緯各算法進(jìn)行逐個(gè)分析計(jì)算。

1.2.1 Bourges算子

Bourges是由Bourges提出的，Bourges算子的公式如下所示。其中ω為弧度參考系數(shù)，t為等效計(jì)日序數(shù)。太陽(yáng)赤緯角的單位為弧度。

δ=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)

弧度參考系數(shù)如下所示：

(2)

計(jì)日序數(shù)的計(jì)算如下所示。其中dn為計(jì)日序數(shù)，即本年度的第幾天。

t=dn-1-n0

(3)

n0的的計(jì)算公式如下，其中n表示年份，如本年度為2022年，則n=2022。

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

(4)

Bourges算法實(shí)現(xiàn)的代碼如下:

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;//等效計(jì)日序數(shù)

double w=2*M_PI/365.2422;//參考弧度系數(shù)

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為該年度中的天數(shù)累計(jì)值，如1月11日，則n=11。太陽(yáng)赤緯角的返回值單位是弧度。

(5)

δ為太陽(yáng)赤緯角。代碼實(shí)現(xiàn)如下：

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提出的，其中太陽(yáng)赤緯角的單位為弧度，θ為日角。算法公式如下所示：

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

日角的計(jì)算公式如下所示：

(6)

Spencer算子的實(shí)現(xiàn)代碼如下：

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提出來(lái)的，是在Spencer算子的基礎(chǔ)上做的簡(jiǎn)化，計(jì)算得到的太陽(yáng)赤緯角單位為弧度。θ為日角，公式如下:

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

(7)

其中的日角的計(jì)算公式如下:

(8)

實(shí)現(xiàn)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提出的，其中太陽(yáng)赤緯角的單位為弧度，其中n為計(jì)日序數(shù)，即本年度的第幾天。

Stine算法的實(shí)現(xiàn)代碼如下：

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的基礎(chǔ)上，針對(duì)計(jì)日系數(shù)，修改了部分系數(shù)，并將起算年份從1969改為了1985，使得其更加貼近真實(shí)值。計(jì)算所得的太陽(yáng)赤緯角的單位為弧度。公式如下所示，與Bourges算子相同。太陽(yáng)赤緯角的單位為弧度。

δ=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)

其中的弧度參考系數(shù)也與Bourges相同。如下所示:

(10)

以及等效計(jì)日序數(shù)的算法也是相同的，如下:

t=dn-1-n0

(11)

與Bourges不同的是對(duì)n0的計(jì)算，如下所示：

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

(12)

Wang算法的實(shí)現(xiàn)代碼如下所示:

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;//等效計(jì)日序數(shù)

double w=2*M_PI/365.2422;//參考弧度系數(shù)

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 數(shù)值擬合法

數(shù)值擬合法是由李文提出的，根據(jù)2015年-2018年的天文年歷，提出使用積日因子β作為多項(xiàng)式的中間系數(shù)，進(jìn)行擬合計(jì)算。得到適2015-2018年的年份的改進(jìn)公式。

太陽(yáng)赤緯角的單位為弧度，其中太陽(yáng)赤緯角的計(jì)算公式如下，其中的β為積日因子，ak為系數(shù)。

(13)

積日因子的求解公式如下所示。其中dn為積日系數(shù)。

(14)

n0的計(jì)算公式如下所示：

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

(15)

數(shù)值擬合算法的代碼如下所示:

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 傅里葉展開(kāi)法

傅里葉展開(kāi)法是由李文提出的，在數(shù)值擬合法的基礎(chǔ)上，進(jìn)行的改進(jìn)算法，通過(guò)觀察連續(xù)的太陽(yáng)赤緯角的變化規(guī)律，發(fā)現(xiàn)太陽(yáng)赤緯角的變化具備周期性，每4年一個(gè)循環(huán)，因此將算法的周期，從1年修改為4年，定義算法中的年份為每個(gè)周期中的第三年，并且將計(jì)日序數(shù)從之前的1年中的第幾天，改成整個(gè)4年周期中的第幾天，即總的計(jì)日序數(shù)。

得到的太陽(yáng)赤緯角的值，單位為弧度，其中β為積日因子。計(jì)算公式如下所示:

(16)

積日因子的計(jì)算公式如下所示，其中dn為積日系數(shù)。

(17)

式中,n0表達(dá)式如下所示:

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

(18)

傅里葉展開(kāi)法的計(jì)算公式如下所示:

double FourierBaseNF(cDateTime date)

{

int n=getDaysOfYear(date);

int year=date.year;

int y=0;//周期內(nèi)的第3個(gè)年份

int dn=0;//周期內(nèi)的總計(jì)日序數(shù)

if(year>=2015)

{

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

}else{

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

}

int temp=(year-3)%4;//這是周期內(nèi)的第幾年

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理論的進(jìn)一步完善，VSOP理論是用來(lái)描述太陽(yáng)系范圍內(nèi)的行星的軌道，在很長(zhǎng)時(shí)間內(nèi)的周期性變化的一種半分析理論。

VSOP形形理論除了可以算出精確的行星軌道參數(shù)，還可以直接計(jì)算出行星的位置，在此所言的位置，包含各種坐標(biāo)系在內(nèi)，如黃道坐標(biāo)系。

計(jì)算VSOP87計(jì)算太陽(yáng)赤緯角的過(guò)程，大致可分為七步。

1)計(jì)算儒略日數(shù)，用J表示。

2)計(jì)算相對(duì)儒略實(shí)際數(shù)，用T表示，如下所示:

(19)

3)計(jì)算太陽(yáng)幾何平黃經(jīng)，用L表示，如下所示:

L=280.466456+36000.76982779×T+0.003032028×

(20)

4)計(jì)算平黃赤交角，用MB表示，如下所示:

(21)

5)計(jì)算太陽(yáng)平近點(diǎn)角，用MSs表示，如下所示:

MSs=357.52191+35999.0503×T-

(22)

6)計(jì)算太陽(yáng)真黃經(jīng)，用SYS表示，如下所示:

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

(23)

7)計(jì)算太陽(yáng)地心赤緯角，用Dec表示，如下所示:

(24)

上式中涉及的radeg為弧度與度的轉(zhuǎn)換量綱，取radeg值為57.295 779 513 07。

用代碼實(shí)現(xiàn)上述公式如下:

double VSOP87(cDateTime date)

{

//儒略日數(shù)

double J=date2JD(date);

//計(jì)算相對(duì)儒略世紀(jì)數(shù)

double T=(J- 2451545.0)/36525;

//太陽(yáng)幾何平黃經(jīng)

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);

//太陽(yáng)平近點(diǎn)角

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

//太陽(yáng)真黃經(jīng)

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);

//太陽(yáng)地心赤緯

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

return Dec;

}

1.3 太陽(yáng)時(shí)角的計(jì)算方法

太陽(yáng)時(shí)角是以太陽(yáng)為參考，其值代表了觀察者所在的子午圈與太陽(yáng)所在的子午圈之間的夾角。太陽(yáng)時(shí)角的單位通?？梢允墙嵌戎?，如度分秒，此時(shí)的取值范圍為0至360度，也可以是時(shí)間值，如時(shí)分秒，此時(shí)的取值范圍為0至24小時(shí)。當(dāng)用角度表示時(shí)，其表示的含義可以理解為兩個(gè)子午圈之間的夾角，當(dāng)使用時(shí)分秒表示時(shí)，其含義可以理解為太陽(yáng)在時(shí)角時(shí)間之前通過(guò)觀察者正上空。如時(shí)角為3.5h，則代表在3.5小時(shí)之前，太陽(yáng)在3.5小時(shí)之前從觀察者所在的正上空通過(guò)。

太陽(yáng)時(shí)角的計(jì)算采用下列公式。其中tS為真太陽(yáng)時(shí)。

(25)

真太陽(yáng)時(shí)的計(jì)算采用下面中的公式計(jì)算，其中t為平太陽(yáng)時(shí)，單位為小時(shí)。在天文觀測(cè)中，很多天文概念都會(huì)加一個(gè)“平”字，其含義往往是代表著未加誤差修正的簡(jiǎn)單值，加上修正以后，通常稱之為“真”。eot即對(duì)時(shí)間的修正值。太陽(yáng)時(shí)角的計(jì)算，其重點(diǎn)就在于誤差修正值eot的計(jì)算。

(26)

下面將介紹6種計(jì)算太陽(yáng)時(shí)角的算子，分別為L(zhǎng)amm算子，Spencer算子，Whillier算子，Wloof算子，Yu算子，以及VSOP87行星理論。其區(qū)別主要便在誤差修正值的計(jì)算上。

1.3.1 Lamm算子

Lamm算子是由Lamm提出的，其特點(diǎn)是考慮了平年和閏年的區(qū)別，其公式如下。以四年為一個(gè)周期，其中n為全周期中的計(jì)日序數(shù)。Ak和Bk為系數(shù)值，在此采用杜春旭提出的參數(shù)值進(jìn)行計(jì)算。計(jì)算所得eot值單位為分鐘值，需要轉(zhuǎn)化為小時(shí)進(jìn)行下一步計(jì)算。

(27)

實(shí)現(xiàn)代碼如下所示:

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;

//轉(zhuǎn)為弧度。1個(gè)小時(shí)對(duì)應(yīng)15°。

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

return w;

}

1.3.2 Spencer算子

Spencer算子除了可以計(jì)算太陽(yáng)赤緯，也提出了計(jì)算太陽(yáng)時(shí)角的公式。其計(jì)算公式如下。其中θ為日角，單位為弧度。

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

(28)

日角的計(jì)算公式如下，周期為1年，n為計(jì)日序數(shù)。

(29)

公式計(jì)算所得的eot值單位為分鐘，需要轉(zhuǎn)換為小時(shí)再進(jìn)行下一步計(jì)算。

實(shí)現(xiàn)代碼如下所示:

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提出，其中θ為日角，單位為弧度。計(jì)算所得eot單位為分鐘，需要轉(zhuǎn)換成小時(shí)。

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

(30)

日角的計(jì)算公式如下，周期為1年，n為計(jì)日序數(shù)。

(31)

實(shí)現(xiàn)代碼如下:

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提出，其中θ為日角，單位為弧度。計(jì)算所得eot單位為分鐘，需要轉(zhuǎn)換成小時(shí)。誤差修正值計(jì)算公式如下:

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

(32)

其中的日角公式如下，周期為1年，n為計(jì)日序數(shù)。

(33)

代碼實(shí)現(xiàn)如下:

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提出，除了可計(jì)算赤緯，可提出了計(jì)算太陽(yáng)時(shí)角的公式，其中θ為日角，單位為弧度。計(jì)算所得eot單位為分鐘，需要轉(zhuǎn)換成小時(shí)。誤差修正值計(jì)算公式如下:

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

其中的日角公式如下，周期為1年，n為計(jì)日序數(shù)。

代碼實(shí)現(xià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行星理論

關(guān)于VSOP87行星理論上文已有介紹，利用該理論計(jì)算的過(guò)程如下。

1)計(jì)算相對(duì)儒略世紀(jì)數(shù)，用T表示:

(34)

2)計(jì)算黃道與月球平軌道升交點(diǎn)黃經(jīng)，用mw表示:

mw=125.04452-1934.136261×T+

(35)

3)計(jì)算太陽(yáng)幾何平黃經(jīng)，用L表示:

L=280.466456+36000.76982779×T+

(36)

4)計(jì)算月球幾何平黃經(jīng)，用L1表示:

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

(37)

5)計(jì)算黃經(jīng)章動(dòng)，用NHJ表示:

(38)

6)計(jì)算平黃赤交角，用MB表示:

(39)

7)計(jì)算任意時(shí)刻格林尼治恒星時(shí)，并加上黃經(jīng)章動(dòng)修正，得到視恒星時(shí)用Q0表示，需要將Q0值轉(zhuǎn)換到0～360°范圍內(nèi)：

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

(40)

8)計(jì)算太陽(yáng)平近點(diǎn)角，用MSs表示:

MSs=357.52191+35999.0503×T-

(41)

9)計(jì)算太陽(yáng)真黃經(jīng)，用SYS表示:

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

(42)

10)計(jì)算太陽(yáng)赤經(jīng)，用RM表示:

(43)

11)計(jì)算當(dāng)?shù)靥?yáng)時(shí)角，用ω，需要將ω值轉(zhuǎn)換到0～360°范圍內(nèi)，計(jì)算結(jié)果為度值，通常需要轉(zhuǎn)換成弧度值，方便后續(xù)計(jì)算:

(44)

ω=Q0+lon-120-RM

(45)

實(shí)現(xiàn)上述公式的代碼如下:

double VSOP87TimeAngle(cDateTime date,double lon)

{

//儒略日數(shù)

double J=date2JD(date);

//計(jì)算相對(duì)儒略世紀(jì)數(shù)

double T=(J- 2451545.0)/36525;

//黃道與月球平軌道升交點(diǎn)黃經(jīng)

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

//太陽(yáng)幾何平黃經(jīng)

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

//月球幾何平黃經(jīng)

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

//黃經(jīng)章動(dòng)

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);

//任意時(shí)刻格林尼治視恒星時(shí)，加黃經(jīng)章動(dòng)修正,得到視恒星時(shí)

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;

//太陽(yáng)平近點(diǎn)角

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

//太陽(yáng)真黃經(jīng)

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);

//太陽(yáng)地心赤經(jīng)

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

//當(dāng)?shù)貢r(shí)角

double w =Q0+lon-120-RM;

int temp=w/360;

w=w-temp*360;

return w/180*M_PI;

}

2 太陽(yáng)高度角及方位角的計(jì)算方法

太陽(yáng)高度角和方位角，即在觀察者位置觀察太陽(yáng)時(shí)的方位俯仰角，與當(dāng)?shù)氐乃矫嫦嚓P(guān)，太陽(yáng)高度角從地表橫切面起算，取值范圍為0～90°，水平角從正北方向起算，取值范圍為0～360°。

計(jì)算太陽(yáng)高度角時(shí)，需要考慮蒙氣差的影響。

2.1 太陽(yáng)高度角的計(jì)算

太陽(yáng)高度角的計(jì)算公式如下，其中ψ為觀察者所處地點(diǎn)的緯度值，δ為太陽(yáng)赤緯角，e為地球橢圓扁率：

(46)

其中地球扁率公式為：

(47)

太陽(yáng)高度角和方位角計(jì)算的代碼如下:

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);//計(jì)算蒙氣差，單位為角秒

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 太陽(yáng)方位角的計(jì)算

太陽(yáng)方位角的計(jì)算在計(jì)算得出高度角的基礎(chǔ)上進(jìn)行，公式如下，其中，α為太陽(yáng)高度角，ψ為觀察者所處緯度值，δ為太陽(yáng)赤緯角。

(48)

2.3 蒙氣差計(jì)算公式

蒙氣差修正，即對(duì)因大氣折射、溫度、氣壓帶來(lái)的誤差進(jìn)行修正。蒙氣差的單位為角秒。在此采用李文提出的計(jì)算公式，以公式適用范圍為依據(jù)，將蒙氣差計(jì)算按俯仰角進(jìn)行劃分。下列公式輸入值α為俯仰角，單位為弧度，輸出值為修正值，輸出為角秒。

2.3.1 蒙氣差對(duì)大氣折射的修正

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

在此范圍內(nèi)，使用傅里葉展開(kāi)法進(jìn)行計(jì)算修正，修正系數(shù)采用真實(shí)記錄值進(jìn)行計(jì)算。計(jì)算公式如下。

(49)

實(shí)現(xiàn)代碼如下：

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°范圍的修正

在此范圍內(nèi)，采用數(shù)值擬合法進(jìn)行修正，修正系數(shù)采用真實(shí)記錄值進(jìn)行計(jì)算，修正公式如下:

(50)

實(shí)現(xiàn)代碼如下:

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°范圍的修正

高俯仰范圍內(nèi)的修正方法較多，常見(jiàn)的有5種算法。

1)算法一的計(jì)算公式如下:

(51)

2)算法二的計(jì)算公式如下:

(51)

3)算法三的計(jì)算公式如下:

(52)

4)算法四的計(jì)算公式如下:

(53)

5)算法五的計(jì)算公式如下:

(54)

實(shí)現(xiàn)上述公式的代碼如下:

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 蒙氣差對(duì)溫度、氣壓的修正

蒙氣差還會(huì)受溫度和濕度的影響，因此需要對(duì)蒙氣差進(jìn)行附加值修正。輸入值為蒙氣差值，單位為角秒，輸出值為蒙氣差修正值，單位為角秒。蒙氣差修正公式如下。其中在高仰角(通常指大于45°)的情況下，大氣折射對(duì)溫度也會(huì)有一個(gè)影響，因此需要對(duì)溫度再進(jìn)行一次修正，其中μ為對(duì)溫度的附加修正。

R=R0(1+μA+B)

(55)

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

(56)

低仰角時(shí)，無(wú)需對(duì)溫度進(jìn)行修正。蒙氣差修正公式如下:

R=R0(1+A+B)

(57)

對(duì)溫度的修正公式如下:

(58)

對(duì)氣壓的修正公式如下:

(59)

上述公式實(shí)現(xiàn)代碼如下:

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 飛行器與太陽(yáng)夾角的計(jì)算方法

在計(jì)算出太陽(yáng)高度角和方位角以后，同時(shí)已知飛行器的方位俯仰角，此時(shí)便能計(jì)算出飛行器與太陽(yáng)之間的夾角。

計(jì)算夾角有兩種方法：一種是向量點(diǎn)乘的方法，另一種是三角函數(shù)推導(dǎo)的方法。

飛行器與太陽(yáng)及觀測(cè)點(diǎn)之間的相互關(guān)系如圖2。F點(diǎn)為飛行器，S點(diǎn)為太陽(yáng)，P點(diǎn)為觀察點(diǎn)，太陽(yáng)夾角即圖中所示∠FPS。

圖2 太陽(yáng)夾角示意圖

計(jì)算時(shí)，為方便計(jì)算，取模的長(zhǎng)度為1，α1為飛行器點(diǎn)的方位角，β1為飛行器的俯仰角，α2為太陽(yáng)的方位角，β2為太陽(yáng)的俯仰角。則F點(diǎn)方向可以假定一個(gè)模長(zhǎng)為1點(diǎn)F′，其坐標(biāo)可以表示為F′(cosβ1sinα1,cosβ1cosα1)，同理，在S點(diǎn)方向可以假定一個(gè)模長(zhǎng)為1的點(diǎn)S′，其坐標(biāo)可以表示為S′(cosβ2sinα2,cosβ2cosα2)。正東方向取為X軸正方向，正北方向取為Y軸正方向。

3.1 向量點(diǎn)乘法

向量點(diǎn)乘的計(jì)算公式如下:

(60)

根據(jù)F′與S′的值，代碼實(shí)現(xiàn)如下:

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 三角函數(shù)推導(dǎo)

在PF與PS兩條線段上，從P點(diǎn)出發(fā)，取一個(gè)長(zhǎng)度為1的線段，記作F′和S′，將F′和S′相連，成一個(gè)等腰三角形，計(jì)算F′S′長(zhǎng)度，從而求出∠F′PS′，即太陽(yáng)夾角。

經(jīng)推導(dǎo)，計(jì)算公式結(jié)果如下:

(61)

代碼實(shí)現(xiàn)如下:

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 飛行器與太陽(yáng)夾角的精度分析

為比較各算法的計(jì)算結(jié)果，進(jìn)行精度分析。

4.1 精度分析方法

用以進(jìn)行數(shù)值精度鑒定的參考數(shù)據(jù)，通常有兩種來(lái)源，一是通過(guò)實(shí)際觀察獲取的實(shí)測(cè)記錄值，另一種來(lái)源是精度更高的計(jì)算值。在本次計(jì)算赤經(jīng)赤緯數(shù)值的精度鑒定中，采用第二種方法，以精度更高的計(jì)算值來(lái)作為參考值的辦法。采用復(fù)雜算法SPA算法的計(jì)算值作為數(shù)據(jù)的參考值，SPA算法是綜合考慮了海拔、壓強(qiáng)、溫度、傾斜度、方位角旋轉(zhuǎn)、大氣折射、時(shí)區(qū)、地球自轉(zhuǎn)時(shí)間與地球時(shí)間之差等各方面因素的算法，其計(jì)算精度可達(dá)±0.000 3°，可以用來(lái)作為參考值對(duì)簡(jiǎn)單算法進(jìn)行精度分析。

因此在對(duì)太陽(yáng)夾角的數(shù)值分析時(shí)，選用SPA算法作為參考值。誤差分析選用均值、方差、均方根三項(xiàng)指標(biāo)進(jìn)行分析。

飛行器飛行軌跡選用樣本，按照每秒一個(gè)點(diǎn)，每秒在方位、俯仰方向上各新增0.000 1°，全時(shí)長(zhǎng)為300秒設(shè)計(jì)。飛行時(shí)間設(shè)定為2018年8月8日上午10點(diǎn)整。起飛地點(diǎn)設(shè)置為東經(jīng)102.241 897 39，北緯27.902 341 42。

數(shù)據(jù)樣本為全飛行軌跡中，每個(gè)點(diǎn)的太陽(yáng)夾角與參考值的差值。

均值即將樣本值求均值。公式如下:

(62)

方差為將樣本值與樣本的均值相減，分別求平方后求和，再將平方和求均值。

(63)

均方根，將所有樣本值分別求平方，再求平方和以后求平方和的均值，最后開(kāi)平方。

(64)

分析對(duì)象為：9種太陽(yáng)赤緯角算法，6種太陽(yáng)時(shí)間算法。

4.2 精度分析結(jié)果

4.2.1 均值分析結(jié)果

均值分析結(jié)果如表1所示。

表1 均值計(jì)算結(jié)果

圖3為折線圖形式。橫坐標(biāo)為赤緯算子，每一條折現(xiàn)代表一種時(shí)角算法。

圖3 均值折線圖

均值能反映統(tǒng)計(jì)樣本的大致的平均水準(zhǔn)，均值越小，說(shuō)明樣本值整體越小。從表1均值計(jì)算結(jié)果分析可得以下結(jié)論。

1)不同的赤緯角算子下，VSOP87的時(shí)角計(jì)算結(jié)果偏差最小。

2)不同的時(shí)角算子下，Wang算法的赤緯角計(jì)算結(jié)果偏差最小。

4.2.2 方差分析結(jié)果

方差分析結(jié)果如表2所示。

表2 方差計(jì)算結(jié)果

圖4為折線圖形式。橫坐標(biāo)為赤緯算子，每一條折現(xiàn)代表一種時(shí)角算法。

圖4 方差折線圖

方差分析反映的是自變量對(duì)數(shù)值的影響。即對(duì)于均值的偏離度的一個(gè)分析，即數(shù)值的穩(wěn)定程度。數(shù)值越小，說(shuō)明數(shù)值越穩(wěn)定，起伏越小。

1)在不同的赤緯角算子，VSOP87的時(shí)角計(jì)算結(jié)果波動(dòng)最小。

2)在不同的時(shí)角算子下，Wang算子的赤緯角計(jì)算結(jié)果波動(dòng)最小。

4.2.3 均方根分析結(jié)果

均方根分析結(jié)果如表3所示。

表3 均方根計(jì)算結(jié)果

圖5為折線圖形式。橫坐標(biāo)為赤緯算子，每一條折現(xiàn)代表一種時(shí)角算法。

圖5 均方根折線圖

均方根和均值一起分析，反映了樣本的一致性情況。均方根和均值越接近，說(shuō)明樣本數(shù)據(jù)一致性越好。

1)在不同的赤緯角算子下，VSOP87的時(shí)角一致性最好。

2)在不同的時(shí)角算子下，Wang算子的赤緯角一致性最好。

5 結(jié)束語(yǔ)

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