基于多角度偏振成像儀數(shù)據(jù)的測試光源變化提取方法研究

葉晨，駱冬根，李怡心，姚萍萍，王羿，向光峰，李揚，李雙，洪津

（1 中國科學院合肥物質(zhì)科學研究院 安徽光學精密機械研究所，合肥 230031）

（2 中國科學技術大學，合肥 230026）

（3 中國科學院通用光學定標與表征技術重點實驗室，合肥 230031）

（4 桂林電子科技大學 信息與通信學院，廣西 桂林 541004）

0 引言

中國科學院安徽光學精密機械研究所研制的星載偏振傳感器多角度偏振成像儀（Directional Polarimetric Camera，DPC）具有低畸變、大視場、畫幅式成像的特性，使其可獲得地氣系統(tǒng)的多角度、多光譜、偏振輻射信息［1-4］，將獲得的偏振信息和大氣特性反演模型相結合，可得到大氣氣溶膠和云的分布狀況及相關應用參數(shù)［5-6］。DPC 載荷在出廠前，需要在實驗室對DPC 響應做穩(wěn)定性能測試［7-8］，實驗要求光源的輻射能盡可能穩(wěn)定。輻射源通常采用鹵鎢燈積分球，周小麗等研制的近紅外光源在穩(wěn)定后10 min 以內(nèi)光強波動約0.1%［9］；袁銀麟等研制的大孔徑可調(diào)光譜積分球光源2 h 內(nèi)的全譜段不穩(wěn)定性約為0.2%［10］；D'AMATO D 闡明了同一鹵鎢燈，其不同波段的相對光譜輻射率不同，且不穩(wěn)定性也不一致［11］。光源使用年限增加導致其性能衰變，結合環(huán)境溫濕度變化和操作規(guī)范性差異等因素影響，光源非穩(wěn)定性［12］參數(shù)會進一步變大。因此有必要對光源進行能量監(jiān)測，以降低多角度偏振成像儀的測量誤差，校正穩(wěn)定性測量參數(shù)［13］。本文通過對多角度偏振成像儀12 個通道的灰度值序列進行小波分解，提取出各通道的光源能量變化量并在數(shù)據(jù)處理中進行校正，以有效提升載荷穩(wěn)定性參數(shù)測量精度。

1 基本原理

1.1 多角度偏振成像儀和陷阱探測器觀測成像模型

搭載在某型號衛(wèi)星上的DPC 在一個探測周期內(nèi)可獲取3 個非偏波段（763 nm、765 nm、910 nm）和3 個偏振波段（490 nm、670 nm、865 nm）共12 個通道的圖像，每個偏振波段設置0°、60°、120°三個檢偏方向。剩余一個本底通道用于測量DPC 的暗電流數(shù)據(jù)，以消除由于探測器自身、電路盒溫度等因素造成的影響。

根據(jù)DPC 探測原理，其圖像數(shù)字量（Digital Number，DN）可表示為

式中，DNDPC為像元的DN 值，GDPC為DPC 的電子增益系數(shù)，λ1、λ2為DPC 探測波段的上下限，Ls(λ)為積分球光源的光譜輻亮度，Apixel為像元面積，F(xiàn)為等效光圈數(shù)，Moptic為光學系統(tǒng)放大倍率，Tatm(λ)為環(huán)境光譜透過率，Toptic-DPC(λ)為DPC 光學系統(tǒng)的光譜透過率，η(λ)為探測器量子效率，tint為積分時間，Ndark為暗電流電子數(shù)總和，Nnoise1為DPC 各種噪聲的總和。因此從式（1）可以看到主要影響DPC 長期穩(wěn)定性的是其中4 個參數(shù)：受積分球輻亮度穩(wěn)定性影響的Ls(λ)、受環(huán)境影響的Τatm(λ)、受探測器溫度影響的暗電流Ndark和探測器的其它噪聲Nnoise1。

陷阱探測器的探測原理可描述為

式中，s為陷阱探測器靈敏度，λ3、λ4為其探測帶寬上下限，對于同一款陷阱探測器來說，λ3、λ4是固定值，MTRAP、ATRAP分別為陷阱探測器光學系統(tǒng)參數(shù)、探測器面積，Noffset、Nnoise2分別為陷阱探測器的固定偏置、噪聲總和。因此陷阱探測器可直接測量獲得光源的穩(wěn)定性。

1.2 載荷非穩(wěn)定性評價

載荷的一項重要指標是像元響應不確定性，即載荷非穩(wěn)定性（Instability，NS），其計算公式為

式中，k為載荷在穩(wěn)定性實驗期間響應的總次數(shù)，y(i)為載荷第i次響應時的輸出信號，為全時間段內(nèi)輸出信號的算術平均值。

1.3 小波分解原理

小波分解可由粗到精地對信號進行多尺度分解，所以能夠很好地提取出信號中的各層低頻成分［14-16］，待分解信號通過低通濾波器得到低頻系數(shù)A，通過高通濾波器得到高頻系數(shù)D，如圖1 所示，信號經(jīng)第一層小波分解后得到A1和D1，然后再將第一層的低頻系數(shù)A1作為第二層的輸入，又得到一組低頻系數(shù)A2和高頻系數(shù)D2。以此類推分解，直到滿足想要的頻率要求。

圖1 信號小波分解過程Fig.1 Signal wavelet decomposition process

待分解的原始信號可看做0 級低頻系數(shù)a0=(f0，f1，...，fn)，那么

式中，G為低通濾波器矩陣，H為高通濾波器矩陣，am、dm分別為當前分解m層時的低頻系數(shù)、高頻系數(shù)，因此

式中，G*、H*為G、H的共軛矩陣。

小波基函數(shù)為

式中，φ(t)為基小波或母小波函數(shù)，經(jīng)過尺度因子a和平移因子b變換后的φa，b(t)統(tǒng)稱小波。

對于離散信號，有

2 各觀測波段光源變化量提取

2.1 實驗數(shù)據(jù)的獲取

以搭載在某型號衛(wèi)星上的多角度偏振儀DPC 穩(wěn)定性測試實驗為例。其穩(wěn)定性實驗測試現(xiàn)場示意圖如圖2 所示，將陷阱探測器和DPC 平行放置于積分球的出光口，控制柜控制積分球的點燈數(shù)和電流，采用直徑為2.5 m 的鹵鎢燈IS2500-1000 型積分球，工作波段是350～2 500 nm，標稱全譜非穩(wěn)定性為0.2%/h，所采用的SI19008 型陷阱探測器的輻射不確定度低于0.65%。先將DPC 和陷阱探測器分別正對積分球輻射源，打開積分球光源及DPC 進行預熱至穩(wěn)定狀態(tài)，設置好相應的儀器參數(shù)，各自記錄響應值，根據(jù)陷阱探測器監(jiān)測的光源穩(wěn)定段所對應的DPC 成像情況，經(jīng)數(shù)據(jù)處理判斷DPC 響應是否穩(wěn)定。

圖2 DPC 穩(wěn)定性實驗示意圖Fig.2 Schematic diagram of DPC stability experiment

DPC 的同一通道在6 595 s 里采了1 024 次，采樣頻率為0.155 Hz，陷阱探測器監(jiān)測頻率是2 Hz，由Nyquist 采樣定律知：DPC 各通道的DN 幀序列能檢測到的最高頻率是0.077 5 Hz，陷阱探測器能識別光源變化最高的頻率是1 Hz。

在穩(wěn)定性實驗中，重點關注的是鹵鎢燈積分球的低頻輻射變化信號［17］，而陷阱探測器在監(jiān)測時引入的噪聲及儀器測量誤差主要是相對高頻噪聲，因此利用小波分解可以從含有高頻噪聲的曲線中提取出低頻的光源緩變信號。同時，在DPC 的CCD 記錄數(shù)字量灰度值DN 值過程中，CCD 的光電響應的噪聲較光源能量變化屬于高頻成分，陷阱探測器的監(jiān)測靈敏度變化等噪聲較光源本身的能量變化也屬于高頻成分，因此光源能量的低頻變化曲線是暗藏在DN 值幀序列和陷阱探測器監(jiān)測曲線中的。又由于陷阱探測器的采樣頻率比CCD 同一通道的采樣頻率高，因此陷阱探測器中蘊含的光源變化信息更多，但是穩(wěn)定性實驗只要求光源能量穩(wěn)定在一定范圍內(nèi)，因此理論上可從DN 值幀序列中提取出光源變化趨勢。

在穩(wěn)定性實驗中，積分球輻射源在DPC 的CCD 上所成的某一幀圖像如圖3 所示，由于觀測位置和條件均保持不變，因此每幀的被輻射區(qū)域位置是固定的，其內(nèi)接正方形尺寸約是160 pixels×160 pixels。考慮用多像元平均法提高數(shù)據(jù)信噪比，因此本次選取的像元監(jiān)測區(qū)域為中央?yún)^(qū)域140 pixels×140 pixels。DPC 采用的是幀轉(zhuǎn)移型面陣CCD，因此圖像的預處理包括扣除當前幀及前后各3 幀本底平均和幀轉(zhuǎn)移效應校正。

圖3 DPC 某一通道接收到的一幀圖像Fig.3 A frame of image received by one channel of DPC

2.2 光源變化曲線的提取

基于像元的光子信號和轉(zhuǎn)換過程中產(chǎn)生的噪聲變換特性，本次小波基選取的是db2 小波，分解層數(shù)為8，從陷阱探測器監(jiān)測曲線中提取出光源的低頻變化趨勢A1～A8，也從各通道的DN 值幀曲線中提取出A1～A8，將兩者之間的低頻曲線做Pearson 相關，取Pearson 相關系數(shù)最大所對應的低頻曲線，認為是各通道的光能量變化趨勢，提取光源變化曲線流程如圖4 所示。Pearson 相關系數(shù)是用來反映兩變量之間相似程度的統(tǒng)計量，衡量兩者之間的線性關聯(lián)性，其對兩組數(shù)據(jù)的異常值較為敏感，即數(shù)據(jù)的異常點會導致兩組數(shù)據(jù)間的Pearson 相關系數(shù)降低。

圖4 提取光源變化曲線流程Fig.4 Process of extracting light source change curve

依次對12 個通道的DN 值幀曲線做8 層小波分解，將每個通道每層的低頻曲線An和陷阱探測器的An計算Pearson 相關系數(shù)。以670 nm（P1）通道為例，此通道的DN 值幀曲線的A1～A8和對應陷阱探測器的A1～A8之間的Pearson 相關系數(shù)如圖5 所示，DPC 的670（P1）、670（P2）、670（P3）之間的Pearson 相關系數(shù)如圖6 所示。

圖5 670 nm（P1）通道的DN 值幀序列的A1～A8和陷阱探測器的A1～A8的相關系數(shù)Fig.5 Correlation coefficient between A1～A8 of trap detector and A1～A8 of DPC 670 nm（P1）

圖6 670 nm（P1）、670 nm（P2）、670 nm（P3）的DN 值幀序列的A1～A8之間的Pearson 相關系數(shù)Fig.6 Pearson correlation coefficient between A1～A8 of DN value frame sequence of 670 nm（P1）、670 nm（P2）、670 nm（P3）channel

由圖5 可看出，對于670 nm（P1）而言，在陷阱探測器的A1～A8曲線中，Pearson 相關系數(shù)最大值均是在DN 值幀序列的A5處取得，且均在0.95 以上，說明670 nm（P1）和光源的670 nm 波段輻射能的變化極度相關。由圖6 可看出670 nm 的三個偏振通道的DN 值幀序列在相同的小波分解層數(shù)下（即對角線方向），其Pearson 相關系數(shù)取得最大值，且在A5時的兩兩之間的Pearson 相關系數(shù)接近0.998，說明三個偏振通道的DN 值幀序列變化趨勢一致，是由光源670 nm 波段的輻射能變化造成的。因此對670 nm 波段而言，其對應的DN 值幀序列的A5最能體現(xiàn)光源此波段的變化趨勢。同理計算得到其余8 個成像通道的光源變化曲線是各自通道DN 值幀序列的A4或A5，說明在整個穩(wěn)定性觀測實驗期間，光源的光輻射能主要以低于0.004 8 Hz的頻率緩慢變化。因此，每個通道的DN 值幀序列小波分解的A4或A5即為實驗期間光源變化量。

以測量的初始時刻的光能量為基準，各通道光源能量變化趨勢和變化量分別如圖7 和表1 所示。光源各波段的能量穩(wěn)定性不一致，其中490 nm、910 nm 和865 nm 波段變化超過0.5%，因此有必要對DPC 的觀測數(shù)據(jù)進行去光源變化處理，以提高穩(wěn)定性參數(shù)的真實性和有效性。

圖7 各通道的光源能量變化Fig.7 Variation of light source energy in all channels

表1 各波段的光能量變化Table 1 Variation of light energy in each wavelength band

利用所提取的各波段光源能量變化情況進行穩(wěn)定性校正，得到利用陷阱探測器監(jiān)測方法和小波分解方法前后的穩(wěn)定性參數(shù)計算結果如表2 所示，可以看到兩種方法獲得的穩(wěn)定性參數(shù)計算結果均有一定程度的提升，而小波分解方法獲得的穩(wěn)定性參數(shù)的提升更為明顯，較好地去除了光源自身的波動、環(huán)境光譜透過率不穩(wěn)定性、DPC 自身溫度變化導致的暗電流不穩(wěn)定性，計算出的不穩(wěn)定性更接近DPC 實際工作時的不穩(wěn)定性。

表2 兩種校正方法的非穩(wěn)定性參數(shù)計算結果對比Table 2 Comparison of calculation results of the instability parameter by two correction methods

3 基于小波分解方法提取穩(wěn)定性參數(shù)的驗證

3.1 實驗室信噪比計算方法

1）多幀圖像數(shù)據(jù)信噪比計算法

CCD 同一通道接收到了k幀圖像，在被輻照區(qū)域內(nèi)選取計算區(qū)域m×n，對于計算區(qū)域內(nèi)的某一個像元(i，j)，在k幀圖像中，該像素點對應k個DN 值，求其平均值μ(i，j)及標準差std(i，j)，定義該像元的信噪比為

以此類推，可求得計算區(qū)域內(nèi)其它像元的信噪比。將m×n個像元的信噪比SNR(i，j)取平均，作為CCD此通道接收到的圖像的信噪比SNR，即

2）相鄰兩幀信噪比計算法

歐洲機器視覺協(xié)會制定的圖像傳感器與相機性能測試標準EMVA Standard 1288［18］中提出用相鄰兩幀計算信噪比，取CCD 相鄰兩幀圖像A、B，假設噪聲源靜止且具有均勻性，則在被輻照區(qū)域內(nèi)選取計算區(qū)域m×n，定義平均灰度值即信號為

式中，yA(i，j)表示圖像A中坐標(i，j)處的像元灰度值，yB(i，j)表示圖像B中坐標(i，j)處的像元灰度值。

兩圖像A、B相減，得到新圖像C，其m×n個像元灰度值的標準差記為stdC，則原圖像A或B的噪聲大小視為

3.2 陷阱探測器修正光源計算信噪比

在穩(wěn)定性實驗中，常用陷阱探測器監(jiān)測光源變化曲線來對DPC 所成的圖像進行去光源變化，將當前圖像的建立時間與所對應的光源曲線的點對應起來，然后再去光源變化，即可計算信噪比。設某一通道的第1幀、第2 幀、...、第k幀圖像建立時對應的陷阱探測器監(jiān)測的光輻射能依次為Q1、Q2、...、Qk，將其光能量歸一化處理，則此通道第i幀的光能量修正系數(shù)為

第i幀圖像修正后的DNTRAP_corrected_i值為

將每幀圖像根據(jù)陷阱探測器的光能量進行DN 值修正后，計算出的信噪比記為陷阱探測器修正后信噪比。

3.3 基于小波分解法修正光源計算信噪比

為證明提出的用小波分解法提取光源的有效性，用基于小波分解法提取出12 個成像通道各自的光源變化曲線，依據(jù)光源變化修正圖像DN 值，得到修正后的幀序列信噪比以及信噪比提升比例分別見圖8 和9，具體提升的信噪比數(shù)值比較見表3。由圖8 可看出，小波分解法修正光源后計算出的信噪比接近相鄰兩幀信噪比，證明提取光源方法的正確性。

表3 陷阱探測器去光源變化法和小波分解法去光源變化法提升的信噪比數(shù)值比較Table 3 Numerical comparison of the improved SNR between trap detector de-illumination variation method and wavelet decomposition de-illumination variation method

圖8 不同的方法計算的信噪比Fig.8 SNR calculated by different methods

陷阱探測器去光源修正后較未修正前信噪比平均提升了0.72%，而小波分解法修正后信噪比平均提升了2.52%。在865 nm 和910 nm 波段，小波分解去光源法將信噪比分別提升了4.34%和12.9%，遠高于陷阱探測器監(jiān)測數(shù)據(jù)修正的結果，反映了實驗期間光源在865 nm 波段和910 nm 波段能量波動較大，也說明用寬譜的陷阱探測器監(jiān)測數(shù)據(jù)來修正865 nm 和910 nm 波段的變化存在局限。在490 nm 波段，小波分解法將信噪比提升了0.54%，但陷阱探測器修正法僅將信噪比提升了0.09%，主要是因為490 nm 波段的光源能量較低，獲取的圖像DN 值低。雖然光源在490 nm 波段的變化最大，但信噪比提升并沒有910 nm 波段顯著。兩種方法在670 nm 波段信噪比提升幅度小則是因為此波段的光源變化小于0.205%。

圖9 陷阱探測器去光源變化法和小波分解法去光源變化法提升的信噪比Fig.9 The improved SNR of the trap detector de-illumination variation method and the wavelet decomposition deillumination variation method

4 結論

積分球光源不同光譜波段的變化趨勢不同，且與光源的使用年限、環(huán)境控制、操作規(guī)范等因素有一定的相關性，但總的來說，大口徑積分球光源的變化趨勢較緩慢。本文提出用小波分解法從DPC 各成像通道的DN 值幀曲線中提取出光源各波段的能量變化趨勢，進而扣除光源變化趨勢，使穩(wěn)定性測試的結果更為客觀和真實。光源變化校正的效果優(yōu)于傳統(tǒng)陷阱探測器監(jiān)測數(shù)據(jù)校正，DPC 的信噪比較未扣除光源波動前平均提升了2.52%，非穩(wěn)定性降低到0.031%，驗證了此方法提取光源變化量的有效性。