Design and optimization analysis of imaging system of polarized skylight pattern of full polarization
-
摘要: 全偏振成像能够获取目标更为丰富的信息, 在目标探测、大气特性研究和医学诊断等领域具有广阔的应用前景. 为了实现大视场天空区域全偏振信息的快速获取, 设计了一套全偏振大气偏振模式成像系统; 针对因系统传输矩阵“性态”的不同使得求解的目标Stokes矢量存在误差的问题, 通过分析传输矩阵的特性并建立目标函数, 将传输矩阵的优化转化为目标函数在多组条件下的求解, 确定了最优系统传输矩阵; 并对系统的四分之一波片的延迟量、偏振片的消光比以及传输矩阵进行标定. 通过开展优化前后偏振信息的对比实验, 结果表明: 优化后偏振角误差较优化前降低了10%以上; 偏振度和线偏振度中最大偏振度带的误差和中性区域的误差较优化前也有不同程度的下降. 在此基础上开展了外场全偏振信息测量实验, 结果表明系统满足设计要求, 能够有效地获取天空全偏振信息.Abstract: Full polarization imaging can obtain more information about target, which has a broad application prospect in the target detection, researches of atmospheric characteristics, and medical diagnosis. This paper develops an imaging system of polarized skylight pattern of full polarization for obtaining the information about full polarization rapidly. Meanwhile, aiming at the problem that the error of the light intensity image obtained by the system due to the different “behavior” of the system transmission matrix is brought into the solution of the target Stokes vector, this paper analyzes the condition number and determinant of the system transmission matrix. Firstly, an objective function is established by combining the three sets of condition numbers and the determinant. Therefore, the problem of solving the optimal transmission matrix is transformed into a multi-condition extremal problem. And then the objective function is minimized to determine the optimal angle of the transmission matrix when the 1 norm condition number, 2 norm condition number and ∞ norm condition number reach the minimum value and the determinant reaches the maximum value. In addition, in order to improve the measurement accuracy, the delay components of quarter wave plate, extinction ratio of polarizer, and the transmission matrix of the system are calibrated. Optimization contrast experiment and outfield experiment are performed. The entropy, mean, and standard deviation are used to quantify the optimized results of the angle of polarization, degree of polarization, and degree of linear polarization. ∆Aop is defined as the difference in absolute value of angle of polarization between the two sides of the symmetry axis to verify the optimization performance of angle of polarization. Experimental results show that the polarization angle error after optimization is reduced by more than 10% compared with that before optimization; the error of the band of maximum polarization and the error of the neutral zone in the degree of polarization and linear polarization also decline to different degrees compared with before optimization. On this basis, an experiment on measuring external field full polarization information is carried out. The results show that the system meets the design requirements and can effectively obtain the sky full polarization information.
-
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所示.
偏振光经过光学器件后发生的变化可以使用Mueller矩阵来描述. 假设入射光的Stokes矢量为
${{{S}}_{{\rm{in}}}} = [{s_0}, {s_1}, {s_2}, {s_3}]$ , QWP的Mueller矩阵为${{{M}}_{{\rm{QWP}}}}$ , 偏振片的Muller矩阵为${{{M}}_{{\rm{Pol}}}}$ , 所以CCD的输出${I_{{\rm{out}}}}$ 为$${I_{{\rm{out}}}}(\beta ) = {{{M}}_{{\rm{Pol}}}} \times {{{M}}_{{\rm{QWP}}}}(\beta ) \times {{S}}_{{\rm{in}}}^{\rm{T}},$$ 1 其中T表示对矩阵或向量求转置; β表示QWP的Fast轴与参考轴的夹角. 将
${{{M}}_{{\rm{QWP}}}}(\beta )$ 和${{{M}}_{{\rm{Pol}}}}$ 带入(1)式有$${I_{{\rm{out}}}}(\beta ) = \frac{1}{2}\left[ {{s_0} + {s_1}{{\cos }^2}2\beta + \frac{1}{2}{s_2}\sin 4\beta - {s_3}\sin 2\beta } \right].$$ 2 当QWP分别旋转
${\beta _0}$ ,${\beta _1}$ ,${\beta _2}$ ,${\beta _3}$ 时, 得到如下方程组:$$\left[ {\begin{array}{*{20}{c}} {{I_{{\rm{out}}}}({\beta _0})} \\ {{I_{{\rm{out}}}}({\beta _1})} \\ {{I_{{\rm{out}}}}({\beta _2})} \\ {{I_{{\rm{out}}}}({\beta _3})} \end{array}} \right] = \frac{1}{2}\left[ {\begin{array}{*{20}{c}} 1&{{{\cos }^2}2{\beta _0}}&{0.5\sin 4{\beta _0}}&{ - \sin2{\beta _0}} \\ 1&{{{\cos }^2}2{\beta _1}}&{0.5\sin 4{\beta _1}}&{ - \sin2{\beta _1}} \\ 1&{{{\cos }^2}2{\beta _2}}&{0.5\sin 4{\beta _2}}&{ - \sin2{\beta _2}} \\ 1&{{{\cos }^2}2{\beta _3}}&{0.5\sin 4{\beta _3}}&{ - \sin2{\beta _3}} \end{array}} \right] \times \left[ {\begin{array}{*{20}{c}} {{s_0}} \\ {{s_1}} \\ {{s_2}} \\ {{s_3}} \end{array}} \right] = {{{M}}_{{\rm{tran}}}} \times {{S}}_{{\rm{in}}}^{\rm{T}}.$$ 3 因此, 当传输矩阵
${{{M}}_{\rm{tran}}}$ 可逆, 就可以通过(3)式反解出${{{S}}_{{\rm{in}}}}$ , 所以目标偏振光的偏振角(angle of polarization, Aop)、偏振度(degree of polarization, Dop)、线偏振度(degree of linear polarization, Dolp)和圆偏振度(degree of circular polarization, Docp)等信息可通过(4)式解得$$\left\{ {\begin{aligned} & {{\rm{Aop}} = \frac{1}{2}\arctan \left( {\frac{{{s_2}}}{{{s_1}}}} \right){{ ,}}} \\ & {{\rm{Dop}} = \frac1{{{s_0}}}{{\sqrt {{s_1^2} + {s_2^2} + {s_3^2}} }}{{,}}} \end{aligned}} \right.\;\;\;\;\;\left\{ {\begin{aligned} & {{\rm{Dolp}} = \frac1{{{s_0}}}{{\sqrt {{s_1^2} + {s_2^2}} }}{{ ,}}} \\ & {{\rm{Docp}} = \frac{{{s_3}}}{{{s_0}}} {{.}}} \end{aligned}} \right.$$ 4 2.2 全偏振测量系统设计
系统的光路传输示意图如图2所示, 系统由鱼眼镜头、中继透镜、滤波片、QWP、偏振片、微距镜头和CCD组成. 鱼眼镜头使得系统具有较大的视场角, 通过使用文献[17]中的测量方法, 得到系统的视场角为142°; 中继透镜则是对鱼眼镜头出射的发散光进行准直; 微距镜头是为了解决系统光斑较小的问题. 此外, 由于可见光范围内488 nm的光波辐射能力强、衰减较少[17], 所以滤波片的中心波长与QWP一致且均为488 nm, 且滤波片的半高宽为(1 ± 0.2) nm. 为了实现全偏振信息的快速获取, 采用压电谐振电机来控制QWP的旋转.
当系统入射光Stokes矢量为
${{S}_{\rm{in}}} = [{s_0}, {s_1}, {s_2}, $ $ {s_3}]$ 时系统的输出${I_{{\rm{out}}}}$ 为$$\begin{split} {I_{{\rm{out}}}} =\;& {{{M}}_{\rm{Macro}}} \times {{{M}}_{\rm{Pol}}} \times {{{M}}_{\rm{QWP}}}(\beta ) \times {{{M}}_{\rm{Filter}}} \times {{{M}}_{{\rm{Lens}}}} \times \\ &{{{M}}_{{\rm{Fishlens}}}} \times {{S}}_{{\rm{in}}}^{\rm{T}} = {{{M}}_{{\rm{System}}}} \times {{S}}_{{\rm{in}}}^{\rm{T}},\\[-10pt]\end{split}$$ 5 其中
${{{M}}_{{\rm{Fishlens}}}}, {{{M}}_{{\rm{Lens}}}}, {{{M}}_{{\rm{Filter}}}}, {{{M}}_{{\rm{Macro}}}}$ 分别表示鱼眼镜头、中继透镜、滤波片、微距镜头的Mueller矩阵. 当QWP旋转到任意四个不同角度${\beta _0}$ ,${\beta _1}$ ,${\beta _2}$ ,${\beta _3}$ 后(3)式变为$$\left[ {\begin{array}{*{20}{c}} {{I_{{\rm{out}}}}({\beta _0})} \\ {{I_{{\rm{out}}}}({\beta _1})} \\ {{I_{{\rm{out}}}}({\beta _2})} \\ {{I_{{\rm{out}}}}({\beta _3})} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{m_{11}}({\beta _0})}&{{m_{12}}({\beta _0})}&{{m_{13}}({\beta _0})}&{{m_{14}}({\beta _0})} \\ {{m_{11}}({\beta _1})}&{{m_{12}}({\beta _1})}&{{m_{13}}({\beta _1})}&{{m_{14}}({\beta _1})} \\ {{m_{11}}({\beta _2})}&{{m_{12}}({\beta _2})}&{{m_{13}}({\beta _2})}&{{m_{14}}({\beta _2})} \\ {{m_{11}}({\beta _3})}&{{m_{12}}({\beta _3})}&{{m_{13}}({\beta _3})}&{{m_{14}}({\beta _3})} \end{array}} \right] \times \left[ {\begin{array}{*{20}{c}} {{s_0}} \\ {{s_1}} \\ {{s_2}} \\ {{s_3}} \end{array}} \right] = {{{M}}_{{\rm{tran}}}} \times {{S}}_{{\rm{in}}}^{\rm{T}},$$ 6 其中
$[{m_{11}}(\beta ), {m_{12}}(\beta ), {m_{13}}(\beta ), {m_{14}}(\beta )]$ 为QWP旋转到任意角度β时${{{M}}_{{\rm{system}}}}$ 的第一行元素, 所以当${{{M}}_{{\rm{tran}}}}$ 可逆时可以解算出${{{S}}_{\rm{in}}}$ , 如(7)式所示:$$\left[ {\begin{array}{*{20}{c}} {{s_0}} \\ {{s_1}} \\ {{s_2}} \\ {{s_3}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{n_{11}}}&{{n_{12}}}&{{n_{13}}}&{{n_{14}}} \\ {{n_{21}}}&{{n_{22}}}&{{n_{23}}}&{{n_{24}}} \\ {{n_{31}}}&{{n_{32}}}&{{n_{33}}}&{{n_{34}}} \\ {{n_{41}}}&{{n_{42}}}&{{n_{33}}}&{{n_{44}}} \end{array}} \right] \times \left[ {\begin{array}{*{20}{c}} {{I_{{\rm{out}}}}({\beta _0})} \\ {{I_{{\rm{out}}}}({\beta _1})} \\ {{I_{{\rm{out}}}}({\beta _2})} \\ {{I_{{\rm{out}}}}({\beta _3})} \end{array}} \right],$$ 7 其中
${n_{ij}}(i, j = 1, \cdot \cdot \cdot, 4)$ 为${{{M}}_{{\rm{tran}}}}$ 的逆矩阵中的元素. 将(7)式中的${{{S}}_{\rm{in}}}$ 带入(4)式即可求得Aop, Dop, Dolp和Docp.3 系统误差分析、优化与标定
3.1 系统误差分析与优化
在大气偏振模式的获取过程中会受到噪声的影响, 假设噪声对
${I_{{\rm{out}}}}$ 产生的扰动为$\Delta {I_{{\rm{out}}}}$ , 此时$\Delta {I_{{\rm{out}}}}$ 对${{{S}}_{\rm{in}}}$ 解算产生的扰动为$\Delta {{{S}}_{\rm{in}}}$ , 所以(6)式变为$$({I_{{\rm{out}}}} + \Delta {I_{{\rm{out}}}}) = {{{M}}_{{\rm{tran}}}} \times ({{S}}_{{\rm{in}}}^{\rm{T}} + \Delta {{S}}_{{\rm{in}}}^{\rm{T}}),$$ 8 所以系统解算出的
${{{S}}_{\rm{in}}}$ 的相对误差可以用(9)式表示[21]:$$\frac{{\left\| {\Delta {{S}}_{{\rm{in}}}^{\rm{T}}} \right\|}}{{\left\| {{{S}}_{{\rm{in}}}^{\rm{T}}} \right\|}} \leqslant \left\| {{{{M}}_{{\rm{tran}}}}} \right\| \times \left\| {{{{M}}_{{\rm{tran}}}}^{ - 1}} \right\| \times \frac{{\left\| {\Delta {I_{{\rm{out}}}}} \right\|}}{{\left\| {{I_{{\rm{out}}}}} \right\|}},$$ 9 其中
$||\cdots||$ 表示矩阵或向量的范数; 由条件数(condition number, Cond)的定义可知:$\left\| {{{{M}}_{{\rm{tran}}}}} \right\| \times $ $ \left\| {{{{M}}_{\rm{tran}}}^{ - 1}} \right\|$ 为${\rm{Cond}}({{{M}}_{{\rm{tran}}}})$ . 由(9)式可以看出: 当${I_{{\rm{out}}}}$ 因噪声产生扰动$\Delta {I_{{\rm{out}}}}$ 时,${{{S}}_{\rm{in}}}$ 的相对误差不超过${I_{{\rm{out}}}}$ 的相对误差乘以${\rm{Cond}}({{{M}}_{\rm{tran}}})$ . 由于$\Delta {I_{{\rm{out}}}}$ 是受噪声影响的不确定的量, 为了使${{{S}}_{\rm{in}}}$ 的相对误差尽量小, 则${\rm{Cond}}({{{M}}_{\rm{tran}}})$ 应该尽量的小. 由于条件数分析法常用于偏振仪的优化中, 所以本文使用条件数来优化传输矩阵[22]. 通过2.2的分析可以知道系统优化的本质是${{{M}}_{{\rm{tran}}}}$ 中β的选择问题, 所以β的选择依据如(10)式所示:$$({\beta _0},{\beta _1},{\beta _2},{\beta _3}) = \min ({\rm{Cond}}({{{M}}_{{\rm{tran}}}}(\beta ))).$$ 10 此外在选择β时还要保证
${{{M}}_{{\rm{tran}}}}$ 是可逆的, 而且选择的β要使得${{{M}}_{{\rm{tran}}}}$ 的行列式(Determinant, Det)尽量的大. 由于Det表示系统对输入信息的缩放程度, 所以选择${\rm{Det}} > 0$ 条件下对应的β.经过分析发现对β以2º, 5º和10º等为间隔进行遍历时, 优化得到的角度均在某些固定的角度左右. 因此, 结合数据量和运算时间选择5º为间隔对β在0º—180º之间进行遍历, 获得36个
$1 \times 4$ 的行向量, 然后求解所有行向量的组合数, 得到58905个${{{M}}_{{\rm{tran}}}}$ , 并求解相应矩阵的1范数Cond(${\rm{Con}}{{\rm{d}}_1}$ )、2范数Cond(${\rm{Con}}{{\rm{d}}_2}$ )、∞范数Cond(${\rm{Con}}{{\rm{d}}_\infty }$ )和Det. 由于矩阵的Cond过大, 所以在图3(a)—(c)中纵坐标为Cond的倒数, 结果如图3所示.图 3${{{M}}_{{\rm{tran}}}}$ 的$1/{\rm{Cond}}$ 与Det变化趋势 (a)$1/{\rm{Con}}{{\rm{d}}_1}$ 变化趋势; (b)$1/{\rm{Con}}{{\rm{d}}_2}$ 变化趋势; (c)$1/{\rm{Con}}{{\rm{d}}_\infty }$ 变化趋势; (d)$Det$ 变化趋势Figure 3. The change trend of$1/{\rm{Cond}}$ and Det of the${{{M}}_{{\rm{tran}}}}$ : (a) Change trend of$1/{\rm{Con}}{{\rm{d}}_1}$ ; (b) change trend of$1/{\rm{Con}}{{\rm{d}}_2}$ ; (c) change trend of$1/{\rm{Con}}{{\rm{d}}_\infty }$ ; (d) change trend of Det.从图3中可以看出, 三个Cond的变化趋势是类似的, 为了直观的体现Cond, Det与β的对应关系, 求解了最小
${\rm{Con}}{{\rm{d}}_1}$ ,${\rm{Con}}{{\rm{d}}_2}$ ,${\rm{Con}}{{\rm{d}}_\infty }$ 和最大Det对应的β, 结果如表1所列.表 1 最大Det、最小Cond与β对应表Table 1. Corresponding table of maximum Det, minimum Cond and β.条件数(Cond) 行列式(Det) 角度(β0, β1,
β2, β3)Cond1 8.6084 0.0923(最大) Det: (5°, 45°, 120°, 155°) Cond2 3.4864 Cond∞ 7.837 Cond1(最小) 7.9448 0.0814 Cond1: (0°, 30°, 115°, 150°) –0.0814 Cond1: (25°, 60°, 90°, 120°) Cond2 4.0092 — — Cond∞ 7.0859 — — Cond1 8.4653 — — Cond2(最小) 3.3381 –0.0919 Cond2: (35°, 70°, 100°, 135°) Cond∞ 7.791 — — Cond1 8.5767 — — Cond2 3.4563 — — Cond∞(最小) 6.2275 0.0911 Cond∞: (0°, 40°, 115°, 150°) 根据上述分析, 从表1中可以得到(5º, 45º, 120º, 155º), (0º, 30º, 115º, 150º)和(0º, 40°, 115º, 150º)三组角度. 因此, 在上述三组约束条件下求最优角度的问题就转化为多条件极值问题, 其中目标函数为
$L({\beta _0}, {\beta _1}, {\beta _2}, {\beta _3})$ 的定义如下式所示:$$L({\beta _0},{\beta _1},{\beta _2},{\beta _3}) = \frac{{{\rm{Con}}{{\rm{d}}_1}({{{M}}_{{\rm{tran}}}}(\beta ))\& {\rm{Con}}{{\rm{d}}_2}({{{M}}_{{\rm{tran}}}}(\beta ))\& {\rm{Con}}{{\rm{d}}_\infty }({{{M}}_{{\rm{tran}}}}(\beta ))}}{{{\rm{Det}}({{{M}}_{{\rm{tran}}}}(\beta ))}}.$$ 11 当
$L({\beta _0}, {\beta _1}, {\beta _2}, {\beta _3})$ 在约束条件(s.t.)下取得最小值时,$({\beta _0}, {\beta _1}, {\beta _2}, {\beta _3})$ 即为所求角度, 从表1得到约束条件(s.t.)如(12)式所示:$${\rm{s}}.{\rm{t}}{{.}}\;\left\{ {\begin{aligned} & {{\rm{7}}.{\rm{9448}} \leqslant {\rm{Con}}{{\rm{d}}_1}({{{M}}_{{\rm{tran}}}}(\beta )) \leqslant {\rm{8}}.{\rm{6048}}}, \\ &{{\rm{3}}.{\rm{3381}} \leqslant {\rm{Con}}{{\rm{d}}_2}({{{M}}_{{\rm{tran}}}}(\beta )) \leqslant {\rm{4}}.{\rm{0092}}}, \\ &{{\rm{6}}.{\rm{2275}} \leqslant {\rm{Con}}{{\rm{d}}_\infty }({{{M}}_{\rm{tran}}}(\beta )) \leqslant {\rm{7}}.{\rm{837}}}, \\ &{{\rm{0}}.{\rm{0814}} \leqslant {\rm{Det}}({{{M}}_{{\rm{tran}}}}(\beta )) \leqslant {\rm{0}}.{\rm{0923}}} . \end{aligned}} \right.$$ 12 由(11)式和(12)式计算得到,
$L({\beta _0}, {\beta _1}, {\beta _2}, {\beta _3})$ 在s.t.下的最优值为(5º, 40º, 120º, 155º)和(10º, 45º, 125º, 160º), 优化前后对比实验采用的角度为(10º, 45º, 125º, 160º).3.2 系统标定
${{{M}}_{{\rm{tran}}}}$ 的标定主要有两部分: 1) QWP的延迟量和偏振片消光比的标定; 2) 不同β下的${{{M}}_{{\rm{tran}}}}$ 的标定. 由于第一部分标定对${{{M}}_{{\rm{tran}}}}$ 的影响是确定的, 所以本文标定的重点是第二部分. 经过标定, 标定后QWP延迟量的误差小于 ± 0.5°, 偏振片消光比的误差小于5%. 第二部分的标定原理如图4所示, 标定系统包括光源、偏振光调制器、Stokes矢量仪和光功率计.本文在定标时采用的是主光线定标, 定标的光源采用的是激光光源, 主要是考虑到激光光源具有一致性较好、输出功率较强且衰减较少等优点, 激光发生器采用Thorlabs公司的高功率LDC4020型激光发生器. 具体定标步骤如下: 1) 选择488 nm波段下的激光光源, 通过偏振光调制器产生不同偏振态的入射光, 使用Stokes矢量仪记录待测系统入射光的偏振态, 并使用光功率计记录待测系统的入射光强值; 2) 把β依次设置为5º, 10º, 40º, 45º, 120º, 125º, 155º, 160°, 使用光功率计记录每个角度下待测系统的出射光强值; 3) 改变入射光的偏振态, 重复(2)得到10组偏振态下待测系统的入射光和每组入射光下8个β对应的输出光强值. 假设10组输入光强为
${I_{\rm{in}}}(i)$ , 归一化的Stokes矢量为${{{S}}_{\rm{in}}}(i)$ , 输出光强为${I_{\beta} }(i)$ , 则${I_{\rm{in}}}(i)$ ,${{{S}}_{\rm{in}}}(i)$ 和${I_{\beta} }(i)$ 满足的关系如下式所示:$$\begin{split} &[{m_{11}}(\beta ),{m_{12}}(\beta ),{m_{13}}(\beta ),{m_{14}}(\beta )] \\ &\times {\left( {{I_{\rm{in}}}(i) \times {{{S}}_{\rm{in}}}(i)} \right)^{\rm{T}}} = {I_{\beta} }(i). \end{split}$$ 13 由于(13)式有32个变量, 80个方程, 故该方程组为超定方程组, 因此使用最小二乘法进行求解, 得到系统传输矩阵的标定结果如表2所列.
表 2 系统传输矩阵标定结果Table 2. Calibration results of system transmission matrixQWP旋转角度 m11 m12 m13 m14 5° 0.500 –0.44 –0.0775 0.225 10° 0.500 –0.377 –0.137 0.299 40° 0.500 0.0102 0.0577 0.497 45° 0.500 0.001 0.144 0.479 120° 0.500 –0.0574 –0.0994 –0.487 125° 0.500 –0.00975 –0.0268 –0.499 155° 0.500 –0.269 0.32 –0.274 160° 0.500 –0.352 0.295 –0.198 4 实验与结果分析
为了对优化和标定结果进行说明, 开展了两项实验: 1) 优化前后偏振信息对比实验; 2) 外场全偏振信息获取实验. 实验地点为合肥工业大学翡翠湖校区操场(E117.214º, N31.778º), 实验时间为2019年12月13日上午9点, 此时太阳的高度角为29.66º, 方位角为152.1º, 天气晴朗.
4.1 优化前后对比实验
在对比实验中, 优化前QWP的旋转角度为(0º, 45º, 60º, 120º), 优化后QWP的旋转角度为(10º, 45º, 125º, 160º°), 本组实验对Aop, Dop和Dolp进行分析, 结果如图5所示. 为了量化优化效果, 本文选用熵、均值、标准差等指标进行说明.
由于Aop在太阳子午线两侧数值的极性是相反的, 因此无法直接使用熵、均值、标准差分析Aop数据的离散度, 但Aop的绝对值关于太阳子午线对称, 所以将Aop对称轴两侧数据的绝对值之差(∆Aop, 如图5(b)和图5(f)所示)作为分析数据, ∆Aop的熵、均值、标准差的求解结果如表3所列.
表 3 优化前后∆Aop数据离散度对比Table 3. Comparisons of data dispersion of the ∆Aop before and after optimization.Index 熵 均值 标准差 优化前 优化后 优化前 优化后 优化前 优化后 ∆Aop 1.996 1.7956 4.9553 4.2633 27.7036 24.709 从表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的差(
$\Delta {\rm{Do}}{{\rm{p}}_{{\rm{BMP}}}}$ 和$\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{BMP}}}}$ )作为分析数据; 而NZ为偏振度最小的区域, 所以NZ的误差分析以NZ中的数据与0的差($\Delta {\rm{Do}}{{\rm{p}}_{{\rm{NZ}}}}$ 和$\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{NZ}}}}$ )作为分析数据, 如下式所示:$$\begin{split} &\Delta {\rm{Do}}{{\rm{p}}_{{\rm{BMP}}}} = 1 - {\rm{Do}}{{\rm{p}}_{{\rm{BMP}}}},\;\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{BMP}}}} = 1 - {\rm{Dol}}{{\rm{p}}_{{\rm{BMP}}}}{{ ,}} \\ &\Delta {\rm{Do}}{{\rm{p}}_{{\rm{NZ}}}} = {\rm{Do}}{{\rm{p}}_{{\rm{NZ}}}},\;\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{NZ}}}} = {\rm{Dol}}{{\rm{p}}_{{\rm{NZ}}}}{{ }}{{.}}\\[-10pt] \end{split} $$ 14 为了保证数据的一致性, 先从Dop和Dolp图中裁剪相同区域的
$\Delta {\rm{Do}}{{\rm{p}}_{{\rm{BMP}}}}$ 和$\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{BMP}}}}$ 以及相同区域的$\Delta {\rm{Do}}{{\rm{p}}_{{\rm{NZ}}}}$ 和$\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{NZ}}}}$ , 然后求解上述区域的熵、均值、标准差, 裁剪结果如图6所示.表 4 优化前后BMP的数据离散度指标对比Table 4. Comparisons of data dispersion index of BMP before and after optimization.Index 熵 均值 标准差 优化前 优化后 优化前 优化后 优化前 优化后 $\Delta {\rm{Do}}{{\rm{p}}_{{\rm{BMP}}}}$ 2.6884 2.5550 0.7496 0.7252 0.0476 0.0359 $\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{BMP}}}}$ 2.6933 2.5564 0.7508 0.7256 0.0481 0.0361 从表4中可以得到优化前后
$\Delta {\rm{Do}}{{\rm{p}}_{{\rm{BMP}}}}$ 的熵、均值、标准差分别降低了4.964%, 3.255%和24.537%; 而$\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{BMP}}}}$ 的熵、均值、标准差分别降低了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 优化前后NZ数据离散度指标对比Table 5. Comparisons of data dispersion index of NZ before and after optimization.Index 熵 均值 标准差 Poa 优化前 优化后 优化前 优化后 优化前 优化后 优化前 优化后 $\Delta {\rm{Do}}{{\rm{p}}_{{\rm{NZ}}}}$ 1.5359 1.5141 0.0446 0.0413 0.0246 0.0222 31.9377 35.3431 $\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{NZ}}}}$ 1.5285 1.5085 0.0441 0.0409 0.0245 0.0223 32.8158 35.9437 从表5可以进一步得到优化前后
$\Delta {\rm{Do}}{{\rm{p}}_{{\rm{NZ}}}}$ 的熵、均值、标准差分别降低了1.423%, 7.414%和9.4%, Poa提高了10.663%; 而$\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{NZ}}}}$ 的熵、均值、标准差分别降低了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分量的分布模式中, 靠近太阳位置区域值较大, 远离太阳位置值较小.
图 7 目标天空区域光强图与偏振模式分布结果 (a) 偏振光强图; (b) Aop; (c) Dop; (d) Dolp; (e) Docp; (f) I分量图; (g) Q分量图; (h) U分量图; (i) V分量图Figure 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°的视场角且满足全偏振信息快速获取的要求; 同时, 针对系统的
${{{M}}_{{\rm{tran}}}}$ 因存在“病态”情况而易受噪声影响的问题, 从Cond和Det的角度对${{{M}}_{{\rm{tran}}}}$ 的进行分析, 得到了最小Cond和最大Det条件下对应的旋转角度, 并建立了多条件极值目标函数, 将旋转角度的求解问题转化为多条件极值问题, 得到最优的四分之一波片旋转角度. 此外, 还对${{{M}}_{{\rm{tran}}}}$ 、四分之一波片的延迟量、偏振片的消光比进行了标定, 并开展了优化前后的对比实验和全偏振大气偏振模式测量实验. 实验结果表明, 优化后系统获取的偏振信息较优化前有了较好的提升, 同时全偏振实验的获取结果也说明天空确实存在全偏振信息. -
图 3
${{{M}}_{{\rm{tran}}}}$ 的$1/{\rm{Cond}}$ 与Det变化趋势 (a)$1/{\rm{Con}}{{\rm{d}}_1}$ 变化趋势; (b)$1/{\rm{Con}}{{\rm{d}}_2}$ 变化趋势; (c)$1/{\rm{Con}}{{\rm{d}}_\infty }$ 变化趋势; (d)$Det$ 变化趋势Figure 3. The change trend of
$1/{\rm{Cond}}$ and Det of the${{{M}}_{{\rm{tran}}}}$ : (a) Change trend of$1/{\rm{Con}}{{\rm{d}}_1}$ ; (b) change trend of$1/{\rm{Con}}{{\rm{d}}_2}$ ; (c) change trend of$1/{\rm{Con}}{{\rm{d}}_\infty }$ ; (d) change trend of Det.图 7 目标天空区域光强图与偏振模式分布结果 (a) 偏振光强图; (b) Aop; (c) Dop; (d) Dolp; (e) Docp; (f) I分量图; (g) Q分量图; (h) U分量图; (i) V分量图
Figure 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.
表 1 最大Det、最小Cond与β对应表
Table 1. Corresponding table of maximum Det, minimum Cond and β.
条件数(Cond) 行列式(Det) 角度(β0, β1,
β2, β3)Cond1 8.6084 0.0923(最大) Det: (5°, 45°, 120°, 155°) Cond2 3.4864 Cond∞ 7.837 Cond1(最小) 7.9448 0.0814 Cond1: (0°, 30°, 115°, 150°) –0.0814 Cond1: (25°, 60°, 90°, 120°) Cond2 4.0092 — — Cond∞ 7.0859 — — Cond1 8.4653 — — Cond2(最小) 3.3381 –0.0919 Cond2: (35°, 70°, 100°, 135°) Cond∞ 7.791 — — Cond1 8.5767 — — Cond2 3.4563 — — Cond∞(最小) 6.2275 0.0911 Cond∞: (0°, 40°, 115°, 150°) 表 2 系统传输矩阵标定结果
Table 2. Calibration results of system transmission matrix
QWP旋转角度 m11 m12 m13 m14 5° 0.500 –0.44 –0.0775 0.225 10° 0.500 –0.377 –0.137 0.299 40° 0.500 0.0102 0.0577 0.497 45° 0.500 0.001 0.144 0.479 120° 0.500 –0.0574 –0.0994 –0.487 125° 0.500 –0.00975 –0.0268 –0.499 155° 0.500 –0.269 0.32 –0.274 160° 0.500 –0.352 0.295 –0.198 表 3 优化前后∆Aop数据离散度对比
Table 3. Comparisons of data dispersion of the ∆Aop before and after optimization.
Index 熵 均值 标准差 优化前 优化后 优化前 优化后 优化前 优化后 ∆Aop 1.996 1.7956 4.9553 4.2633 27.7036 24.709 表 4 优化前后BMP的数据离散度指标对比
Table 4. Comparisons of data dispersion index of BMP before and after optimization.
Index 熵 均值 标准差 优化前 优化后 优化前 优化后 优化前 优化后 $\Delta {\rm{Do}}{{\rm{p}}_{{\rm{BMP}}}}$ 2.6884 2.5550 0.7496 0.7252 0.0476 0.0359 $\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{BMP}}}}$ 2.6933 2.5564 0.7508 0.7256 0.0481 0.0361 表 5 优化前后NZ数据离散度指标对比
Table 5. Comparisons of data dispersion index of NZ before and after optimization.
Index 熵 均值 标准差 Poa 优化前 优化后 优化前 优化后 优化前 优化后 优化前 优化后 $\Delta {\rm{Do}}{{\rm{p}}_{{\rm{NZ}}}}$ 1.5359 1.5141 0.0446 0.0413 0.0246 0.0222 31.9377 35.3431 $\Delta {\rm{Dol}}{{\rm{p}}_{{\rm{NZ}}}}$ 1.5285 1.5085 0.0441 0.0409 0.0245 0.0223 32.8158 35.9437 -
[1] 高隽, 范之国 2014 仿生偏振光导航方法 (北京: 科学出版社) 第2页Gao J, Fan Z G 2014 Bionic polarized light navigation method (Beijing: Science Press) p2 (in Chinese) [2] Wehner R 2003 J. Comp. Physiol. A 189 579 doi: 10.1007/s00359-003-0431-1 [3] Kraft P, Evangelista C, Dacke M, Labhart T, Srinivasan M V 2011 Philos. Trans. R. Soc. London, Ser. B 366 703 doi: 10.1098/rstb.2010.0203 [4] Reppert S M, Zhu H, White R H 2014 Curr. Biol. 14 155 [5] Homberg U 2015 Front. Behav. Neurosci. 9 346 [6] Liang H J, Bai H Y, Liu N, Sui X B 2020 Math. Probl. Eng. 2020 1 [7] Li J S, Chu J K, Zhang R, Chen J H, Wang Y L 2020 Appl. Opt. 59 2955 doi: 10.1364/AO.387770 [8] Dupeyroux J, Viollet S, Serres J R 2019 J. R. Soc. Interface 16 20180878 doi: 10.1098/rsif.2018.0878 [9] 胡帅, 高太长, 李浩, 程天际, 刘磊, 黄威, 江诗阳 2016 物理学报 65 014203 doi: 10.7498/aps.65.014203Hu S, Gao T C, Li H, Cheng T J, Liu L, Huang W, Jiang S Y 2016 Acta Phys. Sin. 65 014203 doi: 10.7498/aps.65.014203 [10] 刘敬, 金伟其, 王霞, 鲁啸天, 温仁杰 2016 物理学报 65 094201 doi: 10.7498/aps.65.094201Liu J, Jin W Q, Wang X, Lu X T, Wen R J 2016 Acta Phys. Sin. 65 094201 doi: 10.7498/aps.65.094201 [11] 王晨光, 张楠, 李大林, 杨江涛, 王飞, 任建斌, 唐军, 刘俊, 薛晨阳 2015 光电工程 42 60 doi: 10.3969/j.issn.1003-501X.2015.03.010Wang C G, Zhang N, Li D L, Yang J T, Wang F, Ren J B, Tang J, Liu J, Xue C Y 2015 Opto-Electron. Eng. 42 60 doi: 10.3969/j.issn.1003-501X.2015.03.010 [12] Horváth G, Barta A, Gál J, Suhai B, Haiman O 2002 Appl. Opt. 41 543 doi: 10.1364/AO.41.000543 [13] Pust N J, Shaw J A 2006 Appl. Opt. 45 22 doi: 10.1364/AO.45.000022 [14] 孙洁, 高隽, 怀宇, 毕冉, 范之国 2016 光电工程 43 45 doi: 10.3969/j.issn.1003-501X.2016.07.008Sun J, Gao J, Huai Y, Fan Z G 2016 Opto-Electron. Eng. 43 45 doi: 10.3969/j.issn.1003-501X.2016.07.008 [15] 戴俊, 高隽, 范之国 2017 中国激光 44 184Dai J, Gao J Fan Z G 2017 Chin. J. Las. 44 184 [16] Hsu W L, Myhre G, Balakrishnan K, Brock N, Ibn-Elhaj M, Pau S. 2014 Opt. Express 22 3063 doi: 10.1364/OE.22.003063 [17] 张忠顺 2014 硕士学位论文 (合肥: 合肥工业大学)Zhang Z S 2014 M. S. Thesis (Hefei: Hefei University of Technology) (in Chinese) [18] 殷玉龙, 孙晓兵, 宋茂新, 陈卫, 陈斐楠 2019 物理学报 68 024203 doi: 10.7498/aps.68.20181553Yin Y L, Sun X B, Song M X, Chen W, Chen F N 2019 Acta Phys. Sin. 68 024203 doi: 10.7498/aps.68.20181553 [19] Kiyohara J, Ueno S, Kitai R, Kurokawa H, Makita M, Ichimoto K 2004 Ground-based Instrumentation for Astronomy (Glasgow: SPIE) p1778 [20] Anan T, Ichimoto K, Oi A, Kimura G, Nakatani Y, Ueno S 2012 Ground-based Instrumentation for Astronomy IV (Glasgow: SPIE) p84461C [21] Isaacson E, Keller H B 2012 Analysis of numerical methods (Massachusetts: Courier) pp54−55 [22] Iniesta J C D T, Collados M 2000 Appl. Opt. 39 1637 doi: 10.1364/AO.39.001637 [23] Worster S, Mouritsen H, Hore P J 2017 J. R. Soc. Interface 14 134 -