Methodology for the Correction of the Spatial Orientation Angles of the Unmanned Aerial Vehicle Using Real Time GNSS, a Shoreline Image and an Electronic Navigational Chart

Abstract

Undoubtedly, Low-Altitude Unmanned Aerial Vehicles (UAVs) are becoming more common in marine applications. Equipped with a Global Navigation Satellite System (GNSS) Real-Time Kinematic (RTK) receiver for highly accurate positioning, they perform camera and Light Detection and Ranging (LiDAR) measurements. Unfortunately, these measurements may still be subject to large errors-mainly due to the inaccuracy of measurement of the optical axis of the camera or LiDAR sensor. Usually, UAVs use a small and light Inertial Navigation System (INS) with an angle measurement error of up to ${0.5}^{\circ }$ (RMSE). The methodology for spatial orientation angle correction presented in the article allows the reduction of this error even to the level of ${0.01}^{\circ }$ (RMSE). It can be successfully used in coastal and port waters. To determine the corrections, only the Electronic Navigational Chart (ENC) and an image of the coastline are needed.

1. Introduction

Unmanned aerial vehicles (UAVs) are increasingly used in numerous marine applications, including intelligence, surveillance and reconnaissance (conducted by the navy or coast guard), observation of the maritime environment (in particular the coastal zone, for example, in the field of environmental protection), supervision and support of maritime traffic (e.g., as part of the Vessel Traffic Service), search and rescue (conducted by SAR services), bathymetric and geodetic measurements (in order to obtain information on the shape of the seabed and for the preparation of cartographic studies, including sea port navigation charts), an inspection of port infrastructure (piers, breakwaters, etc.), and navigation (navigation markings, antenna masts, etc.).

Micro-mini UAVs [1,2], equipped with visible, near-infrared, and multispectral cameras, hyperspectral cameras, thermal cameras and laser scanners, are usually used to perform the indicated tasks [2].

As they rise to a certain height (and most often remain to hover), they enable real-time camera monitoring of large areas of coastal sea waters, including roadways, anchorages, or traffic regulation systems. Navy, coast guard, SAR and VTS can obtain online information about the positions of tracked vessels (also in motion), drifting objects and survivors, supervised navigational markings, or areas of contamination moving on the water surface. It is worth emphasizing that these positions are determined by a camera from considerable distances, reaching even several dozen kilometers. It is also impossible to overestimate the possibility of using USV for bathymetric measurements performed with Light Detection and Ranging (LiDAR) instead of an echosounder or sonar. A measuring system operating in this way allows reducing the workload. Moreover, it is much faster, more effective, and practical to use in increased vessel traffic, for example, ports (where periodic checking of the depth is obligatory). However, it should also be emphasized that the positions of the bottom points are determined using LiDAR from a distance of up to several dozen meters.

An important factor affecting the accuracy of positioning using both the camera and LiDAR is the lack of correlation between the antenna’s position of the positioning system of the mounted UAV and the reflection point of the light wave. Under ideal propagation conditions (ignoring the influence of atmospheric refraction), the light ray can be treated as a straight line with a constant direction to the plumb line and azimuth. Such an assumption makes it possible to link the terrain point (seabed) with the position of the positioning system antenna using constant values of the rotation parameters (most often these are spatial orientation angles: yaw, pitch, roll), which are used to transform the coordinates between the coordinate system of the measuring sensor and the coordinate system connected with the Earth.

However, in real conditions, due to wind or changes in atmospheric pressure and temperature as well as structural imperfections, the values of the UAV’s spatial orientation angles will change during flight. Therefore, it becomes necessary to update them with each new measurement regularly. Only based on their current values can the changed direction of the light beam relative to the plumb line and the azimuth be determined, and then the location of the terrain point (seabed) of the light wave reflection. Currently, most micro-mini UAVs use information about spatial orientation angles obtained from the inertial navigation system (INS) based on data from the inertial measurement unit (IMU) of the type microelectromechanical systems (MEMS) supported by the global navigation satellite system (GNSS) receiver allowing flexible approach for seabed survey through Autonomous Mobile Vehicles [3].

An IMU is typically composed of a 3-axis accelerometer and 3-axis gyroscope. Usually, information from a 3-axis magnetic sensor and a thermal sensor supplements them as well. An IMU combines linear accelerations from an accelerometer and rotations from a gyroscope to deliver navigation parameters and position update information. The accuracy of these parameters is influenced by several errors, which are a function of time. Gyroscope errors are bias stability, angle random walk, and calibration errors. Accelerometer errors are constant bias, velocity random walk, scale factor, and vibration rectification error (VRE). Depending on orientation performance, IMUs belong to different classes; beginning with the most accurate, we have the followings grades: military, marine, navigation, tactical, industrial, automotive, the consumer. The first three are precise, very expensive ($ 100 k–1 M), and not available on the open market. In UAVs, tactical or industrial grade IMUs are used to compromise the price and performance. Typical spatial orientation errors for popular IMUs used in drones are presented in

Table 1. Orientation performance specification of popular IMU used in drones.

Type of IMU *	Price	Roll/Pitch Error (RMS)	Yaw Error (RMS)
Yost Labs 3-Space Sensor	$\$250$	${1.0}^{\circ }$	${2.0}^{\circ }$
Vectornav VN-100	$\$500$	${0.5}^{\circ }$	${2.0}^{\circ }$
XSens MTi-600 series	$\$520$	${0.5}^{\circ }$	${1.0}^{\circ }$
Analog Devices Inc. ADIS16362BMLZ	$\$550$	${0.5}^{\circ }$	${1.0}^{\circ }$
Honeywell Aerospace HG1120CA50	$\$1500$	${0.1}^{\circ }$	${1.0}^{\circ }$
SGB Eclipse 2	$\$4000$	${0.1}^{\circ }$	${0.8}^{\circ }$

* All types use: 3-axis accelerometer, 3-axis gyroscope, 3-axis magnetic sensor.

.

However, the errors presented in Table 1 are only theoretical. In reality, they are higher due to additional errors connected with assembly, such as the IMU coordinate system is not aligned to the axis of the drone, and misalignment between the IMU and camera coordinate systems is also not unknown.

The solution to this problem may be in the methodology of correcting the spatial orientation angles determined by INS, based on the shoreline image generated based on the electronic navigational chart (ENC) as equivalent to the shoreline image (photograph) recorded with a micro-mini UAV camera. ENCs are incredibly widely available on the World-wide ENC Database [4] and cover most of the world’s coastal and port waters [5].

Nevertheless, considering the assumed purpose of applying the methodology, it should be assumed that the calculations should be performed with the smallest possible delay (and preferably in real time) directly on the onboard AUV computer, based on a single photo. Consequently, this excludes the possibility of using the Structure from Motion (SfM) method, commonly used in photogrammetry. On the one hand, it is based on a sequence of overlap photos taken from different positions. On the other hand, it is characterized by data processing in post-processing due to its very high demand for computing power.

In mathematical terms, the proposed methodology uses the existing collinearity of the projection center (whose position is determined by the GNSS RTK receiver), the terrain point (with ENC), and its projection in the photo described by the collinearity equation. By expanding this equation into the Taylor series, limited to the linear part, we obtain correction equations in which the vector of unknowns are the angles: yaw, pitch, roll. Their values are determined by the least squares method due to the iterative process of matching the coastline from the ENC to the shoreline in the photograph.

Considering the proposed concepts of the matching algorithm (based only on changes in orientation angles, taking into account the mean error of the INS angle measurement and the angular resolution of the camera measurement) and the methods of generating the equivalent shoreline image based on ENC (in accordance with the optical parameters of the UAV camera, taking into account the geoid and the added edge points), the methodology can be considered to be unique. The few studies related to the use of ENCs have focused solely on the development of positioning techniques with methods of comparison with respect to the coastline [6,7,8] or methods of geodetic adjustment with respect to navigation marks [9].

On the other hand, those studies concerning optical measurements performed in sea areas focused more on methods of determining spatial orientation angles concerning the line of the horizon [10,11,12,13,14]. In the case of land optical angular measurements based on a single image, most studies focus on the use of specific land signs being encoded [15], in the shape of concentric circles [16], luminous [17], or distinctive terrain points [18].

Considering the proposed methodology from the point of view of conceptual assumptions, the most important of these were determined in the form of a thesis that an equivalent image of the shoreline can be generated only based on exact values of spatial orientation angles corresponding to the photograph taken. The approximate values of the angles of spatial orientation, the exact coordinates of the photograph’s position, the internal parameters of the camera, and the exact terrain coordinates of the points forming the coastline must be known. Due to the fact that the ENC-encoded coordinates of the terrain points forming the coastline are determined by geodetic methods (with millimeter accuracy), the use of the methodology depends primarily on the positioning accuracy of the micro-mini UAV. Until recently, precisely for this reason, the proposed methodology would be useless; however, this is now different, as an increasing number of micro-mini UAVs use GNSS RTK for positioning (with at least centimeter accuracy) [19,20,21,22,23,24,25]. This fact played a crucial role in the decision to conduct research on this methodology and, consequently, present the results in this article.

The remaining part of this article is divided into three sections: The first section presents the methodology of correcting the spatial orientation angles, initially describing the process of preparing a suitable photograph by removing distortions and detecting its edges. The main part of this section is the iterative process of correcting the spatial orientation angles, which consists of generating the shoreline’s image with different spatial orientations and determining the degree of matching of this image to the detected edges in the photograph. Finally, the reader’s attention is focused on the methods of transforming the terrain coordinates to the camera reference system, perspective projection of the edges forming the coastline, the method of fitting the image of the coastline generated based on the ENC to the edges on the photo and the algorithm that automates the correction process.

The second section describes the research on the developed methodology, and is divided into four consecutive stages of research, that is, taking aerial measurements, determining reference angles, determining corrections of the angles, and conducting a precision analysis of the corrected angles in relation to the reference angles in conjunction with the evaluation of the iterative process of determining the corrections of those angles. The third and final section presents generalized conclusions regarding the specifics of the methodology’s operation and application, derived in the light of the results of the conducted research.

2. Method

The proposed methodology for the correction of spatial orientation angles is based on an iterative process of matching an equivalent “computer” generated image of the coastline (based on ENC) to its real image recorded in the form of a photograph taken with a UAV camera. At the outset, to carry out this process, it is necessary to pre-process the photograph. Geometric distortions—caused by the imperfection of the camera’s optical system—have to be removed, and the edges-among which there will be fragments of the coastline, must be detected.

The primary process of correcting the angles of spatial orientation consists of generating (based on the ENC and the known position in which the photograph was taken) an image of the shoreline with different spatial orientations (with added corrections) and then determining the degree of matching this image (according to a given criterion) to the detected edges on the appropriately processed photograph. A single increment or decrement changing the correction value for each of the spatial orientation angles corresponds to the angular value of the measurement resolution, which increases by increments until the nominal resolution of the optical system with which the photograph was taken is reached (this depends primarily on the complementary metal-oxide semiconductor parameters (CMOS) of the matrix and focal length). At the same time, the maximum value of the correction, treated as the absolute value of the sum of all increments, cannot exceed the value of the mean error of the actual measurement of each angle made with the on-board device on the UAV (most often, this is a small INS or IMU).

As a result of the process carried out in this way, the corrected values of spatial orientation angles are finally obtained, corresponding to the best-match coastline image generated based on the ENC for the nominal angular measurement resolution of the optical system, which was used to take a photograph from the UAV.

2.1. Geometric Distortion Removal and Edge Detection

Currently, a number of methods are known to obtain the calibration parameters of the optical system intended to improve the reproducible image. Most often, these parameters are characterized by three main distortions of the optical system being radial, decentering and affine. In situations when distortions caused by radial distortion dominate, the methods described in [26,27,28,29,30,31,32] are proposed. However, in the case of most cameras mounted on UAVs (cheap, with small dimensions, exposed to continuous shocks and vibrations), each of the indicated distortions may occur equally—therefore, the best solution is to use the methods described in publications [33,34,35].

On the other hand, when it comes to detecting the edges in the photograph, the most common methods for this purpose are those that can be divided into two groups, based on either the first-order derivative (gradient) [36,37,38,39] or the second-order derivative (using the Laplace of the Gaussian function) [40,41]. Among them, the method proposed by John Canny [41] was found to be the most suitable for the detection of edges, which are to be fitted in the second place to the coastline generated based on the ENC. This method was chosen mainly due to its three main performance criteria, that is:

• minimizing the number of erroneous detections, where the detection error is both the detection of false edges (false-positive detection) and the omission of the actual edges in the image (false-negative detection),
• ensuring the exact location of the edge—a point classified as an edge point should be as close as possible to the middle point of the real edge,
• generating a single response for each real edge in the image—this is equivalent to one pixel edge detection.

Thus, it is based on the simple assumption that this method should guarantee good detection of edges with a thickness of one pixel, while maintaining their original position in the picture (however, it should be emphasized here that the choice of the edge detection method was not the main subject of research). Figure 1 shows the pier in the original photo:

(a)

taken with a DJI FC6310R camera with the UAV (size: 5472 × 3648 pixels, colorspaces: RGB) [42],

(b)

after removing geometric distortions (camera parameters: $f_x=3670.0$ pixels, $f_y=3663.45$ pixels, $c_x=-2.89$ pixels, $c_y=-0.88$ pixels; distortion coefficients: $k_1=-0.262391$, $k_2=0.111511$, $k_3=-0.0396721$, $p_1=0.000859802$, $p_2=-0.000259255$) [33],

(c)

after edge detection (changing colorspaces: RGB to grayscale; Gaussian blur operation: ${\sigma }_{GBO}=1$; Canny algorithm: hysteresis procedure (150, 225), Sobel operator-two 3 × 3 kernels) [36,41].

2.2. Generating the Coastline Image Based on ENC

Each line or area object is encoded in ENC as a set of nodes using the “chain-node” topology [43]. The location of the node is most often described by ellipsoidal coordinates in the World Geodetic System 1984 (WGS 84) [44], although it should be borne in mind that the “S-57” standard defines 228 of them, each of which can be used [43]. Node pairs form edges that can be utilized, among other things, for the geometric description of the coastline. The “S-57” standard defines, in the “Object Catalog”, several hundred classes of objects, attributes and attribute domains [45]. The object classes and their attributes provided with the standard allow describing of most of the real world objects necessary to be placed on marine navigation charts. For the coding of the land area, the spatial object of the GEO type with the acronym LNDARE (code 71) is used, defined as the solid portion of the Earth’s surface, as opposed to sea or water (IHO Dictionary, S-32, 5th Edition, 2635) [45,46]. Figure 2 shows an image of an electronic navigational chart (display “all other information” [47,48]) with the border of the land area marked.

Based on this, an equivalent image of the coastline can be generated. For this purpose, first of all, the coordinates of the border nodes of the land area should be transformed to the camera reference system, and secondly-a perspective projection on the plane of the edges created between the transformed nodes.

2.2.1. Transforming the Coordinates of Nodes to the Camera Reference System

The transformation of the coordinates of each node from WGS 84 to the camera reference system can be easily described mathematically based on three-dimensional Euclidean space ${\mathbb{E}}^3$ over the field of real numbers $\mathbb{R}$, associated Euclidean space $E^3$,point $\mathbf{O}{\mathbb{E}}^3$, and an orthonormal basis $E^3\left({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3\right)$. The Ortho-Cartesian frame of reference (global benchmark) of space ${\mathbb{E}}^3$ can then be defined as follows:

$$\mathcal{T} = \left\{\mathbf{O},\left({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3\right)\right\},$$

(1)

where the point $\mathbf{O}$ determines its origin (or base point), and versors $\left({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3\right)$ its base. However, the location of each point (node) $\mathbf{P}$ relative to the frame of reference $\mathcal{T}$ can be described as a vector $\overrightarrow{\mathbf{O}\mathbf{P}}$:

$$\mathbf{P}=\overrightarrow{\mathbf{O}\mathbf{P}}={x\mathbf{e}}_1+{y\mathbf{e}}_2+z{\mathbf{e}}_3,$$

(2)

where a set of numbers $\left(x,y,z\right)$ specifies the Ortho-Cartesian coordinates of a point $\mathrm{P}$ relative to the frame of reference $\mathcal{T}$.

Using the above definitions, three Ortho-Cartesian reference frames of space can be established ${\mathbb{E}}^3$: $\mathcal{T} = \left\{\mathbf{O},\left({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3\right)\right\}$, ${\mathcal{T}}^{\mathrm{H}}= \left\{{\mathbf{O}}^{\mathrm{H}},\left({\mathbf{e}}_1^{\mathrm{H}},{\mathbf{e}}_2^{\mathrm{H}},{\mathbf{e}}_3^{\mathrm{H}}\right)\right\}$, and ${\mathcal{T}}^{\mathbf{C}}= \left\{{\mathbf{O}}^{\mathrm{C}},\left({\mathbf{e}}_1^{\mathrm{C}},{\mathbf{e}}_2^{\mathrm{C}},{\mathbf{e}}_3^{\mathrm{C}}\right)\right\}$.

Assuming $\mathcal{T}$ will correspond to the Earth-Centered Earth-Fixed (ECEF) coordinate system [49,50,51], ${\mathcal{T}}^{\mathrm{H}}$ will correspond to the East-North-Up (ENU) coordinate system (Figure 3).

Frame of reference ${\mathcal{T}}^{\mathrm{H}}$ will be obtained as a result of the transformation of $\mathcal{T}$. It will be a point translation $\mathbf{O}$ to the point ${\mathbf{O}}^{\mathrm{H}}$ (which will also be the focal point of the camera ${\mathbf{O}}^{\mathrm{C}}$ located on the UAV) and rotation ${\mathbf{C}}_{\mathrm{ENU}}^{\mathrm{ECEF}}$ base $\left({\mathbf{e}}_1^{\mathrm{H}},{\mathbf{e}}_2^{\mathrm{H}},{\mathbf{e}}_3^{\mathrm{H}}\right)$ relative to the base $\left({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3\right)$ performed according to the dependencies:

$${\mathbf{e}}_1^{\mathrm{H}}=\left[\begin{array}{ccc}{x}_1^{\mathrm{H}}& y_1^{\mathrm{H}}& z_1^{\mathrm{H}}=0\end{array}\right]=-\sin \lambda {\mathbf{e}}_1+\cos \lambda {\mathbf{e}}_2,$$

(3)

$${\mathbf{e}}_2^{\mathrm{H}}=\left[\begin{array}{ccc}{x}_2^{\mathrm{H}}& y_2^{\mathrm{H}}& z_2^{\mathrm{H}}\end{array}\right]=-\sin {\phi }^{{}^{\prime }}·\cos \lambda {\mathbf{e}}_1-\sin {\phi }^{{}^{\prime }}·\sin \lambda {\mathbf{e}}_2+\cos {\phi }^{{}^{\prime }}{\mathbf{e}}_3,$$

(4)

$${\mathbf{e}}_3^{\mathrm{H}}=\left[\begin{array}{ccc}{x}_3^{\mathrm{H}}& y_3^{\mathrm{H}}& z_3^{\mathrm{H}}\end{array}\right]=\cos {\phi }^{{}^{\prime }}·\cos \lambda {\mathbf{e}}_1+\cos {\phi }^{{}^{\prime }}·\sin \lambda {\mathbf{e}}_2+\sin {\phi }^{{}^{\prime }}{\mathbf{e}}_3,$$

(5)

where:

${\phi }^{{}^{\prime }}$—is the geocentric latitude of the point ${\mathbf{O}}^{\mathrm{H}}$,

$\lambda$—means the longitude of the point ${\mathbf{O}}^{\mathrm{H}}$.

Frame of reference ${\mathcal{T}}^{\mathrm{C}}$ will be rigidly related to the projection system (i.e., the projection center and the CMOS projection) of the camera on the UAV in such a way that the direction of the versor ${\mathbf{e}}_3^{\mathrm{C}}$ will be consistent with the direction of the camera’s optical axis. The other two versors, ${\mathbf{e}}_2^{\mathrm{C}},{\mathbf{e}}_3^{\mathrm{C}}$, will lie in a plane parallel to the plane of the camera’s COMS matrix passing through the focal point (projection centre) ${\mathbf{O}}^{\mathrm{H}}={\mathbf{O}}^{\mathrm{C}}$. The versor ${\mathbf{e}}_2^{\mathrm{C}}$ after mapping to the CMOS matrix will lie in the middle of its height, and the ${\mathbf{e}}_1^{\mathrm{C}}$ versor halfway across its width (Figure 4). Rotation ${\mathbf{C}}_{\mathrm{CAM}}^{\mathrm{ENU}}$ of the base $\left({\mathbf{e}}_1^{\mathrm{C}},{\mathbf{e}}_2^{\mathrm{C}},{\mathbf{e}}_3^{\mathrm{C}}\right)$ of the reference system ${\mathcal{T}}^{\mathrm{C}}$ will be performed in relation to the base $\left({\mathbf{e}}_1^{\mathrm{H}},{\mathbf{e}}_2^{\mathrm{H}},{\mathbf{e}}_3^{\mathrm{H}}\right)$ of the reference system ${\mathcal{T}}^{\mathrm{H}}$ according to the dependencies (Figure 4):

$$\begin{array}{c}\hfill {\mathbf{e}}_1^{\mathrm{C}}=\left[\begin{array}{ccc}{x}_1^{\mathrm{C}}& y_1^{\mathrm{C}}& z_1^{\mathrm{C}}\end{array}\right]=\cos \psi ·\cos \varphi +\sin \psi ·\sin \theta ·\sin \varphi {\mathbf{e}}_1^{\mathrm{H}}\\ \hfill +\sin \psi ·\cos \varphi +\cos \psi ·\sin \theta ·\sin \varphi {\mathbf{e}}_2^{\mathrm{H}}-\cos \theta ·\sin \varphi {\mathbf{e}}_3^{\mathrm{H}},\end{array}$$

(6)

$${\mathbf{e}}_2^{\mathrm{C}}=\left[\begin{array}{ccc}{x}_2^{\mathrm{C}}& y_2^{\mathrm{C}}& z_2^{\mathrm{C}}\end{array}\right]=\sin \psi ·\cos \theta {\mathbf{e}}_1^{\mathrm{H}}+\cos \psi ·{\cos \theta \mathbf{e}}_2^{\mathrm{H}}+\sin \theta {\mathbf{e}}_3^{\mathrm{H}},$$

(7)

$$\begin{array}{c}\hfill {\mathbf{e}}_3^{\mathrm{C}}=\left[\begin{array}{ccc}{x}_3^{\mathrm{C}}& y_3^{\mathrm{C}}& z_3^{\mathrm{C}}\end{array}\right]=-\cos \psi ·\sin \varphi +\sin \psi ·\sin \theta ·\cos \varphi {\mathbf{e}}_1^{\mathrm{H}}\\ \hfill +\sin \psi ·\sin \varphi +\cos \psi ·\sin \theta ·\cos \varphi {\mathbf{e}}_2^{\mathrm{H}}-\cos \theta ·\cos \varphi {\mathbf{e}}_3^{\mathrm{H}},\end{array}$$

(8)

where: $\psi$—yaw, $\theta$—pitch, $\varphi$—the roll will be represented by the so-called Tait – Bryan angles [52,53].

Assuming the ellipsoidal coordinates of $\left({\phi }_c,{\lambda }_c,h_c\right)$, point ${\mathbf{O}}^{\mathrm{C}}$ defined in relation to $\mathcal{T}$ (designated GNSS RTK on UAV) are known, and angles ($\psi$, $\theta$, $\varphi$) (corresponding to the yaw, pitch, roll angles obtained from IMU on UAV) bases $\left({\mathbf{e}}_1^{\mathrm{C}},{\mathbf{e}}_2^{\mathrm{C}},{\mathbf{e}}_3^{\mathrm{C}}\right)$ and $\left({\mathbf{e}}_1^{\mathrm{H}},{\mathbf{e}}_2^{\mathrm{H}},{\mathbf{e}}_3^{\mathrm{H}}\right)$ can be determined.

To do this, the ellipsoidal coordinates $\left({\phi }_c,{\lambda }_c,h_c\right)$ to Ortho-Cartesian coordinates $\left(x_c,y_c,z_c\right)$ of the point ${\mathbf{O}}^{\mathrm{C}}$ should be changed according to the dependencies:

$$x_c=\left(\frac{a}{\sqrt{1-e^2·{\sin }^2{\phi }_c}}+h_c\right)·\cos {\phi }_c·\cos {\lambda }_c,$$

(9)

$$y_c=\left(\frac{a}{\sqrt{1-e^2·{\sin }^2{\phi }_c}}+h_c\right)·\cos {\phi }_c·\sin {\lambda }_c,$$

(10)

$$z_c=\left(\frac{a}{\sqrt{1-e^2·{\sin }^2{\phi }_c}}·\left(1+e^2\right)+h_c\right)·\sin {\phi }_c,$$

(11)

where a and b denote the lengths of the major and minor semi-axes of the reference ellipsoid, $e^2=\frac{a^2-b^2}{a^2}$ the square of the first eccentric.

Based on them, the geocentric latitude is determined according to the relationship:

$${\phi }^{{}^{\prime }}=\mathrm{atan}\left(\frac{z_c}{\sqrt{x_c^2+y_c^2}}\right).$$

(12)

Knowing the spherical coordinates $\left({\phi }^{{}^{\prime }},\lambda \right)$ describing the position of the point ${\mathbf{O}}^{\mathrm{C}}={\mathbf{O}}^{\mathrm{H}}$ the base can be set $\mathcal{T}$ in accordance with ${\mathcal{T}}^{\mathrm{H}}$ (3)–(5), which will form the base after rotation by angles ($\psi$, $\theta$, $\varphi$) in accordance with (6)–(8) ${\mathcal{T}}^{\mathrm{C}}$.

Taking the bases so established, $\left({\mathbf{e}}_1^{\mathrm{C}},{\mathbf{e}}_2^{\mathrm{C}},{\mathbf{e}}_3^{\mathrm{C}}\right)$ and $\left({\mathbf{e}}_1^{\mathrm{H}},{\mathbf{e}}_2^{\mathrm{H}},{\mathbf{e}}_3^{\mathrm{H}}\right)$, it is possible to easily transform the Ortho-Cartesian coordinates of each border node ${\mathbf{N}}^{\mathrm{ENC}}=\left(x_{\mathrm{ENC}},y_{\mathrm{ENC}},z_{\mathrm{ENC}}\right)$ (obtained according to (9)–(11) from the source 2D ellipsoidal coordinates $\left({\phi }_{\mathrm{ENC}},{\lambda }_{\mathrm{ENC}}\right)$ encoded in ENC and the ellipsoidal height $h_{\mathrm{ENC}}=\mathrm{geoid}\mathrm{height}$) land area from the frame of reference $\mathcal{T}$ to the frame of reference ${\mathcal{T}}^{\mathrm{C}}$.

This transformation will be a combination of the node coordinate translation ${\mathbf{N}}^{\mathrm{ENC}}$ from $\mathcal{T}$ down ${\mathcal{T}}^{\mathrm{C}}$ and successive coordinate transformations ${\mathbf{C}}_{\mathrm{ENU}}^{\mathrm{ECEF}}$ from the base $\left({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3\right)$ to base $\left({\mathbf{e}}_1^{\mathrm{H}},{\mathbf{e}}_2^{\mathrm{H}},{\mathbf{e}}_3^{\mathrm{H}}\right)$, and then ${\mathbf{C}}_{\mathrm{CAM}}^{\mathrm{ENU}}$ from the base $\left({\mathbf{e}}_1^{\mathrm{H}},{\mathbf{e}}_2^{\mathrm{H}},{\mathbf{e}}_3^{\mathrm{H}}\right)$ to the base $\left({\mathbf{e}}_1^{\mathrm{C}},{\mathbf{e}}_2^{\mathrm{C}},{\mathbf{e}}_3^{\mathrm{C}}\right)$. Relationships for the coordinate transformation carried out in this way, Ortho-Cartesian node ${\mathbf{N}}^{\mathrm{ENC}}=\left(x_{\mathrm{ENC}},y_{\mathrm{ENC}},z_{\mathrm{ENC}}\right)$ in $\mathcal{T}$ to the target Ortho-Cartesian coordinates of the node ${\mathbf{N}}_{\mathrm{ENC}}^{\mathrm{C}}=\left(x_{\mathrm{ENC}}^{\mathrm{C}},{\mathrm{y}}_{\mathrm{ENC}}^{\mathrm{C}},z_{\mathrm{ENC}}^{\mathrm{C}}\right)$ in ${\mathcal{T}}^{\mathrm{C}}$, can be represented in matrix notation as follows:

$${\mathbf{N}}_{\mathrm{ENC}}^{\mathrm{C}}=\left({\mathbf{N}}^{\mathrm{ENC}}-{\mathbf{O}}^{\mathrm{C}}\right)·{\mathbf{C}}_{\mathrm{ENU}}^{\mathrm{ECEF}}·{\mathbf{C}}_{\mathrm{CAM}}^{\mathrm{ENU}},$$

(13)

$$\left[\begin{array}{c}{x}_{\mathrm{ENC}}^{\mathrm{C}}\\ {\mathrm{y}}_{\mathrm{ENC}}^{\mathrm{C}}\\ z_{\mathrm{ENC}}^{\mathrm{C}}\end{array}\right]=\left[\begin{array}{c}{x}_{\mathrm{ENC}}-x_c\\ y_{\mathrm{ENC}}-y_c\\ z_{\mathrm{ENC}}-z_c\end{array}\right]·\left[\begin{array}{ccc}{x}_1^{\mathrm{H}}& y_1^{\mathrm{H}}& z_1^{\mathrm{H}}\\ x_2^{\mathrm{H}}& y_2^{\mathrm{H}}& z_2^{\mathrm{H}}\\ x_3^{\mathrm{H}}& y_3^{\mathrm{H}}& z_3^{\mathrm{H}}\end{array}\right]·\left[\begin{array}{ccc}{x}_1^C& y_1^C& z_1^C\\ x_2^C& y_2^C& z_2^C\\ x_3^C& y_3^C& z_3^C\end{array}\right].$$

(14)

2.2.2. Edge Perspective Projection

Due to the fact that the coastline image should be based on the edges between nodes boundaries of the land area, it becomes necessary to designate points representing them. For this purpose, Bresenham’s algorithm can be successfully applied, assuming a specific distance between consecutive points, for example, equal to the ground sample distance (GSD) in the photograph [54]. Each point of the edge ${\mathbf{N}}_{\mathrm{EDG}}^{\mathrm{C}}=\left(x_{\mathrm{EDG}}^{\mathrm{C}},{\mathrm{y}}_{\mathrm{EDG}}^{\mathrm{C}},z_{\mathrm{EDG}}^{\mathrm{C}}\right)$ so designated can then be subjected to the projection process assuming that the distance of the focal point from the projection, and the projection parameters (i.e., its size and resolution) correspond to the focal length and the CMOS matrix parameters of the UAV camera. The coordinates then obtained as a result of projection ${\mathbf{N}}_{\mathrm{EDG}}^{\mathrm{P}}=\left(x_{\mathrm{P}},y_{\mathrm{P}}\right)$ of each point (forming a single point of the coastline on the equivalent image) can be calculated using the following relationships:

$$x_{\mathrm{P}}=\frac{w_p·f·x_{\mathrm{EDG}}^{\mathrm{C}}}{w_m·z_{\mathrm{EDG}}^{\mathrm{C}}},$$

(15)

$$y_{\mathrm{P}}=\frac{h_p·f·y_{\mathrm{EDG}}^{\mathrm{C}}}{{h_m·z}_{\mathrm{EDG}}^{\mathrm{C}}},$$

(16)

where:

$w_p$—the number of pixels on the horizontal line of the matrix (in the versor axis ${\mathbf{e}}_1^{\mathrm{C}}$),

$h_p$—the number of pixels on the vertical line of the matrix (in the versor axis ${\mathbf{e}}_2^{\mathrm{C}}$),

$w_m$—matrix width,

$h_m$—matrix height,

f—focal length.

Figure 5 shows the coastline’s image obtained as a result of the perspective projection of the edge of the land area from the ENC No. PL5GDYNA.000 corresponding to the DJI FC6310R camera, wherein it was assumed that ${\mathbf{O}}^{\mathrm{C}}=$(54${}^{\circ }$32${}^{\prime }$1.116${}^{″}$ N, 18${}^{\circ }$32${}^{\prime }$48.528${}^{″}$ E, 131.44 m),yaw ($\psi )=$ 128.0${}^{\circ }$, pitch $\left(\theta \right)=-0.1{}^{\circ }$, $\mathrm{roll}\left(\varphi \right)=$ 0${}^{\circ }$, $w_p=$ 5472 pixels, $h_p=$ 3648 pixels, $w_m=$ 12.8333 mm, $h_m=$ 8.55554 mm, and $f=$ 8.5797 mm—these values were read out from the eXtensible Metadata Platform (XMP) tags [55] of the .jpg file of the photograph presented in Figure 1 and the database of camera model parameters in Mapper Pix4D [56].

2.3. Match the Coastline Image to the Edge in the Photo

For determining spatial orientation angle corrections ($\psi$, $\theta$, $\varphi$), an iterative method of matching the boundary of the land area obtained on the basis of the ENC (e.g., Figure 5) to the edges detected in the photo (e.g., Figure 1c) was used, with the evaluation of the sum of squares of the distance “deviations” between the identical points or the border of the area of the close neighborhood (20). At each iteration step i, for each point in the set $\mathbb{N}=\left\{{{\mathbf{N}}_{\mathrm{EDG}}^{\mathrm{P}}}_1,{{\mathbf{N}}_{\mathrm{EDG}}^{\mathrm{P}}}_2,\dots ,{{\mathbf{N}}_{\mathrm{EDG}}^{\mathrm{P}}}_{n_N}\right\}$—representing the boundary of the land area with the ENC, a point “identical” to the set $\mathbb{I}=\left\{{\mathbf{I}}_1,{\mathbf{I}}_2,\dots ,{\mathbf{I}}_{n_I}\right\}$ is assigned, representing the edge points detected in the photograph. The association of the points from set $\mathbb{I}$ with points from set $\mathbb{N}$ is based on the nearest neighbor algorithm (NN) [57]. The area in the vicinity of the point from set $\mathbb{N}$ specifies a circle with radius length $l_i$ (changed in each step of the main iteration i), determined according to the dependencies:

$$l_i=0.75·{\sigma }_i/r_{nom}fori=0,1,\dots ,kand{\sigma }_i⩾r_{nom},$$

(17)

$${\sigma }_{i+1}=\frac{{\sigma }_i}{2},$$

(18)

$$r_{nom}=\mathrm{atan}\frac{\sqrt{{w_m}^2+{h_m}^2}}{f·\sqrt{{w_p}^2+{h_p}^2}},$$

(19)

where:

${{\sigma }_0=\sigma }_{max}$—the highest mean error value of angle measurement: yaw, pitch, roll with navigation devices (e.g., INS) on the UAV,

$r_{nom}$—nominal angular measurement resolution of the optical system (of UAV camera),

k—the number of iterations dependent on the value ${\sigma }_i$ and $r_{nom}$.

The point from set $\mathbb{N}$ has the same point if in its close vicinity there is a point from set $\mathbb{I}$. In the case when the neighborhood affiliation concerns more than one point from set $\mathbb{I}$, the closest of them becomes the counterpart (identical). If, however, identical points are at the same distance, one of them is selected randomly (Figure 6).

Then, all possible combinations of corrections are checked using the brute force method [58]—added to the angles of spatial orientation ($\psi$, $\theta$, $\varphi$), the values of which change step by step $\frac{{\sigma }_i}{2}$ in the range of $\left(-{\sigma }_i,{\sigma }_i\right)$. In order to find the best three among them, the combination for which the sum value is selected S squared distances (“deviations”) $d_1$, $d_2,\cdots$, $d_{n_M}$, calculated according to the formula:

$$S=\sum _{j=1}^{n_M}{d_j}^2,$$

(20)

where:

$d_j=\left\{\begin{array}{c}{l}_i,dmin⩾l_i\\ dmin,dmin<l_i\left(\mathrm{area}\mathrm{of}\mathrm{the}\mathrm{close}\mathrm{neighborhood}\right)\end{array}\right.$,

$$dmin=\underset{\begin{array}{c}m=1,2,\dots ,n_I\hfill \end{array}}{min}∥{{\mathbf{N}}_{EDG}^P}_j-{\mathbf{I}}_m∥,$$

is the smallest.

The best corrections are then added to the angles $\psi$, $\theta$, $\varphi$ which, after changing their values, are corrected again with increasing precision in subsequent iterations. It should be noted here that in each subsequent iterative step, corrections are checked, the values of which are changed with increments and in the interval, reduced by half.

The process of correcting the spatial orientation angles (iteration) is considered complete when the nominal measurement resolution of the optical system $r_{nom}$ is less than the iteratively reduced measurement error ${\sigma }_i$. The algorithm (in the form of pseudocode) of the entire process of correction of spatial orientation angles, implemented in accordance with the proposed methodology, is presented below (see Algorithm 1).

Algorithm 1 Spatial orientation angles correction
Input:				ENC; .jpg photograph; ${\sigma }_{max};$ spatial orientation parameters: ${\phi }_c$, ${\lambda }_c$, $h_c$, $\psi$, $\theta$, $\varphi$; camera parameters: $f_{\mathrm{x}}$, $f_{\mathrm{y}}$, $c_{\mathrm{x}}$, $c_{\mathrm{y}}$, $w_m$, $h_m$, $w_p$, $h_p$, f; distortion coefficients: $k_1$, $k_2$, $k_3$, $p_1$, $p_2$
Output:				Spatial orientation angle corrections: $\Delta {\psi }_{cor}$, $\Delta {\theta }_{cor}$, $\Delta {\varphi }_{cor}$
Initialization:				Removing geometric distortions in the photograph as per [33]. Edge detection: changing color spaces RGB to grayscale, Gaussian blur operation, Canny algorithm.
				Result: a black and white image (one bit) with the resolution of the source image with detected edges and a set of points $\mathbb{I}=\left\{{\mathbf{I}}_1,{\mathbf{I}}_2,\dots ,{\mathbf{I}}_{n_I}\right\}$ Set of ${\mathbf{C}}_{\mathrm{ENU}}^{\mathrm{ECEF}}$ for ${\phi }_c$, ${\lambda }_c$, $h_c$ $r_{nom}=\mathrm{atan}\frac{\sqrt{{w_m}^2+{h_m}^2}}{f·\sqrt{{w_p}^2+{h_p}^2}}$; ${\psi }_{Old}=\psi$; ${\theta }_{Old}=\theta$; ${\varphi }_{Old}=\varphi$
1:	while${\sigma }_{max}⩾r_{nom}$do
2:	$l=0.75·{\sigma }_{max}/r_{nom}$
3:		for${\Delta \psi =-\sigma }_{max}$to${\sigma }_{max}$
4:			for${\Delta \theta =-\sigma }_{max}$to${\sigma }_{max}$
5:				for${\Delta \varphi =-\sigma }_{max}$to${\sigma }_{max}$
6:				Set of ${\mathrm{C}}_{\mathrm{CAM}}^{\mathrm{ENU}}$ for $\psi +\Delta \psi$, $\theta +\Delta \theta$, $\varphi +\Delta \varphi$
7:				Determination of the ellipsoidal height $h_{\mathrm{ENC}}=\mathrm{geoid}\mathrm{height}$ of each land node from the ENC with 2D ellipsoidal coordinates $\left({\phi }_{\mathrm{ENC}},{\lambda }_{\mathrm{ENC}}\right)$
8:				3D ellipsoidal coordinate conversion $\left({\phi }_{\mathrm{ENC}},{\lambda }_{\mathrm{ENC}},h_{\mathrm{ENC}}\right)$ of each knot to Ortho Cartesian ${\mathbf{N}}^{\mathrm{ENC}}=\left(x_{\mathrm{ENC}},y_{\mathrm{ENC}},z_{\mathrm{ENC}}\right)$
9:				Transformation of the coordinates of each node ${\mathbf{N}}^{\mathrm{ENC}}$ on ${\mathbf{N}}_{\mathrm{ENC}}^{\mathrm{C}}=\left({\mathbf{N}}^{\mathrm{ENC}}-{\mathbf{O}}^{\mathrm{C}}\right)·{\mathbf{C}}_{\mathrm{ENU}}^{\mathrm{ECEF}}·{\mathbf{C}}_{\mathrm{CAM}}^{\mathrm{ENU}}$
10:				Determination of the coordinates of each point ${\mathbf{N}}_{\mathrm{EDG}}^{\mathrm{C}}=\left(x_{\mathrm{EDG}}^{\mathrm{C}},{\mathrm{y}}_{\mathrm{EDG}}^{\mathrm{C}},z_{\mathrm{EDG}}^{\mathrm{C}}\right)$ of the edge of Bresenham’s algorithm between successive knots ${\mathbf{N}}_{\mathrm{ENC}}^{\mathrm{C}}$
11:				Dropping of each point ${\mathbf{N}}_{\mathrm{EDG}}^{\mathrm{C}}$ and creating a set of points $\mathbb{N}=\left\{{{\mathbf{N}}_{\mathrm{EDG}}^{\mathrm{P}}}_1,{{\mathbf{N}}_{\mathrm{EDG}}^{\mathrm{P}}}_2,\cdots ,{{\mathrm{N}}_{\mathrm{EDG}}^{\mathrm{P}}}_{n_N}\right\}$
12:				Finding the sum of squared distances $S={\sum }_{j=1}^{n_M}{d_j}^2$ and when it is the smallest, saving the correction values: ${\Delta \psi }_{Best}=\Delta \psi$, ${\Delta \theta }_{Best}=\Delta \theta$, ${\Delta \varphi }_{Best}=\Delta \varphi$
13:				$\Delta \varphi +=\frac{{\sigma }_{max}}{2}$
14:				end for
15:			$\Delta \theta +=\frac{{\sigma }_{max}}{2}$
16:			end for
17:		$\Delta \psi +=\frac{{\sigma }_{max}}{2}$
18:		end for
19:	${\sigma }_{max}$ / = 2
20:	Spatial orientation angles correction: $\psi +={\Delta \psi }_{Best}$; $\Delta {\psi }_{cor}={\psi }_{Old}-\psi$ $\theta +={\Delta \theta }_{Best}$; $\Delta {\theta }_{cor}={\theta }_{Old}-\theta$ $\varphi +={\Delta \varphi }_{Best};$ $\Delta {\varphi }_{cor}={\varphi }_{Old}-\varphi$
21:	end while

3. Research

In the actual research phase, the developed methodology for the correction of spatial orientation angles was planned to be accurately assessed under real conditions. It was decided that it would be based on a comparative analysis carried out against very accurate reference measurements made by the GNSS RTK method and tachometer, recognizing that this type of comparison will fully justify the purposefulness and indicate the right directions for further research phases.

The research was divided into four successive stages:

• Performing aerial measurements-consisting of taking a series of photographs of the coastline from a UAV raid.
• Determining reference angles of spatial orientation—based on the known coordinates of the image position and two characteristic photopoints.
• Determining corrections-in accordance with the methodology presented in point 1, based on the photographs of the coastline and the ENC.
• Conducting an accuracy analysis-based on a comparison of the source and corrected values of the spatial orientation angles in relation to their reference values determined by the GNSS RTK method and the tachometer.

3.1. Taking Measurements

The DJI Phantom 4 RTK UAV was used to carry out the measurements [59]. Before proceeding, the calibration parameters of the optical system of the FC6310R camera mounted on the UAV were established [59]. It was decided to use the DJI Assistant 2 For Phantom software application for this purpose [60,61]. The calibration parameters established in this way could later be encoded as metadata in each photograph using the standard on-board UAV software, enabling the determination of spatial orientation angle corrections according to Algorithm 1 only based on the photo and the ENC.

The port of Gdynia was selected as the test area (i.e., a selected fragment of the coastal zone) for the measurements (Figure 7).

In this area, using the DJI GS RTK App [60], a photogrammetric flight with the DJI Phantom 4 RTK UAV was planned and carried out—taking a series of photographs from a height of 100 m (above mean sea level (MSL)), with fixed values of gimbal rotation angles of the camera: $\theta =0{}^{\circ }$ (in the case of the DJI Phantom 4 RTK UAV, this angle corresponded to the value of $-90{}^{\circ }$), $\varphi =0{}^{\circ }$—so that the optical axis of the camera points vertically downwards. The position coordinates of the DJI Phantom 4 RTK UAV (including those in which the photos were taken) were determined with the on-board GNSS receiver using the RTK method, and using corrections sent from the server of the HxGN SmartNet reference stations network (the mean 3D error of their measurement did not exceed 3.7 cm) [61].

Among the photographs taken showing the coastline, the most useful (representative) for the accuracy assessment were selected. When choosing, firstly, the differentiation of the irregular shape of the shoreline and the quality of its detection depending on the occurring “disturbances” caused by, among others, moored ships or strong reflection of light from the water surface, was considered. Although, at this point it should be mentioned that the final decision on their selection was made only after confirming the possibility of identifying two field points (photopoint—the so-called ground control points, whose positions were determined by the GNSS receiver) in each selected photo: ${\mathrm{Z}}_{{\mathrm{O}}^{\mathrm{P}}}$-corresponding to projection centre on the image (projected point ${\mathbf{O}}^{\mathrm{C}}$) and ${\mathrm{Z}}_{{\mathbf{e}}_2^{\mathrm{P}}}$-the unit vector lying on the axis ${\mathbf{e}}_2^{\mathrm{C}}$ on the image.

3.2. Determination of Reference Angles of Spatial Orientation

3D ellipsoidal coordinates were determined for each of the selected photos $({\phi }_{{\mathbf{O}}^{\mathrm{P}}},{\lambda }_{{\mathbf{O}}^{\mathrm{P}}},$$h_{{\mathbf{O}}^{\mathrm{P}}})$ photopoint position ${\mathbf{Z}}_{{\mathrm{O}}^{\mathrm{P}}}$ with the Leica GS18 T RTK receiver with a CS 20 controller using (the same as the UAV) corrections sent from the HxGN SmartNet reference station network server (their mean error 3D measurement did not exceed 1.7 cm) and the reference angle value ${\psi }_r$ measured with the Leica Nova TS60 total station with mini prism GRS101 [62] from the photopoint ${\mathbf{Z}}_{{\mathrm{O}}^{\mathrm{P}}}$ on photopoint ${\mathbf{Z}}_{{\mathbf{e}}_2^{\mathrm{P}}}$ (the mean error of the angle measurement did not exceed 0.5″).

Based on known coordinates $\left({\phi }_{{\mathbf{e}}_2^{\mathrm{P}}},{\lambda }_{{\mathbf{e}}_2^{\mathrm{P}}},h_{{\mathbf{e}}_2^{\mathrm{P}}}\right)$ photopoint position ${\mathbf{Z}}_{{\mathrm{O}}^{\mathrm{P}}}$ and known coordinates $\left({\phi }_c,{\lambda }_c,h_c\right)$ of position ${\mathbf{O}}^{\mathrm{C}}$(saved as XMP metadata in a .jpg file, tag names: XMP-drone-dji:GPSLatitude, XMP-drone-dji: GPSLongtitude, XMP-drone-dji: AbsoluteAltitude [55] reference values of spatial orientation angles were calculated ${\theta }_r$ and ${\varphi }_r$ (

Table 2. A comparison of five photos with the designated values of three types of spatial orientation angles.

Source	Source Angles			Reference Angles			Correct Angles
Photo	of Spatial			of Spatial			of Spatial
(Photo	Orientation			Orientation			Orientation)
Number)	(deg)			(deg)			(deg)
	$\psi$	$\theta$	$\varphi$	${\psi }_r$	${\theta }_r$	${\varphi }_r$	$\psi +$$\Delta {\psi }_{cor}$	$\theta +$$\Delta {\theta }_{cor}$	$\varphi +$$\Delta {\varphi }_{cor}$
1	126.3	0.1	0	123.72	–0.61	–0.31	123.79	–0.6	–0.3
2	–53.2	0.1	0	–55.76	–0.91	0.92	–55.8	–0.9	0.97
3	85	0.1	0	79.32	–0.31	–0.32	79.35	–0.3	–0.3
4	–52.7	0.1	0	–56.21	–0.74	0.35	–56.2	–0.8	0.34
5	125.9	0	0	122.63	0.53	0.28	121.06	2.15	1.14

). For this purpose, the known WGS 84 ellipsoidal coordinates were first converted into Ortho-Cartesian coordinates according to (9)–(11), describing the position of the points with respect to $\mathcal{T}$ and then transformed to the frame of reference ${\mathcal{T}}^{\mathrm{C}}$ according to (14), assuming that the optical axis of the camera is directed vertically downwards (${\psi =\psi }_r,$$\theta =0{}^{\circ },$$\varphi =0{}^{\circ }$), ultimately thus obtaining the $\left(x_{{\mathrm{O}}^{\mathrm{P}}}^{\mathrm{C}},{\mathrm{y}}_{{\mathrm{O}}^{\mathrm{P}}}^{\mathrm{C}},z_{{\mathrm{O}}^{\mathrm{P}}}^{\mathrm{C}}\right)$ coordinates. Secondly, using simple trigonometric relationships, on their basis, the reference values of the angles of spatial orientation were calculated, according to the dependencies:

$${\theta }_r=\mathrm{atan}\left(\frac{{\mathrm{y}}_{{\mathrm{O}}^{\mathrm{P}}}^{\mathrm{C}}}{\left|z_{{\mathrm{O}}^{\mathrm{P}}}^C\right|}\right),$$

(21)

$${\varphi }_r=-\mathrm{atan}\left(\frac{x_{{\mathrm{O}}^{\mathrm{P}}}^{\mathrm{C}}}{\left|z_{{\mathrm{O}}^{\mathrm{P}}}^{\mathrm{C}}\right|}\right).$$

(22)

The accepted idea of determining the angles of spatial orientation ${\theta }_r$, ${\varphi }_r$ based on the coordinates of the photograph’s position ${\mathbf{O}}^{\mathrm{C}}$ and photopoint ${\mathrm{Z}}_{{\mathrm{O}}^{\mathrm{P}}}$ and the angle ${\psi }_r$ are shown in Figure 8.

The values of the angles of spatial orientation determined in this way, (${\psi }_r$, ${\theta }_r$, ${\varphi }_r$) were used as a reference for the comparative precision analysis presented later in the article.

The main reasons for determining the Leica GS18 T RTK position coordinates with different accuracy in relation to the DJI Phantom 4 RTK UAV can be indicated that these are receivers from different manufacturers, the measurements were carried out in motion, and statically, corrections from the HxGN SmartNet reference station network server were sent via the telephone network T-Mobile cellular priority for transmitting data.

3.3. Determination of Corrections According to Algorithm 1

Determination of corrections was carried out in accordance with the author’s Algorithm 1, prepared in RAD Studio C ++ Builder 10.2 [63]. First, this program, using the Open Source Computer Vision Library (OpenCV-3.4.7) [64], performed the task of image processing with two functions, geometric distortion removal, and edge detection.

The following input (initialization) arguments and data processing options were taken:

• for the geometric distortion removal function:
  • camera parameters: $f_{\mathrm{x}}=3670.0\mathrm{pixels}$, $f_{\mathrm{y}}=3663.45\mathrm{pixels}$, $c_{\mathrm{x}}=-2.89$ pixels, $c_{\mathrm{y}}=-0.88\mathrm{pixels}$, $w_p=5472\mathrm{pixels}$, $h_p=3648\mathrm{pixels}$ (these values were read out from the metadata of XMP photos, tag names: Dewarp Data and Exif, tag names: ExifImageWidth, ExifImageHeight);
  • distortion coefficients: $k_1=-0.262391$, $k_2=0.111511$, $k_3=-0.0396721$,
    $p_1=0.000859802$, $p_2=-0.000259255$ (these values were read out from the metadata of XMP photographs, tag names: Dewarp Data).
• for the edge detection function:
  • changing colorspaces: RGB to grayscale (set arbitrarily);
  • Gaussian blur operation: ${\sigma }_{GBO}=1$ (set arbitrarily);
  • Canny algorithm: hysteresis procedure (150, 225), Sobel operator-two 3 × 3 kernels (set arbitrarily).

Secondly, the program performed the correction process iteratively, calling a function that adjusts the edge line image from the ENC to the edges detected in the image. The following input (initialization) arguments and data processing options were taken:

• to find the rotation matrix ${\mathbf{C}}_{\mathrm{ENU}}^{\mathrm{ECEF}}$ - WGS ellipsoidal coordinates of 84 photo positions $\left({\phi }_c,{\lambda }_c,h_c\right)$ (these values were read out from the XMP metadata of the photos, tag names: XMP-drone-dji: GPSLatitude, XMP-drone-dji: GPSLongtitude, XMP-drone-dji: AbsoluteAltitude);
• to designate $r_{nom}$-camera parameters: $w_p=5472\mathrm{pixels}$, $h_p=3648$ pixels, $w_m=w_p{·p}_s=$ 12.8333 mm, $h_m=h_p{·p}_s=$ 8.55554 mm, $f=p_s\frac{\left(f_x+f_y\right)}{2}=8.5797\mathrm{mm}$ (pixel size $p_s=2.34527\mathrm{m}$ read out from the database of camera model parameters in Mapper Pix4D);
• to find the rotation matrix ${\mathbf{C}}_{\mathrm{CAM}}^{\mathrm{ENU}}$-${\psi }_{Old}=\psi$; ${\theta }_{Old}=90{}^{\circ }+\theta$; ${\varphi }_{Old}=\varphi$, (these values were read out from the XMP metadata of the photographs, tag names: XMP-drone-dji: GimbalYawDegree, XMP-drone-dji: GPSLongtitude, XMP-drone-dji: GimbalPitchDegree, XMP-drone-dji: GimbalRollDegree);
• to determine the ellipsoidal height $h_{\mathrm{ENC}}=\mathrm{geoid}\mathrm{height}\approx 29.27\mathrm{m}$ land nodes with WGS 84 2D ellipsoidal coordinates $\left({\phi }_{\mathrm{ENC}},{\lambda }_{\mathrm{ENC}}\right)$ read out from ENC-EGM2008 geoid model [65,66,67],
• for determining the radius length $l_i$ close neighborhood circles: ${\sigma }_{max}=3{}^{\circ }$ (set arbitrarily).

Figure 9 shows a photo after removing geometric distortions and after edge detection with the boundary of the land area plotted based on the ENC in the 1st, 2nd, and 8th iteration.

3.4. Accuracy Analysis

The matching of the land borderline generated based on the ENC was correct in 97% of the photographs taken (321 in total), showing the edges of the pier. The analysis was carried out based on the results of measurements and calculations made for the five selected images (the most represented, where it was possible to identify photopoints). At the beginning, the following were used: source angles of spatial orientation ($\psi$, $\theta$, $\varphi$), reference spatial orientation angles (${\psi }_r$, ${\theta }_r$, ${\varphi }_r$), and improved spatial orientation angles ($\psi +\Delta {\psi }_{cor}$, $\theta +\Delta {\theta }_{cor}$, $\varphi +\Delta {\varphi }_{cor}$). The obtained values of the above-mentioned types of angles, compiled separately for each photograph in order to facilitate their comparison with each other, are presented in Table 2.

The obtained results allow us to state that the proposed methodology for determining the corrections of spatial orientation angles proved to be very reliable in terms of accuracy in the case of photograph numbers 1–4 (80%). This can be proved by the values of statistical parameters characterizing the accuracy of determining the angles ($\psi +\Delta {\psi }_{cor},\theta +\Delta {\theta }_{cor},\varphi +\Delta {\varphi }_{cor}$) referred to the values of the reference angles corresponding to them (${\psi }_r,{\theta }_r,{\varphi }_r$) presented in

Table 3. Values of statistical parameters characterizing the accuracy of determining angles ($\psi +\Delta {\psi }_{cor},\theta +\Delta {\theta }_{cor},phi+\Delta {\varphi }_{cor}$) with respect to the angles (${\psi }_r$, ${\theta }_r$, ${\varphi }_r$) for photographs 1–4.

Angle	Minimal Absolute Value of the Difference (deg)	Maximal Absolute Value Difference (deg)	Average Value of the Difference (deg)	Standard Deviation of the Difference (deg)
$\psi +\Delta {\psi }_{cor}$	0.03	0.08	0.01	0.08
$\theta +\Delta {\theta }_{cor}$	0.02	0.05	0.01	0.05
$\varphi +\Delta {\varphi }_{cor}$	0.01	0.05	−0.03	0.001

.

The results from Table 3 show that the minimum and maximum values of absolute differences in the measurement results oscillated at the level of hundredths of a degree, which can be considered an excellent result. On the other hand, the mean value of the differences in the measurement results was concentrated around the zero value, only three hundredths of a degree away from it. Nevertheless, the standard deviation of the difference indicates a relatively large scatter of the measurement results around the mean value. Certainly, the cause of this phenomenon may be the achievement of a level of accuracy in determining corrections of spatial orientation angles close to (or even higher than) the accuracy of the reference measurements (despite the fact that the mean error value of the measurements ${\psi }_r$ did not exceed 0.08${}^{\circ }$, however, ${\theta }_r$ and ${\varphi }_r$ 0.03${}^{\circ }$).

Unfortunately, the values of the corrected angles of spatial orientation for photograph 5 were affected by gross errors (1.57${}^{\circ }$ for $\psi +\Delta {\psi }_{cor}$, −1.62${}^{\circ }$ for $\theta +\Delta {\theta }_{cor}$, and −0.68${}^{\circ }$ for $\varphi +\Delta {\varphi }_{cor}$). This was mainly due to the fact that the border of the land area was in the shape of a straight line (it had no distinctive bends). In this case, it is possible to “slide” the straight line of the land area generated from the ENC over the straight edge of the coastline in the image and obtain the optimal level of alignment for many combinations of spatial orientation angles (there is simply more than one solution to this problem). In order to carry out a deeper analysis of extreme cases (for which the best and worst results were obtained), it was decided to follow the course of their iterative process of determining corrected angles of spatial orientation. Figure 10 and Figure 11 show processed photograph No. 1 and No. 5 with the land area boundary generated based on the ENC in the 1st, 2nd, and 8th iterative step.

On the other hand, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16 and Figure 17 present the graph values $dmin$, the number of analyzed points (from sets $\mathbb{N}$ and $\mathbb{I}$), values of the determined corrections of spatial orientation angles ($\Delta {\psi }_{cor},\Delta {\theta }_{cor},\Delta {\varphi }_{cor}$) as a function of successive iterations.

The analysis of the photos presented in Figure 10 and Figure 11 shows that the course of the process of adjusting the boundary of the land area generated based on the ENC in subsequent iterations is similar. A significant match can be seen as early as the third iteration. The visible disturbances in the photographs caused by the reflection of sunlight from the water surface (visible in particular in photograph No. 1) and the hulls and mooring lines of the ships standing at the pier do not affect the selection’s accuracy of the appropriate edges.

In turn, the graphs presented in Figure 10 and Figure 11 show that the mean value as well as the standard deviation of the distance $dmin$ between identical points from sets $\mathbb{N}$ and $\mathbb{I}$ decrease in successive iterations oscillating abruptly within the circle of radius $l_i$ (near neighborhood area). Their stabilization (reducing the amplitude of fluctuations) is much faster in the case of photograph No. 1 than photograph No. 5. This allows us to hypothesize that the achievement of stabilization in later iterations may indicate that the determined corrections of the spatial orientation angles ($\Delta {\psi }_{cor},\Delta {\theta }_{cor},\Delta {\varphi }_{cor}$) will be burdened with a gross error.

The graphs presented in Figure 14 and Figure 15 show that, after the seventh iteration, the number of identical points (found) from set $\mathbb{I}$ clearly decreased. In the case of photograph No. 1, identical points were found for about 20% of all searched points from set $\mathbb{N}$. After the seventh iteration, the spatial orientation angle corrections are also stabilized ($\Delta {\psi }_{cor},\Delta {\theta }_{cor},\Delta {\varphi }_{cor}$). This can be clearly seen in the graphs presented in both Figure 16 and Figure 17. Then, the iteratively reduced value of the measurement error ${\sigma }_i$ is very close to the nominal measurement resolution value of the optical system $r_{nom}$.

The presented analysis of the iterative process of determining the corrections of the spatial orientation angles has confirmed the correctness of the methodology proposed in the article. The only identified main drawback was related to the situation where the land area’s border has a shape similar to a straight line, and its use should be avoided then. Therefore, the decision to use the methodology should be preceded by an analysis of the shape of the coastline generated based on the ENC in terms of the presence of refractions. The decision-making process could then take place, for example, based on the internal energy of the active “snake” contour model [68,69] or the sum of the distance from the line section determined following the Ramer-Douglas-Peucker algorithm [70,71]. Additionally, a requirement for a minimum number of identical points selected from the set may also be considered (e.g., calculated by the method of interval estimation for the assumed confidence level), allowing the determination of the corrections with the error within the specified limits.

4. Conclusions

Taking into account the results of the analysis carried out in the stages of preparation, development, and testing of the methodology for the correction of spatial orientation angles presented in the article, the following generalized conclusions can be drawn:

• The conducted research proves that the presented methodology can guarantee a very high accuracy of measuring the spatial orientation angles of UAVs performed in coastal and port waters. The size of the angle measurement error depends primarily on the measuring accuracy of the camera and the UAV flight altitude. It was confirmed that for a camera with a CMOS matrix resolution of 5472 × 3648 pixels and a UAV flight altitude of about 100 m, the accuracy of measuring spatial orientation angles can be obtained after being corrected at the level of hundredths of a degree.
• The conducted research also confirms the correctness of the methodology. On the one hand, in terms of matching the land boundary line generated based on the ENC, it was 97% accurate (concerning a statistical sample composed of all tested photographs) and on the other hand, in terms of the possibility of determining corrections it was at the level of 80% (concerning a statistical sample composed of the most representative photographs).
• The conducted research has additionally shown that the strategy of searching for the best corrections (exploring the solution space) used within the methodology can be adapted (tuned) to the computing power of the computer or the measuring accuracy of the camera to increase or decrease the calculation time as well as to increase or decrease the accuracy of determining corrections, and this only requires the determination of appropriate values of two parameters: the maximum mean error of the angle measurement: yaw, pitch, roll and the distance between the edge points generated based on the ENC.
• The equivalent image of the land area lines generated from the spatial object of the acronym LNDARE encoded in ENC can be successfully used to match with the image of the shoreline edge from the photo. However, it should be remembered that the use of the methodology should result from the analysis of the shape of the shoreline in terms of the presence of “kinks”. The decision-making process could then take place, for example, based on the internal energy of the active contour model “snake” or the sum of deviations from the line segment determined according to the Ramer-Douglas-Peucker algorithm.
• The proposed methods for both removing geometric distortions and the detection of edges at the water-land interface proved successful in all photos, thus allowing for the efficient implementation of matching the points of the land area line generated based on the ENC. However, it should be borne in mind that it is possible to distort the shape of the land area line developed based on the ENC in accordance with the calibration parameters of the camera and thus resignation from the initial process of image processing, certainly requiring much greater computing power. This problem will be one of the threads of future research carried out by the authors of the publications to optimize the methodology in terms of minimizing the calculation time.
• It should also be emphasized that the purpose of future research related to the implementation of the methodology for use (increasing its level of technological readiness) will also be to adapt the algorithm for the correction of spatial orientation angles to perform parallel calculations with it, for example, using hybrid multi-core systems with Graphics Processing Unit (GPU).

These generalized conclusions give rise to the hypothesis that the proposed methodology has great potential in supporting angular measurements made with UAVs in coastal and port waters. With the ongoing development of vision systems and computing power, it will undoubtedly be able to be used more widely in the future. Perhaps eventually using only a single photo and ENC.