Analysis of Assembly Error Effect on Stability Accuracy of Unmanned Aerial Vehicle Photoelectric Detection System

Featured Application

For photoelectric detection devices based on unmanned aerial vehicles, the analysis results are used as standards of manufacturing process improvements.

Abstract

A photoelectric detection system is a typical type of device widely used for detecting purposes based on unmanned aerial vehicles (UAV). Stability accuracy is the key performance index. Compared to traditional analysis methods aimed at unpredictable error-causing sources, assembly errors can be easily controlled during the manufacturing processes. In this research, an analysis method of assembly error effect on stability accuracy is proposed. First, by using kinematics analysis of homogeneous coordinate transformation, stability accuracy is comprehensively modeled and simulated. Then, by analyzing the manufacturing process, assembly errors of axis perpendicularities, run-outs and gyroscopes are defined and modeled. By simulating different carrier movements, the effects caused by assembly errors under various environments are studied. Finally, error sensitivity is proposed by using standard deviation analysis. Results show that the most sensitive assembly errors are identified, and ranked in order of sensitivity as follows: x-component of pitch axis perpendicularity, y-component of the azimuth gyroscope assembly, and z-component of the pitch gyroscope assembly. In conclusion, the results can be used as standards of manufacturing process improvements, and the proposed methods can be used to provide valuable references for real application scenarios.

1. Introduction

The unmanned aerial vehicle (UAV) photoelectric detection system (hereinafter, the system) is a typical type of device widely used in both military and civilian fields. Specifically, the system can be used to transfer target information and images among the fields, client–server stations (CS), and the satellite situation center (SSC). The typical application scenario is shown in Figure 1. This kind of system is especially used for detecting purposes such as target search, identification, localization, tracking, damage effect feedback, etc. Furthermore, they have been playing the roles of information transfer centers, and they are widely considered to be the most critical device for battlefield scouting based on UAV. The field of sight (FOS) can be considered to be the vision plane chosen by the UAV operator. FOS is chosen to observe the target properly, and it is considered temporally stationary relative to the ground (or to the inertia space). The line of sight (LOS) is defined as the central axis of the system camera.

Ideally, to obtain clear images, LOS should always be kept vertical to FOS. However, there are various disturbances in the working process, and those disturbances will result in angular deviations of LOS in real time. Therefore, the stability of LOS is a fundamental basis for proper working. Over a certain period time of the working process, the intention of stabilizing the LOS or the system is to minimize the angular deviations of LOS relative to FOS. The index to evaluate the stability performance of the system is defined as stability accuracy. It is treated as the focus of this research, and will be defined and analyzed later in detail. On the other hand, those disturbances are considered to be casual elements, and cannot be easily predicted, because they are often caused by movements of UAV and environmental factors such as wind and noises, etc.

From the perspective of product mechanism design, the system is often divided into four major parts: photoelectric and infrared detectors, which are used to capture vision images; multi-axis frame (or gimbal), for supporting physical frames and compensation movements of axial rotations; drive, control panel, and related system interfaces, for controlling the motions of motors on axes and transferring data; and passive vibration damper, for dissolving some of the disturbances. The typical composite of the system is shown in Figure 2.

As mentioned, the stability can be affected by various unpredictable factors. Therefore, the servo-control of the system is the core for realizing stability. The overall structure of the servo-control functions shown in Figure 3 mainly consist of loops of stabilizing, tracking, and positioning. The control loops can correspond to the system’s working modes of stabilizing, tracking, and positioning, respectively. The working modes can be selected and switched by the operator.

Based on the commonly used mechanism structure shown in Figure 2 and the servo-control structure shown in Figure 3, the working modes can be explained as follows: while the system is working on the stabilizing mode, the gyroscopes are used to measure angular velocities coupled to the system in real time. The velocity data will be transferred by the system interfaces. Then, the velocities will be compensated by the motor motions of the stabilizing loop, ensuring the stability of LOS to obtain clear images. In addition, the interface performance of the system is very important, both for data transmission speed and even overall system performance. The reasons for interface failure include damages or wrong assembly of hardware, incompatibilities of data formats, and inferiority of control algorithms. The failures of interfaces will directly lead to losing real-time states, and therefore the overall failures of servo-control functions.

While working on the tracking mode, based on the stabilizing loop, the miss distances of the target are obtained by analyzing the vision images. Then, motor rotations will be controlled to keep the LOS always aligned with the target. Specifically, the target is always kept in the center of FOS. While working on the positioning mode, the real-time positions of the frames are obtained by position sensors, and then the motors will be controlled to rotate to the specified positions.

Obviously, the stabilizing function is the core of the whole system, and the basis of further applications such as tracking and positioning. Therefore, stability accuracy of the stabilizing function is the key performance index, and it is treated as the focus of this research. Compared to unpredictable factors, the assembly errors of system components will also affect that index greatly. Those errors can be controlled relatively easily. The relative analysis will provide valuable guidance for the design, development, and manufacturing process.

In this research, a method to analyze assembly error effect on stability accuracy is proposed. Stability accuracy is comprehensively modeled and simulated. The system’s general kinematics, compensation principle, and major components are studied. Assembly errors of axis perpendicularities, run-outs, and gyroscopes are defined, modeled, and introduced into the overall model. The effects caused by assembly errors under different disturbance frequencies are studied. Error sensitivity is proposed by using standard deviation analysis. The most sensitive assembly errors are identified, and the results can be used to determine standards and priorities of improvements in the future. The modeling and simulations are all carried out in MATLAB R2019b Simulink. The MATLAB codes of the simulation works are available online at Supplementary Materials.

2. Literature Review

Recently, many researchers have highlighted the stabilization accuracy control of photoelectric devices. For instance, Song, et al. [1] studied the stabilization precision control methods of a photoelectric aim-stabilized system. Lin, et al. [2] developed an inertial-stabilized platform for airborne remote sensing using magnetic bearings. Lei, et al. [3] proposed a composite control method for two-axis inertial-stabilized platforms. Liu, et al. [4] compared ANFIS and fuzzy PID control model for performance in a two-axis inertial-stabilized platform. Fei, et al. [5] proposed a dynamic model and control method. Kwak, et al. [6] proposed a dual-stage and digital image-based method for sight stabilization. It can be concluded that the stability control of LOS is the key factor to keep the relevant systems working properly. Stability accuracy is always treated as the key performance index to evaluate the whole system.

On the other hand, the overall stability accuracy of the system can be severely affected by various kinds of disturbances, including carrier movements, winds, noises, etc. Therefore, many researchers have been focusing on the effect caused by different disturbances: Wu, et al. [7] studied the video stabilization with total warping variation model. Spampinato, et al. [8] studied a low-cost rot-translational video stabilization algorithm. Mao, et al. [9] proposed the continuous second-order sliding mode control. Chen, et al. [10] studied the disturbance observer-based control and related methods. Dasgupta, et al. [11] designed the disturbance observer by using a Hirschhorn inverse approach. Those unpredictable disturbances contribute to the stability performance as casual elements, and they can only be compensated for by using real-time mechanism motions.

In contrast, the assembly errors of components can be easily controlled during the manufacturing processes. For instance, Xiang, et al. [12] carried out analysis and evaluation of installation uncertainty in rotary pair measurement. Liu, et al. [13] carried out a machining error analysis of a freeform surface off-axis three-mirror system based on optical performance evaluation. Ren and Wang [14] studied the impacts of installation errors on the calibration accuracy of a gyro accelerometer tested on centrifuge. Wang, et al. [15] studied SINS installation error calibration based on multi-position combinations, including installation errors of gyroscopes and accelerometers. Chen, et al. [16] analyzed a novel method based on a magnetically suspended gyroscope in the presence of installation error. However, it is common that only one kind of assembly error is studied, such as that of axis, gyroscope, or accelerometer, without comprehensive modeling and analysis. Furthermore, these researchers barely considered a two-axis stabilization mechanism, which is commonly adopted by UAV photoelectric detection. In summary, the analysis of assembly error effect on the stability accuracy of a UAV photoelectric detection system deserves further study.

3. Modeling and Simulations of the Assembly Error Effect on the Stability Accuracy

To study the dynamic control performance of the system stabilizing function, each component must be modeled and simulated separately. These will be integrated together in the final step for the simulation experiments.

3.1. Definition and Modeling of Stability Accuracy Based on Pointing Error

Based on the kinematics of space mechanisms, those disturbances can be decomposed into angular motions around three vertical axes of inertial space, and they will all be coupled to LOS. To stabilize LOS, at least two of those motions must be compensated for, and typically they are the rotations around the y- and z-axes of inertial space [17]. The compensations will be realized by the two motor-driven frames, namely the frames of azimuth and pitch, respectively. In this research, they are abbreviated as the azimuth and the pitch hereafter.

UAV photoelectric detection commonly adopts this two-axis stabilization mechanism because of advantages of simplicity and low weight [18]. Based on the above background, an index of stability accuracy $△\theta$ is proposed for evaluation of the angular deviations. $△\theta$ is defined as the root mean square (RMS) of the deviations between ideal angles and real angles. In a statistical sense, the value of $△\theta$ can be estimated by sampling within a certain total sampling time $t_s$ (with sampling period $T_s$ and sampling number $n_s$). In that case, $△\theta$ can be obtained by using Equation (1), where $△\overline{\theta }$ is the mean of the sampling values, $△{\theta }_{si}$ is the sample value of the ith sample, and ${\sigma }_{△\theta }$ is the standard error.

$$\{\begin{array}{c}△\theta \approx \sqrt{△{\overline{\theta }}^2+{\sigma }_{△\theta }^2}\\ △\overline{\theta }=\frac{{{\displaystyle \sum }}_{i=1}^{n_s}△{\theta }_{si}}{n_s}\\ {\sigma }_{△\theta }^2=\frac{{{\displaystyle \sum }}_{i=1}^{n_s}{\left(△{\theta }_{si}-△\overline{\theta }\right)}^2}{n_s-1}\end{array},$$

(1)

As shown in Figure 4, in practice, the value of Δθ_si is calculated by using sampling values, which include the angular deviations along the y-axis Δθ_syi and that of z-axis Δθ_szi. The calculation method is expressed as Equation (2). Δθ_syi and Δθ_szi will be obtained by using the signals of the two gyroscopes assembled on the frames.

$$△{\theta }_{si}=\sqrt{△{\theta }_{syi}^2+△{\theta }_{szi}^2},$$

(2)

3.2. Definitions of Coordinates and Kinematics Analysis of the System

To stabilize LOS, the angular deviations must be compensated for by the two motor-driven frames. To realize that, establishment of coordinates and relevant kinematics analysis is essential. To describe the motion relationships among the components in the system, their coordinates and motions should be defined separately. Those definitions of the commonly used azimuth-pitch two-axis mechanism are as follows, and shown in Figure 5.

• The inertial coordinate

This is denoted as $\left\{i\right\}$ with axes of $\left(o_i{x}_i{y}_i{z}_i\right)$. The center point of the Earth is taken as the coordinate origin $o_i$; $o_i{x}_i$ is along with the Earth’s rotation axis; $o_i{y}_i$ and $o_i{z}_i$ are orthogonal to each other in the equatorial plane of the Earth, pointing to two fixed stars. $o_i{x}_i$ and $o_i{y}_i$ form a right-hand coordinate with $o_i{z}_i$ perpendicularly.

• The coordinate of the system base

This is denoted as $\left\{b\right\}$ with $\left(o_b{x}_b{y}_b{z}_b\right)$. To simplify the modeling, the system is treated as fixed to the UAV, neglecting effects of the vibration damper, which is not the focus of this research. Hence, the coordinates of the UAV and the system base are unified as $\left\{b\right\}$. $o_b$ is the geometric center point of the location where the azimuth frame is assembled on the base; $o_b{x}_b$, $o_b{y}_b$ and $o_b{z}_b$ are along the heading, left-moving, and rising directions of the UAV. Then, the linear and angular velocity vectors of the system base in the inertial coordinate $\left\{i\right\}$ are denoted as ${\overrightarrow{v}}_{ib}={\left[v_{ibx},v_{iby},v_{ibz}\right]}^T$ and ${\overrightarrow{\omega }}_{ib}={\left[{\omega }_{ibx},{\omega }_{iby},{\omega }_{ibz}\right]}^T$.

• The coordinate of the azimuth frame and axis

This is denoted as $\left\{a\right\}$ with $\left(o_a{x}_a{y}_a{z}_a\right)$. The coordinates of the azimuth frame and axis are unified as $\left\{a\right\}$ for simplification. $o_a$ is the geometric center where the pitch frame is connected to the azimuth frame. $o_a{x}_a$, $o_a{y}_a$ and $o_a{z}_a$ form a coordinate, parallel to $\left\{b\right\}$ when the system has no azimuth rotation; the compensation azimuth rotation is around $o_a{z}_a$. The angular velocity vector of the azimuth in $\left\{i\right\}$ is denoted as ${\overrightarrow{\omega }}_{ia}={\left[{\omega }_{iax},{\omega }_{iay},{\omega }_{iaz}\right]}^T$; the rotation angle and angular velocity of azimuth compensation driven by the azimuth motor are denoted as ${\theta }_a$ and ${\dot{\theta }}_a$, which are around $o_a{z}_a$ in $\left\{a\right\}$.

• The coordinate of the pitch frame and axis

This is denoted as $\left\{p\right\}$ with $\left(o_p{x}_p{y}_p{z}_p\right)$. $o_p$ is treated as the start point of the LOS, and $o_p{x}_p$ is treated as the LOS for simplification; $\left\{p\right\}$ is parallel to $\left\{a\right\}$ when the system has no pitch rotation; the compensation pitch rotation is around $o_p{y}_p$; the angular velocity vector of the pitch in $\left\{i\right\}$ is denoted as ${\overrightarrow{\omega }}_{ip}={\left[{\omega }_{ipx},{\omega }_{ipy},{\omega }_{ipz}\right]}^T$; the rotation angle and angular velocity of pitch compensation around $o_p{y}_p$ in $\left\{p\right\}$ are denoted as ${\theta }_p$ and ${\dot{\theta }}_p$.

• The coordinates of the gyroscopes

These include $\left\{ga\right\}$ with $\left(o_{ga}{x}_{ga}{y}_{ga}{z}_{ga}\right)$ for the azimuth gyroscope, and $\left\{gp\right\}$ with $\left(o_{gp}{x}_{gp}{y}_{gp}{z}_{gp}\right)$ for the pitch one. $\left\{ga\right\}$ and $\left\{gp\right\}$ are ideally parallel to $\left\{p\right\}$. $o_{ga}$ and $o_{gp}$ are the original points where the gyroscopes are assembled. The sensitive axis for LOS angular velocity of ${\omega }_{ipz}$ is $o_{ga}{z}_{ga}$ of $\left\{ga\right\}$, and that for ${\omega }_{ipy}$ is $o_{gp}{y}_{gp}$ of $\left\{gp\right\}$.

3.3. Modeling and Simulations of the General Kinematics and Compensation Principle

According to kinematics analysis of homogeneous coordinate transformation [19], the rotation transformation matrixes from $\left\{b\right\}$ to $\left\{a\right\}$ and from $\left\{a\right\}$ to $\left\{p\right\}$ can be expressed as Equation (3), where $c{\theta }_a$ means $cos{\theta }_a$, $s{\theta }_a$ means $sin{\theta }_a$, and likewise hereinafter for simplicity.

$$R_b^a=\left[\begin{array}{ccc}c{\theta }_a& s{\theta }_a& 0\\ -s{\theta }_a& c{\theta }_a& 0\\ 0& 0& 1\end{array}\right],R_a^p=\left[\begin{array}{ccc}c{\theta }_p& 0& -s{\theta }_p\\ 0& 1& 0\\ s{\theta }_p& 0& c{\theta }_p\end{array}\right],$$

(3)

The rotation of the azimuth in $\left\{i\right\}$ can be decomposed into that of the azimuth in $\left\{b\right\}$ and that of the system base in $\left\{i\right\}$, which is expressed as Equation (4), where ${\overrightarrow{\omega }}_{ba}=\left[0,0,\theta a\right]$, representing the angular velocity vector of the azimuth in $\left\{b\right\}$.

$${\overrightarrow{\omega }}_{ia}=R_b^a{\overrightarrow{\omega }}_{ib}+{\overrightarrow{\omega }}_{ba},$$

(4)

Similarly, the rotation of the pitch in $\left\{i\right\}$ can be decomposed into the pitch rotation in $\left\{a\right\}$ and the azimuth rotation in $\left\{i\right\}$, which is expressed as Equation (5), where ${\overrightarrow{\omega }}_{ap}=\left[0,{\dot{\theta }}_p,0\right]$, representing the angular velocity vector of the pitch in $\left\{a\right\}$.

$${\overrightarrow{\omega }}_{ip}=R_a^p{\overrightarrow{\omega }}_{ia}+{\overrightarrow{\omega }}_{ap},$$

(5)

By introducing Equation (3) and Equation (4) into Equation (5) and decomposing as needed, Equation (6) is obtained, where ${\overrightarrow{\omega }}_{ib}^p={\left[{\omega }_{ibx}^p,{\omega }_{iby}^p,{\omega }_{ibz}^p\right]}^T$, representing the projection of the system base angular velocity relative to the inertia coordinate in $\left\{p\right\}$, and ${\overrightarrow{\omega }}_{bp}^p={\left[{\omega }_{bpx}^p,{\omega }_{bpy}^p,{\omega }_{bpz}^p\right]}^T$, representing the projection of the pitch angular velocity relative to the system base in $\left\{p\right\}$.

$${\overrightarrow{\omega }}_{ip}=\left[\begin{array}{c}{\omega }_{ipx}\\ {\omega }_{ipy}\\ {\omega }_{ipz}\end{array}\right]=\left[\begin{array}{c}{\omega }_{ibx}^p\\ {\omega }_{iby}^p\\ {\omega }_{ibz}^p\end{array}\right]+\left[\begin{array}{c}{\omega }_{bpx}^p\\ {\omega }_{bpy}^p\\ {\omega }_{bpz}^p\end{array}\right]=\left[\begin{array}{c}{\omega }_{ibx}c{\theta }_{a}c{\theta }_p+{\omega }_{iby}s{\theta }_{a}c{\theta }_p-{\omega }_{ibz}s{\theta }_p\\ -{\omega }_{ibx}s{\theta }_a+{\omega }_{iby}c{\theta }_a\\ {\omega }_{ibx}c{\theta }_{a}s{\theta }_p+{\omega }_{iby}s{\theta }_{a}s{\theta }_p+{\omega }_{ibz}c{\theta }_p\end{array}\right]+\left[\begin{array}{c}{\dot{\theta }}_{a}s{\theta }_p\\ {\dot{\theta }}_p\\ {\dot{\theta }}_{a}c{\theta }_p\end{array}\right],$$

(6)

As mentioned, axial rotation movements of y- and z-axes of the inertial coordinate coupled to the LOS must be compensated. Those are realized by using the motor-driven azimuth and pitch frames to make ${\omega }_{ipy}=0$ and ${\omega }_{ipz}=0$. In that case, the compensation angular velocities ${\dot{\theta }}_a$ and ${\dot{\theta }}_p$ can be inferred from Equation (6), which is expressed as Equation (7), where $tg{\theta }_p$ represents $tan{\theta }_p$.

$$\{\begin{array}{c}{\dot{\theta }}_a=-\frac{{\omega }_{ibz}^p}{c{\theta }_p}=-\left({\omega }_{ibx}c{\theta }_{a}tg{\theta }_p+{\omega }_{iby}s{\theta }_{a}tg{\theta }_p+{\omega }_{ibz}\right)\\ {\dot{\theta }}_p=-{\omega }_{iby}^p=-{\omega }_{iby}c{\theta }_a+{\omega }_{ibx}s{\theta }_a\end{array},$$

(7)

Practically, the azimuth and pitch gyroscopes are assembled on the pitch frame to detect the real-time value of ${\omega }_{ibz}^p$ and ${\omega }_{iby}^p$. The sensitive velocities of the two gyroscopes are denoted as ${\omega }_{gyro\_z}={\omega }_{ibz}^p$ and ${\omega }_{gyro\_y}={\omega }_{iby}^p$. Therefore, the relationship between the gyroscope sensitive velocities and the compensation velocities are expressed as Equation (8):

$$\{\begin{array}{c}{\dot{\theta }}_a=-{\omega }_{ibz}^p/c{\theta }_p=-{\omega }_{gyro\_z}sec{\theta }_p\\ {\dot{\theta }}_p=-{\omega }_{iby}^p=-{\omega }_{gyro\_y}\end{array},$$

(8)

Based on the analysis above, the general kinematics and compensation principle of the system are simulated in Simulink as shown in Figure 6. In the simulation module of the azimuth kinematics from $\left\{b\right\}$ to $\left\{a\right\}$, e.g., the angular velocity vector of the system base in $\left\{i\right\}$, namely ${\overrightarrow{\omega }}_{ib}$, is set as one of the inputs. The other input is the z-axial component of the azimuth velocity vector ${\overrightarrow{\omega }}_{ia}$, namely ${\omega }_{iaz}$. ${\omega }_{iaz}$ can be obtained by using signals of the motor modules, which will be discussed later. The velocity vector of the azimuth in $\left\{i\right\}$, namely ${\overrightarrow{\omega }}_{ia}$, is set as the output. Based on Equation (4), the transformation matrixes of $R_b^a$ are imbedded in the module. The azimuth motor velocity and the rotation angle, namely ${\dot{\theta }}_a$ and ${\theta }_a$, can also be calculated by using inverse operation of Equation (4). Similarly, the simulation module for the pitch kinematics from $\left\{a\right\}$ to $\left\{b\right\}$ can also be interpreted by using Equation (5). Those velocity vectors and parameters will be used for later analysis and modeling.

3.4. Modeling and Simulation of the Motor Kinematics

All the compensation movements are driven by the motors including the azimuth motor and the pitch motor. The equivalent circuit of the motor is shown in Figure 7, where $U_m$ is the motor armature voltage; $I_m$ is the armature current; $R_m$ is the total resistor of the armature circuit; $L_m$ is the total inductance; $M_m$ is the motor output torque; $J_m$ and $J_L$ are the rotary inertia of the motor and the load; and $D_L$ and $K_L$ are the damping and elastic coefficients between the motor and the load. ${\theta }_m$, $d{\theta }_m/dt$, and $d^2{\theta }_m/d{t}^2$ are the outputs of angle, angular velocity, and angular acceleration of the motor side. ${\theta }_L$, $d{\theta }_L/dt$, and $d^2{\theta }_L/d{t}^2$ are those of the load side, respectively.

The equations of voltage balance, back electro-dynamic force (EMF), and moment for DC torque motor are shown in Equation (9), where $E_m$ is the EMF, $C_e$ is the coefficient of EMF, $C_m$ is the coefficient of the motor torque.

$$\{\begin{array}{c}{U}_m=E_m+I_m{R}_m+L_m\left(d{I}_m/dt\right)\\ E_m=C_e\left(d{\theta }_m/dt\right)\\ M_m=C_m{I}_m\end{array},$$

(9)

By combining the three equations of Equation (9), Equation (10) is obtained:

$$U_m=C_e\frac{d{\theta }_m}{dt}+\frac{M_m}{C_m}{R}_m+\frac{L_m}{C_m}\frac{d{M}_m}{dt},$$

(10)

The kinematics of the motor is analyzed and expressed as Equation (11):

$$\{\begin{array}{c}{J}_m\frac{d^2{\theta }_m}{d{t}^2}+D_L\left(\frac{d{\theta }_m}{dt}-\frac{d{\theta }_L}{dt}\right)+K_L\left({\theta }_m-{\theta }_L\right)=M_m\\ J_L\frac{d^2{\theta }_L}{d{t}^2}+D_L\left(\frac{d{\theta }_L}{dt}-\frac{d{\theta }_m}{dt}\right)+K_L\left({\theta }_L-{\theta }_m\right)=0\end{array},$$

(11)

The servo-control of the motor generally adheres to principles of pulse width modulation (PWM) or automatic drive (AD). The input is the compensation velocity as shown in Equation (7), and denoted as ${\dot{\theta }}_L$ here. The output is the control voltage $U_c$. In this research, the driving circuit of the motor is taken as a lagging amplification model. Therefore, the Laplace transform function can be expressed as Equation (12), where $K_p$ is the amplification coefficient; $T_p$ is the lagging time. Because the working frequency of the circuit is much larger than the cut-off frequency of the motor, the effect of lagging time can be ignored.

$$U_c\left(s\right)=K_p{e}^{-T_{p}s}{\dot{\theta }}_L\left(s\right)\approx K_p{\dot{\theta }}_L\left(s\right),$$

(12)

Using Laplace inverse transformation, Equation (12) can be transformed as Equation (13):

$$U_c\left(t\right)=K_p{\dot{\theta }}_L\left(t\right),$$

(13)

The voltage balance of the servo-control circuit is expressed as (14), where $E_s$ is the EMF of the voltage source, and $R_s$ is the resistor.

$$\{\begin{array}{c}{U}_c=E_s+U_m\\ E_s=R_s{I}_m\end{array},$$

(14)

Based on the analysis above, the kinematics and servo-control of the motors are simulated in Simulink. For example, the simulation of the pitch motor is shown in Figure 8. The real-time motor velocity ${\dot{\theta }}_p$ is set as one of the inputs. ${\dot{\theta }}_p$ can be obtained by using output of the pitch kinematics module in Section 3.4. Gyroscope signal of ${\omega }_{gyro\_y}$ is the other input, and it can be obtained by using signals of the gyroscope module, which will be discussed later.

Based on the above analyses of Equation (10), Equation (11), Equation (13), and Equation (14), the y-axial component of the pitch velocity vector ${\overrightarrow{\omega }}_{ip}$, namely ${\omega }_{ipy}$, is equal to ${\dot{\theta }}_L$ by definition. So ${\omega }_{ipy}$ is set as the final output of this module and it will be used in other modules. The simulation module for the azimuth motor is similar, and the output will be ${\omega }_{iaz}$, which is discussed and used in Section 3.3.

3.5. Modeling and Simulation of the Assembly Error Effects on Stability Accuracy

From Equation (6) and Equation (7), it can be seen how the UAV movements and the system attitude generally affect the stability accuracy. However, considerations of other error-causing sources, including environment disturbances, component noises, or assembly errors, are neglected. Those will deteriorate the overall stability accuracy index once they are coupled to the system movements, which cannot be ignored in practice. In this research, major error sources including assembly errors of axis perpendicularities, run-outs and gyroscope locations are modeled. Those analysis results will provide guidance for future possible improvements.

3.5.1. Modeling and Simulation of Assembly Errors of Axis Perpendicularities

Ideally, the coordinate of the azimuth $\left\{a\right\}$ is parallel to $\left\{b\right\}$ when the system has no azimuth rotation. However, due to the effect of assembly error, the real axis location is not perfectly perpendicular to the base, and the perpendicularity deviation can be denoted as a vector of ${\overrightarrow{\delta }}_{ba}={\left[{\delta }_{bax},{\delta }_{bay},0\right]}^T$.

As a result, the coordinate $\left\{a\right\}$ and especially for the rotation axis of the azimuth compensation $o_a{z}_a$ will also deviate from the ideal direction by ${\overrightarrow{\delta }}_{ba}$. (For simplification, $o_a{z}_a$ is taken as the instantaneous rotation axis for error analysis of perpendicularity, rather than the average one.) For the pitch compensation, the perpendicularity deviation vector is denoted as ${\overrightarrow{\delta }}_{ap}={\left[{\delta }_{apx},0,{\delta }_{apz}\right]}^T$ respectively. The meaning of axis perpendicularity deviations is shown in Figure 9. In this case, the rotation transformation matrix affected by errors of axis perpendicularity can be expressed as Equation (15):

$$\begin{gathered} R_a^{\prime }=\left[\begin{array}{ccc}c{\delta }_{bay}& 0& -s{\delta }_{bay}\\ 0& 1& 0\\ s{\delta }_{bay}& 0& c{\delta }_{bay}\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& c{\delta }_{bax}& s{\delta }_{bax}\\ 0& -s{\delta }_{bax}& c{\delta }_{bax}\end{array}\right], \\ {R}_p^{\prime }=\left[\begin{array}{ccc}c{\delta }_{apz}& s{\delta }_{apz}& 0\\ -s{\delta }_{apz}& c{\delta }_{apz}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& c{\delta }_{apx}& s{\delta }_{apx}\\ 0& -s{\delta }_{apx}& c{\delta }_{apx}\end{array}\right], \end{gathered}$$

(15)

The angular velocity vectors of the azimuth and the pitch in $\left\{i\right\}$ affected by errors of axis perpendicularity are denoted as ${\overrightarrow{\omega }}_{ia}^{\prime }$ and ${\overrightarrow{\omega }}_{ip}^{\prime }$, and they can be calculated by using Equation (16):

$$\begin{gathered} {\overrightarrow{\omega }}_{ia}^{\prime }=R_b^a{R}_a^{\prime }{\overrightarrow{\omega }}_{ib}+{\overrightarrow{\omega }}_{ba}, \\ {\overrightarrow{\omega }}_{ip}^{\prime }=R_a^p{R}_p^{\prime }{\overrightarrow{\omega }}_{ia}^{\prime }+{\overrightarrow{\omega }}_{ap}=R_a^p{R}_p^{\prime }{R}_b^a{R}_a^{\prime }{\overrightarrow{\omega }}_{ib}+R_a^p{R}_p^{\prime }{\overrightarrow{\omega }}_{ba}+{\overrightarrow{\omega }}_{ap}, \end{gathered}$$

(16)

Based on Equation (16) and analysis above, the signals of the gyroscopes affected by ${\overrightarrow{\delta }}_{ba}$ and ${\overrightarrow{\delta }}_{ap}$ are simulated in Simulink as shown in Figure 10. Taking the azimuth, for example, based on Equation (16), ${\overrightarrow{\omega }}_{ib}$ is set as the input; the perpendicularity deviation vector ${\overrightarrow{\delta }}_{ba}$ and the transformation matrix $R_a^{\prime }$ are imbedded in the middle; the output is $R_a^{\prime }{\overrightarrow{\omega }}_{ib}$, which will be used in other modules. The value of ${\overrightarrow{\delta }}_{ba}$ can be set as needed in the future simulation experiments. The simulation module for the pitch frame is similar. Additionally, the value of ${\overrightarrow{\delta }}_{ba}$ and ${\overrightarrow{\delta }}_{ap}$ are determined by the assembly accuracies of the axes and frames. They can be measured by using high-precision instruments [20].

3.5.2. Modeling and Simulation of Assembly Errors of Run-Outs

As mentioned, the instantaneous rotation axes are taken as the rotation center lines for error analysis of perpendicularity, rather than the average one. However, in practice, due to bad assembly of the bearings, the rotation center lines can be dynamically unsteady under working state. The maximum deviations between the instantaneous rotation axes and the average ones are defined as the assembly errors of run-outs.

Specifically, run-outs include axial shifts $△s$, radial run-outs $△{\overrightarrow{R}}_o$, and tilting oscillations $△\overrightarrow{\gamma }$ ($△{\overrightarrow{R}}_o$ and $△\overrightarrow{\gamma }$ can be decomposed along two axes of the coordinate [21]). Taking the azimuth for example, the assembly errors of run-outs are shown in Figure 11. Apparently, displacement errors such as $△s$ and $△{\overrightarrow{R}}_o$ had no effect on the stability accuracy, and angular tilting oscillations $△\overrightarrow{\gamma }$ are what this research focuses on.

The oscillation vector of the azimuth is expressed as $△{\overrightarrow{\gamma }}_a={\left[△{\alpha }_{am},△{\beta }_{am},0\right]}^T$, and that of the pitch is expressed as $△{\overrightarrow{\gamma }}_p={\left[△{\alpha }_{pm},0,△{\beta }_{pm}\right]}^T$. $△{\alpha }_{am},△{\beta }_{am},△{\alpha }_{pm},△{\beta }_{pm}$ are the decomposed components of the maximum tilting oscillations along the coordinates $\left\{b\right\}$ and $\left\{a\right\}$, and they can also be measured [22].

From the perspective of dynamic signal processing, the tilting oscillations are periodic signals, with the angular frequencies of ${\dot{\theta }}_a$ and ${\dot{\theta }}_p$. Taking $△{\alpha }_{am}$ for example, its real-time value expression $△{\alpha }_a\left(t\right)$ is a periodic signal with angular frequency of ${\dot{\theta }}_a$ and period of $T_a=2\pi /{\dot{\theta }}_a$, then $△{\alpha }_{am}$ can be expressed as Equation (17) by using Fourier series:

$$\{\begin{array}{cc}△{\alpha }_a\left(t\right)\hfill & =\frac{△{\alpha }_{a0}}{2}+{{\displaystyle \sum }}_{n=1}^{\infty }\left(a_{n}c\left(n{\dot{\theta }}_{a}t\right)+a_{n}s\left(n{\dot{\theta }}_{a}t\right)\right)\hfill \\ & a_n=\frac{2}{T_a}{{\displaystyle \int }}_{-\frac{T_a}{2}}^{\frac{T_a}{2}}△{\alpha }_a\left(t\right)c\left(n{\dot{\theta }}_{a}t\right)dt\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em},\left(n=0,1,2,\dots \right),\hfill \\ & b_n=\frac{2}{T_a}{{\displaystyle \int }}_{-\frac{T_a}{2}}^{\frac{T_a}{2}}△{\alpha }_a\left(t\right)s\left(n{\dot{\theta }}_{a}t\right)dt\hfill \end{array}$$

(17)

For simplicity, $△{\alpha }_a\left(t\right)$ is taken as an odd function, and high-frequency signals with $n>1$ are ignored. Then $△{\alpha }_a\left(t\right)$ can be expressed as Equation (12), and $△{\alpha }_{am}$ is considered to be the amplitude of this signal. Similarly, other tilting oscillations are also shown in Equation (18).

$$\{\begin{array}{c}△{\alpha }_a\left(t\right)\approx \frac{2}{T_a}{{\displaystyle \int }}_{-\frac{T_a}{2}}^{\frac{T_a}{2}}△{\alpha }_a\left(t\right)s\left({\dot{\theta }}_{a}t\right)dt·s\left({\dot{\theta }}_{a}t\right)=△{\alpha }_{am}s\left({\theta }_a\right)\\ △{\beta }_a\left(t\right)=△{\beta }_{am}s\left({\theta }_a\right)\\ \begin{array}{c}△{\alpha }_p\left(t\right)=△{\alpha }_{pm}s\left({\theta }_p\right)\\ △{\beta }_p\left(t\right)=△{\beta }_{pm}s\left({\theta }_p\right)\end{array}\end{array},$$

(18)

Based on the kinematic analysis of axis perpendicularity errors earlier, the rotation transformation matrixes caused by run-outs $△R_a$ and $△R_p$ can be expressed as Equation (19):

$$\begin{gathered} △R_a=\left[\begin{array}{ccc}c△{\beta }_a\left(t\right)& 0& -s△{\beta }_a\left(t\right)\\ 0& 1& 0\\ s△{\beta }_a\left(t\right)& 0& c△{\beta }_a\left(t\right)\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& c△{\alpha }_a\left(t\right)& s△{\alpha }_a\left(t\right)\\ 0& -s△{\alpha }_a\left(t\right)& c△{\alpha }_a\left(t\right)\end{array}\right], \\ △R_p=\left[\begin{array}{ccc}c△{\beta }_p\left(t\right)& s△{\beta }_p\left(t\right)& 0\\ -s△{\beta }_p\left(t\right)& c△{\beta }_p\left(t\right)& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& c△{\alpha }_p\left(t\right)& s△{\alpha }_p\left(t\right)\\ 0& -s△{\alpha }_p\left(t\right)& c△{\alpha }_p\left(t\right)\end{array}\right], \end{gathered}$$

(19)

The angular velocity vectors of the azimuth and the pitch in $\left\{i\right\}$ affected by errors of run-outs are denoted as $△{\overrightarrow{\omega }}_{ia}$ and $△{\overrightarrow{\omega }}_{ip}$, and can be calculated by using Equation (20):

$$\begin{gathered} △{\overrightarrow{\omega }}_{ia}=R_b^a△R_a{\overrightarrow{\omega }}_{ib}+{\overrightarrow{\omega }}_{ba}, \\ △{\overrightarrow{\omega }}_{ip}=R_a^p△R_p△{\overrightarrow{\omega }}_{ia}+{\overrightarrow{\omega }}_{ap}, \end{gathered}$$

(20)

Based on analysis above, the signals of the gyroscopes affected by $△{\overrightarrow{\gamma }}_a$ and $△{\overrightarrow{\gamma }}_p$ are simulated in Simulink as shown in Figure 12. Taking the azimuth, for example, based on Equation (18) and Equation (20), ${\overrightarrow{\omega }}_{ib}$ and the motor rotation angle ${\theta }_a$ are set as the inputs; the maximum tilting oscillations including $△{\alpha }_{am}$ and $△{\beta }_{am}$, and the transformation matrix $△R_a$ are imbedded in the middle; the output is $△R_a{\overrightarrow{\omega }}_{ib}$, which will be used in other modules. The value of $△{\alpha }_{am}$ and $△{\beta }_{am}$ can be set as needed in the future simulation experiments. The value of ${\theta }_a$ can be obtained from the motor module in Section 3.4. The simulation module for the pitch frame is similar.

3.5.3. Modeling and Simulation of Assembly Errors of Gyroscope Locations

The gyroscopes are the core components for realizing the stabilizing function, and the assembly of them will directly affect the stability accuracy. The coordinates of the gyroscopes $\left\{ga\right\}$ and $\left\{gp\right\}$ are ideally parallel to $\left\{p\right\}$. Due to bad assembly, there are angular deviations between the real and the ideal locations. The deviations are denoted as ${\overrightarrow{\epsilon }}_a={\left[{\epsilon }_{ax},{\epsilon }_{ay},0\right]}^T$ for the azimuth gyroscope and ${\overrightarrow{\epsilon }}_p={\left[{\epsilon }_{px},0,{\epsilon }_{pz}\right]}^T$ for the pitch one. Based on the earlier analysis, it is easy to infer the rotation transformation matrixes caused by assembly errors of gyroscope locations $R_{ag}$ and $R_{pg}$, which are expressed as Equation (21).

$$\begin{gathered} R_{ag}=\left[\begin{array}{ccc}c{\epsilon }_{ay}& 0& -s{\epsilon }_{ay}\\ 0& 1& 0\\ s{\epsilon }_{ay}& 0& c{\epsilon }_{ay}\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& c{\epsilon }_{ax}& s{\epsilon }_{ax}\\ 0& -s{\epsilon }_{ax}& c{\epsilon }_{ax}\end{array}\right], \\ {R}_{pg}=\left[\begin{array}{ccc}c{\epsilon }_{pz}& s{\epsilon }_{pz}& 0\\ -s{\epsilon }_{pz}& c{\epsilon }_{pz}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& c{\epsilon }_{px}& s{\epsilon }_{px}\\ 0& -s{\epsilon }_{px}& c{\epsilon }_{px}\end{array}\right], \end{gathered}$$

(21)

Based on Equation (8) and Equation (21), the sensitive velocities of the two gyroscopes in $\left\{i\right\}$ affected by $R_{ag}$ and $R_{pg}$ are denoted as ${\omega }_{gyro\_z}^{\prime }$ and ${\omega }_{gyro\_y}^{\prime }$, and they could be calculated by using Equation (22).

$$\begin{gathered} {\omega }_{gyro\_z}^{\prime }=\left[0,0,1\right]R_{ag}{\overrightarrow{\omega }}_{ib}^p, \\ {\omega }_{gyro\_y}^{\prime }=\left[0,1,0\right]R_{pg}{\overrightarrow{\omega }}_{ib}^p, \end{gathered}$$

(22)

By using a typical second-order model, the transformation function between the input and output signals of gyroscopes $Tran{s}_g\left(s\right)$ could be expressed as Equation (23), where $K_g$, ${\omega }_g=2\pi f_g$ and ${\zeta }_g$ represent the amplification coefficient, natural frequency, and damping coefficient of the gyroscopes.

$$Tran{s}_g\left(s\right)=\frac{K_g{\omega }_g^2}{s^2+2{\zeta }_g{\omega }_{g}s+{\omega }_g^2},$$

(23)

Based on Equation (23) and analysis above, the signals and assembly errors of the gyroscopes are simulated in Simulink as shown in Figure 13. Based on Equation (22), ${\overrightarrow{\omega }}_{ib}^p$ is set as the input. Furthermore, the value of ${\overrightarrow{\omega }}_{ib}^p$ can be considered equal to the real-time value of ${\overrightarrow{\omega }}_{ip}$, which can be obtained from the kinematics module in Section 3.3. The assembly error vectors of ${\overrightarrow{\epsilon }}_a$ and ${\overrightarrow{\epsilon }}_p$, the transformation matrixes $R_{ag}$ and $R_{pg}$, and $Tran{s}_g\left(s\right)$ are imbedded in the middle. The final sensitive velocities of the two gyroscopes ${\omega }_{gyro\_z}^{\prime }$ and ${\omega }_{gyro\_y}^{\prime }$ are set as the output, and they will be used in other modules. The value of ${\overrightarrow{\epsilon }}_a$ and ${\overrightarrow{\epsilon }}_p$ can be set as needed in the future simulation experiments, and also can be measured by using high-precision instruments and relevant methods [23,24].

3.6. Integrating Simulation of the Overall System Stability Accuracy

By integrating all the above modules, the overall system stability accuracy is modeled in Simulink, which is shown in Figure 14. The original system base angular velocity ${\overrightarrow{\omega }}_{ib}={\left[{\omega }_{ibx},{\omega }_{iby},{\omega }_{ibz}\right]}^T$ is set as the input. All the modules of the kinematics, components, and assembly errors are imbedded in the middle. The final index of stability accuracy $△\theta$ is set as the output by referring to Equation (1) and Equation (2).

4. Simulation Experiments and Analysis Based on Error Sensitivity

In the real application process of the system, those disturbances caused by UAV, environment factors such as wind, noises, etc., will all couple to the system pedestal base in the end. They are all integrated and defined as the base disturbances (BD) for simplicity. Apparently, in the normal working process of the system, BD is the strongest factor affecting the stability accuracy. That effect will be deteriorated by those modeled assembly errors. According to previous definitions in Section 3.2, the angular velocities of BD disturbances are denoted as ${\overrightarrow{\omega }}_{ib}={\left[{\omega }_{ibx},{\omega }_{iby},{\omega }_{ibz}\right]}^T$.

By carrying out dynamic simulations, the effects caused by BD and the modeled assembly errors can be analyzed. The analysis results will provide valuable guidance for future possible improvements. For simulation study, according to Equation (1), the sampling period $T_s$ is equal to the Simulink step size, and set as $1ms$. For other parameter settings: Total sampling time: $t_s=30s$; Sampling number: $n_s=30000$. By referring to real measurement data of a newly developing system, the necessary parameters are predefined before the simulation experiments. Those parameters are listed in

Table 1. Necessary measured parameters before simulation experiments.

Parameters	Descriptions	Values	Units
$J_{ma}$	Rotary inertia of the azimuth motor	0.02	$\mathrm{kg}·{\mathrm{m}}^2$
$J_{La}$	Rotary inertia of the azimuth load	0.04	$\mathrm{kg}·{\mathrm{m}}^2$
$J_{mp}$	Rotary inertia of the pitch motor	0.01	$\mathrm{kg}·{\mathrm{m}}^2$
$J_{Lp}$	Rotary inertia of the pitch load	0.02	$\mathrm{kg}·{\mathrm{m}}^2$
$R_m$	Total resistor of the two armature circuits	8.60	$\Omega$
$L_m$	Total inductance of the armature circuits	0.01	H
$D_L$	Damping coefficient of the two motors	10	/
$K_L$	Elastic coefficient of the two motors	≈100000	/
$C_e$	EMF coefficient of the two motors	0.33	/
$C_m$	Motor torque coefficient of the two motors	0.33	/
$K_p$	Amplification coefficient of the driving circuit	≈20	/
$R_s$	Resistor of the driving circuit	0.4	$\Omega$
$K_g$	Amplification coefficient of the gyroscopes	15/π	/
${\omega }_g$	Natural frequency of the gyroscopes	≈200π	rad/s
${\zeta }_g$	Damping coefficient of the gyroscopes	0.7	/

, and the prototype of the product is shown in Figure 15.

To study the sensitivity of the assembly error effect on the stability accuracy, the simulation experiments affected by each kind error source was carried out separately. Specifically, the simulations of each type of error source are carried out with other errors are equal to zero.

The variables of each group of simulations will be the levels of each error source, and the disturbance frequencies of BD. The two variables are determined by referring to principles as follows: according to engineering experiences and the real conditions of manufacturing accuracies, the error source levels are set to range from $1m^{°}$ to $50m^{°}$. By referring to [25,26], the angular velocities of the three axes of BD disturbances, specifically ${\overrightarrow{\omega }}_{ib}={\left[{\omega }_{ibx},{\omega }_{iby},{\omega }_{ibz}\right]}^T$, are all simulated as sinusoidal signals of low frequencies, and they range from ${1.5}^{°}/0.5Hz$ to ${1.5}^{°}/4Hz$. The simulation results of stability accuracy obtained $△\theta$ are presented with three decimal places, and listed in the central area of

Table 2. Simulation results under BD Dis of ${1.5}^{°}/0.5Hz$.

Stability Accuracy $△\mathit{\theta }\left[\mathit{\mu }\mathit{r}\mathit{a}\mathit{d}\right]$ (Data in the Central Area with Three Decimal Places)		$\mathbf{Error}\mathbf{Levels}\left[{\mathit{m}}^{°}\right]$ (For Each Column, the First Data Is the Error Level of This Column, and for the Following Each Line of This Column, Each Relevant Error Type on the Left Is Equal to the Error Level, with All of the Other Errors Are Equal To 0 Simultaneously.)
		1	5	10	20	50
Error Types (For each line, the error level is increasing corresponding to the error level on the top, and other error types are all equal to 0 simultaneously)	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{x}}$	112.117	112.143	112.056	112.127	112.261
	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{y}}$	112.074	112.024	112.257	112.103	112.108
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{x}}$	112.359	113.085	113.733	115.457	120.446
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{z}}$	112.222	112.148	112.153	112.138	112.104
	$△{\mathit{\alpha }}_{\mathit{a}\mathit{m}}$	112.155	112.205	112.065	112.123	112.111
	$△{\mathit{\beta }}_{\mathit{a}\mathit{m}}$	112.191	112.085	112.096	112.145	112.069
	$△{\mathit{\alpha }}_{\mathit{p}\mathit{m}}$	112.101	112.178	112.235	112.044	112.066
	$△{\mathit{\beta }}_{\mathit{p}\mathit{m}}$	112.146	112.155	112.097	112.150	112.109
	${\mathit{\epsilon }}_{\mathit{a}\mathit{x}}$	112.028	112.183	112.094	112.123	112.144
	${\mathit{\epsilon }}_{\mathit{a}\mathit{y}}$	111.922	111.302	110.480	108.808	103.845
	${\mathit{\epsilon }}_{\mathit{p}\mathit{x}}$	112.057	112.212	112.121	112.100	112.077
	${\mathit{\epsilon }}_{\mathit{p}\mathit{z}}$	112.216	112.331	112.689	113.233	115.184

,

Table 3. Simulation results under BD Dis of ${1.5}^{°}/1Hz$.

Stability Accuracy $△\mathit{\theta }\left[\mathit{\mu }\mathit{r}\mathit{a}\mathit{d}\right]$		Error Levels $\left[{\mathit{m}}^{°}\right]$
		1	5	10	20	50
Error Types	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{x}}$	39.544	39.552	39.527	39.558	39.526
	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{y}}$	39.535	39.492	39.499	39.534	39.524
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{x}}$	39.614	39.929	40.316	41.128	43.521
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{z}}$	39.546	39.535	39.539	39.531	39.520
	$△{\mathit{\alpha }}_{\mathit{a}\mathit{m}}$	39.540	39.516	39.515	39.518	39.505
	$△{\mathit{\beta }}_{\mathit{a}\mathit{m}}$	39.496	39.528	39.525	39.571	39.564
	$△{\mathit{\alpha }}_{\mathit{p}\mathit{m}}$	39.511	39.548	39.526	39.530	39.508
	$△{\mathit{\beta }}_{\mathit{p}\mathit{m}}$	39.536	39.538	39.548	39.523	39.555
	${\mathit{\epsilon }}_{\mathit{a}\mathit{x}}$	39.519	39.531	39.530	39.532	39.542
	${\mathit{\epsilon }}_{\mathit{a}\mathit{y}}$	39.434	39.093	38.755	37.982	35.669
	${\mathit{\epsilon }}_{\mathit{p}\mathit{x}}$	39.521	39.511	39.518	39.542	39.505
	${\mathit{\epsilon }}_{\mathit{p}\mathit{z}}$	39.567	39.727	39.922	40.366	41.700

,

Table 4. Simulation results under BD Dis of ${1.5}^{°}/2Hz$.

Stability Accuracy $△\mathit{\theta }\left[\mathit{\mu }\mathit{r}\mathit{a}\mathit{d}\right]$		Error Levels $\left[{\mathit{m}}^{°}\right]$
		1	5	10	20	50
Error Types	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{x}}$	15.892	15.898	15.885	15.892	15.897
	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{y}}$	15.884	15.896	15.901	15.898	15.899
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{x}}$	15.933	16.075	16.267	16.626	17.780
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{z}}$	15.890	15.884	15.890	15.890	15.892
	$△{\mathit{\alpha }}_{\mathit{a}\mathit{m}}$	15.903	15.894	15.901	15.895	15.891
	$△{\mathit{\beta }}_{\mathit{a}\mathit{m}}$	15.889	15.907	15.896	15.904	15.899
	$△{\mathit{\alpha }}_{\mathit{p}\mathit{m}}$	15.886	15.890	15.889	15.902	15.893
	$△{\mathit{\beta }}_{\mathit{p}\mathit{m}}$	15.903	15.886	15.902	15.891	15.905
	${\mathit{\epsilon }}_{\mathit{a}\mathit{x}}$	15.897	15.881	15.896	15.904	15.898
	${\mathit{\epsilon }}_{\mathit{a}\mathit{y}}$	15.856	15.708	15.531	15.166	14.105
	${\mathit{\epsilon }}_{\mathit{p}\mathit{x}}$	15.896	15.880	15.904	15.888	15.896
	${\mathit{\epsilon }}_{\mathit{p}\mathit{z}}$	15.904	16.026	16.141	16.413	17.237

,

Table 5. Simulation results under BD Dis of ${1.5}^{°}/3Hz$.

Stability Accuracy $△\mathit{\theta }\left[\mathit{\mu }\mathit{r}\mathit{a}\mathit{d}\right]$		Error Levels $\left[{\mathit{m}}^{°}\right]$
		1	5	10	20	50
Error Types	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{x}}$	9.770	9.779	9.777	9.780	9.783
	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{y}}$	9.774	9.780	9.776	9.777	9.780
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{x}}$	9.800	9.891	10.014	10.248	10.975
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{z}}$	9.778	9.774	9.777	9.771	9.774
	$△{\mathit{\alpha }}_{\mathit{a}\mathit{m}}$	9.771	9.777	9.777	9.776	9.777
	$△{\mathit{\beta }}_{\mathit{a}\mathit{m}}$	9.779	9.776	9.778	9.774	9.779
	$△{\mathit{\alpha }}_{\mathit{p}\mathit{m}}$	9.781	9.778	9.772	9.782	9.771
	$△{\mathit{\beta }}_{\mathit{p}\mathit{m}}$	9.774	9.781	9.780	9.778	9.771
	${\mathit{\epsilon }}_{\mathit{a}\mathit{x}}$	9.780	9.777	9.770	9.774	9.782
	${\mathit{\epsilon }}_{\mathit{a}\mathit{y}}$	9.765	9.659	9.551	9.315	8.649
	${\mathit{\epsilon }}_{\mathit{p}\mathit{x}}$	9.771	9.776	9.778	9.774	9.776
	${\mathit{\epsilon }}_{\mathit{p}\mathit{z}}$	9.792	9.870	9.971	10.154	10.744

and

Table 6. Simulation results under BD Dis of ${1.5}^{°}/4Hz$.

Stability Accuracy $△\mathit{\theta }\left[\mathit{\mu }\mathit{r}\mathit{a}\mathit{d}\right]$		Error Levels $\left[{\mathit{m}}^{°}\right]$
		1	5	10	20	50
Error Types	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{x}}$	7.032	7.028	7.032	7.029	7.024
	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{y}}$	7.026	7.027	7.026	7.032	7.028
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{x}}$	7.045	7.116	7.197	7.371	7.913
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{z}}$	7.023	7.030	7.026	7.029	7.026
	$△{\mathit{\alpha }}_{\mathit{a}\mathit{m}}$	7.030	7.032	7.028	7.027	7.023
	$△{\mathit{\beta }}_{\mathit{a}\mathit{m}}$	7.029	7.028	7.025	7.028	7.030
	$△{\mathit{\alpha }}_{\mathit{p}\mathit{m}}$	7.028	7.030	7.032	7.024	7.028
	$△{\mathit{\beta }}_{\mathit{p}\mathit{m}}$	7.026	7.026	7.027	7.026	7.023
	${\mathit{\epsilon }}_{\mathit{a}\mathit{x}}$	7.027	7.024	7.025	7.028	7.026
	${\mathit{\epsilon }}_{\mathit{a}\mathit{y}}$	7.012	6.948	6.857	6.696	6.217
	${\mathit{\epsilon }}_{\mathit{p}\mathit{x}}$	7.031	7.026	7.023	7.029	7.026
	${\mathit{\epsilon }}_{\mathit{p}\mathit{z}}$	7.042	7.096	7.170	7.324	7.792

. The results of each kind of BD disturbance are listed in a separate table.

In those tables, for each column, the first piece of data is the error level of this column, and for the rest of the grid, each relevant error type on the left is equal to the error level, with all of the other errors equal to 0 simultaneously. For each line, the first grid is the error type of this line, and for the rest of the grid, the error increases corresponding to the error level on the top of the column, and other error types are still all equal to 0 simultaneously.

The standard deviations of the experimental samples under different error levels are defined as error sensitivity, which are denoted as $S\left(x\right)$. $x$ represents the objective error source. $S\left(x\right)$ is an index which is used to evaluate the fluctuations changing along with the error levels. Therefore, $S\left(x\right)$ can be used to evaluate the sensitivity of those error sources.

$S\left(x\right)$ is calculated by using Equation (24), where $n_l$ is the number of error levels; $△{\theta }_i\left(x\right)$ is the accuracy sample value of the $ith$ error level of $x$; $\overline{△\theta \left(x\right)}$ is the average value of $△{\theta }_i\left(x\right)$. The calculation results of $S\left(x\right)$ are listed in

Table 7. Error sensitivity of different error sources.

Error Sensitivity $\mathit{S}\left(\mathit{x}\right)$ (Data in the Center Area with Three Decimal Places)		$\mathbf{BD}\mathbf{Frequency}\left[\mathit{H}\mathit{z}\right]$ (For Each Column, the First Data Is the Value of the BD Frequency, and the Following Are the $\mathit{S}\left(\mathit{x}\right)$ of Each Error Type on the Left Side)
		0.5	1	2	3	4
Error Types (Each line contains the $S\left(x\right)$ of the error type under different BD frequency on the top)	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{x}}$	7.483	1.448	0.517	0.487	0.332
	${\mathit{\delta }}_{\mathit{b}\mathit{a}\mathit{y}}$	8.704	2.007	0.673	0.261	0.249
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{x}}$	324.476	156.998	74.231	47.222	34.866
	${\mathit{\delta }}_{\mathit{a}\mathit{p}\mathit{z}}$	4.305	0.968	0.303	0.277	0.277
	$△{\mathit{\alpha }}_{\mathit{a}\mathit{m}}$	5.213	1.287	0.502	0.261	0.339
	$△{\mathit{\beta }}_{\mathit{a}\mathit{m}}$	5.007	3.078	0.704	0.217	0.187
	$△{\mathit{\alpha }}_{\mathit{p}\mathit{m}}$	7.989	1.612	0.612	0.507	0.297
	$△{\mathit{\beta }}_{\mathit{p}\mathit{m}}$	2.646	1.223	0.838	0.421	0.152
	${\mathit{\epsilon }}_{\mathit{a}\mathit{x}}$	5.817	0.817	0.853	0.477	0.158
	${\mathit{\epsilon }}_{\mathit{a}\mathit{y}}$	325.078	150.672	70.202	44.542	31.868
	${\mathit{\epsilon }}_{\mathit{p}\mathit{x}}$	6.014	1.408	0.912	0.265	0.308
	${\mathit{\epsilon }}_{\mathit{p}\mathit{z}}$	121.427	86.072	53.342	38.140	30.272

. The error sensitivities $S\left(x\right)$ of different error sources changing along with the BD disturbance frequencies are shown in Figure 16.

$$S\left(x\right)=100\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^{n_l}{\left[△{\theta }_i\left(x\right)-\overline{△\theta \left(x\right)}\right]}^2}{n_l-1}},\overline{△\theta \left(x\right)}=\frac{{{\displaystyle \sum }}_{i=1}^{n_l}△{\theta }_i\left(x\right)}{n_l},$$

(24)

From the above tables and Figure 16, the most sensitive error sources are identified as x-axial component of pitch axis perpendicularity deviation ${\delta }_{apx}$, y-axial component of the azimuth gyroscope assembly deviation ${\epsilon }_{ay}$, and z-axial component of the pitch gyroscope assembly deviation ${\epsilon }_{pz}$. Error sensitivities of the identified three errors are ranked as the order of $S\left({\delta }_{apx}\right)>S\left({\epsilon }_{ay}\right)>S\left({\epsilon }_{pz}\right)$. Additionally, with the rising of disturbance frequency, the stability accuracies $△\theta$ can be improved, and the error sensitivities $S\left(x\right)$ decrease.

5. Discussions, Conclusions, and Future Work

The major contributions of this research are listed as follows: (1) Stability accuracy of UAV detection system is comprehensively modeled and simulated. The system’s general kinematics, compensation principle, and major components including servo-motors and gyroscopes are all covered. (2) Assembly errors of axis perpendicularities, run-outs, and gyroscopes are defined, modeled, and introduced into the overall model. The effects caused by various assembly error sources can be comprehensively analyzed. (3) The effects caused by those assembly errors under different disturbance frequencies are studied. The simulation experiment results provide valuable data for real application scenarios. (4) Error sensitivity is proposed by carrying out standard deviation analysis. Based on the analysis, the most sensitive assembly errors are identified, and ranked in order of sensitivity as follows: x-axial component of pitch axis perpendicularity deviation, y-axial component of the azimuth gyroscope assembly deviation, and z-axial component of the pitch gyroscope assembly deviation. In real processes of manufacturing and assembly, the analysis results can be used to determine standards and priorities for possible improvements.

The limits of this research are listed as follows: (1) So far, only assembly errors and UAV angular disturbances are studied. Other error sources and situations to further verify the proposed modeling and analysis method should be considered. (2) Simulation results provide only limited information, and further validations based on real experiments are needed. (3) The limitations, constraints used in this research, and the possibility of uncertainty will affect the results, such as the changes of error distributions, UAV disturbances, and target locations. (4) The effects of series versus parallel (or hybrid) component arrangement are not considered. (5) The effects of interface failures are not covered.

As the next step, validations based on real experiments will be carried out. Other error sources and situations will be investigated, such as frictions, structural eccentricity, and inertia coupling. The linear disturbances of UAV and target movements may have great influence on the results. The effects of series versus parallel (or hybrid) component arrangement will be further studied. The interface performance of the system is very important. The reasons for interface failure include damage or wrong assembly of hardware, incompatibilities of data formats, and inferiority of control algorithms, and will be studied in detail. A similar control mechanism of a more complex system with more dimensions or new materials will be studied, such as the structures of two-axis–four-frame and three-axis–four-frame used on larger carriers. How to incorporate the proposed method to other working models of the system will be studied, such as target identification, tracking, positioning, etc.