ASO-S/HXI太陽指向鏡算法研究?

喻福 蘇 楊 張哲 黃 宇

(1中國科學(xué)院紫金山天文臺南京210023)

(2中國科學(xué)技術(shù)大學(xué)天文與空間科學(xué)學(xué)院合肥230026)

(3中國科學(xué)院暗物質(zhì)與空間天文重點實驗室南京210023)

1 引言

先進天基太陽天文臺(Advanced Space-based Solar Observatory,ASO-S)是中國首顆太陽專用觀測衛(wèi)星[1],硬X射線成像儀(Hard X-ray Imager,HXI)作為其3臺載荷之一,主要負(fù)責(zé)在30–200 keV能段對耀斑源區(qū)進行成像、能譜和光變觀測,以研究耀斑磁重聯(lián)中的能量釋放和高能電子加速等物理過程[2].由于在該能段上僅能看到耀斑區(qū)域,而無法獲得全日面像,因此無法獲知源區(qū)在日面上的位置.因此太陽指向鏡作為HXI的組成部分,可實現(xiàn)對太陽中心位置的精確測量,提供指向中心坐標(biāo),并結(jié)合平臺的旋轉(zhuǎn)指向,以便對耀斑X射線源進行定位[2].指向鏡的本質(zhì)在于使用CCD(Charge-coupled Device)/CMOS(Complementary Metal Oxide Semiconductor)探測器對太陽進行可見光成像,再將得到的太陽光斑像近似為圓形,進而對圓的中心和半徑參數(shù)進行計算.太陽在可見光的輻射主要來自光球,其邊緣附近的輻射強度隨日心角距增加變化非常陡峭,形成銳利的太陽邊緣[3],這有利于提取太陽像的圓形輪廓特征.

HXI太陽指向鏡的原理與太陽導(dǎo)行鏡[4–5]相似(但無需反饋環(huán)節(jié)),都需計算日面中心在探測器視場中的坐標(biāo)位置,或者相對偏移量.而太陽望遠鏡的導(dǎo)行模式主要經(jīng)過了點、線、面的發(fā)展過程[6–8].在20世紀(jì)90年代以前多采用以太陽邊緣對稱的4個點來導(dǎo)行跟蹤的四象限探測器法,如早期的懷柔多通道望遠鏡[6];到80年代以后線陣探測器技術(shù)逐漸完善,德國GCT[9](Gregory Coud′e Telescope)、日本Yohkoh[10]和美國RHESSI[11–12](Reuven Ramaty High Energy Solar Spectroscopic Imager)均采用此方法.特別需要指出的是,正是由于RHESSI對指向(日面邊緣)的高精度測量,為了解太陽的實際形狀提供了最精確的結(jié)果[13].而今隨著計算機處理能力的提高,面陣探測器得到廣泛使用,如SOHO[14](Solar and Heliospheric Observatory).面陣探測器的優(yōu)點在于能夠獲取更多太陽的信息以及利用圖像處理技術(shù),更精確地反演日面中心位置.例如文獻[7]提到的全日面像相關(guān)算法和質(zhì)心算法,測試顯示其在x和y方向檢測精度的均方差分別約為0.0271′′和0.0101′′與0.1142′′和0.0828′′.

根據(jù)設(shè)計要求(來自《先進天基太陽天文臺衛(wèi)星(ASO-S)硬X射線成像儀初樣設(shè)計報告》,2019年,文件編號KX-07-HXI00-JB-04),HXI太陽指向鏡中心位置的測量精度要優(yōu)于2′′.由于定位精度越高,越有利于多波段的數(shù)據(jù)分析,因此我們在本文中設(shè)定測試目標(biāo)精度為1′′.由于數(shù)據(jù)量、采樣率等方面的要求和限制條件,我們無法下傳所有的圖像.為了節(jié)省在軌計算資源和減少傳輸?shù)臄?shù)據(jù)量,提出了四象限法和邊緣擬合法兩種獨立又可以結(jié)合的測量方法.為此,CMOS探測器每次曝光完成后,太陽指向鏡將對圖像進行在軌二值化處理,再采集兩種數(shù)據(jù)下傳,分別是14行數(shù)據(jù)和4個邊角的面積,而中心坐標(biāo)的計算在地面完成.其中,指向鏡在軌記錄的14行像元值用于獲得太陽邊緣,再利用最小二乘法進行擬合,得到太陽中心參數(shù).而CMOS 4個正方形邊角區(qū)域記錄的像元值為1的數(shù)量,用于四象限法的計算.本文中,我們利用MATLAB軟件和SDO(Solar Dynamic Observatory)/AIA(Atmospheric Imaging Assembly)[15]4500?A Level 1.0的數(shù)據(jù)主要針對算法本身進行了測試和評估.

2 太陽指向鏡測量原理

太陽指向鏡的鏡頭部分安裝在HXI載荷準(zhǔn)直器框架的前基板上,焦距f≈1200 mm[2].當(dāng)HXI對太陽觀測時,太陽光斑成像如圖1所示.其中,正方形ABCD為尺寸2048×2048 pixel的CMOS探測器,太陽光斑像的尺寸與CMOS的尺寸相當(dāng),如紅色圓所示.然而地球繞日公轉(zhuǎn)軌道為橢圓,近日點為1.471×108km、遠日點為1.521×108km[3],且太陽半徑為696342 km[16].因此,太陽視直徑隨地球繞日公轉(zhuǎn)軌道變化最大為32.55′,最小為31.48′,這會導(dǎo)致在CMOS上的成像隨之變化.再考慮到HXI指向鏡光軸相對于平臺指向的最大偏差(約為142 pixel)等因素,對應(yīng)在CMOS上成像光斑的邊緣被限制在綠色和黃色圓之間,對應(yīng)的半徑分別為Rmax=1174.6 pixel和Rmin=857.1 pixel.平行于x軸的14條紅色虛線與太陽像邊緣的交點作為邊緣擬合法的數(shù)據(jù)點.4個邊角上正方形的邊長被選為L=445 pixel,從而保證太陽像無論如何移動(在設(shè)計安裝誤差范圍內(nèi)),太陽邊緣都與4個區(qū)域的內(nèi)邊長有兩個交點.正方形陰影區(qū)被二值化并分別求和后,理論上對應(yīng)于圖1中深色區(qū)域的面積S1、S2、S3、S4,以此作為四象限法解算的基礎(chǔ).

圖1 太陽指向鏡算法原理圖.正方形ABCD為2048×2048 pixel CMOS探測器,紅色圓為太陽像光斑,綠色圓和黃色圓之間為儀器允許的光斑邊緣移動的區(qū)域.其中,在H(≈278 pixel)范圍內(nèi)選取了5行數(shù)據(jù);其他各符號的定義見2.2節(jié).Fig.1 Principle diagram of solar aspect system algorithm.The square AB CD is a CMOS detector with size of 2048×2048 pixels.The red circle represents the image of the solar disk with its limb limited between the green and yellow circles.And f ive rows of data are selected in the range of H(≈278 pixel).See Sec.2.2 for the def initions of more symbols.

2.1 邊緣擬合法

最小二乘法(Least Square Method,LSM)是用于曲線擬合最常用的方法,擬合過程快速簡便,且擬合精度相對較高.文獻[17]中利用最小二乘圓擬合法對激光光斑輪廓進行擬合,從而定出了其中心坐標(biāo)和半徑.該方法同樣可以用于指向鏡成像的日面邊緣擬合.設(shè)太陽像光斑的中心坐標(biāo)為(x0,y0),半徑為r0,則滿足的圓方程為:

設(shè)(xi,yi)為從太陽像光斑邊緣提取的第i個坐標(biāo)點,將其代入(1)式并取殘差為:

再求殘差平方和函數(shù):

根據(jù)最小二乘原理,有:

最后,解出滿足(4)式的方程組,可以得到圓的中心坐標(biāo)和半徑參數(shù)[17]:

2.2 四象限法

此處的方法同傳統(tǒng)四象限法不同,但由于算法本身的特征仍然稱之為四象限法(Four Quadrant Method,FQM),其原理如圖1所示.太陽光斑像(紅色圓面)分別與CMOS探測器4個角的正方形陰影區(qū)域相交,對應(yīng)的面積記為S1、S2、S3、S4,而該面積又可由紅色圓的參數(shù):圓心坐標(biāo)O(x0,y0)和半徑r0決定.以面積S1的計算為例,其面積等于扇形OA1B1的面積減去三角形OP1A1和OP1B1的面積,即

基于此,且以CMOS探測器陣列的頂點A為坐標(biāo)原點,向右和向下為正方向,建立直角坐標(biāo)系,則4個陰影區(qū)域的面積S1、S2、S3、S4滿足

其中,帶有下標(biāo)的θ、h、a、b為中間參量,分別表示扇形對應(yīng)的圓心角、三角形的高和底邊長度.而面積S1、S2、S3、S4分別對應(yīng)下面4組關(guān)系:

再將各中間參量的表達式(7)式分別代入面積S1、S2、S3、S4對應(yīng)的表達式(6)式中,則可以得到以太陽中心位置(x0,y0)和半徑r0為參數(shù)的方程組(8)式,即太陽的中心位置由面積S1、S2、S3、S4可以定出.

在解算時,S1、S2、S3、S4為CMOS探測器測量的已知量,由于方程組含有4個方程和3個未知數(shù),且是含有反三角函數(shù)的復(fù)雜方程,所以考慮利用MATLAB的fsolve1fsolve,https://www.mathworks.com/help/optim/ug/fsolve.html函數(shù)進行數(shù)值求解.fsolve可用于求解非線性方程組,需要給定初值,求解時可使用邊緣擬合法得到的(x0,y0,r0)作為數(shù)值求解的初值.此外,可以考慮兩種做法:FQM1,將(8)式和已知的面積S1、S2、S3、S4開方后再數(shù)值求解;FQM2,考慮降低噪聲的情況,將(8)式按如下方式改寫構(gòu)成新的方程組:(S1+S2)?(S3+S4),(S1+S3)?(S2+S4),S1?S4,S2?S3.

3 算法實現(xiàn)與誤差評估

3.1 數(shù)據(jù)來源與測試方法

(1)由于指向鏡數(shù)據(jù)來自太陽可見光,為使測試盡量接近真實情況,選擇了SDO/AIA 4500?A Level 1.0的光球數(shù)據(jù),該數(shù)據(jù)沒有作despike2despike,https://www.lmsal.com/sdodocs/doc/rep/sdod180/fid258/zip/entry/index.html處理.此外4500?A的數(shù)據(jù)尺寸大小為4096×4096 pixel,可以使用imresize函數(shù)將其縮放到2048×2048 pixel.若在縮放前先進行裁剪,則可以控制生成的太陽像半徑大小r0;

(2)為了評定指向鏡算法的精度,首先直接讀取f its文件為2維數(shù)組像元亮度值(Digital Number,DN),經(jīng)裁剪和縮放作為第1張AIA圖像(日面中心在圖像中心附近),并計算其日面中心為(X0,Y0).其次對圖像沿x和y方向平移(δx,δy)以產(chǎn)生大量測試圖像,其理論中心為(XT,YT)=(X0+δx,Y0+δy).其中δx、δy從?150 pixel到150 pixel按3 pixel等間隔選取,且滿足,即指向鏡光軸相對于平臺指向的最大偏差條件;

(3)另一方面,通過算法計算出測試圖像的日面中心坐標(biāo)(X,Y),再用計算得到的坐標(biāo)與理論坐標(biāo)相減得到x和y方向的誤差(X?XT,Y?YT),則算法的誤差可以表示為[6]

此外,一般認(rèn)為圖像是由一個個方格狀的像素單元組成,方格的數(shù)量與像素的數(shù)量對應(yīng).為了反映這種認(rèn)知,可以將CMOS圖像的坐標(biāo)原點(即A點)改為(0.5,0.5)pixel.于是對于2048×2048 pixel的CMOS探測器,其中心在(1024.5,1024.5)pixel,故在解算后需要對原中心坐標(biāo)作+0.5的修正;

(4)根據(jù)太陽在CMOS探測器上成像光斑大小的可能變化,最小直徑為1997.8 pixel,最大直徑為2065.7 pixel,故設(shè)置了8組不同半徑大小的AIA圖像進行測試,r0分別約為:998、1004、1009、1014、1019、1024、1029、1035 pixel.此外考慮了圖像二值化時的閾值和噪聲影響,噪聲是在圖像二值化之后添加的.關(guān)于添加噪聲,首先隨機產(chǎn)生2048×2048的只含有0和1的二值化矩陣.其中1表示噪點,其數(shù)量與0和1總數(shù)量的比值p定義為噪聲比例.然后再將噪聲矩陣與二值化后的AIA圖像做“或”邏輯運算,得到具有一定噪聲水平p的測試圖像.

3.2 測試結(jié)果與分析

為了探索隨機噪聲對指向鏡算法精度的影響,圖2中給出了四象限法的兩種解法FQM1和FQM2以及LSM的誤差隨噪聲的變化規(guī)律.其中橫坐標(biāo)的噪聲比例在0%–1%間等間距取10個點,在2%–5%間等間距取4個點.圖中每一個噪聲比例下對應(yīng)的藍色*點來自821張偏移測試圖像.這些測試圖像從?150 pixel到150 pixel按9 pixel等間隔產(chǎn)生且滿足指向鏡光軸最大偏差條件.這些誤差點總體上對誤差的上限有一個限制.而紫紅色的線條是測試圖像在約70 pixel固定偏移下,其平均誤差隨噪聲比例的變化曲線.這說明隨著噪聲的增加,FQM算法的誤差呈增大趨勢.此外,圖3也給出了FQM和LSM的誤差隨噪聲比例和圖像偏移量的分布圖.

圖2 不同比例隨機噪聲下的指向鏡算法誤差(閾值=1500 DN,r 0=1014 pixel).每個噪聲比例下的藍色*點來自于821張測試圖像,紫紅色線條是偏移約為70 pixel測試圖像的平均誤差曲線.(a)和(b)中的子圖是滿足誤差小于1 pixel(紅色直線以下)的局部放大圖.Fig.2 The error of solar aspect system algorithm under different proportions of random noise(threshold=1500 DN,r 0=1014 pixel).The blue points*under each noise ratio are the errors of 821 test images,and the magenta curve is the average error of the test image with 70 pixel of fset.The subgraphs in(a)and(b)are the local enlarged ones satisfying the error less than 1 pixel(under the red line).

由圖2可知,FQM1和FQM2滿足誤差小于1 pixel的噪聲比例分別約為0.1%和0.85%(對于偏移量為70 pixel的圖像,FQM能忍受的噪聲上限分別約為0.4%和1.6%).而LSM的誤差隨噪聲的變化相對穩(wěn)定且小于0.4 pixel.此外,由圖2(a)的子圖可知FQM1的誤差是隨噪聲逐漸增加的,在大約小于0.02%的噪聲水平時,此方法的精度優(yōu)于FQM2和LSM.而LSM對噪聲的適應(yīng)性以及精度總體上好于FQM,故后面的測試主要基于FQM1和LSM,當(dāng)然3種計算方法可以用來相互驗證.

圖3 FQM和LSM的誤差隨噪聲比例和圖像偏移量的分布圖,對應(yīng)于圖2的結(jié)果.Fig.3 The error distribution of FQM and LSM with noise percentage and image of fset in Fig.2.

基于圖2和圖3的初步分析,接著做了更細致的測試.圖4給出了對邊緣擬合法和四象限法在固定閾值、隨機閾值和附加隨機噪聲情況下的6種測試結(jié)果.其中橫坐標(biāo)為8組不同半徑r0測試圖像的統(tǒng)一編號,大約56000張.縱坐標(biāo)為算法誤差?T的浮動范圍.表1也給出了相應(yīng)誤差的統(tǒng)計結(jié)果,包括均方根(Root Mean Square,RMS)、平均值、極值以及99.5%分位值.

圖4 太陽指向鏡算法在不同條件下的誤差分布:(a)和(b)為固定閾值1500 DN,(c)和(d)在1200?1900 DN之間隨機取閾值,(e)和(f)除了隨機閾值外再附加0.05%的隨機噪聲.Fig.4 The error distribution of solar aspect system algorithm under different conditions:(a)and(b)with fixed thresholds of 1500 DN,(c)and(d)with random thresholds of 1200?1900 DN,(e)and(f)with noise percentage of 0.05%in addition to random thresholds.

表1 太陽指向鏡算法的誤差評估:LSM和FQM 1Table 1 Error evaluation of solar aspect system algorithm:LSM and FQM 1

結(jié)合圖4和表1可知,在固定閾值下四象限法的誤差最小,誤差為0.05′′以內(nèi)的數(shù)據(jù)占99.5%,最大誤差不超過0.06′′,優(yōu)于邊緣擬合法的結(jié)果.在單獨考慮隨機閾值以及再附加0.01%的隨機噪聲情況下,FQM1的各項指標(biāo)依然優(yōu)于LSM.當(dāng)對測試圖像添加0.05%的噪聲時,FQM1的極值誤差接近0.8′′,相對于LSM增加更為顯著.添加隨機噪聲后的圖4(f)的誤差分布相對于無噪聲的圖4(d)變化更明顯.不過根據(jù)圖4(f),極值誤差與太陽像的半徑和偏移量有關(guān)系,故誤差的評估也需要借助其他統(tǒng)計量,如RMS誤差為0.23′′.由于四象限法的精度取決于CMOS上太陽像經(jīng)二值化后4個邊角的面積,可以預(yù)計隨著噪聲的繼續(xù)增加,FQM對面積的分辨變差會導(dǎo)致誤差增大,這也與圖2(a)和(b)的結(jié)果相符.而對于LSM,比較圖4(a)、(c)和(e)可知,其誤差波動都比較隨機均勻,再結(jié)合表1中LSM的統(tǒng)計結(jié)果,有無噪聲相差不大;也符合圖2(c)的預(yù)期.

圖5選取了圖4(e)和(f)中半徑r0=1014 pixel的數(shù)據(jù)點,畫出了LSM和FQM1的誤差隨圖像偏移量的分布.其中橫、縱坐標(biāo)分別為沿x和y方向的圖像偏移量,左圖為邊緣擬合法,其誤差分布相對均勻.右圖為四象限法,根據(jù)色值可知通常在太陽圖像的偏移量較小的時候其誤差較小.該結(jié)果也可與圖2和圖3的分析相互印證.

圖5 太陽指向鏡算法誤差隨圖像偏移量的分布.對應(yīng)圖4(e)和4(f)中r 0=1014 pixel的測試數(shù)據(jù).Fig.5 The error distribution of the solar aspect system algorithm with the image of fset and corresp onding to the test data of r 0=1014 pixel in Fig.4(e)and 4(f).

4 總結(jié)與討論

我們針對HXI的耀斑源區(qū)定位需求以及HXI太陽鏡指向精度優(yōu)于2′′的設(shè)計要求,在查閱和借鑒太陽導(dǎo)行鏡相關(guān)資料的基礎(chǔ)上,對兩種基礎(chǔ)算法,即最小二乘法和四象限法的測量精度及其受噪點水平的影響進行了測試.測試時利用了SDO/AIA 4500?A的Level 1.0數(shù)據(jù).

總體而言,邊緣擬合法和四象限法的精度(在0.05%噪聲下RMS誤差分別約為0.11′′和0.23′′)均優(yōu)于2′′的設(shè)計要求和1′′的測試目標(biāo),且可提供獨立測量結(jié)果,用于交叉驗證.四象限法依賴于4個角的面積對太陽中心坐標(biāo)的約束,所以其精度對面積的變化比較敏感.前面提到的兩種解法,FQM1可以達到最高精度時的噪聲容忍約為0.02%,FQM2雖然能夠容忍約0.85%的噪聲,但是精度總體不如LSM好.而邊緣擬合法的定位精度雖然相對略低,但是由于該算法是基于擬合太陽像的邊緣點數(shù)據(jù),其所選取的14行像素相對于總的2048行像素而言,在概率上能保證具有較強的抗隨機噪聲能力.

但需要指出的是HXI太陽指向鏡在軌提取太陽圖像時,將無法對原始數(shù)據(jù)作復(fù)雜的矯正處理,故在二值化后如果還存在較多噪點,則邊緣擬合法會更適合作為指向鏡算法.就測試使用的SDO/AIA 4500?A Level 1.0數(shù)據(jù)而言,四象限法的FQM1解法在低噪聲時的指向精度更高.但Level 1.0數(shù)據(jù)是Level 0數(shù)據(jù)經(jīng)過一定處理得到的,當(dāng)前測試中,我們尚未考慮這兩種數(shù)據(jù)的具體差異,所以測試會有局限性,特別是對于四象限算法.在本文所得結(jié)果的基礎(chǔ)上,我們將進一步優(yōu)化算法,如采用窮舉矩陣方法以及改進四象限法方程組的求解方法等,以改進太陽中心位置的測量精度.