Solar vector solving method based on photovoltaic array and its application in attitude estimation of navigation system

Abstract

Solar-powered unmanned aerial vehicles usually fly at high altitudes, and they are mainly powered by the photocells covering the body of unmanned aerial vehicles. Considering that the solar vector cannot be affected by the disturbing magnetic field and harmful acceleration, a unit solar vector solving method based on photovoltaic array is proposed in this paper. The photocells with different installation angles are selected to form the photovoltaic array. The solar vector is solved by the least-squares method on the basis of normalization by using the output currents of the photovoltaic array. For eliminating the influence of faults and reflected light on the solving of the solar vector, an adaptive least-squares unit solar vector solving method is proposed. In addition, a solar vector measuring device is designed in order to verify the effectiveness of the proposed methods. By employing the structural advantages of the device, the current generated by the reflected light of the sky can be solved according to the currents generated by all photocells of this device. Thus, the solved current generated by the reflected light of the sky is more accurate. Moreover, strapdown inertial navigation system/solar vector/global positioning system integrated navigation Kalman filtering algorithm is proposed, in which the adaptive least-squares unit solar vector solving method is applied to the measurement update of the filter. The effectiveness of the methods proposed in this paper is illustrated by some numerical and physical simulations.

Keywords

Photovoltaic array photocell current solar vector attitude estimation Kalman filtering solar-powered UAV

Introduction

For most unmanned aerial vehicles (UAVs), improving the accuracy and stability of the navigation system using relatively low-cost sensors is one of the crucial issues to be considered. The MARG (magnetic, angular rate, and gravity) sensors are widely used for attitude determination of UAVs or other vehicles.^1–4 An inertial measurement unit (IMU) is usually composed of gyroscope and accelerometer which are used to measure the angular rate and acceleration of a UAV, respectively. In addition, the accelerometer is also used to sense the gravity field simultaneously. The magnetometer is used to measure the geomagnetic field. However, accelerometers and magnetometers are sensitive to harmful acceleration and disturbing magnetic field in the process of measuring the gravity field and geomagnetic field, respectively. Therefore, it is necessary to develop sensors which are free from harmful acceleration and disturbing magnetic field or correspondinganti-disturbance algorithm for attitude determination of UAVs.

In the paper by Hemerly et al.,⁵ an extended Kalman filter for implementing an attitude and heading reference system (AHRS) with acceleration compensation by the global positioning system (GPS) was introduced. In the filter, the derivative of the velocity obtained by GPS is removed from the acceleration measured by the accelerometer which is also used to sense gravity field. For exploring anti-disturbance attitude measurement methods, an improved antenna array–based communication system⁶ was applied for attitude estimation by computing the line-of-sight (LOS) path between the base station and the UAV. In addition, a complete framework for attitude estimation by exploiting three-dimensional (3D) LOS vector obtained from the antenna array system was proposed. In order to improve the capability of resisting harmful acceleration, an attitude estimation method which is robust for harmful accelerative environment was designed using MARG sensors and GPS receiver in the work of No et al.⁷ The navigation system excludes accelerometer and estimates pseudo-attitude information from the velocity measured by the GPS receiver. In addition, a magnetic-pitch angle determination method which extracts pitch angle and yaw angle information from magnetic vector was proposed in order to complement the deviation of pseudo-attitude. The paper by Tong et al.⁸ proposed a novel adaptive extended Kalman filter which includes a multiplicative extended Kalman filter (MEKF) and a hidden Markov model (HMM) recognizer. Moreover, the HMM recognizer can obtain precise attitude information using MARG sensors in the circumstance of 3D motion tracking.

For the problem of attitude estimation of UAV using Kalman filter, the disturbance of harmful acceleration on the gravity field measured by accelerometer can be reduced by constructing an adaptive observation noise variance matrix. The adaptive adjustment of the observation noise variance matrix requires the estimation of disturbance. Disturbance-observer-based control (DOBC) and related methods have been researched and applied in various industrial sectors in the last four decades. An observation mechanism is employed to estimate the total disturbance in the DOBC method. The paper proposed by Chen et al.⁹ reviews many widely used linear and nonlinear disturbance/uncertainty estimation techniques and discusses various compensation techniques. For enhancing the performance of disturbance rejection, the modeling error is actively estimated by a Kalman filter using an active-model-based control scheme.¹⁰

Celestial navigation system (CNS) has been widely used for attitude and position estimation in the navigation system of UAVs.^11,12 In recent years, research works on the solar vector measurement using atmospheric polarization modes for attitude determination have achieved fruitful results.^13–15 Information on atmospheric polarization modes can also be obtained in the field of views obscured by barriers. For this reason, using the atmospheric polarization mode to obtain the solar vector has become a research direction. As shown in Figure 1, solar-powered UAV usually flies above clouds. With the continuous advancement of photovoltaic technology, solar-powered UAV can obtain more energy from the sun for flight. In addition, solar-powered UAVs have broad application prospects. Hence, it is necessary to develop the techniques of acquiring the solar vector through photocells for attitude determination of UAVs. The paper by Yang et al.¹⁶ proposed a method of solving the angle of the solar ray using the alternating short-circuit currents generated by two rotating photocells. This method demands a rotating device of photocells and has no fault tolerance for photocell failure.

Figure 1.

Solar-powered UAV flying at high altitudes.

Considering that the solar vector is free from harmful acceleration and disturbing magnetic field, a unit solar vector solving method based on photovoltaic array is proposed in this paper. The solar vector is solved by the least-squares method on the basis of normalization using the output currents of the photovoltaic array, which significantly simplifies the solving process of the solar vector. Moreover, an adaptive least-squares unit solar vector solving method is also proposed, which greatly improves the fault tolerance of the solar vector solving. In order to verify the effectiveness of the proposed methods, a physical device with photovoltaic array is designed in this paper. The current generated by the reflected light is solved according to the currents generated by all photocells of this device. Thus, the solved current generated by the reflected light is more reliable. In simulations of the proposed strapdown inertial navigation systems (SINS)/solar vector/GPS integrated navigation algorithm, the effectiveness of the proposed solar vector solving method and its auxiliary effect for attitude estimation are verified.

Mathematical model of solar vector

The relationship between two coordinate frames which refer to “B” and “N” is used to represent the attitude of vehicles:

“B” represents the body-fixed frame which is a sensor frame with each axis pointing along each of the IMU axes.

“N” represents the local earth-fixed reference coordinate frame taken from the East-North-Up (ENU) reference system. The XYZ axes are positive to geographic East, North, and Up directions, respectively.

Hence, $R_{N}^{B}$ and $R_{B}^{N}$ are the rotation matrixes constructed by the quaternions of the real attitudes of the UAV dynamic model in the numerical simulation. Correspondingly, the rotation matrixes $C_{N}^{B}$ and $C_{B}^{N}$ are constructed by estimated quaternions of the filter.

The model of solar vector in the “N” coordinate frame

The angle between the line which connects the earth and the sun and the equatorial plane of the earth will vary when the earth rotates around the sun. The rotation of the earth causes the point of direct sunray to move around the surface of the earth. In order to construct the model of the solar vector in the navigation coordinate frame, assume that the sun rotates around the earth. According to the results proposed by Chong and Wong,¹⁷ the model of the solar vector in the “N” coordinate frame can be derived.

The solar vector in geocentric inertial coordinate frame $O - X_{I} Y_{I} Z_{I}$ is depicted in Figure 2, and the model of the solar vector in geocentric inertial coordinate frame is given by

s_{I} = \frac{1}{1 - e \cos E} R_{x} (- i_{s}) R_{z} (- ω_{s}) s_{o}

(1)

where $s_{o} = (\begin{matrix} \cos E - e & \sin E \sqrt{(1 - e^{2})} & 0 \end{matrix})^{T}$ , $i_{s}$ is the obliquity of the ecliptic, $ω_{s}$ is the argument of perihelion of the earth revolution orbit, $E$ is the eccentric anomaly of the earth revolution orbit, $e$ is the eccentricity of the earth revolution orbit, and $R ()$ is the rotation matrix of rotating the corresponding angle around the corresponding axis.

Figure 2.

The solar vector in geocentric inertial coordinate frame.

According to Kepler’s equation

E - e \sin E = M = Ω_{s} (t - t_{s})

(2)

the real-time solution of the eccentric anomaly $E$ can be obtained. In equation (2), $M$ is the mean anomaly of the earth revolution orbit, $t_{s}$ is the moment when the earth passes perihelion, and $Ω_{s}$ is the earth revolution angular velocity.

Assume that the latitude and longitude of the local position are $l_{a}$ and $l_{o}$ , respectively. Consequently, the model of unit solar vector in the “N” coordinate frame is

s_{N} = R_{y} (\frac{π}{2}) R_{x} (\frac{π}{2}) R_{y} (- l_{a}) R_{z} (l_{o} + Ω_{e} t) s_{I}

(3)

where $Ω_{e}$ is the earth rotation angular velocity. Thus, as long as the current date, time, and the latitude and longitude can be determined, the unit solar vector in the navigation frame can be obtained.

The model of the unit solar vector in the “B” coordinate frame

The azimuth angle and elevation angle of the unit solar vector $s_{B}$ in the “B” coordinate frame are defined as $ξ$ and $ζ$ . The projection of the unit solar vector in the “B” coordinate frame is shown in Figure 3.

Figure 3.

The projection of the unit solar vector in the “B” coordinate frame.

According to the projection of the unit solar vector shown in Figure 3, the unit solar vector in the “B” coordinate frame¹⁷ is given by

s_{B} = [\begin{matrix} s_{x} \\ s_{y} \\ s_{z} \end{matrix}] = [\begin{matrix} \cos ζ \cos ξ \\ \cos ζ \sin ξ \\ \sin ζ \end{matrix}]

(4)

The unit solar vector is usually obtained by solving the azimuth angle $ξ$ and elevation angle $ζ$ of the solar vector in the “B” coordinate frame.¹⁶ This paper proposes a method of directly solving the unit solar vector by using the output currents of the photovoltaic array. The solving progress of the solar vector is significantly simplified.

Design of the unit solar vector solving method

Photocell is a special semiconductor diode which can convert light energy into electrical energy. In order to improve the anti-disturbance capability of attitude estimation, a unit solar vector solving method is proposed using the structural advantages of photovoltaic array and the principle of energy conversion. The short-circuit current output by the photovoltaic array is collected through the corresponding circuit. Usually, solar-powered UAVs are covered with photocells having the same area. Each of these photocells has a fixed installation angle. Thus, there is a fixed and known rotation matrix $R_{i}$ from the coordinate frame of photocell $i$ to the “B” coordinate frame. In order to reduce the influence of wing vibration on solar vector solving, the photocells with different installation angles at the wing root and airframe are selected to form the photovoltaic array. In addition, the photocells are labeled as 1 – N, and the rotation matrix $R_{i}$ is calculated according to the installation angle of photocell $i$ .

Least-squares unit solar vector

The rotation matrix $R_{i}$ from the coordinate frame of photocell $i$ to the “B” coordinate frame is known. Thus, the normal unit vector of photocell $i$ in the “B” coordinate frame can be presented as

n_{i} = R_{i} z_{i}

(5)

where $z_{i} = [\begin{matrix} 0 & 0 & 1 \end{matrix}]^{T}$ is the normal unit vector of photocell $i$ in the coordinate frame.

In the “B” coordinate frame, the angle between the normal unit vector $n_{i}$ of photocell $i$ and the unit solar vector $s_{B}$ is defined as $φ_{i}$ . Thus, the direction cosine of the vector $n_{i}$ and $s_{B}$ can be given by

\cos φ_{i} = 〈 n_{i}, s_{B} 〉 = 〈 R_{i} (3), s_{B} 〉

(6)

where the symbol 〈〉 represents the inner product of two vectors, and $R_{i} (3)$ is the third column of rotation matrix $R_{i}$ .

The current output by photocell $i$ contains the currents generated by the reflected light and the sunray, respectively. Assume that all photocells of the photovoltaic array are exposed to the sky. Thus, define $I_{i}^{F}$ and $I_{i}$ as the currents generated by the reflected light of the sky and the sunray, respectively. In addition, the currents generated by the reflected light of the sky can be treated as the same for different photocells. Thus, the current¹⁸ output by photocell $i$ can be represented as

I_{i}^{S} = I_{F} + I_{i}

(7)

where $I_{F}$ is the current generated by the reflected light of the sky.

The expression of current generated by the sunray¹⁸ is approximately

I_{i} = I_{⊥} \cos φ_{i}

(8)

where $I_{⊥}$ is the current generated by the sunray when the sunray is incident vertically to the photocell of a certain area. Its value is mainly related to the light intensity. For each photocell in the photovoltaic array, the value of $I_{⊥}$ is treated as the same.

According to equations (6)–(8), the current output by photocell $i$ can be expressed as

\begin{matrix} I_{i}^{S} & = I_{F} + I_{⊥} 〈 R_{i} (3), s_{B} 〉 \\ = f_{i} (I_{F}, I_{⊥}, s_{B} (1), s_{B} (2), s_{B} (3)) \end{matrix}

(9)

where constants $I_{F}$ and $I_{⊥}$ , and the vector $s_{B}$ are all unknown quantities, and $f_{i} ()$ is a function with five variables. It takes only five photocells in the photovoltaic array that the variable $I_{F}$ can be solved. Thus, the current generated by the sunray of photocell $i$ is known, and the value is

I_{i} = I_{i}^{S} - I_{F}

(10)

According to equations (6) and (8), the expression of current generated by the sunray of photocell $i$ is

\begin{matrix} I_{i} & = I_{⊥} 〈 R_{i} (3), s_{B} 〉 \\ = I_{⊥} [\begin{matrix} R_{i} (1, 3) & R_{i} (2, 3) & R_{i} (3, 3) \end{matrix}] s_{B} \end{matrix}

(11)

Assume that the photocells in the photovoltaic array are all exposed to the sunray and generate currents normally. Thus, the current vector output by the photovoltaic array is

\begin{matrix} I_{Σ} = {[\begin{matrix} I_{1} & I_{2} & \dots & I_{N} \end{matrix}]}^{T} = I_{⊥} R_{Σ} s_{B} \\ = I_{⊥} [\begin{matrix} R_{1} (1, 3) & R_{1} (2, 3) & R_{1} (3, 3) \\ R_{2} (1, 3) & R_{2} (2, 3) & R_{2} (3, 3) \\ ⋮ & ⋮ & ⋮ \\ R_{N} (1, 3) & R_{N} (2, 3) & R_{N} (3, 3) \end{matrix}] s_{B} \end{matrix}

(12)

where $R_{Σ}$ is the rotation matrix of the photovoltaic array.

There is a linear relationship between the current vector $I_{Σ}$ and the unit solar vector $s_{B}$ . However, matrix $R_{Σ}$ is irreversible. The unit solar vector $s_{B}$ cannot be directly solved. The installation angles of photocells in the photovoltaic array are different from each other. It can be inferred that $R_{Σ}$ must be a matrix with full column rank. In addition, $I_{⊥}$ cannot be directly measured. Hence, the least-squares unit solar vector can be solved according to $R_{Σ}$ and $I_{Σ}$ using the normalized idea. The least-squares unit solar vector is

s_{B} = \frac{{(R_{Σ}^{T} R_{Σ})}^{- 1} R_{Σ}^{T} I_{Σ}}{I_{⊥}} = \frac{{(R_{Σ}^{T} R_{Σ})}^{- 1} R_{Σ}^{T} I_{Σ}}{‖ {(R_{Σ}^{T} R_{Σ})}^{- 1} R_{Σ}^{T} I_{Σ} ‖}

(13)

According to equation (13), the relation of $I_{⊥}$ and $| | (R_{Σ}^{T} R_{Σ})^{- 1} R_{Σ}^{T} I_{Σ} | |$ satisfies

I_{⊥} = ‖ {(R_{Σ}^{T} R_{Σ})}^{- 1} R_{Σ}^{T} I_{Σ} ‖

(14)

Adaptive least-squares unit solar vector

During the flight of UAV, the angles between the normal vectors of some photocells and the solar vector may be greater than 90°, which prevents the sunray from reaching the photocells. Moreover, each photocell still outputs a certain amount of current due to the irradiation of the reflected light. Thus, it is necessary to determine whether the photocell is exposed to the sunray and adaptively construct the current vector $I_{Σ}$ and the rotation matrix $R_{Σ}$ of the photovoltaic array. In addition, the current $I_{F}$ generated by the reflected light of the sky can be obtained only under the condition that the photocells are exposed to sunray.

Suppose that the number $N$ of photocells constituting the photovoltaic array is sufficiently large. The current $I_{i}^{S} (i = 1, 2, \dots, N)$ output by the photocell array ranges from large to small using a sorting algorithm. Thus, the ordered current vector of the photovoltaic array is

I_{seq}^{S} = {[\begin{matrix} I_{m_{1}}^{S} & I_{m_{2}}^{S} & \dots & I_{m_{N}}^{S} \end{matrix}]}^{T}

(15)

where the serial numbers $m_{1}, m_{2}, \dots, m_{N}$ correspond to the labels of the photocells. Similarly, the ordered rotation matrix of the photovoltaic array is given as

R_{seq, Σ} = [\begin{matrix} R_{m_{1}} (1, 3) & R_{m_{1}} (2, 3) & R_{m_{1}} (3, 3) \\ R_{m_{2}} (1, 3) & R_{m_{2}} (2, 3) & R_{m_{2}} (3, 3) \\ ⋮ & ⋮ & ⋮ \\ R_{m_{N}} (1, 3) & R_{m_{N}} (2, 3) & R_{m_{N}} (3, 3) \end{matrix}]

(16)

Take the first five large currents from current vector $I_{seq}^{S}$ to form the current vector $I_{seq, 5}^{S}$ as follows

I_{seq, 5}^{S} = {[\begin{matrix} I_{m_{1}}^{S} & I_{m_{2}}^{S} & \dots & I_{m_{5}}^{S} \end{matrix}]}^{T}

(17)

Under the condition that the photocells work normally and the number of photocells $N$ is sufficiently large, the photocells outputting the current vector $I_{seq, 5}^{S}$ must be exposed to the sunray. Thus, the current $I_{F}$ generated by the reflected light of the sky can be solved according to equations (9) and (17). Moreover, the current vectors $I_{seq} = [\begin{matrix} I_{m_{1}} & I_{m_{2}} & \dots & I_{m_{N}} \end{matrix}]^{T}$ and $I_{seq, 5} = [\begin{matrix} I_{m_{1}} & I_{m_{2}} & \dots & I_{m_{5}} \end{matrix}]^{T}$ are obtained after the current $I_{F}$ is subtracted from the components of current vectors $I_{seq}^{S}$ and $I_{seq, 5}^{S}$ .

According to equation (13), the least-squares unit solar vector satisfies

s_{B}^{5} = \frac{{({(R_{seq, 5})}^{T} R_{seq, 5})}^{- 1} {(R_{seq, 5})}^{T} I_{seq, 5}}{‖ {({(R_{seq, 5})}^{T} R_{seq, 5})}^{- 1} {(R_{seq, 5})}^{T} I_{seq, 5} ‖}

(18)

where

R_{seq, 5} = [\begin{matrix} R_{m_{1}} (1, 3) & R_{m_{1}} (2, 3) & R_{m_{1}} (3, 3) \\ R_{m_{2}} (1, 3) & R_{m_{2}} (2, 3) & R_{m_{2}} (3, 3) \\ R_{m_{3}} (1, 3) & R_{m_{3}} (2, 3) & R_{m_{3}} (3, 3) \\ R_{m_{4}} (1, 3) & R_{m_{4}} (2, 3) & R_{m_{4}} (3, 3) \\ R_{m_{5}} (1, 3) & R_{m_{5}} (2, 3) & R_{m_{5}} (3, 3) \end{matrix}]

$s_{B}^{5}$ is the least-squares unit solar vector solved by $I_{seq, 5}$ and $R_{seq, 5}$ .

According to equations (5) and (18), the angle $φ_{m_{i}}$ $(i > 5)$ between the normal unit vector $n_{m_{i}}$ of photocell $m_{i}$ and the solved unit solar vector $s_{B}^{5}$ can be determined. Thus, according to equations (11) and (14), the theoretical value of the current generated by photocell $m_{i}$ with $φ_{m_{i}} < 90 °$ is given as

\begin{matrix} I_{m_{i}}^{p} = ‖ {({(R_{s e q, 5})}^{T} R_{s e q, 5})}^{- 1} {(R_{s e q, 5})}^{T} I_{s e q, 5} ‖ \\ \times [\begin{matrix} R_{m_{i}} (1, 3) & R_{m_{i}} (2, 3) & R_{m_{i}} (3, 3) \end{matrix}] s_{B}^{5} \end{matrix}

(19)

In addition, the criterion of determining whether the photocell labeled m_i (i > 5) is available is

| I_{m_{i}}^{p} - I_{m_{i}} | \leq τ

(20)

where $τ$ is an empirically given threshold. When the current generated by photocell $m_{i}$ does not satisfy equation (20), the serial number of photocell $m_{i - 1}$ is saved as $m_{η}$ .

Therefore, an adaptive photovoltaic array can be established according to the photocells labeled $m_{1}, \dots, m_{η}$ . In addition, the current vector generated by the photocells is defined as $I_{seq, η}$ . According to equation (13), the adaptive least-squares unit solar vector is

s_{B}^{a} = \frac{{({(R_{seq, η})}^{T} R_{seq, η})}^{- 1} {(R_{seq, η})}^{T} I_{seq, η}}{‖ {({(R_{seq, η})}^{T} R_{seq, η})}^{- 1} {(R_{seq, η})}^{T} I_{seq, η} ‖}

(21)

where the rotation matrix $R_{seq, η}$ of the photovoltaic array is given as

R_{seq, η} = [\begin{matrix} R_{m_{1}} (1, 3) & R_{m_{1}} (2, 3) & R_{m_{1}} (3, 3) \\ R_{m_{2}} (1, 3) & R_{m_{2}} (2, 3) & R_{m_{2}} (3, 3) \\ ⋮ & ⋮ & ⋮ \\ R_{m_{η}} (1, 3) & R_{m_{η}} (2, 3) & R_{m_{η}} (3, 3) \end{matrix}]

(22)

Solar vector measuring device and system model

The solar vector and the gravity field are used to assist SINS for attitude estimation of UAV. GPS is used to assist SINS for position and velocity estimation of UAV.

Solar vector measuring device and its model

The temperature drift of each photocell and the difference between each photocell are ignored. In order to verify the validity of the designed unit solar vector solving method, a solar vector measuring device shown in Figure 4 is designed in this paper. In addition, the device is placed on the roof, and all photocells are exposed to the clear sky. Seven photocells with the same area are uniformly laid on seven surfaces of the six frustum planes and are labeled as 1–7, respectively. The angle between each side face and the bottom surface is $α = 40 °$ . The angle between each side of the bottom surface is $β = 60 °$ . Moreover, the acquisition circuit of short-circuit current is designed for the photocells. Since the light is too strong during the test, the short-circuit currents are somewhat distorted. Thus, half of each photocell is blocked with white paper.

Figure 4.

Solar vector measuring device placed on the roof.

Note that $\cos β = - \cos 2 β$ and $\sin β = \sin 2 β$ . Under the condition that the sunray can reach all the photocells of the solar vector measuring device, the rotation matrix of the photovoltaic array is constructed as

\begin{matrix} R_{Σ, 7} & = [\begin{matrix} R_{1} (1, 3) & R_{1} (2, 3) & R_{1} (3, 3) \\ R_{2} (1, 3) & R_{2} (2, 3) & R_{2} (3, 3) \\ ⋮ & ⋮ & ⋮ \\ R_{7} (1, 3) & R_{7} (2, 3) & R_{7} (3, 3) \end{matrix}] \\ = [\begin{matrix} 0 & - \sin α & \cos α \\ \sin α \sin β & - \sin α \cos β & \cos α \\ \sin α \sin β & \sin α \cos β & \cos α \\ 0 & \sin α & \cos α \\ - \sin α \sin β & \sin α \cos β & \cos α \\ - \sin α \sin β & - \sin α \cos β & \cos α \\ 0 & 0 & 1 \end{matrix}] \end{matrix}

(23)

According to equation (6), the direction cosine of the normal unit vector $n_{i} (i = 1, 2, \dots, 7)$ and $s_{B}$ is given as

\begin{matrix} {\begin{matrix} \cos φ_{1} = - s_{B} (2) \sin α + s_{B} (3) \cos α \\ \cos φ_{2} = s_{B} (1) \sin α \sin β - s_{B} (2) \sin α \cos β \\ + s_{B} (3) \cos α \\ \cos φ_{3} = s_{B} (1) \sin α \sin β + s_{B} (2) \sin α \cos β \\ + s_{B} (3) \cos α \\ \cos φ_{4} = s_{B} (2) \sin α + s_{B} (3) \cos α \\ \cos φ_{5} = - s_{B} (1) \sin α \sin β + s_{B} (2) \sin α \cos β \\ + s_{B} (3) \cos α \\ \cos φ_{6} = - s_{B} (1) \sin α \sin β - s_{B} (2) \sin α \cos β \\ + s_{B} (3) \cos α \\ \cos φ_{7} = s_{B} (3) \end{matrix} \end{matrix}

(24)

According to equation (24), the following equation can be obtained

{\begin{matrix} \cos φ_{1} + \cos φ_{4} = 2 \cos α \cos φ_{7} \\ \cos φ_{2} + \cos φ_{5} = 2 \cos α \cos φ_{7} \\ \cos φ_{3} + \cos φ_{6} = 2 \cos α \cos φ_{7} \end{matrix}

(25)

According to equations (7), (8), and (25), the following equation can be obtained

{\begin{matrix} I_{1}^{S} + I_{4}^{S} - 2 I_{F} = 2 (I_{7}^{S} - I_{F}) \cos α \\ I_{2}^{S} + I_{5}^{S} - 2 I_{F} = 2 (I_{7}^{S} - I_{F}) \cos α \\ I_{3}^{S} + I_{6}^{S} - 2 I_{F} = 2 (I_{7}^{S} - I_{F}) \cos α \end{matrix}

(26)

The average value of $I_{F}$ obtained from the three expressions in equation (26) is

I_{F} = \frac{I_{1}^{S} + I_{2}^{S} + I_{3}^{S} + I_{4}^{S} + I_{5}^{S} + I_{6}^{S} - 6 I_{7}^{S} \cos α}{6 (1 - \cos α)}

(27)

In addition, the model of the solar vector measuring device can be simply constructed as

I_{Σ, 7} = {[\begin{matrix} I_{1} & I_{2} & \dots & I_{7} \end{matrix}]}^{T} = I_{⊥} R_{Σ, 7} s_{B, N} + n_{I}

(28)

where $s_{B, N} = R_{N}^{B} s_{N}$ ; $n_{I}$ represents the random white noise drift vector with variance $σ_{I}^{2}$ .

The model of IMU sensors and GPS

It is assumed that the IMU sensors have been calibrated. Thus, the models of gyroscope and accelerometer can be established as follows

{\begin{matrix} ω_{B} = ω_{t} + ε_{r}^{ω} + n_{ω} \\ a_{B} = R_{N}^{B} g_{N} + a_{v} + b_{a} + ε_{r}^{a} + n_{a} \end{matrix}

(29)

where $ω_{B}$ and $a_{B}$ are measuring outputs of gyroscope and accelerometer in the “B” coordinate frame, respectively; $ω_{t}$ is the actual angular velocity of the UAV dynamic model; $g_{N}$ is the gravity field in the “N” coordinate frame; $a_{v}$ and $b_{a}$ are vehicle linear acceleration and bias of accelerometer obtained by the accelerometer calibration, respectively; $n_{ω}$ and $n_{a}$ are the random white noise drift vectors of gyroscope and accelerometer, respectively; and $ε_{r}^{ω}$ and $ε_{r}^{a}$ are first-order Markov process drifts of the gyroscope and accelerometer, respectively. The dynamic models of $ε_{r}^{ω}$ and $ε_{r}^{a}$ are as follows

{\begin{matrix} {\overset{\cdot}{ε}}_{r}^{ω} = - \frac{1}{T_{ω}} ε_{r}^{ω} + ω_{r} \\ {\overset{\cdot}{ε}}_{r}^{a} = - \frac{1}{T_{a}} ε_{r}^{a} + ω_{a} \end{matrix}

(30)

where $T_{ω}$ and $T_{a}$ are the correlation time, and $ω_{r}$ and $ω_{a}$ are the random white noise drifts of the first-order Markov process.

The position and velocity model in the “N” coordinate frame obtained by GPS is established as follows

{\begin{matrix} p^{G} = p_{t}^{G} + n_{p} \\ v^{G} = v_{t}^{G} + n_{v} \end{matrix}

(31)

where $p_{t}^{G}$ and $v_{t}^{G}$ are the actual position and velocity of the UAV dynamic model in the “N” coordinate frame, respectively, and $n_{p}$ and $n_{v}$ are the random white noise drift vectors with variances of $σ_{p}^{2}$ and $σ_{v}^{2}$ , respectively.

State equation of UAV dynamic model

The motion states of UAV usually include attitude quaternions, positions, and velocities in 3D space. Thus, the solar vector of the UAV dynamic model is constructed as

x = {[\begin{matrix} q^{T} & p^{T} & v^{T} \end{matrix}]}_{10 \times 1}^{T}

(32)

where $q$ , $p$ , and $v$ are attitude quaternions, positions, and velocities of the UAV, respectively.

In order to make the numerical simulation more close to the actual situation of hardware, the dynamic model of UAV is established according to the measurement data of sensor models which have been established in equation (29). Thus, the dynamic model of attitude of UAV is given by

\overset{\cdot}{q} = \frac{1}{2} Ω q

(33)

where

Ω = [\begin{matrix} 0 & - ω_{B}^{x} & - ω_{B}^{y} & - ω_{B}^{z} \\ ω_{B}^{x} & 0 & ω_{B}^{z} & - ω_{B}^{y} \\ ω_{B}^{y} & - ω_{B}^{z} & 0 & ω_{B}^{x} \\ ω_{B}^{z} & ω_{B}^{y} & - ω_{B}^{x} & 0 \end{matrix}]

(34)

The linear acceleration of UAV measured by the accelerometer model in the “N” coordinate frame is given as

a_{N} = C_{B}^{N} (a_{B} - b_{a}) - g_{N}

(35)

which is taken as the input of the second-order integral dynamic model of UAV. The second-order integral dynamic model is given as

{\begin{matrix} \overset{\cdot}{p} = v \\ \overset{\cdot}{v} = a_{N} \end{matrix}

(36)

According to equations (33), (35), and (36), the state equation of UAV based on sensor models can be obtained as

\overset{\cdot}{x} = Ax + B a_{N}

(37)

where

{\begin{matrix} A = {[\begin{matrix} 1 / 2 Ω & 0_{4 \times 3} & 0_{4 \times 3} \\ 0_{3 \times 4} & 0_{3 \times 3} & I_{3 \times 3} \\ 0_{3 \times 4} & 0_{3 \times 3} & 0_{3 \times 3} \end{matrix}]}_{10 \times 10} \\ B = {[\begin{matrix} 0_{7 \times 3} \\ I_{3 \times 3} \end{matrix}]}_{10 \times 3} \end{matrix}

(38)

Observation matrix of the state

The observed unit of the solar vector and gravity field is

y_{k} = [\begin{matrix} C_{N}^{B} \\ C_{N}^{B} \end{matrix}] [\begin{matrix} s_{N} \\ g_{N} \end{matrix}] = h (q)

(39)

Expand $h (q)$ by Taylor at the state of quaternions. The observation matrix of the solar vector and gravity field is

H_{sg, k} = {\frac{\partial h (q)}{\partial q} |}_{q = {\hat{q}}_{k | k - 1}}

(40)

where ${\hat{q}}_{k | k - 1}$ is the estimated quaternion after time update.

Since the filter also observes the position and velocity information obtained by GPS, the observation matrix of the filter is constructed as

H_{k} = {[\begin{matrix} H_{sg, k} & 0_{6 \times 6} \\ 0_{6 \times 4} & I_{6 \times 6} \end{matrix}]}_{12 \times 10}

(41)

SINS/solar vector/GPS integrated navigation Kalman filtering algorithm

Before the measurement update of the Kalman filtering algorithm, the unit solar vector needs to be solved in real time through the output current vector of the photovoltaic array. The detailed steps of the algorithm are given as follows:

Time update

Since the state equation (37) is a continuous time system, equation (37) is discretized by fourth-order Runge–Kutta method in order to update the state in discrete time. Then the time update of the state can be obtained

{\begin{matrix} S_{1} = A {\hat{x}}_{k - 1} + B a_{N} \\ S_{2} = A ({\hat{x}}_{k - 1} + S_{1} T_{h} / 2) + B a_{N} \\ S_{3} = A ({\hat{x}}_{k - 1} + S_{2} T_{h} / 2) + B a_{N} \\ S_{4} = A ({\hat{x}}_{k - 1} + S_{3} T_{h}) + B a_{N} \end{matrix}

(42)

{\begin{matrix} dx = (S_{1} + 2 S_{2} + 2 S_{3} + S_{4}) T_{h} / 6 \\ {\hat{x}}_{k / k - 1} = {\hat{x}}_{k - 1} + dx \end{matrix}

(43)

where $T_{h}$ is the sampling period, ${\hat{x}}_{k - 1}$ is the state estimation at the moment $k - 1$ , and ${\hat{x}}_{k / k - 1}$ is the state estimation after the time update of the filter.

The time update of state estimation covariance matrix is given as

{\begin{matrix} dP = (A P_{k - 1} + P_{k - 1} A^{T} + B R_{u} B^{T} + Q) T_{h} \\ P_{k / k - 1} = P_{k - 1} + dP \end{matrix}

(44)

where $R_{u}$ is the input noise variance matrix, $Q$ is process noise variance matrix of state equation (37), $P_{k - 1}$ is the state estimation covariance matrix at the moment $k - 1$ , and $P_{k / k - 1}$ is the covariance matrix after the time update of the filter.

Solving of the unit solar vector

According to equations (15)–(20), the ordered current vector $I_{seq, η}$ and the ordered rotation matrix $R_{seq, η}$ of the photovoltaic array can be obtained. Thus, according to equation (21), the adaptive unit solar vector is given as

s_{B}^{a} = \frac{{({(R_{seq, η})}^{T} R_{seq, η})}^{- 1} {(R_{seq, η})}^{T} I_{seq, η}}{‖ {({(R_{seq, η})}^{T} R_{seq, η})}^{- 1} {(R_{seq, η})}^{T} I_{seq, η} ‖}

(45)

Measurement update

The measurement update is given as

{\begin{matrix} K = P_{k / k - 1} H_{k}^{T} {(H_{k} P_{k / k - 1} H_{k}^{T} + R)}^{- 1} \\ {\hat{x}}_{k} = {\hat{x}}_{k / k - 1} + K (Z_{k} - H_{k} {\hat{x}}_{k / k - 1}) \\ P_{k} = (I - K H_{k}) P_{k / k - 1} \end{matrix}

(46)

where $K$ is the gain matrix of the Kalman filtering algorithm and $R$ is the measurement noise variance matrix. In addition, the measurement information $Z_{k}$ is

Z_{k} = [{\begin{matrix} {(s_{B}^{a})}^{T} & {(a_{B} - b_{a} - a_{v})}^{T} & {(p^{G})}^{T} & {(v^{G})}^{T}] \end{matrix}}^{T}

(47)

where the motion acceleration $a_{v}$ is solved by the measured velocity differentiation of UAV according to the method proposed by Hemerly et al.,⁵ and $s_{B}^{a}$ is the adaptive least-squares unit solar vector in the “B” coordinate frame.

Numerical simulation and physical simulation

Numerical simulation

The algorithms designed in this paper are simulated by MATLAB. The simulation step size is given as $Δ t = 10^{- 2} s$ . The simulation time is given as $t = 200 s$ . The unit solar vector and the unit gravity field vector in the “N” coordinate frame are given as $s_{N} = [\begin{matrix} 1 / 2 & 0 & \sqrt{3} / 2 \end{matrix}]^{T}$ and $g_{N} = [\begin{matrix} 0 & 0 & - 1 \end{matrix}]^{T}$ , respectively. The value of $I_{⊥}$ is given as 10 mA. The standard deviation of the random white noise of the solar vector measuring device is given as $σ_{I} = 0.03 mA$ . The random white noise and first-order Markov process drift of the gyroscope are given as 0.002 and 0.001 rad/s, respectively, and the correlation time of the Markov process drift is given as $T_{w} = 7200 / s$ . The random white noise and first-order Markov process drift of the accelerometer are given as 0.8mg and 2mg, respectively, and the correlation time of the Markov process drift is given as $T_{a} = 3600 s$ . The random white noise drifts of position and velocity obtained by GPS are given as 2 m and 0.3 m/s. Assume that the UAV is a multirotor, and the multirotor flies at a fixed altitude. For the convenience of simulation, suppose that the motion accelerations in the x-direction and the y-direction are proportional to the pitch and roll angles of UAV, respectively. Assume that the yaw angle of UAV is fixed; the simulation curves of pitch and roll angles of UVA are given in Figure 5.

Figure 5.

The simulation curves of pitch and roll angles of UAV.

The position and velocity estimation errors of the filter are shown in Figures 6 and 7; the attitude estimation errors of the filter assisted by the gravity field and the solar vector are shown in Figure 8; and the attitude estimation errors of the filter assisted only by the gravity field are shown in Figure 9. According to the results shown in Figures 6 and7, the estimation errors of position and velocity converge to a relatively small range. As shown in Figure 8, the attitude estimation errors of the filter are within 0.5°. As shown in Figure 9, the attitude estimation errors of the pitch and roll angles are all within 0.5°, but the estimation errors of the yaw angle are divergent. Thus, according to the simulation results of Figures 8 and 9, we can see the assistant effect of the solar vector on attitude estimation.

Figure 6.

Position estimation errors of the filter in numerical simulation.

Figure 7.

Velocity estimation errors of the filter in numerical simulation.

Figure 8.

Attitude estimation errors of the filter assisted by the gravity field and the solar vector in numerical simulation.

Figure 9.

Attitude estimation errors assisted only by the gravity field in numerical simulation.

Physical simulation

The simulation step size and simulation time are also given as $Δ t = 10^{- 2} s$ and $t = 200 s$ , respectively. On a bright afternoon, the solar vector measuring device shown in Figure 4 is placed horizontally on the roof. Moreover, the photocells of the device are all exposed to the sky and the sunray. The collected currents of photocells are also filtered by a digital low-pass filter, and then the currents are sent to MATLAB. The collected currents of photocells through the solar vector measuring device are shown in Figure 10. The current generated by the reflected light of the sky is solved according to equation (27). The currents of photocells generated by the sunray are shown in Figure 11. Figures 10 and 11 show that the seven photocells of the solar vector measuring device have the same characteristics according to the variation trend of the current curves. Therefore, it can be inferred that the design of the solar vector measuring device is reasonable to some extent.

Figure 10.

The collected currents of photocells through the solar vector measuring device.

Figure 11.

The currents of photocells which have removed the current $I_{F}$ .

Since the IMU sensors and magnetometer in Pixhawk flight controller of an UAV shown in Figure 12 have been corrected, export data of these sensors for physical simulation. The UAV with Pixhawk flight controller is placed on the ground, and it is assumed that the UAV model in numerical simulation is also in a static state. The transmission speed of IMU and magnetometer data in the MAVLink module of the PX4 source code of the Pixhawk flight controller is changed to 100 Hz. The data of IMU and magnetometer are collected to log file through ground control station for more than 200 s, and then the log file is read through MATLAB for 200 s. The collected IMU data are simulated together with the current vector output by the solar vector measuring device.

Figure 12.

The UAV with Pixhawk flight controller.

Since various errors exist in the solar vector measuring device, the projection in the “N” coordinate system of the solved unit solar vector by the solar vector measuring device at the initial moment is used as the unit solar vector in the “N” coordinate frame. In addition, the Euler angles of UAV can be solved according to the collected data of accelerometer and magnetometer when the UAV is still. Hence, the data of IMU and magnetometer are collected once at test 1 when the yaw, pitch, and roll angles of UAV are $ψ_{1} = 0.507 °$ , $θ_{1} = - 0.718 °$ , and $γ_{1} = 1.401 °$ , respectively, and the data of IMU and magnetometer are collected once again at test 2 when the yaw, pitch, and roll angles of UAV are $ψ_{2} = 17.321 °$ , $θ_{2} = - 22.023 °$ , and $γ_{2} = 14.248 °$ , respectively.

For determining the assistant effect of solar vector on attitude estimation, the following experiments are carried out. The attitude estimation errors obtained by gyroscope integration are shown in Figures 13 and 14, when the collected IMU data come from tests 1 and 2, respectively. The attitude estimation errors of the filter assisted by the gravity field and the solar vector are shown in Figures 15 and 16, when the collected IMU data come from tests 1 and 2, respectively. The attitude estimation errors of the filter assisted only by the solar vector are given in Figures 17 and 18, when the collected IMU data come from tests 1 and 2, respectively. The attitude estimation errors of the filter assisted only by the gravity field are presented in Figures 19 and 20, when the collected IMU data come from tests 1 and 2, respectively.

Figure 13.

Attitude estimation errors obtained by gyroscope integration with IMU data coming from test 1.

Figure 14.

Attitude estimation errors obtained by gyroscope integration with IMU data coming from test 2.

Figure 15.

Attitude estimation errors assisted by the gravity field and the solar vector with IMU data coming from test 1.

Figure 16.

Attitude estimation errors assisted by the gravity field and the solar vector with IMU data coming from test 2.

Figure 17.

Attitude estimation errors assisted only by the solar vector with IMU data coming from test 1.

Figure 18.

Attitude estimation errors assisted only by the solar vector with IMU data coming from test 2.

Figure 19.

Attitude estimation errors assisted only by the gravity field with IMU data coming from test 1.

Figure 20.

Attitude estimation errors assisted only by the gravity field with IMU data coming from test 2.

It can be seen from Figures 17 and 18 that the cumulative errors of attitude estimation assisted by the solar vector are smaller than those of gyroscope integration given in Figures 13 and 14. As presented in Figures 15 and 16, the cumulative errors of yaw angles are smaller than those shown in Figures 19 and 20. Therefore, according to the comparison of attitude estimation errors illustrated in Figures 13 – and 20, we can see that the solar vector has an assistant effect on attitude estimation. Consequently, the numerical and physical simulations have verified that the solar vector solving method is effective and has an assistant effect on attitude estimation.

Conclusion

MARG sensors in navigation system of UAV are vulnerable to disturbing magnetic field and harmful acceleration. Considering that the solar-powered UAVs usually fly at high altitudes, the unit solar vector and adaptive unit solar vector solving methods are proposed in this paper. The unit solar vector is solved by the least-squares method. Thus, the obtained unit solar vector has higher precision and fault tolerance when there exist errors in photovoltaic array and acquisition circuit of short-circuit currents. For solar-powered UAVs flying at high altitudes, photocells selected from the surface of the body can be used to assist attitude estimation of the UAV. This method has an assistant effect on attitude estimation without adding too much extra weight while the photocells are providing energy to the UAV. Moreover, the solar vector cannot be affected by the disturbing magnetic field and harmful acceleration. Thus, it can be concluded that the unit solar vector solving methods proposed in this paper have certain value in engineering application.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is partially supported by the National Natural Science Foundation of China (61703286) and the series project of Liaoning Province Education Department (L201711).

ORCID iD

Huaqing Zhang

References

Wang

Effective adaptive Kalman filter for MEMS-IMU/magnetometers integrated attitude and heading reference systems. J Navig 2013; 66(1): 99–113.

Zhang

Jia

Chen

, et al. Attitude solution for strapdown inertial navigation system based on quaternion algorithm. Optic Precis Eng 2008; 16(10): 1963–1970 (in Chinese).

Grivon

Vezzetti

Violante

MG.

Development of an innovative low-cost MARG sensors alignment and distortion compensation methodology for 3D scanning applications. Robot Auton Syst 2013; 61(12): 1710–1716.

Liu

, et al. Attitude measurement and optimization for horizontal directional drilling based on MARG sensor. Adv Mat Res 2013; 616–618: 455–460.

Hemerly

Maciel

BCO

Milhan

, et al. Attitude and heading reference system with acceleration compensation. Aircr Eng Aerosp Tec 2012; 84(2): 87–93.

Cordeiro

TFK

Da Costa

JPCL

de Sousa

, et al. Improved Kalman-based attitude estimation framework for UAVs via an antenna array. Digit Signal Process 2016; 59: 49–65.

Cho

Kee

Attitude estimation method for small UAV under accelerative environment. GPS Solut 2015; 19(3): 343–355.

Tong

Han

, et al. Adaptive EKF based on HMM recognizer for attitude estimation using MEMS MARG sensors. IEEE Sens J 2018; 18(8): 3299–3310.

Chen

Yang

Guo

, et al. Disturbance—observer—based control and related methods—an overview. IEEE Trans.Ind Electron 2016; 63(2): 1083–1095.

10.

Liang

, et al. Active-model-based control for the quadrotor carrying a changed slung load. Electronics 2019; 8(4): 461.

11.

Ushaq

Cheng

Jamshaid

A fault tolerant SINS/GPS/CNS integrated navigation scheme realized through federated Kalman filter. Appl Mech Mater 2013; 332: 104–110.

12.

Min

Xiong

Xing

, et al. An improved SINS/GNSS/CNS federal filter based on dual quaternions. Acta Armament 2018; 39(2): 315–324 (in Chinese).

13.

Wang

Chu

Zhang

, et al. A novel autonomous real-time position method based on polarized light and geomagnetic field. Sci Rep 2015; 5: 9725.

14.

Ren

Jun

Tang

, et al. Polarized skylight pattern-based approach to attitude determination. IEEE Sensors J 2015; 15(9): 4917–4927.

15.

Fan

Chen

Wang

, et al. Three-dimensional attitude information obtained by the skylight polarization pattern. Optic Precis Eng 2016; 24(6): 1248–1256 (in Chinese).

16.

Yang

Shi

Liu

, et al. Research on measurement method of the large range sunray incident angle. Acta Energ Solar Sin 2011; 32(2): 216–219 (in Chinese).

17.

Chong

Wong

CW.

General formula for on-axis sun-tracking system and its application in improving tracking accuracy of solar collector. Sol Energy 2009; 83(3): 298–305.

18.

Santoni

Bolotti

Attitude determination of small spinning spacecraft using three axis magnetometer and solar panels data. In: Proceedings of the 2000 IEEE aerospace conference, Big Sky, MT, 25 March 2000, pp. 127–133. Piscataway, NJ: IEEE.