應用于低對比度擴展目標觀測的大型陣列太陽自適應光學電子系統的設計與實現

費瑋瑋 王長清 劉濮鯤 蔡惠智 姜愛民 唐清善(中國科學院電子學研究所 北京 100190)(中國科學院聲學研究所 北京 100190)(中國科學院國家天文臺 北京 100012)

1 引言

自適應光學由于可以有效補償光的大氣傳輸效應，達到或接近光學系統的衍射極限，被廣泛地應用于天文觀測中。目前已有的自適應光學系統主要是用于夜間對恒星的觀測[1-4]。與這一相對成熟的夜天文自適應光學技術相比，白天對太陽觀測的太陽自適應光學更具挑戰性，成為自適應光學研究的難點。

針對太陽黑子或太陽表面米粒結構這一類低對比度擴展目標，通常對恒星采用的運算量較小的質心算法不再適用，一般采用協方差相關算法或者差絕對值和法[5]，這對電子系統提出了大計算量、高實時性的要求。目前基于FPGA的系統設計研究[6]，取得一定的成就，但隨著系統規模的擴大，單獨使用FPGA尚不能滿足實際應用需求。DSP因具有強大的并行數據處理能力，被引入電子系統的設計[7]。如今，將DSP與FPGA結合起來，采用合適的優化算法，使得構建更大規模、更好實時性的自適應光學系統成為可能。

本文針對太陽表面低對比度擴展目標的波前校正的應用，設計了一種基于FFT協方差相關算法的電子系統。該系統采用陣列DSP與FPGA協同合作的架構，由低階相關跟蹤系統和高階波前處理系統組成。系統具有較好的靈活性、通用性和擴展性。下文先對采用算法進行簡單介紹，接著分析系統需求，最后給出系統的整體設計方案。該電子系統目前已初步應用于76子孔徑，109單元變形鏡的自適應光學驗證系統，實驗結果驗證方案可行。

2 采用算法介紹

系統的控制帶寬，是衡量自適應光學系統好壞的主要指標。在系統采樣頻率一定的情況下，電子系統的計算延時會直接影響到自適應光學系統的控制帶寬。因此必須特別地對電子系統采用的算法進行分析和優化。

目前的相關跟蹤算法中，差絕對值和法與協方差相關算法均能使用。本文選用FFT協方差相關算法。

整個電子系統算法應用過程如圖1所示。先確定一個參考圖像進行FFT運算后作為模板，每個實時采集的圖像經過預處理后進行FFT運算，復共軛后與參考圖像FFT點乘，結果再進行逆FFT變換得到相關函數值，找到匹配點后，進行5×5的拋物線擬合，所得到的相對位移與波前復原互作用矩陣進行部分積運算，得到電壓控制量。

圖1 相關算法框圖

具體分析如下：

2.1 圖像預處理

在目標為太陽米粒組織圖像或太陽黑子的情況下，采用圖像相關處理算法，其本質是尋找圖像的最佳匹配點，從而得到圖像的運動估計。對于N×N的兩幅圖像，根據相關的定義：

其中CLR(x,y)為相關函數，IR(x,y)是參考圖像，IL(x,y)是目標活動圖像。目標圖像和參考圖像的重疊度和相似度越大，則相關值越大，因此通過計算相關函數最大值的位置，就可以得到圖像的相對移動量，也就是目標相對于參考圖像的位移。

協方差相關算法來源于標準化協方差相關，就是對相關函數進行去均值和歸一化處理。

2.2 基于FFT協方差相關算法

圖像在預處理之后，就可以進行相關計算來求相對偏移量，但是直接計算運算量過于巨大，必須將互相關函數首先轉化成FFT變換來求解。

即使直接采用FFT運算量也不小，必須對算法進行進一步優化。由于圖像是純實數，可以通過變換處理使得圖像序列變成一個復數序列的實部和虛部，圖像的傅里葉變換是 2維的，2維離散傅里葉變換(DFT)可以利用1維基2FFT來計算優化[9]。相關算法按照步驟主要包括一次2維FFT，一次頻譜點乘，一次2維逆FFT。

2.3 拋物面擬合

相關函數求得的峰值位置，只能計算出位置偏差的整像元數，為了提高計算精度，還需對相關峰進行曲面擬合，求得最接近的亞像元偏移位置。在這里取相關函數最大值(設為C(0,0))附近5×5的子矩陣作拋物面擬合，根據最小二乘法可以推得相移坐標X，Y：

其中xmax，ymax是相關函數最大點的行列坐標。

2.4 直接斜率法[10]計算波前復原矩陣

采用直接斜率法校正波前時，首先在系統開環的情況下，對變形鏡的驅動器分別施加單位控制電壓，并由哈特曼波前傳感器分別測量相應的波前斜率，根據控制電壓V與波前斜率矩陣S間的關系式

求得斜率響應矩陣H的值，H維數為2n×m(n為變形鏡校正器驅動單元的數目，m為哈特曼傳感器的子孔徑數)。在畸變波前進行校正閉環工作時，根據實時動態測量出的畸變波前斜率S，可根據式(7)求得所需的波前復原電壓V

式中H*為H的廣義逆矩陣，亦稱波前復原矩陣。當2n＞m時，它們之間滿足最小二乘關系，能夠保證對各校正驅動單元施加的控制電壓是以最小方差擬合波前探測的相位斜率。本系統中n=76，m=109。

將2.3節得到的偏移量與H*進行部分積運算，結果運用相應PID控制算法就可以求出高階波前校正系統需要的反饋電壓控制量(低階相關跟蹤系統只用偏移量即可計算控制電壓)。

3 系統設計及實現

3.1 系統需求分析

參考第2節的算法，根據系統需求可以估算出系統的計算量，進行電子系統初步設計。

驗證系統采用10×10的哈特曼微透鏡陣列。高階波前處理系統波前探測采用高幀頻CMOS相機，窗口輸出幀頻最高可達2500 Hz。一共100個子孔徑(又稱為子窗口)，有效子孔徑為76個，子孔徑的劃分如圖2所示。根據系統需求，要求高階系統處理能力必須達到20.52 Gflops，完成每個子孔徑計算耗時＜0.04 ms。

低階相關跟蹤器探測相機采用CCD為DALSA系列，輸出32 pixel×32 pixel，輸出幀頻＞5000 Hz，要求處理能力達到4.4 Gflops，完成全部計算在主頻500 MHz的情況下，耗時＜0.18 ms。

圖2 子孔徑陣列排布

3.2 系統整體設計

根據算法特點和系統要求，本文設計了基于DSP與FPGA協同合作的，由低階相關跟蹤系統和高階波前校正系統組成的自適應光學電子系統，其結構如圖3所示。其中，高階波前處理系統其工作流程為：高幀頻CMOS相機采集的高速數字圖像信號通過被采集進入圖像采集分發板，經過圖像分割處理，分割后的圖像數據通過 Link口(單向速度為450 MB/s)傳輸給m塊圖像計算處理板(每塊板4片DSP)進行并行計算，計算的結果再通過 Link口匯總給控制輸出板，再將校正計算的結果輸出給n塊32路輸出的D/A板，進而產生所需控制電壓對變形鏡進行反饋控制。低階相關跟蹤系統工作流程與高階系統類似，但需要在一塊板上完成從圖像采集、計算控制到D/A輸出的所有功能。高階低階兩個子系統均需通過PCI總線，將計算結果、控制數據存儲并傳至控制計算機，在計算機上實時顯示監控圖像。低階高階之間設計有連線，方便兩個子系統交互數據。

圖3 天文圖像信號處理實現方案

3.3 圖像處理部分的設計

圖像處理設計是整個系統設計的核心部分，為實現對所觀測的天文圖像自適應地進行高速、實時的處理，本文采用了基于DSP+FPGA的圖像處理設計，包括圖像采集單元、圖像分割策略以及圖像計算單元的設計。

3.3.1圖像采集單元 在3.1節里面介紹過，電子系統通過CMOS和CCD相機來獲取圖像信息，高階和低階波前探測單元圖像接口分別采用 Camera Link和RS422接口。其中Camera Link接口方式的數據傳輸速率高達660 MB/s，因此在設計中需要利用FPGA內部的高速、大容量Block RAM作為存儲單元，對該方式的數據進行緩存，并以 Gray碼對存儲單元的讀寫地址進行編碼，自行構建高速、大容量的異步FIFO。通過該FIFO的數據緩存后，系統能夠實現通過多路DSP-TS201的 Link口同步傳輸數據給所有圖像計算單元。該方式的優點是，任意的計算單元能夠按時鐘同步獲取相同的數據，處理完畢后同步進行圖像波前校正，避免出現不同的計算單元在同一時刻處理不同幀圖像的情況。

3.3.2圖像分割策略 在波前處理系統中，傳統的策略是處理器獲取一幀完整圖像并儲存，再通過總線尋址來獲取所需數據，這對于低幀頻圖像比較易于實現，對本系統不再適用。對于圖2所示的子孔徑分布，需運用塊處理的策略，將輸入數據分割成塊，發揮DSP并行計算的優勢。在采用流水線處理的方式下，所有的DSP處于任務級同步工作，以實現高速實時處理天文圖像。

為實現這一目的，本文采用基于FPGA的準實時流水線進行圖像分割策略。第1步，按照上文設計，建立第1級FIFO，通過Camera Link接口取得原始圖像，同時產生與數據同步的地址；第2步，由目標DSP產生圖像分割參數表，具體而言就是將所選取分割區域的每個邊界拐點坐標按照“包頭+坐標地址”的格式寫入一個參數表，由Link口發送給 FPGA，表的容量根據 Link口的要求，必須為256的倍數，實際系統中選取為1 kbyte；第3步，FPGA將每一行的數據讀入第2級FIFO，再根據分割參數表得到的有效地址，將所需的數據寫入第3級FIFO，通過Link口發送給對應的DSP進行處理。

3.3.3圖像計算單元 針對整個系統的特點，選用ADI高性能DSP芯片TS201。該芯片具有高處理器時鐘頻率，合理的結構和新的存儲器技術，以及高IO接口帶寬，使得基于該芯片的系統性能提升，大大降低了開發成本、功耗和芯片尺寸。其主要的優點為：提供高性能靜態超標量處理器，專對大的信號處理任務和通信結構進行優化；4條相互獨立的128位寬的內部數據總線連接到片內 6個 4 MbitDRAM，片內帶寬33.6 GB/s，運行在600 MHz時，ADSP-TS201S可以提供48億次40位MAC運算或者12億次80位MAC運算；14通道的DMA控制器支持優先級中斷和嵌套中斷；4個全雙工的Link口支持單向最高500 MB/s的傳輸速度[11]。

圖像計算單元DSP板的核心組成部分是由4片TS201、大量存儲器(128MB的SDRAM)及PCI接口構成的DSP簇。圖4描述了計算版上DSP簇內的Link口連接結構。

圖4 TS201的Link口連接結構

特別地，針對76子孔徑，109變形鏡的實驗系統，根據圖2所示的子孔徑分布，系統計算的瓶頸在中間4-7行每行10個子孔徑的并行計算。可以直接使用DSP板代替控制輸出板，直接采用DSP板上固定的Link口連接，起到了簡化系統的作用。如圖5所示，將76個子孔徑分為左右兩部分，分配給DSP板1和DSP板2，在4-7行需計算一行10個子孔徑時，DSP板3上的DSP0和DSP2也參與計算。這樣分配的結果在不影響系統傳輸時延的情況下，可以最大發揮DSP并行計算的優勢，減少系統復雜度。根據圖像分割策略，以左半邊子孔徑為例，每塊DSP預計承擔的計算任務為：DSP板1中的DSP0-4以及DSP板3中的DSP0需計算子孔徑數目分別為6,10,8,10,4個。

3.4 通用性和擴展性

本系統具備強大的通用性和擴展能力，既能滿足小型系統的需要，又能滿足大型系統的需要。采用的DSP計算板和D/A板數目可以根據實際需求調節。系統改動時，僅需更改FPGA中的圖像分割算法以及Link口的連接，硬件架構并不需要特別改動，具有較高的實用價值。

4 實驗測試結果

圖5 實際數據流圖

圖6為DSP運行在500 MHz，系統調試測試時的示波器截圖。圖 6(a)為計算板 0上的其中兩個DSP工作波形，CH1是DSP0，CH2是DSP2。第1個周期(即第1幀圖像)為參考模板的計算，后面每個周期為活動圖像和參考圖像的相關計算。圖6(b)，6(c)為參考模板的計算，在一幀圖像中，每個子孔徑計算時間約為24 μs(參考模板只需要計算FFT即可)，圖6(d)，6(e)為活動圖像與參考圖像的相關計算，每個子孔徑計算時間約為35 μs。從圖上可以得到周期為 400 μs，反映出系統的控制帶寬為 2500 Hz。

圖6 示波器截圖

每一幀圖像由10個子孔徑行組成，波形圖上矩形脈沖數目的不同反應了計算分配量的不同。從圖上看出，在一幀周期內，DSP0中間有 6個矩形脈沖，DSP1中間有8個矩形脈沖，這表明DSP0只有在第3-8行的子孔徑參與了計算，DSP2則在第2-9行參與了計算。同理可以推斷，DSP1與DSP3承擔了每一行的相關運算，原因是他們幾乎不參與數據在DSP板上的分配和集中，因而可以承擔更多的計算任務。當圖像傳輸處于第4-7行子孔徑時，DSP板3的DSP0和DSP2也需要參加運算。這樣分配的結果符合預期的設計，在采用流水線處理的方式下，所有的DSP處于任務級同步工作，系統總延時完全滿足實時要求。

5 結論

本文針對低對比度擴展目標觀測，設計了基于DSP與FPGA協同架構的自適應光學電子系統，選取FFT協方差相關算法，以并行流水的方式完成對相關函數的實時波前計算，達到了預先設計的目標。該系統具有較強的通用性和擴展能力，可以根據需要對計算單元和 D/A單元進行擴展，硬件實現簡單，具有較高的實用價值。下一步工作是通過光、機、電系統聯調，使該系統應用到太陽望遠鏡中。

致謝感謝戴妍峰博士，王景宇老師在論文工作期間給予的幫助。

[1] 饒長輝, 姜文漢, 張雨東等. 云南天文臺1.2 m望遠鏡61單元自適應光學系統[J]. 量子電子學報, 2006, 23(3): 295-302.Rao Chang-hui, Jiang Wen-han, and Zhang Yu-dong,et al..61-element adaptive optical system for 1.2 m telescope of Yunnan Observatory.Chinese Journal of Quantum Electronics, 2006, 23(3): 295-302.

[2] 鄭文佳, 王春鴻, 姜文漢等. 基于脈動陣列的自適應光學實時波前處理機設計[J]. 光電工程, 2008, 34(5): 44-49.Zheng Wen-jia, Wang Chun-hong, and Jiang Wen-han,et al..Design and analysis of real-time adaptive optics wavefront processor based on systolic array.Opto-Electronic Engineering, 2008, 34(5): 44-49.

[3] 閆光輝, 王春鴻, 黃奎. 基于雙ADSP-TS201的波前信號處理試驗平臺設計[J]. 儀器儀表用戶, 2008, 15(5): 89-91.Yan Guang-hui, Wang Chun-hong, and Huang Kui. The design of wavefront processing test platform based on double ADSP-TS201.Electronic Instrumentation Customer, 2008,15(5): 89-91.

[4] 高越, 趙丹培, 姜志國. 基于小波變換的空間目標圖像去噪方法[J]. 電子器件, 2009, 32(3): 716-720.Gao Yue, Zhao Dan-pei, and Jiang Zhi-guo. A denoising method based on wavelet for image of space.Chinese Journalof Electron Devices, 2009, 32(3): 716-720.

[5] Rao Chang-hui, Jiang Wen-han, and Ling Nin. Tracking algorithm for low contrast extened object[J].Acta Astronmica Sinica, 2001, 42(3): 329-338.

[6] 彭曉峰, 李梅, 饒長輝. 基于絕對差分算法的相關 HS波前處理機設計[J]. 光電工程, 2008, 35(12): 18-22.Peng Xiao-feng, Li Mei, and Rao Chang-hui. Design of correlating Hartmann-Shack wavefront processor based on absolute difference algorithm.Opto-Electronic Engineering,2008, 35(12): 18-22.

[7] Rimmele T, Richards K, and Hegwer S,et al.. First results from the NSO/NJIT solar adaptive optics system[J].SPIE,2004, 5171: 179-186.

[8] Smithson R and Tarbell T D. Correlation tracking study for meter-class solar telescope on space shuttle. Lockheed Palo Alto Res.Lab.1977.

[9] 姜愛民. 相關跟蹤器研制[D]. [博士論文], 中國科學院國家天文臺, 2004.Jiang Ai-min. Studies on correlation tracker[D]. [Ph.D.dissertation], National Astronomical Observatories, Chinese Academy of Sciences, 2004.

[10] Jiang W H and Li H G. Hartmann-Shack wavefront sensing and wavefront control algorithm[J].SPIE, 1990, 1271: 82-93.

[11] 劉書明, 羅勇江. ADSP-TS20XS系列 DSP原理與應用設計[M]. 北京: 電子工業出版社, 2007: 6-9.