蔡 弘，張 然，關(guān) 樂(lè)，萬(wàn)振華，褚金奎※

（1.大連理工大學(xué)機(jī)械工程學(xué)院，遼寧大連116023；2.大連大學(xué)機(jī)械工程學(xué)院，遼寧大連116622）

0 引言

非偏振的太陽(yáng)光傳輸?shù)酱髿庵?，?huì)被大氣中的塵埃、分子以及氣溶膠等微粒散射和吸收，從而變成部分偏振光[1]。晴朗天空的偏振模式符合瑞利單次散射理論，在某一時(shí)刻某一位置對(duì)應(yīng)的天空偏振模式是穩(wěn)定的，太陽(yáng)的位置隨著時(shí)間及觀測(cè)地點(diǎn)的變化而有規(guī)律的移動(dòng)，因此在天空偏振模式下獲取載體坐標(biāo)系下太陽(yáng)方位角可以用于導(dǎo)航定向[2-3]。

當(dāng)前天空偏振導(dǎo)航系統(tǒng)主要分為點(diǎn)源式和圖像式。本課題組較早使用點(diǎn)源式偏振導(dǎo)航系統(tǒng)，并基于瑞利散射理論獲取載體坐標(biāo)系下太陽(yáng)方位角[4-7]。但由于點(diǎn)源式偏振導(dǎo)航系統(tǒng)局限于測(cè)量天空某幾個(gè)點(diǎn)的偏振信息，不能很好地測(cè)量全天空域的偏振信息，易受環(huán)境噪聲的影響。于是開(kāi)發(fā)了圖像式偏振導(dǎo)航系統(tǒng)[8-9]。國(guó)防科技大學(xué)使用分通道式的圖像偏振導(dǎo)航系統(tǒng)在晴朗天氣下最大定向誤差為0.5°[10]，由于多個(gè)通道不是經(jīng)過(guò)同一鏡頭入射到CMOS焦平面上，存在瞬時(shí)視場(chǎng)誤差。

當(dāng)前在晴朗天氣下利用天空偏振場(chǎng)圖獲取太陽(yáng)方向的算法主要為識(shí)別太陽(yáng)子午線，因?yàn)樵谇缋侍鞖庀缕穹轿唤顷P(guān)于太陽(yáng)子午線具有很好的對(duì)稱性。清華大學(xué)提出基于霍夫線變換提取太陽(yáng)子午線的方法，在晴朗無(wú)云天氣下精度為0.313°[11]。中北大學(xué)使用分通道式的圖像偏振導(dǎo)航系統(tǒng)基于最小二乘法擬合太陽(yáng)子午線[12]，本課題組提出分箱法擬合太陽(yáng)子午線[13]。

由于太陽(yáng)子午線上偏振光信息并沒(méi)有天空?qǐng)D像中全局偏振光信息豐富，利用擬合太陽(yáng)子午線方法比基于瑞利散射理論獲取太陽(yáng)方向方法所使用的天空偏振光束少，擬合太陽(yáng)子午線方法在完全被遮擋情況下可能會(huì)失效，同時(shí)研究表明，生物中昆蟲(chóng)的復(fù)眼也是利用多個(gè)方向的偏振小眼感知大視場(chǎng)天空區(qū)域的偏振光[14-15]。本文為了獲取更大視場(chǎng)的天空偏振光，使用魚(yú)眼鏡頭結(jié)合偏振相機(jī)搭建圖像式偏振定向測(cè)量系統(tǒng)，提出一種以信標(biāo)發(fā)光二極管為基準(zhǔn)的魚(yú)眼鏡頭天頂角標(biāo)定方法，利用瑞利散射理論和太陽(yáng)矢量約束計(jì)算載體坐標(biāo)系下的太陽(yáng)方位角。室外實(shí)驗(yàn)表明本文方法比擬合太陽(yáng)子午線的兩種方法（最小二乘法、分箱法）誤差更小。

1 理論

1.1 獲取天空偏振場(chǎng)圖

圖1 所示為實(shí)驗(yàn)所用設(shè)備，實(shí)驗(yàn)時(shí)主要使用三腳架支撐偏振探測(cè)系統(tǒng)，通過(guò)網(wǎng)口連接電腦采集天空?qǐng)D像。偏振探測(cè)系統(tǒng)包含魚(yú)眼鏡頭和偏振相機(jī)，偏振相機(jī)主要由微透鏡陣列、偏振片陣列以及像素陣列組成，其中偏振片陣列包含0°、45°、90°和135°共4個(gè)偏振方向。

圖1 實(shí)驗(yàn)設(shè)備

斯托克斯矢量表達(dá)式如下：

I、Q、U三個(gè)變量的獲取如下：

偏振度DOP表征天空偏振態(tài)，計(jì)算如下：

偏振角AOP是載體坐標(biāo)系中天空觀測(cè)點(diǎn)的偏振方向矢量與載體（相機(jī)）體軸的夾角，計(jì)算如下：

偏振方位角AOE是天空觀測(cè)坐標(biāo)系中天空觀測(cè)點(diǎn)的偏振方向矢量與當(dāng)?shù)刈游缇€切線方向的夾角，計(jì)算如下：

式中：（i,j）為像素點(diǎn)坐標(biāo)；h為圖像的高；w為圖像的寬。

1.2 獲取載體坐標(biāo)系下太陽(yáng)方位角

圖2 所示為晴朗天空偏振分布模式，以地面觀測(cè)點(diǎn)為原點(diǎn)建立載體坐標(biāo)系（即相機(jī)坐標(biāo)系XcYcZc），在載體坐標(biāo)系中太陽(yáng)方向矢量OS與軸Zc的夾角γs為太陽(yáng)天頂角，OS與軸Xc的夾角a′s為載體坐標(biāo)系下太陽(yáng)方位角。天空觀測(cè)方向矢量OP與軸Zc的夾角γ為天空觀測(cè)點(diǎn)天頂角，OP與軸Xc的夾角a為天空觀測(cè)點(diǎn)方位角，天頂角的角度范圍為閉區(qū)間[0°，90°]，方位角以軸Xc正方向?yàn)?°，逆時(shí)針?lè)较驗(yàn)檎?，順時(shí)針?lè)较驗(yàn)樨?fù)，角度范圍為閉區(qū)間[-180°，+180°]。

圖2 晴朗天空偏振分布模型

圖中的偏振方向矢量E在載體坐標(biāo)系中表示為cE，在天空觀測(cè)點(diǎn)坐標(biāo)系中表示為iE。在相機(jī)坐標(biāo)系中天空觀測(cè)點(diǎn)P位置的偏振方向矢量cE表示為：

以天空觀測(cè)點(diǎn)P為原點(diǎn)建立天空觀測(cè)坐標(biāo)系XiYiZi，P位置的偏振方向矢量iE與軸Xi正方向的夾角φ為偏振方位角即AOE，該矢量表示如下：

cRi表示天空觀測(cè)坐標(biāo)系相對(duì)載體坐標(biāo)系的旋轉(zhuǎn)矩陣，包含了天空觀測(cè)點(diǎn)方位角a和天空觀測(cè)點(diǎn)天頂角γ的具體展開(kāi)形式為：

由于魚(yú)眼鏡頭采集的圖像存在較大的徑向畸變，徑向畸變主要影響天空觀測(cè)點(diǎn)的天頂角，因此為了獲取準(zhǔn)確的天空觀測(cè)點(diǎn)的天頂角，對(duì)像素距離（圖像中像素點(diǎn)與圖像中心的距離）與天空觀測(cè)點(diǎn)的天頂角進(jìn)行標(biāo)定。圖3所示為魚(yú)眼鏡頭的天頂角標(biāo)定原理，固定偏振相機(jī)與魚(yú)眼鏡頭，將發(fā)光二極管與魚(yú)眼鏡頭的光軸對(duì)齊，并繞鏡頭的光心逆時(shí)針旋轉(zhuǎn)，每旋轉(zhuǎn)一次對(duì)應(yīng)的觀測(cè)點(diǎn)的天頂角γ累加10°并采集圖像，共旋轉(zhuǎn)9次，旋轉(zhuǎn)角度范圍為閉區(qū)間[0°，90°]，模擬了半球狀的天空分布。每一位置的發(fā)光二極管在圖像中對(duì)應(yīng)一個(gè)像素點(diǎn)，例如位置A處的發(fā)光二極管在圖像中對(duì)應(yīng)著像素點(diǎn)A′，同理B對(duì)應(yīng)著B(niǎo)′。要標(biāo)定觀測(cè)點(diǎn)的天頂角γ與像素距離（圖像中像素點(diǎn)與圖像中心A′的距離）的映射關(guān)系。將實(shí)際得到的觀測(cè)點(diǎn)的天頂角γ作為因變量，對(duì)應(yīng)的像素距離作為自變量，可以得到9個(gè)點(diǎn)：（36.0,10），（73.0,20），（110.0,30），（146.0,40），（182.0,50），（219.0,60），（255.0,70），（291.1,80），（327.1,90）。用Matlab曲線擬合工具箱擬合觀測(cè)點(diǎn)的天頂角γ與像素距離（dis）之間的映射關(guān)系，當(dāng)γ=0.2688dis∧1.004時(shí)，均方根誤差最小，為0.105。

圖3 魚(yú)眼鏡頭天頂角標(biāo)定原理

根據(jù)瑞利單次散射理論，任意天空觀測(cè)點(diǎn)P的偏振方向矢量cE垂直于由天空觀測(cè)點(diǎn)P、太陽(yáng)S及地面觀測(cè)點(diǎn)O組成的平面，則cE垂直于OS。天空中所有符合瑞利散射的點(diǎn)在相機(jī)坐標(biāo)系中的偏振方向矢量組成矩陣?yán)硐肭闆r下：

但實(shí)際由于誤差存在，需對(duì)太陽(yáng)方向矢量作最優(yōu)估計(jì)：

式（13）為對(duì)式（12）中函數(shù)f(OS)求導(dǎo)，為了獲取該函數(shù)的最小值，令導(dǎo)數(shù)等于0。由于OS不為零向量，則則f(OS)min=λmin，則太陽(yáng)方向矢量OS的最優(yōu)估計(jì)為矩陣的最小特征值λmin對(duì)應(yīng)的特征向量。

2 實(shí)驗(yàn)與結(jié)果

實(shí)驗(yàn)時(shí)間為2020年12月19日，天空晴朗無(wú)云，實(shí)驗(yàn)地點(diǎn)為大連理工大學(xué)機(jī)械學(xué)院2號(hào)樓樓頂，坐標(biāo)為（121.538 857，38.886 469）。實(shí)驗(yàn)時(shí)用圖1所示設(shè)備采集11：25—14：16這一時(shí)間段的天空強(qiáng)度圖。大概是每隔10 min采集一次，如表1所示，共采集了14組數(shù)據(jù)。實(shí)驗(yàn)過(guò)程以11：25時(shí)采集的天空?qǐng)D像獲取載體坐標(biāo)系下太陽(yáng)方位角為例，如圖3所示。

表1 實(shí)驗(yàn)數(shù)據(jù)（°）

圖4 包含4張圖像，均裁剪為400×400大小并做偽彩色處理。第一張為天空強(qiáng)度（Image）圖，為斯托克斯矢量中的I變量，第二張為偏振度（DOP）圖，第三張為偏振角（AOP）圖，第四張為偏振方位角（AOE）圖，AOE圖中Xc軸為載體體軸，圖中標(biāo)記的黑線表示太陽(yáng)子午線的位置，夾角a′s為載體坐標(biāo)系下太陽(yáng)方位角。圖4中Image圖顯示了存在太陽(yáng)過(guò)曝區(qū)域，該區(qū)域?qū)?yīng)的DOP=0，實(shí)驗(yàn)時(shí)將DOP=0的像素點(diǎn)過(guò)濾，然后用本文提出的算法獲取載體坐標(biāo)下太陽(yáng)方位角。

圖4 獲取太陽(yáng)方向?qū)嶒?yàn)過(guò)程

表1 所示為14組實(shí)驗(yàn)數(shù)據(jù)，每一組標(biāo)注了地理坐標(biāo)系下太陽(yáng)方位角和太陽(yáng)天頂角，通過(guò)查閱天文年歷可以獲得。并列出了最小二乘法擬合太陽(yáng)子午線[12]、分箱法擬合太陽(yáng)子午線[13]以及本文算法這3種方法獲取的載體坐標(biāo)系下太陽(yáng)方位角。

圖5 所示為14組數(shù)據(jù)的太陽(yáng)相對(duì)轉(zhuǎn)角的曲線，其中▇曲線表示理想太陽(yáng)相對(duì)轉(zhuǎn)角，通過(guò)將每組的地理坐標(biāo)系下太陽(yáng)方位角與第一組作差得到，●曲線表示用最小二乘法得到的太陽(yáng)相對(duì)轉(zhuǎn)角組成的曲線，▲曲線表示用分箱法得到的太陽(yáng)相對(duì)轉(zhuǎn)角組成的曲線，▼曲線表示用本文算法得到的太陽(yáng)相對(duì)轉(zhuǎn)角組成的曲線，這3種方法獲得的太陽(yáng)相對(duì)轉(zhuǎn)角均是將每組的載體坐標(biāo)系下太陽(yáng)方位角與第一組作差得到?？梢钥闯鰣D中的4條曲線具有相同的變化趨勢(shì)，但其中●曲線與▇曲線（理想太陽(yáng)相對(duì)轉(zhuǎn)角）每組數(shù)據(jù)有較大偏差，而▼曲線相比于▲曲線，其與▇曲線每組數(shù)據(jù)更較為接近，說(shuō)明本文算法獲取的太陽(yáng)相對(duì)轉(zhuǎn)角更加準(zhǔn)確。

圖5 太陽(yáng)相對(duì)轉(zhuǎn)角

為了更好地說(shuō)明本文算法的優(yōu)越性，對(duì)上述3種算法分別與理想太陽(yáng)相對(duì)轉(zhuǎn)角作差并擬合了3條曲線，如圖6所示。由誤差曲線圖可以看出▼曲線波動(dòng)最小且一直在0附近變化。又分別計(jì)算3種算法的均方根誤差（RMSE），最小二乘法為3.12°，分箱法為0.9°，本文算法為0.28°。本文算法的均方根誤差最小，因此，本文算法要優(yōu)于另兩種算法，在晴朗無(wú)云的天氣下具有可行性。

圖6 誤差曲線

3 結(jié)束語(yǔ)

本文首先基于斯托克斯矢量獲取天空偏振場(chǎng)圖，然后提出一種以信標(biāo)發(fā)光二極管為基準(zhǔn)的魚(yú)眼鏡頭天頂角標(biāo)定方法，利用瑞利散射理論和太陽(yáng)矢量約束計(jì)算載體坐標(biāo)系下的太陽(yáng)方位角，最后進(jìn)行實(shí)驗(yàn)分析。通過(guò)與擬合太陽(yáng)子午線的兩種算法（最小二乘法和分箱法）比較，本文算法具有最小的均方根誤差，在實(shí)際的晴朗無(wú)云天氣下進(jìn)行導(dǎo)航定向具有可行性。