全偏振大氣偏振模式成像系統的設計與優化分析*

王成 范之國 金海紅3)? 汪先球 華豆

1) (合肥工業大學計算機與信息學院, 合肥 230601)

2) (合肥工業大學, 教育部大數據知識工程重點實驗室, 合肥 230601)

3) (安徽建筑大學電子與信息工程學院, 合肥 230601)

全偏振成像能夠獲取目標更為豐富的信息, 在目標探測、大氣特性研究和醫學診斷等領域具有廣闊的應用前景.為了實現大視場天空區域全偏振信息的快速獲取, 設計了一套全偏振大氣偏振模式成像系統; 針對因系統傳輸矩陣“性態”的不同使得求解的目標Stokes矢量存在誤差的問題, 通過分析傳輸矩陣的特性并建立目標函數, 將傳輸矩陣的優化轉化為目標函數在多組條件下的求解, 確定了最優系統傳輸矩陣; 并對系統的四分之一波片的延遲量、偏振片的消光比以及傳輸矩陣進行標定.通過開展優化前后偏振信息的對比實驗, 結果表明: 優化后偏振角誤差較優化前降低了10%以上; 偏振度和線偏振度中最大偏振度帶的誤差和中性區域的誤差較優化前也有不同程度的下降.在此基礎上開展了外場全偏振信息測量實驗, 結果表明系統滿足設計要求, 能夠有效地獲取天空全偏振信息.

1 引 言

大氣偏振模式是太陽光在大氣傳輸過程中由于大氣的散射、輻射和吸收等作用而產生的偏振光所形成動態地、穩定地偏振態分布, 是地球重要的自然屬性之一[1].自從Arago于1809年首次發現天空存在偏振現象以來, 很多學者對大氣偏振模式開展了大量研究, 其中動物實驗研究表明諸如沙蟻、蜜蜂、蝴蝶、蝗蟲等[2?5]動物能利用大氣偏振模式進行導航和定向.因此, 大氣偏振模式的研究對動物導航行為、大氣特性、地表環境特征以及仿生偏振光導航[6?11]具有重要的理論與應用價值.

在進行大氣偏振模式的研究中, 國內外許多學者設計了多種天空偏振成像系統.如2002年Horváth等[12]開發一套三相機全天空偏振測試儀,并在紅、綠、藍等波段下進行了多云天氣下的天空偏振信息獲取實驗.2006年, Pust和Shaw[13]利用LCVR設計了全天域觀測設備, 開展了晴天和多云天氣下偏振測量實驗.2016年孫潔等[14]設計了一套偏振測量裝置, 并開展了多波段、不同天氣以及不同地表環境下的大氣偏振模式獲取實驗.雖然天空中的圓偏振分量較少, 但學者們發現圓偏振比線偏振具有更好的特性, 比如具有較強的后項散射保持能力[15], 因此全偏振測量系統變得越來越重要.2014年Hsu等[16]研制了一套全Stokes偏振成像儀, 成像儀獲取的線偏振度和圓偏振度的誤差小于7%, 標準差小于5%.同年張忠順[17]設計了全偏振測量成像系統, 開展了全天域外場實驗并研究了中性點的變化.2019年殷玉龍等[18]研制一套全Stokes同時偏振成像系統并進行標定, 研究了1/2波片和1/4波片相位延遲對系統測量誤差的影響.

為了在獲取全偏振信息的同時提高系統獲取全偏振信息的精度, 本文在以往學者研究的基礎上使用魚眼鏡頭、中繼透鏡、濾波片、四分之一波片、偏振片以及微距鏡頭等研制一套全偏振大氣偏振模式成像系統, 并針對因系統傳輸矩陣“性態”的不同而使得系統獲取的光強圖的誤差對目標Stokes矢量求解帶來誤差的問題進行了研究.本文對傳輸矩陣的條件數、行列式等進行分析, 建立目標函數將傳輸矩陣的優化問題轉化為多條件極值問題, 獲得兩組傳輸矩陣的最優角度.為了說明優化效果和系統的有效性, 還開展了優化前后對比實驗與外場全偏振大氣偏振模式獲取實驗.

2 全偏振測量原理與系統設計

全偏振獲取方法主要有兩類: 1) 分振幅法; 2)偏振光調制法.分振幅法是使用分光棱鏡分割入射光的光強, 再經過偏振片或四分之一波片(quarter wave plate, QWP)后由多個CCD同時成像的方法.該方法不需旋轉器件, 但系統結構復雜, 還涉及CCD響應一致性配準的問題.偏振光調制法包含“LCVR + 偏振片”[13]和“127°波片 + 偏振片”[19,20]以及“QWP + 偏振片”[17]三種方式.在綜合考慮系統復雜性和變量的個數的基礎上, 選擇“QWP + 偏振片”方式.

2.1 全偏振測量原理

根據對上述測量方法的分析, 本系統采用的測量原理如圖1所示.

圖1 全偏振測量原理圖Fig.1.Principle diagram of full polarization measurement.

偏振光經過光學器件后發生的變化可以使用Mueller矩陣來描述.假設入射光的Stokes矢量為Sin=[s0,s1,s2,s3] , QWP的Mueller矩陣為MQWP,偏振片的Muller矩陣為MPol, 所以CCD的輸出Iout為

其中T表示對矩陣或向量求轉置;b表示QWP的Fast軸與參考軸的夾角.將MQWP(β) 和MPol帶入(1)式有

當QWP分別旋轉β0,β1,β2,β3時, 得到如下方程組:

因此, 當傳輸矩陣Mtran可逆, 就可以通過(3)式反解出Sin, 所以目標偏振光的偏振角(angle of polarization, Aop)、偏振度(degree of polarization,Dop)、線偏振度(degree of linear polarization,Dolp)和圓偏振度(degree of circular polarization,Docp)等信息可通過(4)式解得

2.2 全偏振測量系統設計

系統的光路傳輸示意圖如圖2所示, 系統由魚眼鏡頭、中繼透鏡、濾波片、QWP、偏振片、微距鏡頭和CCD組成.魚眼鏡頭使得系統具有較大的視場角, 通過使用文獻[17]中的測量方法, 得到系統的視場角為142°; 中繼透鏡則是對魚眼鏡頭出射的發散光進行準直; 微距鏡頭是為了解決系統光斑較小的問題.此外, 由于可見光范圍內488 nm的光波輻射能力強、衰減較少[17], 所以濾波片的中心波長與QWP一致且均為488 nm, 且濾波片的半高寬為(1 ± 0.2) nm.為了實現全偏振信息的快速獲取, 采用壓電諧振電機來控制QWP的旋轉.

圖2 成像系統光路示意圖Fig.2.Schematic diagram of the optical path of the imaging system.

當系統入射光Stokes矢量為Sin=[s0,s1,s2,s3]時 系統的輸出Iout為

其中MFishlens,MLens,MFilter,MMacro分別表示魚眼鏡頭、中繼透鏡、濾波片、微距鏡頭的Mueller矩陣.當QWP旋轉到任意四個不同角度β0,β1,β2,β3后(3)式變為

其 中 [m11(β),m12(β),m13(β),m14(β)] 為QWP旋轉到任意角度b時Msystem的第一行元素, 所以當Mtran可逆時可以解算出Sin, 如(7)式所示:

其中||···||表示矩陣或向量的范數; 由條件數(condition number, Cond)的定義可知:由(9)式可以看出: 當

其中nij(i,j=1,···,4) 為Mtran的逆矩陣中的元素.將(7)式中的Sin帶入(4)式即可求得Aop, Dop,Dolp和Docp.Iout因噪聲產生擾動 ?Iout時,Sin的相對誤差不超過Iout的相對誤差乘以 C ond(Mtran).由于 ?Iout是受噪聲影響的不確定的量, 為了使Sin的相對誤差盡量小, 則 C ond(Mtran) 應該盡量的小.由于條件數分析法常用于偏振儀的優化中, 所以本文使用條件數優化的本質是Mtran中b的選擇問題, 所以b的選擇依據如(10)式所示:

3 系統誤差分析、優化與標定

來優化傳輸矩陣[22].通過2.2的分析可以知道系統

3.1 系統誤差分析與優化

在大氣偏振模式的獲取過程中會受到噪聲的影響, 假設噪聲對Iout產生的擾動為 ?Iout, 此時 ?Iout對Sin解算產生的擾動為 ?Sin, 所以(6)式變為

此外在選擇b時還要保證Mtran是可逆的, 而且選擇的b要使得Mtran的行列式(Determinant,Det)盡量的大.由于Det表示系統對輸入信息的縮放程度, 所以選擇 D et>0 條件下對應的b.

所以系統解算出的Sin的相對誤差可以用(9)式表示[21]:

圖3 M tran 的 1 /Cond 與Det變化趨勢 (a) 1 /Cond1 變化趨勢; (b) 1 /Cond2 變化趨勢; (c) 1 /Cond∞ 變化趨勢; (d) D et 變化趨勢Fig.3.The change trend of 1 /Cond and Det of the M tran : (a) Change trend of 1 /Cond1 ; (b) change trend of 1 /Cond2 ; (c) change trend of 1 /Cond∞ ; (d) change trend of Det.

經過分析發現對b以2o, 5o和10o等為間隔進行遍歷時, 優化得到的角度均在某些固定的角度左右.因此, 結合數據量和運算時間選擇5o為間隔對b在0o—180o之間進行遍歷, 獲得36個 1×4 的行向量, 然后求解所有行向量的組合數, 得到58905個Mtran, 并求解相應矩陣的1范數Cond( C ond1)、2范數Cond( C ond2)、∞范數Cond( C ond∞)和Det.由于矩陣的Cond過大, 所以在圖3(a)—(c)中縱坐標為Cond的倒數, 結果如圖3所示.

從圖3中可以看出, 三個Cond的變化趨勢是類似的, 為了直觀的體現Cond, Det與b的對應關 系, 求解了最小 C ond1, C ond2, C ond∞和 最 大Det對應的b, 結果如表1所列.

根據上述分析, 從表1中可以得到(5o, 45o,120o, 155o), (0o, 30o, 115o, 150o)和(0o, 40°, 115o,150o)三組角度.因此, 在上述三組約束條件下求最優角度的問題就轉化為多條件極值問題, 其中目標函數為L(β0,β1,β2,β3) 的定義如下式所示:

當L(β0,β1,β2,β3) 在約束條件(s.t.)下取得最小值時, (β0,β1,β2,β3) 即為所求角度, 從表1得到約束條件(s.t.)如(12)式所示:

由(11)式和(12)式計算得到,L(β0,β1,β2,β3)在s.t.下的最優值為(5o, 40o, 120o, 155o)和(10o,45o, 125o, 160o), 優化前后對比實驗采用的角度為(10o, 45o, 125o, 160o).

3.2 系統標定

Mtran的標定主要有兩部分: 1) QWP的延遲量和偏振片消光比的標定; 2) 不同b下的Mtran的標定.由于第一部分標定對Mtran的影響是確定的,所以本文標定的重點是第二部分.經過標定, 標定后QWP延遲量的誤差小于 ± 0.5°, 偏振片消光比的誤差小于5%.第二部分的標定原理如圖4所示, 標定系統包括光源、偏振光調制器、Stokes矢量儀和光功率計.

表1 最大Det、最小Cond與b對應表Table 1.Corresponding table of maximum Det,minimum Cond and b.

圖4 系統標定原理框圖Fig.4.Block diagram of system calibration principle.

本文在定標時采用的是主光線定標, 定標的光源采用的是激光光源, 主要是考慮到激光光源具有一致性較好、輸出功率較強且衰減較少等優點, 激光發生器采用Thorlabs公司的高功率LDC4020型激光發生器.具體定標步驟如下: 1) 選擇488 nm波段下的激光光源, 通過偏振光調制器產生不同偏振態的入射光, 使用Stokes矢量儀記錄待測系統入射光的偏振態, 并使用光功率計記錄待測系統的入射光強值; 2) 把b依次設置為5o, 10o, 40o, 45o,120o, 125o, 155o, 160°, 使用光功率計記錄每個角度下待測系統的出射光強值; 3) 改變入射光的偏振態, 重復(2)得到10組偏振態下待測系統的入射光和每組入射光下8個b對應的輸出光強值.假設10組輸入光強為Iin(i) , 歸一化的Stokes矢量為Sin(i) , 輸出光強為Iβ(i) , 則Iin(i) ,Sin(i) 和Iβ(i) 滿足的關系如下式所示:

由于(13)式有32個變量, 80個方程, 故該方程組為超定方程組, 因此使用最小二乘法進行求解, 得到系統傳輸矩陣的標定結果如表2所列.

表2 系統傳輸矩陣標定結果Table 2.Calibration results of system transmission matrix.

4 實驗與結果分析

為了對優化和標定結果進行說明, 開展了兩項實驗: 1) 優化前后偏振信息對比實驗; 2) 外場全偏振信息獲取實驗.實驗地點為合肥工業大學翡翠湖校區操場(E117.214o, N31.778o), 實驗時間為2019年12月13日上午9點, 此時太陽的高度角為29.66o, 方位角為152.1o, 天氣晴朗.

4.1 優化前后對比實驗

在對比實驗中, 優化前QWP的旋轉角度為(0o, 45o, 60o, 120o), 優化后QWP的旋轉角度為(10o, 45o, 125o, 160o°), 本組實驗對Aop, Dop和Dolp進行分析, 結果如圖5所示.為了量化優化效果, 本文選用熵、均值、標準差等指標進行說明.

由于Aop在太陽子午線兩側數值的極性是相反的, 因此無法直接使用熵、均值、標準差分析Aop數據的離散度, 但Aop的絕對值關于太陽子午線對稱, 所以將Aop對稱軸兩側數據的絕對值之差(?Aop, 如圖5(b)和圖5(f)所示)作為分析數據, ?Aop的熵、均值、標準差的求解結果如表3所列.

圖5 優化前后偏振信息對比 (a)?(d) 優化前Aop, ?Aop, Dop和Dolp; (e)?(h) 優化后Aop, ?Aop, Dop和DolpFig.5.Comparisons of polarization information before and after optimization: (a)?(c) Aop, ?Aop, Dop and Dolp before optimization; (d)?(f) Aop, ?Aop, Dop and Dolp after optimization.

表3 優化前后?Aop數據離散度對比Table 3.Comparisons of data dispersion of the ?Aop before and after optimization.

從表3中可以進一步得到優化前后?Aop的熵、均值、標準差分別降低了10.04%, 13.96%和10.81%, 所以優化后?Aop的三個數據離散度指標下降均在10%以上, 同時從圖5(b)和圖5(f)看出噪聲對Aop的對稱性附近的數據有一定的影響,區別主要在反太陽位置處, 但整體保持對稱.

對于優化前后Dop和Dolp的對比, 主要分析Dop和Dolp的最大偏振度帶(band of maximum polarization, BMP)[23]和中性點區域(neutral zone, NZ), 因為這兩部分是Dop和Dolp圖中的固有特征.由于Dop和Dolp的數據均分布在0—1之間, BMP為偏振度最大的區域, 所以BMP的誤差分析以BMP中的數據與1的差( ? DopBMP和?DolpBMP)作為分析數據; 而NZ為偏振度最小的區域, 所以NZ的誤差分析以NZ中的數據與0的差( ? DopNZ和 ? DolpNZ)作為分析數據, 如下式所示:

為了保證數據的一致性, 先從Dop和Dolp圖中裁剪 相同區域的 ? DopBMP和 ? DolpBMP以及相同區域 的 ? DopNZ和 ? DolpNZ, 然后求解上 述區域的熵、均值、標準差, 裁剪結果如圖6所示.

對圖6中BMP區域的熵、均值、標準差的求解結果如表4所列.

圖6 優化前后Dop, Dolp的BMP和NZ對比 (a)?(d) 優化前Dop, Dolp的BMP和NZ; (e)?(h) 優化后Dop,Dolp的BMP和NZFig.6.Comparisons of BMP and NZ of Dop and Dolp before and after optimization: (a)?(d) BMP and NZ of Dop and Dolp before optimization; (e)?(h) BMP and NZ of Dop and Dolp after optimization.

表4 優化前后BMP的數據離散度指標對比Table 4.Comparisons of data dispersion index of BMP before and after optimization.

從表4中可以得到優化前后 ? DopBMP的熵、均值、標準差分別降低了4.964%, 3.255%和24.537%;而 ? DolpBMP的熵、均值、標準差分別降低了5.084%,3.357%和25.058%.從Dop和Dolp的圖中可以發現, 優化之后的Dop和Dolp的BMP較優化前明顯變大, 如圖5(c)、圖5(d)、圖5(g)和圖5(h)中的紅色圓圈所示.

對于NZ的分析, 使用熵、均值、標準差以及面積比(probility of area, Poa)等指標.其中由于NZ的值較小, 所以在計算Poa時, 將閾值設置為0.02.因此, 當Dop和Dolp的值小于0.02時, 則Poa為中性區域; 否則, 為非中性區域.對圖6中的NZ的熵、均值、標準差、Poa對比結果如表5所列.

從表5可以進一步得到優化前后 ? DopNZ的熵、均值、標準差分別降低了1.423%, 7.414%和9.4%, Poa提高了10.663%; 而 ? DolpNZ的熵、均值、標準差分別降低了1.287%, 7.177%和8.938%,Poa提高了9.532%.NZ的Poa指標說明優化后的NZ要比優化前更為集中.

4.2 外場實驗

經過計算得到目標天空區域偏振光的I,Q,U和V四個分量以及Dop, Aop, Dolp和Docp分布結果如圖7所示.圖7中除Aop, Dop以及Dolp外均做了歸一化處理.其中,I分量的數據歸一化到0—1之間,Q,U和V歸一化到–1—1之間.從圖7(b)可以看出, Aop圖的分布是符合“∞”字型分布的, 這與大氣偏振模型的Rayleigh仿真結果是相吻合的, 并且Aop圖的對稱性很好.從圖7(c)和圖7(d)可以發現Dop和Dolp在分布模式上非常相似, 通過計算得到Dop的均值為0.0914,Dolp的均值為0.0911.在I分量的分布模式中, 靠近太陽位置區域值較大, 遠離太陽位置值較小.

表5 優化前后NZ數據離散度指標對比Table 5.Comparisons of data dispersion index of NZ before and after optimization.

圖7 目標天空區域光強圖與偏振模式分布結果 (a) 偏振光強圖; (b) Aop; (c) Dop; (d) Dolp; (e) Docp; (f) I分量圖; (g) Q分量圖; (h) U分量圖; (i) V分量圖Fig.7.The light intensity map and polarization mode distribution results of the target sky area: (a) Polarization intensity diagram;(b) Aop; (c) Dop; (d) Dolp; (e) Docp; (f) I component diagram; (g) Q component diagram; (h) U component diagram; (i) V component diagram.

值得注意的是, 目標天空區域存在Docp, 但是此時數值較小并且有正有負, 分布范圍在–0.04—0.06之間, Docp取值較大的區域位于BMP附近.對比圖7(e)和圖7(i)可以發現:V分量和Docp的分布模式類似, 其數值也較小同樣有正有負, 區別在于V分量中的負值區域占比接近100%, 而Docp中的負值區域的占比為41.89%, 同時發現Docp的分布和V分量的分布類似.

5 結 論

本文使用魚眼鏡頭、中繼透鏡、濾波片、四分之一波片、偏振片以及微距鏡頭研制了一套全偏振大氣偏振模式成像系統, 該系統具有142°的視場角且滿足全偏振信息快速獲取的要求; 同時, 針對系統的Mtran因存在“病態”情況而易受噪聲影響的問題, 從Cond和Det的角度對Mtran的進行分析,得到了最小Cond和最大Det條件下對應的旋轉角度, 并建立了多條件極值目標函數, 將旋轉角度的求解問題轉化為多條件極值問題, 得到最優的四分之一波片旋轉角度.此外, 還對Mtran、四分之一波片的延遲量、偏振片的消光比進行了標定, 并開展了優化前后的對比實驗和全偏振大氣偏振模式測量實驗.實驗結果表明, 優化后系統獲取的偏振信息較優化前有了較好的提升, 同時全偏振實驗的獲取結果也說明天空確實存在全偏振信息.