2. TRANSFORMATION FROM IMAGE SPACE TO OBJECT SPACE

Figure 2 shows the transformation mechanism between image space and object space when taking observations using the fisheye camera. In Figure 2, P is the point in image space, P′ is the point in object space, O′ is the intersection between the principal optic axis and image plane, named the principal point. O – X _c Y _c Z _c is defined as the camera coordinate system, in which the origin O represents the optical centre of the camera, axis Z _c coincides with the principal optic axis, axis X _c is perpendicular to axis Z _c and parallel to axis x of the image plane coordinate system and finally axis Y _c is determined by the right-hand rule. As is demonstrated, A _c is defined as the intersection angle between the vertical plane through P ^′ and the plane X _c OZ _c is named azimuth, and θ_c is defined as the intersection angle between the vector through P ^′ and the axis Z _c, named the semi-angular field. Theoretically, point P ^′ could be determined with A _c and θ_c.

Figure 2. Principle of lens projection.

To achieve precise conversion between coordinates in object space and image space, a projection model and a distortion model of the fisheye lens must be constructed. Normally, fisheye projection models include an equidistance projection model, an equisolid projection model, a stereoscopic projection model, and an orthogonal projection model, etc. (Wang, Reference Wang2006). In this paper, the lens chosen in experiments complies with the equisolid projection model. In addition, unlike in Gaussian optics, the “non-similar” theory is applied to the fisheye lens to meet the requirements of a super-wide FOV. Therefore, a distortion model is needed. Normally, the optical distortion of the fisheye lens could be divided into radial distortion, eccentric distortion and in-plane distortion. Among these, the eccentric and in-plane distortions are in the magnitude of 10⁻⁵ ~ 10⁻⁷ when compared with radial distortion, and thus they can be neglected. In order to describe radial distortion, a fourth order polynomial is used in this paper, which is shown as follows (Yuan, Reference Yuan2012):

(1)$$\theta_{c} =2\arcsin {r \over 2f}+k_{1} \left(\arcsin {r \over 2f}\right)^{2}+k_{2} \left(\arcsin {r \over 2f}\right)^{3}+k_{3} \left(\arcsin {r \over 2f}\right)^{4}$$

Where, θ_c denotes the semi-angular field of object point P ^′,  f denotes the focal length of the fisheye camera, (k ₁, k ₂, k ₃) denotes the radial distortion parameters, r denotes the distance between image point P and principal point O ^′ in the image plane, which is defined as

(2)$$r=\sqrt{(x-x_{0})^{2}+(y-y_{0})^{2}}$$

(x, y) denotes image point P’s coordinates in the image plane, and (x ₀, y ₀) denotes principal point O ^′’s coordinates in the image plane.

In the image plane, the intersection angle between vector P O ^′ and axis x is equivalent to the azimuth of P ^′ in the camera coordinate system. Therefore, A _c can be calculated as

(3)$$A_{c} =\left\{ \matrix{\arctan{y-y_{0} \over x_{0} -x} \hfill& x_{0} -x>0,y-y_{0} >0 \cr \arctan{y-y_{0} \over x_{0} -x}+\rpi \hfill& x_{0} -x< 0 \hfill\cr \arctan{y-y_{0} \over x_{0} -x}+2\rpi & x_{0} -x>0,y-y_{0} < 0 }\right.$$

(A _c, θ _c) could be converted to Three-Dimensional (3D) Cartesian coordinates as follows:

(4)$${\bf S}_{c} =\left[\matrix{ x_{c} \cr y_{c} \cr z_{c} } \right]=\left[\matrix{ \sin \theta_{c} \cos A_{c} \cr \sin \theta_{c} \sin A_{c} \cr \cos \theta_{c} } \right]$$

3. EXTRACTION OF THE MOON CENTRE IN THE IMAGE PLANE

The average flattening of the Moon is approximately 1/3,476, and the lunar geometric centre does not coincide with its centroid because the centroid is biased about 2 km towards the earth (Ouyang, Reference Ouyang2005; Li et al., Reference Li, Jing and Huang2012b). However, due to the long distance between the Moon and the Earth, the Moon disk observed on the Earth could be approximately considered as circular, with an angular radius of about 15′.

3.1. Algorithm to fit the Moon centre

Generally, a circular target will present an approximate ellipse on the image plane when shot by a fisheye camera. Furthermore, with the increase of the semi-angular field, the target's image will get closer to the boundary of the FOV, and the ellipticity will become more obvious. Similarly, the lunar edge line is certain to form an approximate ellipse on the image plane. To fit the centre of the Moon, an algorithm—Direct Least Squares Fitting of Ellipses (DLS), is introduced, whose fundamental idea is briefly conveyed as follows (Fitzgibbon et al., Reference Fitzgibbon, Pilu and Fischer1996).

The general Equation of an ellipse is

(5)$$F(x,y)=ax^{2}+bxy+cy^{2}+dx+ey+f=0,b^{2}-4ac< 0$$

The inequality constraint b ² − 4ac<0 could be changed as an equality constraint 4ac − b ² = 1.

By introducing vectors

(6)$${\bf A}=(a,b,c,d,e,f)^{T}\quad {\bf P}=(x^{2},xy,y^{2},x,y,1)$$

F(x, y) could be rewritten as the vector form

(7)$$F_{{\bf A}} ({\bf P})={\bf PA}$$

| F _A (P_i)| is defined as the algebraic distance of the point (x _i, y _i) to the fitting ellipse F(x, y) = 0, which reflects the fitting residual error at the point (x _i, y _i). Then a vector formed by algebraic distances of all the fitting participants can be created, whose norm can be regarded as a fitting accuracy index. The specific vector A minimising the norm mentioned above is the direct least squares fitting value of elliptic coefficients.

Suppose the number of points participating in elliptical fitting is M, the coordinates of the points are (x ₁, y ₁)…(x _M, y _M). By introducing matrices

(8)$${\bf D}=\left(\matrix{ x_{1}^{2} & x_{1} y_{1} & y_{1}^{2} & x_{1} & y_{1} & 1 \cr \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \cr x_{i}^{2} & x_{i} y_{i} & y_{i}^{2} & x_{i} & y_{i} & 1 \cr \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \cr x_{M}^{2} & x_{M} y_{M} & y_{M}^{2} & x_{M} & y_{M} & 1 }\right)\quad {\bf C}=\left(\matrix{ 0 & 0 & 2 & 0 & 0 & 0 \cr 0 & {-1} & 0 & 0 & 0 & 0 \cr 2 & 0 & 0 & 0 & 0 & 0 \cr 0 & 0 & 0 & 0 & 0 & 0 \cr 0 & 0 & 0 & 0 & 0 & 0 \cr 0 & 0 & 0 & 0 & 0 & 0 \cr }\right)$$

the fitting assignment can be rewritten as (Halif and Flusser, Reference Halif and Flusser1998)

(9)$$\min\limits_{{\bf A}} \left\| {\bf DA} \right\|\;\;st\;\;{\bf A}^{{\bf T}}{\bf CA}=1$$

which can be solved by generalised eigenvectors

(10)$$\left\{ \matrix{{\bf D}^{{\bf T}}{\bf DA}=\lambda {\bf CA} \cr {\bf A}^{{\bf T}}{\bf CA}=1 \hfill}\right.$$

It can be proved that Equation (10) has a sole positive eigenvalue λ, whose corresponding eigenvector is exactly the coefficient vector of the best-fit ellipse for the given set of points (Gander, Reference Gander1981).

3.2. Method to vote in the real edge line of the Moon phase

The Moon has a periodic phase change when observed from the Earth. The Moon phase is composed of a twilight demarcation line and a real edge line, in which only the real edge line can be used to extract the lunar centre. As a consequence, it is necessary to distinguish the two respective lines.

Normally, such problems tend to be solved by means of robust estimation, for example, the method named Random Sample Consensus (RANSAC), whose basic idea is to search for a best-fit model utilising inliers mixed up with outliers (Fischler and Bolles, Reference Fischler and Bolles1981; Hast et al., Reference Hast, Nysjö and Marchetti2013). However, for the Moon, the outliers (the twilight demarcation line) appear to frequently vote for a regular curve, therefore, this method might occasionally become invalid. Furthermore, for better results, the setting of thresholds might have to be artificially changed with the Moon phase cycle, which is considered to be a serious obstacle in engineering practice. In order to overcome these problems and vote in the real edge line accurately and adequately, a two-step scheme is proposed, which includes a half-intercepted search for preliminary voting and a cyclic search for refinement voting. Making use of the feature of the Moon phase, this method demonstrates good effectiveness and stability.

3.2.1. Half-intercepted Search

Figure 3 is a sketch of the Moon's phases, in which the solid line denotes the real edge line and the dotted line denotes the twilight demarcation line, whose length is always shorter than that of the former. For that reason, if intercepting a half quantity of total Moon edge points as the search step, taking the norm of the algebraic distance vector (ǁDAǁ in Equation (9)) as the ellipse fitting accuracy index, the specific half point set minimising the norm will definitely belong to the lunar real edge line. The implementation steps are carried out as follows:

1. (1) Suppose the total number of edge points detected from the Moon phase is N, a half-rounded continuous quantity ([N/2]) of points are expected to be intercepted for ellipse fitting;

2. (2) N/2 is designated as the length of the search step, during which the search should be kept underway until the step coincides with the first one chosen. In this way, ellipse fitting is repeated N times;

3. (3) Calculate the norms of the algebraic distance vectors;

4. (4) The set of points corresponding to the minimal norm are the ones preliminarily voted in.

Figure 3. Moon phase.

3.2.2. Cyclic Search

The majority of points on the real edge line of the Moon phase can be voted in utilising the method above. However, if we fit the Moon centre merely with them, part of the precision may be lost. Considering this, a further search for the points left should be implemented. To achieve this, a cyclic search method is proposed based on the average of algebraic distances. The detailed implementation steps are carried out as follows:

1. (1) Suppose the total number of Moon phase edge points is N, the real edge points preliminarily voted in are (x ₁, y ₁), (x ₂, y ₂), ……, (x _i, y _i), ……, (x _[N/2], y _[N/2]), and the coefficient vector of the fitting ellipse is A ₁. By introducing vector

   (11)$${\bf P}_{i} =(x_{i}^{2},x_{i} y_{i} ,y_{i}^{2},x_{i} ,y_{i} ,1)$$
   the average of algebraic distances from these points to the fitting ellipse can be expressed as
   (12)$$\overline{d}_{1} ={1 \over [N/2]}\sum_{i=1}^{[N/2]} \vert {\bf P}_{{\bf i}} {\bf A}_{{\bf 1}} \vert $$

2. (2) Suppose one specific point left after the preliminary search process is (x _t, y _t). By introducing vector

   (13)$${\bf P}_{t} =(x_{t}^{2},x_{t} y_{t} ,y_{t}^{2},x_{t} ,y_{t} ,1)$$
   the algebraic distance from this point to the fitting ellipse in step (1) is
   (14)$$d_{t} =\vert F_{{\bf A}_{{\bf 1}}} ({\bf P}_{t} )\vert =\vert {\bf P}_{{\bf t}} {\bf A}_{{\bf 1}}\vert $$

3. (3) If $d_{t} <\overline d_{1} $, point (x _t, y _t) can be treated as one of the Moon's real edge points.

4. (4) Suppose there are W qualified points voted in from step (3). Arriving here, [N/2 +W real edge points can be obtained, including (x ₁, y ₁), (x ₂, y ₂), ……, (x _i, y _i), ……, (x _[N/2]+W, y _[N/2]+W). Using these points for the second elliptical fit, the coefficient vector received is A ₂. Similarly, the average of algebraic distances from these points to the latest fitting ellipse are expressed as

   (15)$$\overline{d}_{2} ={1 \over [N/2]+W}\sum_{i=1}^{[N/2]+W} \vert {\bf P}_{{\bf i}} {\bf A}_{{\bf 2}} \vert $$

5. (5) Calculate the algebraic distances from the points still left to the latest fitting ellipse and determine whether they belong to the Moon's real edge line, just as step (2) and step (3) have conducted. Steps (2)-(4) should be repeated until no point could be voted in from step (3). Arriving here, the real edge line of the Moon could be considered to be completely discovered.

Figure 4 is a flowchart of the process described above.

Figure 4. Flowchart of cyclic search.

4. VESSEL ATTITUDE DETERMINATION BASED ON OBSERVATION OF THE HORIZON

According to the method described above, the coordinates of the Moon centre on the image plane can be extracted and can be further converted to camera coordinates. In order to achieve a Moon-observing heading determination, the coordinates of the Moon centre should be revised to the horizontal plane. The horizon is a natural intersection between the sky and the water at sea, which supplies a natural horizontal datum perpendicular to the local plumb line (Li et al., Reference Li, Luo, Zheng and Zhang2016). Given this beneficial condition, an algorithm to solve the camera's attitude based on horizon observation is proposed.

In Figure 5, O is the camera optical centre, OA and OB are two lines of view for the observation of the horizon, which are at a tangent to the sea at points A and B. All the tangent points constitute a horizontal circle, where O′ is the centre, r is the radius, and E is the geocentre. OS is perpendicular to the sea surface, and point S is located on the sea surface. The length of OS is h.

Figure 5. Principle of horizon observation.

In Δ OAE, ∟ OAE is a right angle, then

(16)$$\vert {OA} \vert =\sqrt{\vert {OE} \vert ^{2}-\vert {AE} \vert ^{2}} =\sqrt{(R+h)^{2}-R^{2}} =\sqrt{h^{2}+2Rh}$$

∟ O O′A is right angle, then

(17)$$r=\vert O'A \vert ={\vert {OA} \vert \cdot \vert {AE} \vert \over \vert {OE} \vert }={R\sqrt{h^{2}+2Rh} \over R+h}$$

Figure 6 is a sketch for the camera shooting the horizon, in which O – X _c Y _c Z _c is the camera coordinate system. O – X _ch Y _ch Z _ch is defined as the camera horizontal coordinate system. In the camera horizontal coordinate system, origin O is the optical centre of the camera, axis Z _ch points to the zenith, axis X _ch is axis X _c ’s projection on the horizontal plane and axis Y _ch is determined by the right-hand rule. W is a point on the horizon and H _W is the elevation angle of W in the camera horizontal coordinate system.

Figure 6. Relationship between coordinate systems.

Two steps are required for the transformation from O – X _c Y _c Z _c to O – X _ch Y _ch Z _ch:

1. (1) Rotate γ radians around axis X_c, and the rotation matrix can be expressed as R_X (γ);

2. (2) Rotate ψ radians around axis Y_c, and the rotation matrix can be expressed as R_Y (ψ);

${\bf R}_{c}^{ch}$ is defined as the rotation matrix from O – X _c Y _c Z _c to O – X _ch Y _ch Z _ch, which can be calculated as

(18)$${\bf R}_{c}^{ch} ={\bf R}_{{\rm Y}} (\psi ){\bf R}_{{\rm X}} (\gamma)$$

Its holonomic form is

(19)$${\bf R}_{c}^{ch} =\left[\matrix{ \cos \psi &\sin \psi \sin \gamma &-\sin \psi \cos \gamma \cr 0 &\cos \gamma &\sin \gamma \cr \sin \psi &-\cos \psi \sin \gamma &\cos \psi \cos \gamma }\right]$$

S_c denotes the 3D vector of point W in the camera coordinate system and it can be transformed into the camera horizontal coordinate system as:

(20)$${\bf S}_{ch} =\left[\matrix{ x_{ch} \cr y_{ch} \cr z_{ch} }\right]={\bf R}_{c}^{ch} \cdot {\bf S}_{c} =\left[\matrix{ x_{c} \cos \psi +y_{c} \sin \psi \sin \gamma -z_{c} \sin \psi \cos \gamma \cr y_{c} \cos \gamma +z_{c} \sin \gamma \cr x_{c} \sin \psi -y_{c} \cos \psi \sin \gamma +z_{c} \cos \psi \cos \gamma }\right]$$

The tangent of point W’s elevation angle is

(21)$$\tan H_{W} ={z_{ch} \over \sqrt{x_{ch}^{2}+y_{ch}^{2}}}$$

In fact, H _W is equivalent to the value of ∟ OAO′. Furthermore, in Figure 5, it is clear that $\angle OAO'=\angle OEA$. So the following equation can be deduced

(22)$$\tan H_{W} ={\vert {OA} \vert \over \vert {AE} \vert }={\sqrt{h^{2}+2Rh} \over R}$$

By combining Equations (21) and (22), the following relationship can be established:

(23)$${z_{ch}^{2} \over x_{ch}^{2}+y_{ch}^{2}}={h^{2}+2Rh \over R^{2}}$$

Then the error equation for least square estimation can be written as

(24)$$v={z_{ch}^{2} \over x_{ch}^{2}+y_{ch}^{2}}-{h^{2}+2Rh \over R^{2}}$$

which includes only two unknown parameters, namely, the pitch angle ψ and the roll angle γ. The holonomic form of Equation (24) is

(25)$$v=f(\psi ,\gamma )={(x_{c} \sin \psi -y_{c} \cos \psi \sin \gamma +z_{c} \cos \psi \cos \gamma )^{2} \over (x_{c} \cos \psi +y_{c} \sin \psi \sin \gamma -z_{c} \sin \psi \cos \gamma )^{2}+(y_{c} \cos \gamma +z_{c} \sin \gamma )^{2}}-{h^{2}+2Rh \over R^{2}}$$

Linearizing Equation (25), and taking the initial value ${\bf X}_{0} =\left[\matrix{ \psi_{0} &\gamma_{0} } \right]^{T}$ for the two unknown parameters, the error equation can be represented as

(26)$${\bf V}={\bf A}\bf{\rdelta}\mathop{{\bf X}}\limits^{\wedge} +{\bf L}$$

In Equation (26), V is the residual vector, A is the coefficient matrix, $\bf{\rdelta}\mathop{{\bf X}}\limits^{\wedge } $ is the correction vector for unknown parameters, and L is the free vector for the error equation. A detailed interpretation of Equation (26) is presented in the Appendix.

Based on the least squares criterion, the estimated vector of the attitude parameters can be calculated as

(27)$$\mathop{{\bf X}}\limits^{\wedge} =-({\bf A}^{T}{\bf A})^{-1}{\bf A}^{T}{\bf L}+{\bf X}_{0}$$

The process above requires iterative calculations. Generally, the two attitude parameters do not fluctuate more than ±10° and they are sensitive to observations. Therefore, in the normal calculation process, the initial values of the parameters can be taken conveniently as 0, and the residual error can be controlled below 0 · 1^ʹʹ by iterating 4~5 times.

5. METHOD TO RESOLVE THE HEADING ANGLE

The Moon centre position on the image plane can be transformed into the camera coordinate system utilising the projection model and distortion model presented in Section 2. Then utilising the camera attitude parameters calculated by the method presented in Section 4, the Moon centre position will be further transformed into the camera horizontal coordinate system. In order to establish contact with the horizon coordinate system (defined as a left-handed system here), the axis Y _ch of the camera horizontal coordinate system should be reversed for conversion to a left-handed system. Figure 7 is a sketch reflecting the key principle of heading determination by lunar observation, where O denotes the camera's optical centre, O – X _h Y _h Z denotes the horizon coordinate system, and O-X_chʹ Y_chʹ Z denotes the converted camera horizontal coordinate system (left-handed system), which share axis Z with the horizon coordinate system. At the moment of camera exposure, L is the horizontal angle of the Moon centre in O-X_chʹ Y_chʹ Z, and A is its azimuth in O – X _h Y _h Z. ξ is the heading angle—the intersection angle between X_chʹ and X _h. Also, ξ is equal to the intersection angle between the projection of axis X _c on the horizontal plane and the south direction.

Figure 7. Principle of heading determination by lunar observation.

The formula to calculate heading angle is

(28)$$\xi =A-L$$

We now briefly introduce the calculation process for horizontal angle L.

Suppose S_chʹ is the vector of the Moon centre in O-X_chʹ Y_chʹ Z, which can be calculated as

(29)$${\bf S}_{c{h}^{\prime}} =\left[\matrix{ x_{c{h}^{\prime}} \cr y_{c{h}^{\prime}} \cr z_{c{h}^{\prime}} }\right]={\bf P}_{Y} {\bf R}_{c}^{ch} {\bf S}_{c}$$

In Equation (29), S_c can be obtained from Equation (4), ${\bf R}_{c}^{ch} $ can be obtained from Equation (19), and P_Y is the reversal matrix of axis Y. Therefore, L can be calculated as

(30)$$L=\left\{ \matrix{ \arctan {y_{c{h}^{\prime}} \over x_{c{h}^{\prime}}} \hfill& x_{c{h}^{\prime}} >0,y_{c{h}^{\prime}} >0 \cr \arctan {y_{c{h}^{\prime}} \over x_{c{h}^{\prime}} }+\pi \hfill& x_{c{h}^{\prime}} < 0 \hfill\cr \arctan {y_{c{h}^{\prime}} \over x_{c{h}^{\prime}} }+2\pi \hfill& x_{c{h}^{\prime}} >0,y_{c{h}^{\prime}} < 0 }\right.$$

As for the azimuth A, it can be calculated taking advantage of the exposure moment, the lunar ephemeris and the vessel's location. Generally, the hour angle method is used to solve this (Adams, Reference Adams1968; Robbins, Reference Robbins2013).

(31)$$\tan A={\cos \delta \sin t \over \sin \varphi \cos \delta \cos t-\sin \delta \cos \varphi}$$

In Equation (31), φ is the astronomical latitude of the vessel, δ is the lunar declination, and t is the lunar hour angle, which can be calculated as

(32)$$t=S_{0} +T_{{\rm ut1}} +\lambda -\alpha$$

In Equation (32), α is the lunar right ascension and λ is the astronomical longitude of the vessel. The AVP written as (ᵩ, λ) can be obtained from sequential shooting of the Moon and information from the inertial navigation system about relative movement during the shooting. T _ut1 is the Greenwich Universal Time (UT1) at the camera exposure moment, which should be calculated from Coordinated Universal Time (UTC) and the Earth orientation parameter —UT1-UTC. The former is recorded by a high-precision timer and the latter is received from the bulletin communique issued regularly by International Earth Rotation Service (IERS). S ₀ is the Greenwich sidereal time when the UT1 becomes zero.

6. EXPERIMENT ANALYSIS

To test the reliability of the heading determination method, two experiments were carried out. The first was conducted in Huanghai coastal waters on the night of 31 October 2012, while the second was in Bohai Bay on the night of 13 July 2017.

The Charge Coupled Devices (CCD) used in the two experiments were different, notably in the array size (3056 × 3056 pixel and 2048 × 2048 pixel respectively) and the quantum efficiency (64% and 96% respectively). The weather conditions during the two experiments met the demand - good visibility of both the Moon and the stars. However, there was a little mist on the sea surface. Fortunately, under moonlight illumination, the observation of the horizon was almost unaffected. Since the ship was always in a state of roll, the heading angle varied from time to time. Thus, it was difficult to obtain the true value of heading angle at the camera exposure moment. To solve this, the average heading result from hundreds of stars was taken as a reference to evaluate the accuracy of Moon-observing heading determination.

During the 31 October 2012 experiment, 15 images were collected, while the number was 50 for the 13 July 2017 experiment. Figure 8 demonstrates the mean square error of heading angle for each image calculated from multiple stars, and the fluctuation ranges are 2·0 × 10⁻³~ 6·6 × 10⁻³ ° and 2·1 × 10⁻³~ 4·9 × 10⁻³ ° respectively. Among the former 15 images, the average mean square error is 4·6 × 10⁻³ °, while the magnitude of error in the later 50 images is 3·4 × 10⁻³ °. Both experiments fully meet the accuracy demand for a vessel's heading and can act as a reference for Moon-observing heading determination. The accuracy of the 13 July 2017 experiment seems to be a little higher, which could be attributed to the CCD's better photosensitive performance, that is, the higher quantum efficiency may help the camera alleviate the influence of fog and clouds, and thus it can capture a clearer horizon and more stars, contributing to more stable and reliable results.

Figure 8. Results of star-observing heading determination.

The data in Figure 9 show the deviation of heading results between the Moon and the stars for each image. As is indicated in the upper figure, the deviations of the 31 October 2012 experiment vary from − 0·045° to 0·026°. Similarly, the situation for the 13 July 2017 experiment in the lower figure is − 0·047° to 0·04°. In addition, the average deviations of the two experiments are 0·021° and 0·019°, respectively. As can be seen from the two graphs, in general, the results of Moon-observing heading determination have stable randomness, mainly fluctuating in the range of − 0·04° ~ 0·04° and − 0·05° ~ 0·05°. Following this, the results of the first experiment tended to behave more steadily, and this is probably because the higher resolution ratio of the CCD (3056 × 3056 pixel) can give a larger Moon image in the image plane, which is helpful for the extraction of the Moon's centre.

Figure 9. Results of Moon-observing heading determination.

In general, the heading errors can be generated by methodological, algorithmic, environmental or hardware causes. Normally, the main possible error causes are:

1. (1) The elevation angle of the Moon changes over time, and the magnitude of influence on the horizontal angle L from the tangential component of the lunar centre extraction error may change with elevation angle.

2. (2) The horizon extraction error may lead to a systematic deviation of camera attitude parameters, further contributing to the L error.

3. (3) The inaccurate camera parameters (especially the principal point's coordinate (x ₀, y ₀)), stemming from the camera calibration error or parameter drift caused by changes in environmental factors (temperature, humidity, air pressure, etc.), may lead to the L error.