日暈光度計圖像日心自動定位方法*

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

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

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

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

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

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

表1 日暈光度計的4個波段的觀測曝光時間、通透波長、通透帶寬以及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)于日暈光度計儀器參數(shù)的詳細(xì)資料見文［1］，詳細(xì)數(shù)據(jù)見文［2］和文［3］。

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

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

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

3 日面中心自動定位方法1:日面總強度法

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

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

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

圖1 (a)經(jīng)過日面的一條線上的強度變化。(b)日面區(qū)域的等高線。在觀測條件好的時候，觀測得到的日面非常接近圓形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)和日面半徑，則所有在這個半徑內(nèi)的點的強度和會達(dá)到最大值，因為只要這個中心點稍微移動，就會使圓環(huán)偏離日面，使得內(nèi)部強度和迅速降低。這一點可以從日面區(qū)域的等高線上看出，如圖1(b)，靠近日面邊緣等高線變化快，說明強度變化快，相比之下日面中心強度變化要慢，圖2為日面大小圓環(huán)內(nèi)強度和橢圓環(huán)中心變化，偏離實際日心，強度和迅速降低，也說明了這一點。搜尋方法是選取一定范圍內(nèi)的所有點，每個點計算所有距該點為固定半徑的圓環(huán)內(nèi)的點的強度和，找到強度和最大的那個點作為日心。這個半徑盡量選取日面邊緣強度變化最大時對應(yīng)的值。具體計算步驟如下。

(1)自動尋找日面中心

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

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

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

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

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

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

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

(3)計算離日心一定距離處的日暈強度

圖2 自動尋找日心的算法圖示每個白色圓環(huán)有一個中心，相應(yīng)的也有在這個圓環(huán)內(nèi)部所有的點的強度總和。圖示曲線為經(jī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

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

圖3 (a)去除衍射光以及支架影子部分后的數(shù)據(jù)。(b)Isky(r)隨r變化的曲線，曲線中為0的部分是因為去除了衍射光后留下的空隙Fig.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)日面邊緣強度變化最劇烈的半徑的確定

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

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

圖4 (a)Isky(r)－Isky(r+1)隨r變化的曲線，在r=47處強度變化最大，取R=47。(b)以點(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)計算機程序完整步驟:

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

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

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

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

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

該方法原理簡要介紹如下，首先定義一個與真實數(shù)據(jù)大小相同(765 pixel×510 pixel)的數(shù)組。然后以ND4濾光片投影中心(x0，y0)(即上述計算機程序第1步的成像系統(tǒng)中心)為初始日心坐標(biāo)值，采用半徑R=47(由方法1中的第4步計算得到)個像素的圓環(huán)為日面。數(shù)組的日面內(nèi)部點賦值為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)進行快速傅里葉相關(guān)分析，就得到mask數(shù)據(jù)與目標(biāo)數(shù)據(jù)在相關(guān)最大時的Δx與Δy值，于是得到日面中心的實際坐標(biāo)(x0+Δx，y0+Δy)。其中傅里葉變換的相關(guān)定理表達(dá)式是:

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

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

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

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

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

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

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

圖5 4個波段所有數(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 兩種方法計算得的各波段數(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ù)用兩種方法求得的日心強度差別的百分比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)后日暈強度的差別百分比Fig.7 Differences between the calibrated solar-halo intensity values from the two Methods(in percentages)

6 結(jié)論

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

圖8 由方法1計算的4波段日面中心強度(左圖)以及定標(biāo)后日暈強度(右圖)隨觀測時間的變化曲線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個參量(日暈強度以及定標(biāo)后的日暈強度)的差別情況看，這種坐標(biāo)間的略微差別對取得的這些歸算參量影響極小。因此，這兩種方法均可用于實際計算，用它們計算得到的正確日心坐標(biāo)奠定了進一步數(shù)據(jù)處理的基礎(chǔ)(日心強度隨時間變化的曲線和定標(biāo)后的日暈強度隨時間變化曲線如圖8)，2個參量隨觀測時間的變化顯示出觀測地大氣對太陽光的散射情況，通過計算可以驗證大氣對太陽光的散射規(guī)律，并求得觀測地大氣的氣溶膠、塵埃和水汽等大氣成分。

［1］劉順慶，段輯，張雪飛，等.現(xiàn)代日暈光度計多波段測光系統(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］.天文學(xué)報，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.