基於瑞利模型的E矢量叉乘太陽方位角解算方法

        太陽光在其傳播過程中由於受到大氣層中空氣分子、氣溶膠粒子等的散射作用，以及地表植被、水面等的反射作用而產生了偏振光。雖然人類的視覺無法直接感知偏振光，但許多生物如沙蟻、蜜蜂、蝗蟲等昆蟲，甚至大耳蝙蝠，卻能夠利用其獨特的視覺結構，感知並利用光的偏振現象獲取信息，從而進行導航、覓食、交流、遷徙等活動。一系列的解剖生物學研究證明，部分生物具有偏振光導航能力是依靠其復眼中特殊的偏振神經感光結構能夠對快速變化的大氣偏振模式信息進行檢測和處理，從而提取可靠的羅盤信息進行導航。這一研究結果為仿生偏振光導航提供了生物基礎，開辟了基於仿生的偏振光導航的研究方向。

        大氣偏振模式是天空中的偏振光最終形成的具有一定規律的偏振態分布，其中蘊含重要的方向信息。生物研究表明，沙蟻、蜜蜂等昆蟲的導航不是應用天空單點的大氣偏振模式信息，而是利用整個或局部具有穩定分布的大氣偏振模式信息，並從中提取體軸與太陽子午線的夾角，從而實現導航功能。因此，利用大氣偏振模式有效獲取作為天空顯著點的太陽位置信息，就可以為導航提供參考基准，是實現偏振光導航的前提和關鍵技術，同時，也對太陽位置追蹤和利用有着重要的意義。

        大氣偏振模式的瑞利表征方法：晴朗天氣情況下，大氣粒子對光的散射過程，滿足瑞利散射條件，因此可通過瑞利散射理論得到大氣偏振模式的表征。瑞利表征認為太陽光在天球表面上的各點位置同大氣粒子發生一次散射后，即被觀測者接收，這能准確描述在理想條件下大氣中光纖呈現的分布特征，是應用最為廣泛的一種表征方法。而大氣偏振模式則可看成是大氣中各個點一次散射作用的集合效應。瑞利散射模型是現有描述理想條件下大氣偏振模式最為經典的表征方法，與實際大氣偏振模式具有較高的相似性，通常使用偏振度（DOP）和偏振角度（AOP，也稱E-矢量方向）來描述，下圖給出了大氣偏振模式的瑞利表征建模坐標系。圖中，觀測者為原點，太陽位置為，天定點為，向天頂方向為軸，正東方向為軸，正北方向為軸。球面上任意一點的天頂角為，方位角為，以正北方向為0，正北往東為正值。太陽空間位置表示為，其中，表示太陽天頂角，表示太陽方位角，太陽的高度角為，太陽子午線在地理坐標系中的位置以方位角來表示。

        偏振化方向角定義為被測點處E-矢量方向同所在子午線夾角，根據瑞利散射定律，任一被測點處的E-矢量為：

 其中，，。

        當太陽位置坐標在三維坐標下表示為，其高度角為，方位角為時，可以獲得天空中任意一點，高度角為，方位角為時的偏振化方向角的理論值為：

        下圖表示瑞利模型下的幾種坐標系之間的關系，包括地理坐標系、相機坐標系、入射光坐標系。圖中，表示地理北， 表示地理坐標系，表示相機坐標系，表示入射光坐標系，與分別與圖像的橫縱坐標軸平行，與相機平面垂直，指向天定點；表示太陽位置，表示觀測點，與分別為太陽位置在地理坐標系中的天頂角和方位角，與分別為觀測點在相機坐標系中的天頂角和方位角；為在入射光坐標系中的E-矢量，該矢量與、、三點構成的平面垂直，由矢量與矢量叉乘得到；入射光坐標系的軸與E-矢量之間的夾角即為該處的偏振化方向角（簡稱偏振角）。

        入射光坐標系中的E-矢量為關於偏振角的函數，為求解E-矢量則首先應求解偏振角，求解偏振角的步驟如下：首先獲得晴朗天氣下的天空圖像，使用相機NiKon-D800拍攝、、和四個偏振角度的大氣偏振圖像，相機設備如下圖所示，在相機鏡頭前安裝魚眼鏡頭，並在魚眼鏡頭前放置偏振片，旋轉偏振片至不同角度並拍攝以獲取不同偏振角度的大氣的偏振圖像。

         拍攝得到的四幅不同角度的大氣偏振圖像如下，拍攝時間為2019年6月27日19點04分，從左至右的偏振圖像的角度分別為、、和。

                                    

       設四幅偏振圖像的光強分別為、、和，由此計算該偏振態下的斯托克斯矢量分量、、：

        由、、計算每一個像素處的偏振角為：

 

        重繪圖，使其沿着太陽子午線中心反對稱且呈∞字形分布，計算公式為：

        設天空中不同觀測點位置對應像素中的兩個坐標分別為和，則兩個坐標對應的天頂角和方位角分別為：

其中為圖像的中心位置的坐標，為相機的焦距。

       偏振角圖像中，處的像素灰度值，為該兩觀測點處的散射光的偏振角的值，則在入射光坐標系中表示E-矢量為：

       在三維空間的坐標系為，繞三個坐標軸其中一軸旋轉得到新的坐標系，則原先坐標系中的點用旋轉后的坐標系表示時，點的坐標轉換關系如下：

       為了計算太陽矢量，需要將E-矢量從入射光坐標系轉換到相機坐標系。若要將入射光坐標系轉化到相機坐標系，需要將入射光坐標系旋轉兩次：將入射光坐標系繞軸向上旋轉角度，得到新的坐標系；將坐標系繞軸（軸與軸重合）向右旋轉角度，得到相機坐標系。每旋轉一次坐標系，左乘對應的旋轉矩陣，旋轉矩陣中的角度對應於坐標軸旋轉的角度，由此可得入射光坐標系到相機坐標系的坐標轉換矩陣為：

則將，處的E-矢量在相機坐標系中表示為：

       由任意兩個E-矢量叉乘得到太陽矢量的規律，可表示太陽矢量為：

       上述計算太陽矢量時，選取一對像素處的E-矢量叉乘得到太陽矢量，這樣只進行一次運算得到的結果誤差太大，需要選取多組E-矢量叉乘得到一組太陽矢量，對這組太陽矢量利用最小二乘法求解方位角，這樣得到的結果更加准確。

       利用最小二乘法求解方位角：設為所求太陽矢量，設矩陣由向量構成（）。構造損失函數為：

其中為3x3的對稱矩陣。

      利用拉格朗日乘數法求解上述損失函數，定義函數為：

      對函數關於求導，可得：

其中為3x3的單位矩陣，為矩陣的特征值。由於為非零向量，故當為0時，有。

       將代入到中，可得：

所以損失函數的最優解對應於矩陣的最大特征值的特征向量，該特征向量即為所求的太陽矢量，則有方位角為：

      本例中選取相機坐標系的軸和地理北極重合，所以航向角滿足：或者。

      上述算法的C++程序如下給出，在VS2017搭配opencv3.4.1環境下運行，另需配置Eigen數學庫：

#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <cmath>
#include <algorithm>
#include <opencv2\opencv.hpp>
#include <opencv2\core.hpp>
#include <opencv2\highgui.hpp>
#include <opencv2\imgproc.hpp>

#include <Eigen\core>
#include <Eigen\Dense>
#include <opencv2\core\eigen.hpp>

using namespace cv;
using namespace std;

typedef std::numeric_limits<double> Info;
double const NAN_d = Info::quiet_NaN();

void Vector2Mat(vector<vector<double>>src, Mat dst);
void cv_to_eigen(const Mat& input, Eigen::Matrix3d& output);

int main()
{
    string loc1 = "E:\\sunvector\\sunvector\\0.JPG";
    string loc2 = "E:\\sunvector\\sunvector\\45.JPG";
    string loc3 = "E:\\sunvector\\sunvector\\90.JPG";
    string loc4 = "E:\\sunvector\\sunvector\\135.JPG";

    Mat src1 = imread(loc1);
    Mat src2 = imread(loc2);
    Mat src3 = imread(loc3);
    Mat src4 = imread(loc4);

    src1.convertTo(src1, CV_64FC3, 1.0, 0);
    src2.convertTo(src2, CV_64FC3, 1.0, 0);
    src3.convertTo(src3, CV_64FC3, 1.0, 0);
    src4.convertTo(src4, CV_64FC3, 1.0, 0);

    int x1 = 1175; int x2 = x1 + 401;
    int y1 = 1863; int y2 = y1 + 401;

    int rr = 200; int oo = 200; int dd = 400;

    Mat image0 = src1(Rect(Point(x1, y1), Point(x2, y2)));
    Mat image45 = src2(Rect(Point(x1, y1), Point(x2, y2)));
    Mat image90 = src3(Rect(Point(x1, y1), Point(x2, y2)));
    Mat image135 = src4(Rect(Point(x1, y1), Point(x2, y2)));

    vector<Mat> planes0;
    split(image0, planes0);
    Mat I0 = (planes0[0] + planes0[1] + planes0[2]) / 3;

    vector<Mat> planes45;
    split(image45, planes45);
    Mat I45 = (planes45[0] + planes45[1] + planes45[2]) / 3;

    vector<Mat> planes90;
    split(image90, planes90);
    Mat I90 = (planes90[0] + planes90[1] + planes90[2]) / 3;

    vector<Mat> planes135;
    split(image135, planes135);
    Mat I135 = (planes135[0] + planes135[1] + planes135[2]) / 3;

    Mat I = (I0 + I45 + I90 + I135) / 2;
    Mat Q = I0 - I90;
    Mat U = I45 - I135;

    int rows = I.rows; int cols = I.cols;
    double PI = 3.1415926;
    Mat aop(rows, cols, CV_64FC1);

    for (unsigned int i = 0; i < rows; i++)
    {
        const double* Qptr = Q.ptr<double>(i);
        const double* Uptr = U.ptr<double>(i);
        double* aopPtr = aop.ptr<double>(i);
        for (unsigned int j = 0; j < cols; j++)
        {
            *aopPtr = (0.5*atan2(*Uptr, *Qptr)) * (180 / PI);
            aopPtr++;
            Uptr++; Qptr++;
        }
    }

    Mat aop_last = Mat::zeros(rows, cols, CV_64FC1);

    for (unsigned int i = 0; i < rows; i++)
    {
        const double* aopPtr = aop.ptr<double>(i);
        double* aop_lastPtr = aop_last.ptr<double>(i);
        for (unsigned int j = 0; j < cols; j++)
        {
            double x = i; double y = j;
            if (j == oo)
            {
                *aop_lastPtr = *aopPtr + (PI / 2) * (180 / PI);
            }
            else
            {
                *aop_lastPtr = *aopPtr - (atan((x - oo) / (oo - y))) * (180 / PI);
                if (*aop_lastPtr < -90)
                {
                    *aop_lastPtr = 180 + *aop_lastPtr;
                }
                else if (*aop_lastPtr > 90)
                {
                    *aop_lastPtr = *aop_lastPtr - 180;
                }
            }
            aopPtr++; aop_lastPtr++;
        }
    }

    Mat aop_circle = Mat::zeros(rows, cols, CV_64FC1);

    for (unsigned int i = 0; i < rows; i++)
    {
        const double* aop_lastPtr = aop_last.ptr<double>(i);
        double* aop_circlePtr = aop_circle.ptr<double>(i);
        for (unsigned int j = 0; j < cols; j++)
        {
            if ((i - oo)*(i - oo) + (j - oo)*(j - oo) < rr * rr)
            {
                *aop_circlePtr = *aop_lastPtr;
            }
            else
            {
                *aop_circlePtr = NAN_d;
            }
            aop_lastPtr++; aop_circlePtr++;
        }
    }

    Mat aop_color; Mat aop_show = Mat::zeros(rows, cols, CV_64FC3);
    double min_aop, max_aop, alpha_aop;

    minMaxLoc(aop_circle, &min_aop, &max_aop);
    Mat aop_circle1 = aop_circle;
    alpha_aop = 255.0 / (max_aop - min_aop);
    aop_circle1.convertTo(aop_circle1, CV_8U, alpha_aop, -min_aop * alpha_aop);
    applyColorMap(aop_circle1, aop_color, COLORMAP_JET);

    imwrite("aop.bmp", aop_color);
    string aop_map1 = "E:\\sunvector\\sunvector\\aop.bmp";

    Mat aop_map = imread(aop_map1);
    aop_map.convertTo(aop_map, CV_64FC3, 1 / 255.0, 0);
    //imshow("aop",aop_map);

    vector<Mat> channels(3);
    split(aop_map, channels);

    for (unsigned int i = 0; i < rows; i++)
    {
        const double* channel0Ptr = channels[0].ptr<double>(i);
        const double* channel1Ptr = channels[1].ptr<double>(i);
        const double* channel2Ptr = channels[2].ptr<double>(i);
        double* aop_showPtr = aop_show.ptr<double>(i);
        for (unsigned int j = 0; j < cols; j++)
        {
            if ((i - oo)*(i - oo) + (j - oo)*(j - oo) < rr * rr)
            {
                *aop_showPtr++ = *channel0Ptr; *aop_showPtr++ = *channel1Ptr; *aop_showPtr++ = *channel2Ptr;
            }
            else
            {
                *aop_showPtr++ = 1; *aop_showPtr++ = 1; *aop_showPtr++ = 1;
            }
            channel0Ptr++; channel1Ptr++; channel2Ptr++;
        }
    }

    imshow("aop_show",aop_show);
    aop_show.convertTo(aop_show, CV_8UC3, 255.0, 0);
    imwrite("aop_save.bmp", aop_show);

    int amount = 10000; int q = 0;
    vector<vector<double>> Su;

    for (unsigned int inter = 0; inter < amount; inter++)
    {
        int i1 = rand() % (rows - 1);
        int j1 = rand() % (cols - 1);
        int i2 = rand() % (rows - 1);
        int j2 = rand() % (cols - 1);
        if (((i1 - oo)*(i1 - oo) + (j1 - oo)*(j1 - oo) < rr * rr) && ((i2 - oo)*(i2 - oo) + (j2 - oo)*(j2 - oo) < rr * rr))
        {
            q++;
            double x1 = i1 - oo; double y1 = j1 - oo;
            double x2 = i2 - oo; double y2 = j2 - oo;

            double k1 = (aop_circle.at<double>(i1, j1))*(PI / 180);
            double k2 = (aop_circle.at<double>(i2, j2))*(PI / 180);

            double b1 = atan(y1 / x1);
            double b2 = atan(y2 / x2);

            double a1 = atan(sqrt(x1*x1 + y1 * y1)*(0.0087 / 8));
            double a2 = atan(sqrt(x2*x2 + y2 * y2)*(0.0087 / 8));

            Eigen::Matrix<double, 3, 3> C11;
            Eigen::Matrix<double, 3, 3> C12;
            C11 << cos(a1), 0, -sin(a1), 0, 1, 0, sin(a1), 0, cos(a1);
            C12 << cos(b1), sin(b1), 0, -sin(b1), cos(b1), 0, 0, 0, 1;
            Eigen::Matrix<double, 3, 3> Cli1 = C11 * C12;

            Eigen::Matrix<double, 3, 3> C21;
            Eigen::Matrix<double, 3, 3> C22;
            C21 << cos(a2), 0, -sin(a2), 0, 1, 0, sin(a2), 0, cos(a2);
            C22 << cos(b2), sin(b2), 0, -sin(b2), cos(b2), 0, 0, 0, 1;
            Eigen::Matrix<double, 3, 3> Cli2 = C21 * C22;

            Eigen::Matrix<double, 1, 3> PEi1;
            Eigen::Matrix<double, 1, 3> PEi2;
            PEi1 << cos(k1), sin(k1), 0;
            PEi2 << cos(k2), sin(k2), 0;

            Eigen::Matrix<double, 1, 3> e1 = PEi1 * Cli1;
            Eigen::Matrix<double, 1, 3> e2 = PEi2 * Cli2;

            Eigen::Vector3d Evec1 = Eigen::Vector3d(e1(0, 0), e1(0, 1), e1(0, 2));
            Eigen::Vector3d Evec2 = Eigen::Vector3d(e2(0, 0), e2(0, 1), e2(0, 2));
            Eigen::Vector3d S0 = Evec1.cross(Evec2);
            double normS0 = sqrt(S0[0] * S0[0] + S0[1] * S0[1] + S0[2] * S0[2]);
            Eigen::Vector3d S = S0 / normS0;

            vector<double> s(3);
            s[0] = S[0]; s[1] = S[1]; s[2] = S[2];

            Su.push_back(s);
            if (q > 999)
            {
                break;
            }
        }
    }

    Mat Smat(Su.size(), 3, CV_64FC1);
    Vector2Mat(Su, Smat);
    Mat S1;
    S1 = (Smat.t())*Smat;

    Eigen::Matrix3d SunMat;
    cv_to_eigen(S1, SunMat);
    Eigen::JacobiSVD<Eigen::Matrix3d> svd(SunMat, Eigen::ComputeFullU | Eigen::ComputeFullV);
    Eigen::Vector3d SigularVal = svd.singularValues();
    Eigen::Matrix3d V = svd.matrixV();

    Eigen::Vector3d SunVector = Eigen::Vector3d(V(0, 0), V(1, 0), V(2, 0));
    double posAngle = atan2(SunVector(1), SunVector(0)) * (180 / PI);

    cout << "The orientation angle is:" << posAngle << " degree" << endl;

    waitKey(0);
    return 0;
}

void Vector2Mat(vector<vector<double>>src, Mat dst)
{
    assert(dst.rows == src.size());
    MatIterator_<double> it = dst.begin<double>();
    for (int i = 0; i < src.size(); i++)
    {
        for (int j = 0; j < src[0].size(); j++)
        {
            *it = src[i][j];
            it++;
        }
    }
}

void cv_to_eigen(const Mat& input, Eigen::Matrix3d& output)
{
    cv2eigen(input, output);
}

      由於瑞利模型是光進行一次散射的理想化模型，而真實的大氣偏振模式是多次散射的結果，所以為了提高計算精度，選取圖像的中心區域（半徑為200個像素）作為像素取樣板。在所選中心圖像區域中，任意選取兩個像素對應的E-矢量叉乘得到太陽矢量，這樣取點次數越多，計算精度越高，本次實驗選取1000對像素。本次實驗運行得到的AOP圖如下圖所示：

       計算得到的方位角為，該夾角為AOP圖中太陽子午線與地理北極之間的夾角，知道了這個夾角，就可以知道所在位置的方位，從而進行導航。

       本次實驗中，計算得到的太陽矢量不僅可以計算得到太陽位置的方位角，也可計算得到太陽位置的高度角。但是本例計算得到的高度角誤差較大，因此在不需要空間導航，僅在二維水平面中導航的情況下，本例基於E-矢量叉乘計算太陽矢量，從而計算航向角的方法是一種很好的基於物理意義的大氣偏振光導航方案。