許偉維 廖新浩 周永宏 許雪晴

(1中國(guó)科學(xué)院上海天文臺(tái) 上海 200030) (2上?？萍即髮W(xué)物質(zhì)學(xué)院 上海 201210) (3中國(guó)科學(xué)院行星科學(xué)重點(diǎn)實(shí)驗(yàn)室 上海 200030) (4中國(guó)科學(xué)院大學(xué) 北京 100049)

天體測(cè)量法探測(cè)系外行星?

許偉維1,2,3,4?廖新浩1,3?周永宏1,3許雪晴1,3

(1中國(guó)科學(xué)院上海天文臺(tái) 上海 200030) (2上海科技大學(xué)物質(zhì)學(xué)院 上海 201210) (3中國(guó)科學(xué)院行星科學(xué)重點(diǎn)實(shí)驗(yàn)室 上海 200030) (4中國(guó)科學(xué)院大學(xué) 北京 100049)

在目前已發(fā)現(xiàn)的系外行星中,絕大多數(shù)是由視向速度法和凌星法探測(cè)得到的,天體測(cè)量法僅發(fā)現(xiàn)了1顆.gaia衛(wèi)星數(shù)據(jù)即將發(fā)布,天體測(cè)量法將逐步成為系外行星探測(cè)的重要方法之一.基于天體測(cè)量法給出的恒星位置參數(shù)序列,討論了在求解行星質(zhì)量和軌道參數(shù)時(shí)涉及的動(dòng)力學(xué)條件方程計(jì)算問題,給出了具體微分改正公式,同時(shí)也進(jìn)行了必要的仿真模擬計(jì)算.建立的方法可以較容易地推廣到多行星系統(tǒng).

系外行星,天體測(cè)量學(xué),天體力學(xué):軌道計(jì)算與定軌,方法:數(shù)據(jù)分析

1 引言

系外行星探測(cè)和研究是當(dāng)今國(guó)際天文學(xué)研究的熱點(diǎn).1995年,M ayor和Queloz采用視向速度法在飛馬座51附近發(fā)現(xiàn)了第1顆圍繞著主序恒星運(yùn)動(dòng)的木星質(zhì)量量級(jí)的系外行星飛馬座51 b[1].近年來,隨著天文探測(cè)技術(shù)的進(jìn)步,特別是Kep ler空間望遠(yuǎn)鏡的成功發(fā)射和觀測(cè),新的系外行星不斷被發(fā)現(xiàn),其數(shù)量明顯增多.根據(jù)系外行星網(wǎng)站(exop lanets.eu)統(tǒng)計(jì),截止2015年12月1日,已確認(rèn)發(fā)現(xiàn)2004顆系外行星,隸屬于1269個(gè)行星系統(tǒng),包括498個(gè)多行星系統(tǒng);此外,400多顆系外行星候選體需進(jìn)一步認(rèn)證.在這些已發(fā)現(xiàn)的系外行星中,絕大多數(shù)為氣態(tài)巨行星和熱木星(約占70%),部分為超級(jí)地球(20%)和地球質(zhì)量大小的行星(2%),但尚未確認(rèn)存在生命的類地行星系統(tǒng).

系外行星的探測(cè)方法分為間接探測(cè)和直接成像,間接探測(cè)有視向速度法、天體測(cè)量法、脈沖星計(jì)時(shí)法、凌星法、微引力透鏡法等[2],目前實(shí)際觀測(cè)中采用最多的是視向速度法和凌星法.天體測(cè)量法至今只發(fā)現(xiàn)了一顆系外行星—HD 176051 b[3].2009年以前,由于HARPS(High Accuracy Radial velocity Planet Searcher)的應(yīng)用[4],視向速度法為主要的探測(cè)方法;2009年美國(guó)國(guó)家航空航天局(National Aeronautics and Space Adm inistration,NASA)建造的開普勒太空望遠(yuǎn)鏡發(fā)射升空,凌星法逐漸成為主流方法,目前該方法已發(fā)現(xiàn)超過1200顆系外行星(exop lanets.eu).2013年12月19日,歐洲空間局(European Space Agency,ESA)研制的gaia空間望遠(yuǎn)鏡升空,其目的是以前所未有的精度對(duì)銀河系內(nèi)數(shù)以十億計(jì)的恒星進(jìn)行觀測(cè),測(cè)量他們的位置、距離和運(yùn)動(dòng),其觀測(cè)效率將比同樣由ESA發(fā)射的依巴谷衛(wèi)星高出數(shù)百萬倍[5].因此,人們對(duì)于利用gaia開展太陽系與系外行星的探測(cè)寄予厚望.

天體測(cè)量法是通過測(cè)量恒星空間位置變化來反演其附近行星的質(zhì)量和軌道參數(shù)的,相比于目前通常采用的視向速度法和凌星法,天體測(cè)量法有其獨(dú)特的優(yōu)勢(shì):(1)直接確定行星質(zhì)量.視向速度法也是通過測(cè)量恒星空間位置變化來反演其附近行星的軌道參數(shù)的,但因其觀測(cè)方程中僅包含了行星質(zhì)量m與行星軌道傾角的正弦sin i的乘積項(xiàng),所以只能推算出m sin i,不能直接確定行星的質(zhì)量.天體測(cè)量法對(duì)應(yīng)的觀測(cè)方程不僅包含了m cos i項(xiàng),而且包含m項(xiàng),故可以通過觀測(cè)資料同時(shí)直接推算出行星的質(zhì)量與軌道傾角[6].(2)尋找長(zhǎng)周期行星.根據(jù)開普勒第三定律可知,行星繞恒星運(yùn)動(dòng)的軌道半長(zhǎng)徑3次冪與行星軌道周期的平方成正比;同樣質(zhì)量的行星,行星軌道半長(zhǎng)徑越長(zhǎng),對(duì)于恒星繞系統(tǒng)質(zhì)心運(yùn)動(dòng)的擾動(dòng)就越大,因此越容易被天體測(cè)量法分辨出來.所以天體測(cè)量法對(duì)長(zhǎng)周期行星探測(cè)較為敏感.(3)檢測(cè)多星系統(tǒng)是否共面.在已知的系外行星中,大多數(shù)具有較大的軌道偏心率,其原因可能是由于受到一顆或多顆大質(zhì)量天體的引力作用[7],該作用力使得行星的軌道面不與恒星的赤道面共面,若在多星系統(tǒng)中,該作用力將使得各行星軌道不共面[8].研究行星質(zhì)量和相互夾角對(duì)于多星系統(tǒng)的動(dòng)力學(xué)模型、動(dòng)力學(xué)穩(wěn)定以及軌道共振具有重要意義.(4)太陽系附近尋找類地行星.利用空間干涉的天體測(cè)量能達(dá)到亞微角秒級(jí)的精度,對(duì)于造成恒星擾動(dòng)的行星質(zhì)量檢測(cè)可達(dá)到m⊕(地球質(zhì)量)量級(jí),滿足在太陽系附近尋找類地行星的要求[9].

對(duì)于由1顆恒星和1顆行星組成的二體系外行星系統(tǒng),W right和Howard于2009年建立了基于軌道根數(shù)的天體測(cè)量法對(duì)應(yīng)的觀測(cè)方程[10],其在赤道坐標(biāo)系下,包含了視差、自行以及行星質(zhì)量和軌道根數(shù)等參數(shù).通過分析二體運(yùn)動(dòng)的解和參數(shù)估計(jì)理論,可以給出這些參數(shù)的最優(yōu)解.稍加分析即可發(fā)現(xiàn),基于軌道根數(shù)的天體測(cè)量方法,當(dāng)行星軌道偏心率e=0時(shí),因狀態(tài)轉(zhuǎn)移矩陣奇異而失效,特別是對(duì)于多行星系統(tǒng),基于軌道根數(shù)的狀態(tài)方程和狀態(tài)轉(zhuǎn)移矩陣計(jì)算都非常復(fù)雜,不便于實(shí)際應(yīng)用,因此需要建立簡(jiǎn)潔可行的確定系外行星軌道的天體測(cè)量方法.本文針對(duì)多體系統(tǒng)動(dòng)力學(xué)特性,提出采用坐標(biāo)速度描述,不僅可以避免e=0對(duì)應(yīng)的狀態(tài)轉(zhuǎn)移矩陣奇異問題,而且可以顯著降低整個(gè)計(jì)算的復(fù)雜性.本文第2節(jié)對(duì)gaia測(cè)量數(shù)據(jù)用于系外行星的確定進(jìn)行了必要的精度分析;第3節(jié)較全面地總結(jié)了基于軌道根數(shù)的天體測(cè)量方法涉及的理論基礎(chǔ);第4節(jié)給出了基于坐標(biāo)速度的天體測(cè)量方法對(duì)應(yīng)的觀測(cè)方程、運(yùn)動(dòng)方程和狀態(tài)轉(zhuǎn)移方程;第5節(jié)列出了仿真模擬計(jì)算結(jié)果并作了一些分析;第6節(jié)推導(dǎo)了行星軌道確定時(shí)初始估計(jì)滿足的約束條件,由此可以顯著提高軌道確定微分改正的收斂效率;第7節(jié)推導(dǎo)了三體問題(1顆恒星2顆行星)的計(jì)算公式;最后1節(jié)給出本文的結(jié)論,同時(shí)也討論存在的問題及其解決途徑.

2 天體測(cè)量觀測(cè)精度分析

當(dāng)恒星附近存在行星時(shí),由于引力作用,使得恒星與行星繞著它們的質(zhì)心以相同的角速度作橢圓運(yùn)動(dòng).假設(shè)恒星和行星圍繞系統(tǒng)質(zhì)心運(yùn)動(dòng)軌道半長(zhǎng)徑分別為as和ap,則在地面上觀察恒星軌道半長(zhǎng)徑時(shí),可產(chǎn)生一個(gè)角度2θ[11](如圖1).

其中,Ms和mp分別為恒星和行星的質(zhì)量,G為萬有引力常數(shù),D為恒星與觀測(cè)者之間的距離,p為二體運(yùn)動(dòng)的軌道周期.角半徑θ也稱作恒星的天體測(cè)量特征角,表征了行星對(duì)恒星的擾動(dòng)大小,除行星質(zhì)量外,還與恒星的質(zhì)量、軌道周期及觀測(cè)距離有關(guān).根據(jù)目前已發(fā)現(xiàn)的系外行星,表1列出了它們的質(zhì)量與軌道周期分布范圍,以及所繞恒星的質(zhì)量和與觀測(cè)距離的分布范圍(參見網(wǎng)站exop lanets.eu),其中行星質(zhì)量單位為MJ(木星質(zhì)量),恒星質(zhì)量單位為M⊙(太陽質(zhì)量).針對(duì)上述參數(shù)分布范圍,本文選取了5種模型,模型中恒星質(zhì)量皆取為1 M⊙,行星質(zhì)量、軌道周期以及觀測(cè)距離不等,計(jì)算5種模型中恒星的天體測(cè)量特征角,見表2.

圖1 二體運(yùn)動(dòng)的軌道示意圖Fig.1 The orbita l sketch of the tw o-body m otion

表1 已知系外行星的部分參數(shù)范圍Tab le 1 The partial param eters’range of the know n exop lanets

表2 幾種模型下恒星天體測(cè)量特征角Tab le 2 The astrom etric signatu re of the star in d ifferen t m odels

由于gaia對(duì)于不同星等的恒星觀測(cè)精度不同,所以還需分析gaia對(duì)于系外行星的可觀測(cè)性.表3給出了對(duì)于太陽質(zhì)量級(jí)的恒星,當(dāng)其觀測(cè)距離在1000 pc內(nèi)時(shí)的星等及gaia對(duì)應(yīng)的位置測(cè)量精度,由表可知:在1000 pc內(nèi),恒星的視星等mV=14.75mag,gaia的觀測(cè)精度高于10μas,其值均小于表2中各模型的天體測(cè)量特征角,因此gaia可進(jìn)行有效觀測(cè).

表3 太陽質(zhì)量級(jí)恒星在不同距離的視星等及gaia對(duì)其的觀測(cè)精度Tab le 3 The ap paren t m agn itude of the sun like stars at d ifferen t d istances and the accu racies of the gaia

3 基于軌道根數(shù)的觀測(cè)方程

利用天體測(cè)量法間接探測(cè)系外行星時(shí),通過觀測(cè)恒星的周期性位置變化來判斷其是否存在.對(duì)于1顆孤立的恒星,可通過5個(gè)天體測(cè)量觀測(cè)值來描述其在天球上的角位置,分別是赤道坐標(biāo)系下的赤經(jīng)α0、赤緯δ0(在某一特定的坐標(biāo)系下某一歷元時(shí)的值,現(xiàn)在一般都取在國(guó)際天球參考系下J2000.0的值),恒星自行相應(yīng)的正交分量μα、μδ,以及視差?;繞著恒星作橢圓運(yùn)動(dòng)的行星軌道可由6個(gè)軌道根數(shù)描述,即軌道半長(zhǎng)徑a、傾角i、升交點(diǎn)經(jīng)度?、近星點(diǎn)經(jīng)度ω、偏心率e及平近點(diǎn)角M;行星質(zhì)量m可由p和a通過Kep ler第三定律確定.因此,對(duì)一個(gè)含有n顆行星的系統(tǒng)進(jìn)行軌道擬合時(shí),需要5+7n個(gè)參數(shù).

建立如圖2所示觀測(cè)示意圖,恒星的運(yùn)動(dòng)軌道處于O-X YZ坐標(biāo)系內(nèi),O-xyz系為觀測(cè)面所在坐標(biāo)系,z軸表示觀測(cè)者視線方向,x軸指向北極,y方向由右手法則確定,i、?、ω為3個(gè)方位角.

圖2 二體運(yùn)動(dòng)的觀測(cè)示意圖Fig.2 T he sketch of the observation of the tw o-body m otion

恒星在空間位置的變化反映了其橢圓軌道運(yùn)動(dòng)在觀測(cè)面上的投影,由赤經(jīng)和赤緯兩個(gè)位置參數(shù)決定,分別用?η(t)、?ξ(t)表示,對(duì)應(yīng)的觀測(cè)方程為

等式左邊為觀測(cè)值,右邊各項(xiàng)依次為初始值、視差改正、自行改正以及天體擾動(dòng)項(xiàng),其中,?α、?δ表示視差改正參數(shù);?xs、?ys表示恒星運(yùn)動(dòng)在觀測(cè)面上兩個(gè)方向的投影,單位為角度;為了計(jì)算和公式表達(dá)的方便,選取適當(dāng)?shù)挠?jì)算單位是必要的,本文分別采取太陽質(zhì)量M⊙、太陽赤道半徑A⊙作為基本質(zhì)量和長(zhǎng)度單位,取引力常數(shù)G=1,則相應(yīng)的時(shí)間單位為

故可通過觀測(cè)距離D將?xs、?ys歸算到系統(tǒng)長(zhǎng)度單位下,分別用xs、ys表示

考慮恒星附近僅存在1顆行星情況,即n=1,由天體力學(xué)二體問題分析解可將xs、ys表示為

其中,

A、B、F、G稱為Thiele-Innes常數(shù)[13],可由軌道根數(shù)描述,

恒星視差與自行改正可通過一定量天體測(cè)量數(shù)據(jù)計(jì)算出來,故在后面的計(jì)算中將視差與自行視為已知量.至此,我們可將觀測(cè)方程改寫為

這里,η′(t)和ξ′(t)包含了自行和視差改正項(xiàng),已歸算為系統(tǒng)長(zhǎng)度單位.假設(shè)初始時(shí)刻恒星軌道根數(shù)和行星質(zhì)量的估計(jì)值為a?、i?、??、ω?、e?、M?及m?,為方便起見,用矢量將其表示為σ=(σ1,σ2,σ3,σ4,σ5,σ6,σ7).對(duì)上式在恒星初始運(yùn)動(dòng)位置處展開,有:

其中,Pi(t)、Qi(t)依次表示對(duì)7個(gè)參數(shù)的微分結(jié)果.如果有k次觀測(cè),那么2k個(gè)觀測(cè)方程可以寫成矩陣形式

其中,

上述方程(10)包含7個(gè)未知量,因此問題可解至少需要觀測(cè)4次.

將由估計(jì)理論首次解得的?σ(0)疊加到初始時(shí)刻近似軌道上,進(jìn)行迭代計(jì)算,當(dāng)?shù)趌次迭代給出的|?σ(l)|小于設(shè)定的收斂標(biāo)準(zhǔn)ε時(shí),求解完成.本文統(tǒng)一設(shè)定誤差限ε= 10?6.

由上述獲得的恒星軌道根數(shù)與行星質(zhì)量的估計(jì)值,可以求出恒星初始時(shí)刻的坐標(biāo)→rs與速度˙→rs,由此根據(jù)二體關(guān)系式,可以求出行星此時(shí)對(duì)應(yīng)的坐標(biāo)→rp與速度˙→rp,

對(duì)于多行星系統(tǒng),即n≥2,對(duì)應(yīng)的運(yùn)動(dòng)沒有分析解,故上述方法不再適用.

4 基于坐標(biāo)速度的觀測(cè)方程

觀測(cè)方程與軌道根數(shù)法對(duì)應(yīng)的恒星觀測(cè)方程相同,即(2)式.假定t=t0時(shí)刻恒星位置與速度和行星質(zhì)量的估計(jì)值分別為(t0)、(t0)、mp,其近似值(即初始估計(jì))為、、m?,兩者相差為?、?、?m,將觀測(cè)方程在、、m?處展開有

當(dāng)t=t0時(shí),狀態(tài)轉(zhuǎn)移方程涉及的初始條件為

這里,H為3階單位矩陣.

5 軌道特征分析與仿真計(jì)算

目前,gaia的巡天數(shù)據(jù)仍未釋放,因此為了驗(yàn)證上述方法的有效性,只能進(jìn)行仿真模擬計(jì)算.假定t0時(shí)刻1組恒星的軌道根數(shù)和行星質(zhì)量σ(t0),由該組參數(shù)可模擬出1組恒星位置的觀測(cè)數(shù)據(jù),然后用其反解σ(t0),由此可以驗(yàn)證和分析評(píng)估軌道確定的方法和精度.

5.1 軌道特征的影響

5.1.1 小偏心率

對(duì)于小偏心率軌道,軌道根數(shù)法不能進(jìn)行軌道確定,但坐標(biāo)速度法仍有效,其原因是e=0時(shí)近星點(diǎn)經(jīng)度不能確定所致.具體體現(xiàn)在當(dāng)e=0時(shí),xs、ys可表示為:

由上式可以看出

即狀態(tài)轉(zhuǎn)移矩陣中有關(guān)ω和M的兩列相關(guān),導(dǎo)致狀態(tài)轉(zhuǎn)移矩陣奇異,故對(duì)于小偏心率的情況,軌道根數(shù)法失效.

5.1.2 軌道傾角

對(duì)于傾角i=0°或180°的軌道,兩種方法都不適用.

軌道根數(shù)法中,由于xs、ys中包含的i項(xiàng)都以cos i出現(xiàn),故當(dāng)i=0°或180°時(shí),xs、ys對(duì)i的偏導(dǎo)數(shù)都為0,具體體現(xiàn)在

坐標(biāo)速度法中,當(dāng)i=0°時(shí),中不顯含z分量,即z?和分量.故有

兩種情況都導(dǎo)致線性方程組的系數(shù)矩陣奇異,無唯一解.

此外,天體測(cè)量法描述的是恒星X、Y方向的運(yùn)動(dòng)在觀測(cè)面上的投影,對(duì)于Z方向沒有直接約束,這就造成了當(dāng)i?π?i時(shí),σ(t0)一樣能滿足觀測(cè)方程.但此時(shí)除傾角項(xiàng)外,升交點(diǎn)的經(jīng)度?和近地點(diǎn)的經(jīng)度ω也存在相應(yīng)的變化,這一點(diǎn)我們可從Thiele-Innes常數(shù)來理解.我們將其對(duì)應(yīng)關(guān)系整理如下:

5.2 觀測(cè)誤差的影響

觀測(cè)即有誤差,因此不同的仿真模型應(yīng)增加相應(yīng)的觀測(cè)誤差.以第2節(jié)中模型A為例,mV=5.63mag＜10mag,此時(shí)gaia的觀測(cè)精度?θ=4μas,若將其視為天體測(cè)量觀測(cè)誤差,對(duì)于觀測(cè)距離D=15 pc,該誤差s的大小為

將其歸算到系統(tǒng)長(zhǎng)度單位A⊙,即得該模型的觀測(cè)誤差?s為

通過M atlab軟件可計(jì)算出任意兩組滿足均值為0、方差為0.013的高斯白噪聲,將其分別增加到恒星的x和y方向的觀測(cè)數(shù)據(jù)上,則可獲得含有觀測(cè)誤差的模擬觀測(cè)數(shù)據(jù).

由含有誤差觀測(cè)數(shù)據(jù)可確定初始根數(shù)σ(t0)的估計(jì)值ˉσ(t0),兩者之差為?σ.表4列出了幾種模型對(duì)于不同觀測(cè)誤差和不同觀測(cè)次數(shù)時(shí)確定的初始軌道估計(jì)值誤差情況.模型A表示一個(gè)軌道周期為2 yr的木星繞一個(gè)距離地球15 pc的太陽作橢圓運(yùn)動(dòng)(見表2),此太陽在天球上位置變化最大角度(特征角)為100μas,考慮不同觀測(cè)誤差與觀測(cè)次數(shù)的組合,形成了3組模擬觀測(cè)數(shù)據(jù):A-1、A-2和A-3.模型A的初始軌道根數(shù)取為:a=1.587 au,e=0.4000,i=40.00°,?=80.00°,ω=50.00°,M=50.00°, m=9.552×10?4Ms;軌道確定誤差統(tǒng)計(jì)分為x、y方向的殘差及軌道根數(shù)的相對(duì)誤差.以A模型為例,對(duì)于特征角為100μas的系外行星系統(tǒng),考慮觀測(cè)隨機(jī)誤差分別為2μas和4μas以及觀測(cè)點(diǎn)數(shù)分別為20和40的情況.模型A-1對(duì)應(yīng)觀測(cè)誤差2μas,觀測(cè)點(diǎn)數(shù)20,由坐標(biāo)速度法確定的軌道殘差的x和y分量分別為23.75μas和24.20μas,由軌道根數(shù)法確定的軌道半長(zhǎng)徑和軌道偏心率的相對(duì)誤差分別為1.65%和0.92%,其他根數(shù)的誤差均小于10%,質(zhì)量確定的精度較高,可好于千分之一;模型A-2對(duì)應(yīng)觀測(cè)誤差2μas,觀測(cè)點(diǎn)數(shù)40,相較于A-1各根數(shù)的相對(duì)誤差均顯著減小,即擬合結(jié)果更優(yōu);模型A-3對(duì)應(yīng)觀測(cè)誤差4μas,觀測(cè)點(diǎn)數(shù)40,相較于A-2大部分根數(shù)的相對(duì)誤差有所增大.通過比較不同模型對(duì)應(yīng)的計(jì)算結(jié)果表明,觀測(cè)數(shù)據(jù)越多,相對(duì)誤差越小,軌道確定的擬合結(jié)果越可靠.

表4 模型A-1、A-2和A-3的計(jì)算結(jié)果Tab le 4 T he ca lcu lation s of the m od els A-1,A-2,an d A-3

5.3 初值的選取

由于觀測(cè)方程(2)中關(guān)于恒星初始軌道根數(shù)或恒星初始坐標(biāo)速度函數(shù)的表示均是非線性的,故在軌道確定中需要采用微分迭代方法求解軌道估計(jì)值,因此需要選定恰當(dāng)?shù)牡踔?以使求解過程收斂.對(duì)于軌道初值的選取,最簡(jiǎn)單的方法是在一定區(qū)間內(nèi)對(duì)每個(gè)待解參數(shù)進(jìn)行掃描計(jì)算,當(dāng)?shù)踔当贿x取在參數(shù)收斂區(qū)間時(shí),便可給出估計(jì)值.本文討論的問題涉及的參數(shù)收斂區(qū)間是7維的,因此要獲得收斂的估計(jì)值,需要較大的掃描計(jì)算量.為了方便計(jì)算,采用固定?、ω、M、m進(jìn)行a、i、e 3維檢索,再固定a、i、e、m進(jìn)行?、ω、M 3維檢索.計(jì)算結(jié)果列于表5,表中列出的是相對(duì)誤差百分比值,需要說明的是其中角度量的誤差相對(duì)量取180,偏心率的誤差相對(duì)量取為1,主要是為了避免小偏心率作為相對(duì)量導(dǎo)致相對(duì)誤差失真.模型A-1和A-3的軌道半長(zhǎng)徑收斂區(qū)間約為5%,其他根數(shù)的收斂區(qū)間相對(duì)較大,一般大于10%,有些甚至可達(dá)20%.從初步的掃描計(jì)算結(jié)果看,軌道半長(zhǎng)徑a的迭代初值選取是影響計(jì)算收斂與否的關(guān)鍵因素,需要盡量靠近估計(jì)值.

表5 模型A-1和A-3的收斂區(qū)間Tab le 5 The convergen ce in terval of the m od els A-1 and A-3

6 關(guān)于初始估計(jì)取值的約束

在上節(jié)討論中,為減少計(jì)算時(shí)間,對(duì)7維參數(shù)空間僅進(jìn)行了兩個(gè)3維參數(shù)空間掃描,且行星質(zhì)量均是固定的,沒有參與掃描計(jì)算.事實(shí)上,可通過恒星空間位置觀測(cè)資料時(shí)間序列對(duì)上述7維參數(shù)空間進(jìn)行適當(dāng)?shù)募s束,降低參數(shù)空間的維數(shù).具體做法如下:已知由天體測(cè)量方法給出的某顆恒星的位置時(shí)間序列xs(tj),ys(tj),j=1,2,···,k,其包含了運(yùn)動(dòng)周期信息,通過對(duì)其傅里葉譜分析(FFT)可獲得軌道周期p,再根據(jù)開普勒第三定律可以給出恒星軌道半長(zhǎng)徑a和行星質(zhì)量mp之間的關(guān)系式,

一般認(rèn)為Ms的值可由天體物理方法獲得.

傅里葉譜分析通常針對(duì)等間距數(shù)據(jù)(均勻采樣)進(jìn)行,故需要先對(duì)恒星的觀測(cè)數(shù)據(jù)進(jìn)行插值,再利用插值數(shù)據(jù)進(jìn)行譜分析.此外,還采用了另一種方法——LSQ(最小二乘擬合),將插值后的數(shù)據(jù)進(jìn)行最小二乘,擬合成一條正弦曲線,該正弦曲線的周期即為恒星運(yùn)動(dòng)周期.

已知gaia的觀測(cè)年限為5 yr,故進(jìn)行傅里葉譜分析所選取的模型周期不應(yīng)該大于5 yr,否則無法觀測(cè)到一個(gè)完整周期的數(shù)據(jù);又根據(jù)奈奎斯特取樣定理,所選模型周期應(yīng)限制在2.5 yr內(nèi).此處,我們選取5.2節(jié)中所建立的模型進(jìn)行計(jì)算,圖3給出了A-3模型(任意選取)x、y方向誤差數(shù)據(jù)的分布及進(jìn)行傅里葉分析的結(jié)果,并將計(jì)算結(jié)果列于表6中.

表6 模型A-1、A-2、A-3進(jìn)行最小二乘擬合和傅里葉分析的結(jié)果Tab le 6 T he LSQ an d FFT resu lts of the m od els A-1,A-2,and A-3

圖3 模型A-3的x、y方向數(shù)據(jù)進(jìn)行傅里葉分析的結(jié)果Fig.3 T he resu lts of x and y com p onen ts of the m odel A-3 with FFT

根據(jù)統(tǒng)計(jì),LSQ和FFT均能實(shí)現(xiàn)觀測(cè)數(shù)據(jù)的周期分析.模型A的精確周期為730 d,3個(gè)子模型進(jìn)行FFT的結(jié)果分別為720 d、716.67 d和716.67 d,明顯要好于LSQ的結(jié)果.

對(duì)序列{xs(tj),ys(tj),j=1,2,···,k}在時(shí)間域上作平均,可以獲得如下關(guān)系式,假定已知行星質(zhì)量mp、恒星的平近點(diǎn)角M、恒星的升交點(diǎn)經(jīng)度?和恒星的近星點(diǎn)經(jīng)度ω,可以由上述3個(gè)約束關(guān)系式(21)–(23)式,近似給出恒星的軌道半長(zhǎng)徑、偏心率和軌道傾角,這樣將對(duì)7維參數(shù)空間掃描問題降低到4維參數(shù)空間掃描,可以顯著提高計(jì)算效率.當(dāng)然,在由觀測(cè)時(shí)間序列計(jì)算上述平均值時(shí),會(huì)產(chǎn)生計(jì)算誤差,將會(huì)直接影響到上述約束條件的精度.

7 多行星系統(tǒng)的計(jì)算

坐標(biāo)速度法可以作為天體測(cè)量法探測(cè)系外行星的基本方法,不僅是因?yàn)樗梢杂糜谛∑穆受壍烙?jì)算,而且較容易將其推廣到多行星系統(tǒng),求解各個(gè)行星的參數(shù)(位置、速度及質(zhì)量).下面僅推導(dǎo)由1顆恒星和2顆行星組成的三體系統(tǒng)中恒星觀測(cè)方程計(jì)算具體關(guān)系式,n體問題的推導(dǎo)思路類似.

建立如圖4所示慣性坐標(biāo)系,恒星與兩顆行星質(zhì)量分別為Ms、m1和m2,O為系統(tǒng)的質(zhì)量中心,空間位置矢量如圖所示.

圖4 三體運(yùn)動(dòng)坐標(biāo)示意圖Fig.4 T he coord inate’s sketch of the th ree-body m otion

3個(gè)天體對(duì)應(yīng)的運(yùn)動(dòng)方程為

因系統(tǒng)質(zhì)心為坐標(biāo)原點(diǎn),由質(zhì)心運(yùn)動(dòng)定理知:

方程(25)左邊分別用η′(t)、ξ′(t)表示,可將其改寫為

由(24)–(25)式可求出行星運(yùn)動(dòng)方程為

狀態(tài)轉(zhuǎn)移方程涉及到的初始條件為:當(dāng)t=t0時(shí),

這里,H為3階單位矩陣.

8 結(jié)論與討論

本文基于天體測(cè)量法給出的恒星空間位置參數(shù)序列,討論了在求解行星質(zhì)量和軌道參數(shù)時(shí)涉及的動(dòng)力學(xué)條件方程計(jì)算問題;對(duì)于采用Kepler軌道根數(shù)和采用坐標(biāo)速度描述形式分別給出了具體微分改正公式,同時(shí)也進(jìn)行了初步必要的仿真模擬計(jì)算;通過天體力學(xué)二體問題運(yùn)動(dòng)關(guān)系式,對(duì)恒星空間位置時(shí)間序列做平均分析,推導(dǎo)出了待估參數(shù)間存在的3個(gè)約束條件關(guān)系式,可以顯著提升計(jì)算效率.

本文沒有對(duì)恒星自行和視差等天體測(cè)量參數(shù)做解算研究,以后若有實(shí)測(cè)數(shù)據(jù),可根據(jù)實(shí)際情況作具體處理:或單獨(dú)求解,或與行星軌道和質(zhì)量一起求解.

[1]M ayor M,Q ueloz D.Natu re,1995,368:355

[2]M ichael P.T he Exop lanet Handbook.New York:Cam b ridge Un iversity P ress,2011:3-5

[3]Luisa G J,Lu is A,Barbara P.M NRAS,2014,443:260

[4]姜明達(dá),肖東,朱永田.天文學(xué)進(jìn)展,2012,30:214

[5]Dou J P,Zhu Y T,Ren D Q.Ch inese Jou rna l of Natu re,2014,36:124

[6]張牛,季江徽.天文學(xué)進(jìn)展,2009,27:14

[7]Lin D N C,Ida S.A p J,1997,477:781

[8]Papaloizou J C B,Terquem C.M NRAS,2001,325:221

[9]金文敬.天文學(xué)進(jìn)展,2013,31:375

[10]W right J T,How ard A W.Ap JS,2009,182:205

[11]Sara S.Exop lanets.A rizona:Un iversity of A rizona P ress,2010:327-329

[12]M ichael P,Joel H,G aspar B,et a l.A p J,2014,797:14

[13]Th iele T N.A stronom ische Nach rich ten,1883,104:245

Exop lanet Detection by A strom etric M ethod

XUWei-wei1,2,3,4LIAO Xin-hao1,3ZHOU Yong-hong1,3XU Xue-qing1,3

(1 Shanghai A stronom ica l O bserva to ry,Chinese A cadem y of Scien ces,Shanghai 200030) (2 Schoo l of Physica l Scien ce an d Techno logy,ShanghaiTech Un iversity,Shanghai 201210) (3 K ey Labo ra to ry of P laneta ry Scien ces,Chinese A cadem y of Scien ces,Shanghai 200030) (4 Un iversity of Chinese A cadem y of Scien ces,Beijing 100049)

Aswe known,the exoplanetsaremostly detected by themethodsof radial velocity and transit,only one is found by the astrometric method.As the data of the gaia to be released,astrom etry w ill becom e one of themost im portantmethod for detecting exoplanets gradually.Based on the position sequence of stars,this paper discusses the calculation of the equations of dynam ics conditions involved in solving themass and the orbit parameters of the planet.Due to the deficiency of the available theory(orbital elem ent m ethod),we put forward a new m ethod(coordinate velocity method).The differential correction formulae of the twomethodsare presented,aswell as the necessary simulation.In addition,themethod established in this paper can be app lied to themulti-planet system easily.

exoplanet,astrometry,celestialmechanics:orbit calculation and determ ination,m ethods:data analysis

P135;

A

10.15940/j.cnki.0001-5245.2016.04.004

2015-12-08收到原稿,2016-03-04收到修改稿

?國(guó)家自然科學(xué)基金項(xiàng)目(11133004)和中國(guó)科學(xué)院先導(dǎo)B項(xiàng)目(XDBO9000000)資助

?1143880832@qq.com

?xhliao@shao.ac.cn