溫 驍，劉 煜，宋騰飛，張雪飛，劉念平

(中國科學(xué)院國家天文臺/云南天文臺，云南 昆明 650011)

現(xiàn)代日暈光度計(jì)是用于西部太陽選址的重要儀器，它能觀測并記錄多波段的大氣參量。隨著觀測的進(jìn)行，觀測到的數(shù)據(jù)量不斷龐大，自動(dòng)處理觀測資料也因此成為必要。在對觀測數(shù)據(jù)的處理中，日面中心坐標(biāo)的確定對數(shù)據(jù)處理結(jié)果具有決定性影響，這是自動(dòng)處理數(shù)據(jù)的第一步，也是從觀測資料中獲取信息最重要的一步。發(fā)展了兩種方法用于日面中心的自動(dòng)定位并在文中逐一介紹，最后對比了兩種方法所得的日心坐標(biāo)的差別及其對測量結(jié)果的影響。

1 數(shù)據(jù)介紹

日暈光度計(jì)的數(shù)據(jù)為4個(gè)波段，每幀圖像大小為765 pixel×510 pixel。在實(shí)際采集時(shí)一般是每個(gè)數(shù)據(jù)之間間隔15 s，每一分鐘采集一套4波段數(shù)據(jù)。在天氣好的時(shí)候，一天可觀測時(shí)間約為10 h，共2 400個(gè)數(shù)據(jù)，大小約1.8 GB。

文中使用的數(shù)據(jù)是在云南天文臺麗江高美古觀測站2011年2月25日觀測得到的，選取的時(shí)間是從8:30到13:30(UT為0:30至5:30)，數(shù)據(jù)質(zhì)量非常好，4個(gè)波段共1044幀數(shù)據(jù)。此段時(shí)間內(nèi)天氣非常好，天空萬里無云。4個(gè)波段的觀測曝光時(shí)間、濾光片通透中心波長、帶寬以及透過率列于表1。

表1 日暈光度計(jì)的4個(gè)波段的觀測曝光時(shí)間、通透波長、通透帶寬以及ND4濾光片在各波段的透過率Table 1 The exposure times，central wavelengths，widths，and throughputs for the four bands of the ND4 filter on the solar-halo photometer

關(guān)于日暈光度計(jì)儀器參數(shù)的詳細(xì)資料見文［1］，詳細(xì)數(shù)據(jù)見文［2］和文［3］。

2 數(shù)據(jù)處理目的

觀測數(shù)據(jù)中包括日心強(qiáng)度和日暈強(qiáng)度，計(jì)算的均是這些參量在CCD上的讀數(shù)，在下文用日心強(qiáng)度表示日心強(qiáng)度在CCD上的讀數(shù)，日暈強(qiáng)度表示日暈強(qiáng)度在CCD上的讀數(shù)。這些數(shù)據(jù)作為后期處理則需要考慮ND4濾光片的透過率，但這并不是本文考慮的范圍，因此并未修正。將實(shí)際圖像中的日面像記為日面，其中心在圖像中的坐標(biāo)記為日心坐標(biāo)。

由于數(shù)據(jù)中含有儀器衍射光，和ND4濾光片支架的影子，這些區(qū)域的數(shù)據(jù)是不需要的，應(yīng)當(dāng)去掉;進(jìn)而分離出日暈區(qū)域，用于計(jì)算距離日心一定半徑處的日暈強(qiáng)度。由于需要求取日心強(qiáng)度以及距離日心一定半徑內(nèi)的日暈強(qiáng)度平均值、日心強(qiáng)度以及定標(biāo)后的日暈強(qiáng)度(日暈強(qiáng)度與日心強(qiáng)度的比值)，這些都與日心坐標(biāo)直接相關(guān)，如果日心坐標(biāo)選取錯(cuò)誤，將直接導(dǎo)致數(shù)據(jù)結(jié)果出現(xiàn)偏差。因此發(fā)展一套自動(dòng)并且準(zhǔn)確的日心坐標(biāo)尋找方法非常必要。

3 日面中心自動(dòng)定位方法1:日面總強(qiáng)度法

要找到日面中心，需討論選取日面以及求取日面中心的標(biāo)準(zhǔn)。

在可見光下太陽極度近似圓形，在理想的觀測條件下，日面應(yīng)該是標(biāo)準(zhǔn)圓形，其中心就是這個(gè)圓形的質(zhì)心(幾何中心)?？紤]使用日面輪廓的質(zhì)心為日面坐標(biāo)，但是在實(shí)際觀測中有諸多限制:觀測條件限制，如系統(tǒng)略微的雜散光，濾光片或CCD上的灰塵，在陽光路徑上飄過的塵土，或是看不見的薄云，薄云會影響日面的亮度梯度，進(jìn)而影響邊緣輪廓的選取，ND4之間的多次反射成像。這些都可能使日面非理想成像，造成日心坐標(biāo)的較大差異。在選取日面輪廓的時(shí)候需要自定義一個(gè)邊界數(shù)值，這樣加入了人工選擇因素。因此使用日面輪廓質(zhì)心定日心坐標(biāo)的方法并不可取。

參看圖1(a)，太陽日面強(qiáng)度隨半徑成梯度變化，離太陽中心越遠(yuǎn)強(qiáng)度越小，在靠近日面邊緣處強(qiáng)度變化極度劇烈。

圖1 (a)經(jīng)過日面的一條線上的強(qiáng)度變化。(b)日面區(qū)域的等高線。在觀測條件好的時(shí)候，觀測得到的日面非常接近圓形Fig.1 (a)The brightness across the solar disk.(b)Contours on the solar disk.The solar disk in the image is nearly a perfect circle under good condition

如果選取正確的日心坐標(biāo)和日面半徑，則所有在這個(gè)半徑內(nèi)的點(diǎn)的強(qiáng)度和會達(dá)到最大值，因?yàn)橹灰@個(gè)中心點(diǎn)稍微移動(dòng)，就會使圓環(huán)偏離日面，使得內(nèi)部強(qiáng)度和迅速降低。這一點(diǎn)可以從日面區(qū)域的等高線上看出，如圖1(b)，靠近日面邊緣等高線變化快，說明強(qiáng)度變化快，相比之下日面中心強(qiáng)度變化要慢，圖2為日面大小圓環(huán)內(nèi)強(qiáng)度和橢圓環(huán)中心變化，偏離實(shí)際日心，強(qiáng)度和迅速降低，也說明了這一點(diǎn)。搜尋方法是選取一定范圍內(nèi)的所有點(diǎn)，每個(gè)點(diǎn)計(jì)算所有距該點(diǎn)為固定半徑的圓環(huán)內(nèi)的點(diǎn)的強(qiáng)度和，找到強(qiáng)度和最大的那個(gè)點(diǎn)作為日心。這個(gè)半徑盡量選取日面邊緣強(qiáng)度變化最大時(shí)對應(yīng)的值。具體計(jì)算步驟如下。

(1)自動(dòng)尋找日面中心

觀測中日面在圖像上經(jīng)常變化。但是日暈光度計(jì)自動(dòng)跟蹤系統(tǒng)和人為調(diào)整可以保證日面始終處在減光片遮體的衍射光環(huán)內(nèi)。

在程序中選取半徑為R=47 pixel，下文討論選取不同半徑的日心識別情況，以及是如何選取這個(gè)半徑的。

由于日面并未超出減光片遮體的衍射光環(huán)，為了減少不必要的計(jì)算時(shí)間，日心坐標(biāo)的搜尋就在這個(gè)范圍內(nèi)進(jìn)行。將得到的日面中心記為點(diǎn)(Cx，Cy)。

(2)程序中R選取不同值時(shí)的日心識別情況

在不能確定R的取值時(shí)，可以先定幾個(gè)值分別計(jì)算日心坐標(biāo)，后面將有算法用以確定R的取值，取不同R得到的日心坐標(biāo)見表2。

表2 取不同R得到的日心坐標(biāo)Table2 The solar-disk center coordinates for different R values

因此選擇點(diǎn)(389，263)作為例圖最終的日心。

(3)計(jì)算離日心一定距離處的日暈強(qiáng)度

圖2 自動(dòng)尋找日心的算法圖示每個(gè)白色圓環(huán)有一個(gè)中心，相應(yīng)的也有在這個(gè)圓環(huán)內(nèi)部所有的點(diǎn)的強(qiáng)度總和。圖示曲線為經(jīng)過日面中心的一條線上的點(diǎn)，以及按上述方法求得的強(qiáng)度分布Fig.2 Illustration for the algorithm of calculating the solar central positionEach white circle has a center and a total enclosed intensity.The curve shows the variation of the total enclosed intensity of a circle when its center is moving along the white line through the solar-disk center

將原圖去除衍射光以及支架影子部分后，計(jì)算數(shù)據(jù)中所有距日心為r的點(diǎn)強(qiáng)度的算術(shù)平均為距離日心r處的日暈強(qiáng)度，記為Isky(r)。去除衍射光和支架影子后的數(shù)據(jù)見圖3(a)，日暈強(qiáng)度隨離日心距離的變化曲線見圖3(b)。

圖3 (a)去除衍射光以及支架影子部分后的數(shù)據(jù)。(b)Isky(r)隨r變化的曲線，曲線中為0的部分是因?yàn)槿コ搜苌涔夂罅粝碌目障禙ig.3 (a)The image where the system scattered light and the shadow of the three trestles have been removed.(b)Isky(r)versus r.The gap in this curve is due to the removal of the scattered light

(4)日面邊緣強(qiáng)度變化最劇烈的半徑的確定

使用Isky(r)分析程序中R的選取，在日面邊緣處有最大的強(qiáng)度變化。這里做了類似于函數(shù)中的一次求導(dǎo)的處理，定義函數(shù):

該函數(shù)隨r變化曲線見圖4(a)，取得最大值時(shí)對應(yīng)的r為程序中所需的R值。此半徑確認(rèn)后，對于處理未變動(dòng)儀器所取得的數(shù)據(jù)時(shí)均不變。取日心坐標(biāo)為 (Cx，Cy)和R=47為半徑的圓環(huán)畫在原圖上，顯示此日心坐標(biāo)和R的選取很準(zhǔn)確。

圖4 (a)Isky(r)－Isky(r+1)隨r變化的曲線，在r=47處強(qiáng)度變化最大，取R=47。(b)以點(diǎn)(cx，cy)為中心，以r=47為半徑，標(biāo)在原圖上Fig.4 (a)Isky(r)－ Isky(r+1)versus r，with the maximum intensity change at r=47.(b)Image with a white circle centered at(cx，cy)and of radius 47

(5)計(jì)算機(jī)程序完整步驟:

①找到成像系統(tǒng)的中心，去除濾光片支架以及鏡筒的散射光;

②找到支架位置，去除其影子在CCD上成像的那部分?jǐn)?shù)據(jù);

③在減光片遮體衍射光環(huán)內(nèi)部，選擇半徑R=47的圓環(huán)內(nèi)強(qiáng)度和最大值的點(diǎn)為日心，記為(Cx，Cy);

④將去除濾光片支架和鏡筒散射光以及支架影子后的數(shù)據(jù)，計(jì)算距日心為r的日暈強(qiáng)度分布Isky(r);計(jì)算所有距日心為4.5實(shí)際太陽半徑和7.8實(shí)際太陽半徑內(nèi)的點(diǎn)平均值為Ihalo;計(jì)算以點(diǎn)(Cx，Cy)為中心，加上下左右4個(gè)點(diǎn)共5個(gè)點(diǎn)取算術(shù)平均，為日心強(qiáng)度，記為Io。

4 日面中心自動(dòng)定位方法2:傅里葉變換相關(guān)法

該方法原理簡要介紹如下，首先定義一個(gè)與真實(shí)數(shù)據(jù)大小相同(765 pixel×510 pixel)的數(shù)組。然后以ND4濾光片投影中心(x0，y0)(即上述計(jì)算機(jī)程序第1步的成像系統(tǒng)中心)為初始日心坐標(biāo)值，采用半徑R=47(由方法1中的第4步計(jì)算得到)個(gè)像素的圓環(huán)為日面。數(shù)組的日面內(nèi)部點(diǎn)賦值為1，數(shù)組其它像素賦值為0，生成mask數(shù)據(jù)g(x，y)。其次，類似于方法1，將原圖的鏡筒口套圈衍射光以及ND4支架影子部分去除(圖3a)，使之成為目標(biāo)數(shù)組f(x，y)。最后，利用傅里葉變換中的相關(guān)定理，對函數(shù)f(x，y)和g(x，y)進(jìn)行快速傅里葉相關(guān)分析，就得到mask數(shù)據(jù)與目標(biāo)數(shù)據(jù)在相關(guān)最大時(shí)的Δx與Δy值，于是得到日面中心的實(shí)際坐標(biāo)(x0+Δx，y0+Δy)。其中傅里葉變換的相關(guān)定理表達(dá)式是:

兩個(gè)程序可以使用C、MATLAB、IDL等語言編寫。程序能夠得到觀測數(shù)據(jù)中的有效數(shù)據(jù)信息，經(jīng)處理后的數(shù)據(jù)可用于進(jìn)一步分析。

5 兩種方法識別的日心坐標(biāo)對比

計(jì)算了兩種方法在相同參數(shù)條件下的日心坐標(biāo)，由于兩種方法的判斷標(biāo)準(zhǔn)不一樣，需要判定識別日心坐標(biāo)的準(zhǔn)確性和兩方法識別坐標(biāo)的差別，以及日心坐標(biāo)的差別對所獲取的數(shù)據(jù)的影響，對比的指標(biāo)有日心坐標(biāo)、日暈強(qiáng)度、日心強(qiáng)度以及定標(biāo)后的日暈強(qiáng)度。日暈強(qiáng)度和日心強(qiáng)度均按照方法1的算法計(jì)算。

對每組數(shù)據(jù)均對比了以上指標(biāo)，對比結(jié)果見圖5～7。

從圖5可以看出兩種方法求得的日心坐標(biāo)并不完全一致。為了準(zhǔn)確地描述兩種方法的差異，將方法1得到的坐標(biāo)減方法2得到的坐標(biāo)，對坐標(biāo)的差異進(jìn)行了統(tǒng)計(jì)，以占總數(shù)據(jù)的百分比列于表3。

從表3可以看出兩種方法得到的結(jié)果相同的較多，相差1個(gè)像素的非常少，沒有相差超過1個(gè)像素的情況。經(jīng)仔細(xì)檢查程序，認(rèn)為這種差別是由于兩種方法的差別造成的。

從圖6、7可以看出由于坐標(biāo)的差別造成數(shù)據(jù)結(jié)果的差別:日心強(qiáng)度的差別均在1%以下，定標(biāo)后的日暈強(qiáng)度(日暈強(qiáng)度與日心強(qiáng)度的比值)差別多數(shù)在0.5%以下，整體未超出1%?？偟膩碚f，兩種方法得到的坐標(biāo)差別對結(jié)果的影響非常小。

圖5 4個(gè)波段所有數(shù)據(jù)的日心坐標(biāo)識別對比上面為x坐標(biāo)的差別，下面為y坐標(biāo)的差別。橫軸是方法1得到的坐標(biāo)，縱軸是方法2得到的坐標(biāo)Fig.5 Comparison between the solar-disk center coordinates in the four wavelength bands measured with the two methods The top row is for the x coordinate，and the bottom row is for the y coordinate

表3 兩種方法計(jì)算得的各波段數(shù)據(jù)的日心坐標(biāo)x坐標(biāo)和y坐標(biāo)的相差情況Table 3 The discrepancies of the x and y coordinates from the two methods(in percentages)

圖6 各波段數(shù)據(jù)用兩種方法求得的日心強(qiáng)度差別的百分比Fig.6 Differences between the solar-center intensity values from the two methods for all four wavelength bands(in percentages)

圖7 各波段數(shù)據(jù)用兩種方法求得的定標(biāo)后日暈強(qiáng)度的差別百分比Fig.7 Differences between the calibrated solar-halo intensity values from the two Methods(in percentages)

6 結(jié)論

本文采用兩種方法對日暈光度計(jì)觀測數(shù)據(jù)的日心坐標(biāo)位置進(jìn)行自動(dòng)計(jì)算。第1種方法是基于日面總強(qiáng)度流量最大值原理，第2種方法是基于實(shí)際和理論日面模式識別原理。兩種日心識別方法雖然本質(zhì)相近，但實(shí)際采用的判斷標(biāo)準(zhǔn)不一樣。利用它們對同一時(shí)段內(nèi)的數(shù)據(jù)進(jìn)行計(jì)算，發(fā)現(xiàn)取得的日心坐標(biāo)結(jié)果在4個(gè)波段均相差甚微(一個(gè)像素之內(nèi))。在實(shí)際情況下得到100%正確的日心坐標(biāo)是很困難的，另一方面，這兩種方法得到的日心坐標(biāo)的微小差別也是在計(jì)算方法本身的理論誤差范圍內(nèi)，即不大于一個(gè)像素。因此認(rèn)為這兩種方法均可提供有效的日心坐標(biāo)信息，使用者可以根據(jù)自己的偏好選擇其中的任何一種方法進(jìn)行日心坐標(biāo)的計(jì)算。

圖8 由方法1計(jì)算的4波段日面中心強(qiáng)度(左圖)以及定標(biāo)后日暈強(qiáng)度(右圖)隨觀測時(shí)間的變化曲線Fig.8 Variations of solar-center intensity(left)and calibrated solar-halo intensity with observations(right)for the four wavelength bands as calculated by Method 1

另外，從對比2個(gè)參量(日暈強(qiáng)度以及定標(biāo)后的日暈強(qiáng)度)的差別情況看，這種坐標(biāo)間的略微差別對取得的這些歸算參量影響極小。因此，這兩種方法均可用于實(shí)際計(jì)算，用它們計(jì)算得到的正確日心坐標(biāo)奠定了進(jìn)一步數(shù)據(jù)處理的基礎(chǔ)(日心強(qiáng)度隨時(shí)間變化的曲線和定標(biāo)后的日暈強(qiáng)度隨時(shí)間變化曲線如圖8)，2個(gè)參量隨觀測時(shí)間的變化顯示出觀測地大氣對太陽光的散射情況，通過計(jì)算可以驗(yàn)證大氣對太陽光的散射規(guī)律，并求得觀測地大氣的氣溶膠、塵埃和水汽等大氣成分。

［1］劉順慶，段輯，張雪飛，等.現(xiàn)代日暈光度計(jì)多波段測光系統(tǒng)初步測試結(jié)果分析 ［J］.天文研究與技術(shù)——國家天文臺臺刊，2012，9(2):168－175.Liu Shunqing，Duan Ji，Zhang Xuefei，et al.Analysis of the Preliminary Measurements Taken by the Multi-color Photometric System of the Sky Brightness Monitor［J］.Astronomical Research＆ Technology——Publications of National Astronomical Observatories of China，2012，9(2):168－175.

［2］Lin Haosheng，Penn Matthew J.The Advanced Technology Solar Telescope site survey Sky Brightness Monitor［J］.The Publications of the Astronomical Society of Pacific，2004，116(821):652－666.

［3］劉念平，劉煜，申遠(yuǎn)燈，等.日暈測量與日暈光度計(jì)外緣雜散光抑制試驗(yàn) ［J］.天文學(xué)報(bào)，2011，52(2):160－170.Liu Nianping，Liu Yu，Shen Yuandeng，et al.Measurement of Sky Brightness and Suppression of Scattering in Sky Brightness Monitor［J］.Acta Astronomica Sinica，2011，52(2):160－170.