Real-time onboard propeller fault diagnosis of autonomous delivery drones through vibration analysis

Abstract

Delivery drones have become increasingly important in recent years. It is advantageous for commercialization that suppliers are able to deliver orders autonomously and directly to their customers via air transport. However, the safety aspect must be considered. Real-time inspection of delivery drones during operation helps preventing accidents and a threat to civilians. For continuous monitoring, the sensors must be installed on the drone throughout the flight. A promising approach uses the inherent excitation of the servomotors for vibration-based Structural Health Monitoring (SHM). Vibrations can be recorded using triaxial acceleration sensors and analyzed using suitable methods such as stochastic subspace-based fault detection or histogram difference. In comparison to nondestructive testing, reference measurements of the intact structure are necessary in SHM. As soon as laboratory conditions no longer exist and environmental parameters, for example, wind, influence the vibration spectrum, classification methods are necessary for compensation. The recording of comprehensive reference datasets with different environmental conditions is limited by the battery life. This work focuses on diagnosing irreversible rotor blade damages of an 8.5 kg delivery drone. A parametric analysis taking into account systematic fusion of damage indicators and a specific number of considered references were determined for this purpose in order to assess the severity of the existing damage. Onboard SHM to evaluate the airworthiness in real time for linear and hovering flights was achieved.

Introduction

The importance of drones for private and industrial purposes has grown rapidly in recent years. In Europe alone, it is estimated that a total of 400,000 commercial drones and around 7 million drones for hobby purposes will be authorized in the next 30 years. Commercial use in particular focuses on the autonomous operation of drones without pilots, as this is economically favorable. EU Implementing Regulation 2019/947 of May 24, 2019 makes it clear that the operation of unmanned aircraft should be just as safe as manned aviation. However, drones (both airplane and multi-rotor types) that have been flying for some time are already showing mechanical reliability problems and related failures. ¹ An estimated 9.2% of all commercial drone failures are due to a malfunction of the mainframe. ² The causes are mainly due to material fatigue in aluminum, collision damage in composite material, or cracks in 3D-printed plastic parts. Sixty-four percent $64$ of accidents and incidents are caused by equipment problems, for example, actuator faults, ³ while $34$ % are due to human errors. ⁴ In general, the failure rate of drones is approximately $25$ %. ⁵ At this point, Structural Health Monitoring (SHM) becomes important.

Vibration-based methods are often used in the literature for global damage identification in complex material systems and structures. ⁶ Related works dealt with the vibration-based monitoring of the railway bridge KW51 between Leuven and Brussels (Belgium), ⁷ the prototype of a 5 MW wind turbine on the North Sea (Germany), ⁸ the Warren-truss bridge in Japan, ⁹ a bridge on the SS114 in Silicy (Italy), ¹⁰ the Z24 bridge between Zurich and Bern (Switzerland), ¹¹ the Tianjin-Yonghe bridge in China, ¹² and the ARTEMIS satellite antenna. ¹³ Output-only methods are particularly relevant here, as usually only the measured output signals are available for SHM. ¹⁴ Damage indicators (DIs) are proposed as dimensionless scalars that are calculated to assess the severity of the existing structural condition. ⁹

In recent years, algorithms based on modal parameters have become increasingly influential for vibration-based SHM. Modal parameters generally include natural frequencies, damping ratios and mode shapes. ⁹ The automatic, data-based calculation of modal parameters in the operating state of a material system is referred as operational modal analysis. ⁸

The challenge is to discriminate structural changes, that is, damage, from changing environmental and operational influences such as temperature variations. ¹⁵ For special applications, these parameters may also be unknown. In order to be able to observe and analyze various weather-related influences outside of a laboratory in all seasons, a measurement period of at least 1 year is worthwhile. The use of machine learning could be helpful here to train a neural network without application-specific physical models. However, it should be noted that especially microcontroller units (MCU) do not have a large working memory and therefore the complexity must be kept to a minimum, although this leads to a loss of accuracy. Nevertheless, Zonzini et al. were able to prove that a one class classifier neural network on a MCU reach an accuracy of $94$ % with high performance. ¹¹

One part of current research is to design vibration-based SHM systems on drones for their damage identification during operation. Most current publications report on offline SHM of drones. ¹⁶ This means that the evaluation of the raw data is only carried out on an external computer after the measurement procedure. Ghazali and Rahiman developed an initial approach to an online monitoring system of drones to assess the sensitivity of various vibration-based sensors with respect to different rotor blade damage on a Storm Drone 8. Experimental studies were carried out under laboratory conditions in the “warm-up stage” (propellers rotated with the drone on ground). Propeller damage was detected in real-time when a threshold of the vibration-based time signals was exceeded. The data were transmitted to a smartphone via Bluetooth for further processing. ¹⁷ Ghalamchi et al. used a Kalman filter to take into account the nonuniform distribution of forces in real time that occur from the servomotors due to rotor blade damage. Two quadrocopters and one hexacopter of different masses with the same flight distances were tested. ¹⁸ The analysis of the vibration behavior due to mass increases was performed on a T50 helicopter by Bektash and la Cour-Harbo ¹⁹ and on a Q8000 delivery drone by Ibrahim et al. ²⁰ An experiment that followed the complete delivery process of a package via a ferromagnetic unloading station with a DJI F550 hexacopter was set up and carried out by Carroll et al. ²¹

Damage sizes in the range of centimeters on propellers are reported in the literature. Liu et al. ²² shortened up to $3.6$  cm, Gomez et al. ²³ to $2$  cm and Ghazali and Rahiman ¹⁷ up to $3.5$  cm. Compared to the healthy structure of each propeller, the percentages of shortening are $14$ %, $17$ %, and $50$ %.

The novelty in this article is to operate and demonstrate the entire onboard SHM system for autonomous delivery drones in real-time during linear and hovering flights on the airfield outside. In this work, we considered a shortening of the propeller to demonstrate the overall damage detection performance of the proposed SHM system. It should be noted that this type of propeller damage is a simplification compared, for example, to fatigue damage. In addition to the onboard sensors, the results are transmitted and displayed through a wireless communication link to a Rasperry Pi and an Android smartphone. The algorithms for calculating DI have been extended by a classification method to extract information on wind conditions and compensation maneuvers of the control technology via the motor control signals. Another new feature is the combination of DI, calculated by the stochastic subspace-based fault detection (SSFD) and histogram difference (HD) method to improve damage detection performance. Two case studies have been performed and the results compared, that is, linear and hovering flights. Linear flights are characterized by a straight-line movement at a constant speed, so that the resulting acceleration of the drone system is zero due to the principle of inertia. In hovering flights, the drone maintains its position above the ground.

The remainder of the article is organized in the following way: “Experimental setup” section describes the experimental setup, which is divided into the sensors on the drone and the ground station. “Damage identification techniques” section gives an overview of implemented algorithms for damage detection: The SSFD and HD method, their combination and the classification method. “Experimental results: onboard propeller fault diagnosis during flight in real-time” section presents the experimental results for the linear and hovering flights separately. Finally, a short summary and an outlook for further research is given in “Conclusions” section.

Experimental setup

The delivery drone to be monitored for structural changes in the present work is a Q8000 octocopter manufactured by Emqopter GmbH in Würzburg (Germany). The battery for the power supply is a 6S LiPo battery that provides a voltage of $25$  V and a flight time of $12$  min on average. According to the manufacturer, the weight of the drone is approximately $8.5$  kg, whereby a total weight of $12$  kg should not be exceeded either with or without an additional transport box.

An STM32 NUCLEO-F767ZI MCU from STMicroelectronics is attached above the Emqopter circuit boards, to which the sensors are connected. The MCU itself is based on a $32$ bit architecture with a CPU operation frequency of $216\mathrm{MHz}$ . The working memory includes a $2\mathrm{MB}$ flash memory and a $512$ kbit SRAM. One hundred forty-four $144$ free assignable pins can be connected to arbitrary sensors with jumper cables. The photographs of the drone, the MCU, and the acceleration sensor on the airfield in Ober-Mörlen (Germany) are shown in Figure 1. The complete sensor network, including the functionalities and interfaces that are used, is depicted in Figure 2 as a flow chart.

Figure 1.

Photo of the drone showing the microcontroller board under the dome and the position of the accelerometer.

Figure 2.

System concept for real-time onboard SHM of autonomous drones through vibration analysis.

The used triaxial ADXL357 accelerometer from ANALOG DEVICES measures accelerations along three axes with a sampling rate of $4$  kHz and with the chosen $\pm 20g$ setting, where $g$  = 9.81 ms⁻² is the gravitational acceleration. In addition, temperatures from $-{40}^{°}\mathrm{C}$ to ${125}^{°}\mathrm{C}$ can be measured. The sensor is connected to a SPI bus interface of the MCU. The acceleration data are temporarily locally stored in a FIFO that is a memory unit with $96$ registers for $21$  bit data, where $20$ bit per register are reserved for acceleration data and $1$ bit for marking the x-axis. If $60$ registers of the FIFO are filled, an interrupt will be triggered for reading with a pointer, starting with the oldest recorded bit. The frequency of the external clock generator must correspond to $1024$  kHz and is provided by the MCU for synchronization purposes.

The bytes of the FIFO are stored in two buffers implemented on the MCU. A buffer is a memory unit for storing acceleration data of a complete measurement interval. If the storage capacity of buffer A is exhausted, the other buffer B is used for storage. If the storage capacity of buffer B is also exhausted, buffer A is overwritten with current data to avoid overloading the working memory with further buffers. After reading, the FIFO can be overwritten with further data. The baud rate specifies here how many characters per second are transferred from the FIFO to the MCU via jumper cables. It must be well adjusted to the cable length so that the bits are not detected with the next clock signal.

A DS3231 real-time clock (RTC) from AZ-Delivery is connected to an I²C bus. It must be set correctly during initialization and remains in operation after the MCU is switched off due to the $3$ V button cell. The task of the RTC is to record the time stamp before, and after the evaluation of a measurement and to check whether real-time operation is maintained.

The micro SD card reader from pzsmocn is used to write data to the inserted 32 GB Samsung Evo+ micro SD card by a maximum of four channels. Supported interfaces are SPI or SDIO interfaces (for this MCU, the interface is called SDMMC). The file system module is FatFs with extensive functions that can be activated in the IDE.

UART interfaces are used to transmit and receive data packages. They are required for the WiFi board module ESP8266-01S from AZ-Delivery and for retrieving the flight parameters with a sampling rate of $10$  Hz from the flight controller that is called Pixhawk. For smooth communication, the baud rate between the transmitter and the receiver must be the same, as there are no clock signals with UART. It is 460,800 for the Wi-Fi board and 115,200 for the Pixhawk, respectively.

To display the monitoring results, data are sent to a Raspberry Pi 4 model B by the Wi-Fi module from the MCU and forwarded to an evaluation app on an Android smartphone via the transmission control protocol in small data packages. Both devices form the main part of the ground station. Furthermore, a remote control and a notebook are used to fly the drone manually and automatically by an autopilot in the MissionPlanner software. In this approach, the Raspberry Pi acts as a server and uses a listen call to the socket interface to check whether the MCU or the smartphone would like to connect as separate clients via a connect call. Before the data exchange between two instances can begin, the server must accept the connection request. The connection can theoretically exists for any length of time or be terminated by the server via a close call.

For the structural analysis, two series of measurements were recorded on two different days under different environmental conditions and aircraft movements. In order to simulate damage, a part of a rotor blade was clipped off with pliers in each series of measurements so that a total of three structural states could be compared with each other (see Figure 3). The measurement plan can be found in Table 1.

Figure 3.

Documentation of propeller damage during (a) hovering and (b) linear flight.

Table 1.

Measurement plan, divided into two series of measurement.

Series of measurement	Flight mode	Type of measurement	Duration
No. 1: Rotor blade shortening at servomotor no. 6 (August 10, 2023)	Linear flight	Baseline 2.6 cm shortening 5.2 cm shortening	4 × 12 s 3 × 12 s 3 × 12 s
No. 2: Rotor blade shortening at servomotor no. 3 (February 14, 2023)	Hovering flight	Baseline 2.3 cm shortening 4.8 cm shortening	3 × 20 s 2 × 20 s 3 × 20 s

Damage identification techniques

This section gives an overview about the procedure of real-time damage identification in three phases, as shown in Figure 4. In the first phase, the vibrations of the drone are measured using an acceleration sensor with a sampling rate of f _s  = 4 kHz. In parallel, the motor control signals of the eight servomotors from the Pixhawk and the temperature data from the acceleration sensor are retrieved with a sampling rate of f _n  = 10 Hz. The raw data are stored on the micro SD card so that it is available for the SHM algorithms and possible offline evaluation.

Figure 4.

Workflow of the SHM system including data acquisition, training, and monitoring phase.

A measurement interval of $1$  s was selected in which the MCU is able to perform the calculation of DI with the limited working memory. The number of acceleration data and motor control signals within the measurement interval is therefore $N_{\mathit{s}}=4000$ and $N_{\mathit{n}}=100$ , respectively. Real-time capability is characterized by the fact that the SHM algorithms are run immediately after recording a 1 s measurement interval, while the next 1 s measurement dataset is being recorded. The evaluation must therefore only require 1 s so that the next evaluation can start immediately afterwards.

In the second phase, the SHM system is trained on references in order to monitor the drone for structural integrity in the third phase. References are datasets which are recorded under healthy rotor blade conditions and collected during the training phase. SSFD and HD were chosen as SHM algorithms, in which the acceleration signals of the reference structure and the damaged structure are linked by mathematical operations to calculate DI. While the SSFD is based on the cross-correlations of the acceleration signals in time domain, the HD subtracts the frequency distributions of two structural states for each bin.

The choice of algorithms was related to real-time capability and flexible parameter setting. Modal methods were unstable relating to the short measurement duration of $1$  s so that the focus was on time domain methods. The SSFD was also used successfully in many applications for vibration-based SHM.^13,20,24,25 The strength of HD lies in its time-efficient evaluation.

For both algorithms, a different number of reference datasets $N$ can be considered for statistical certainty. Three cases are separated here, whether only the optimal ( $N=1$ ) reference, all references (N = N_R) or an intermediate form ( $1<N<N_{\mathrm{R}}$ ) are considered. In the latter case, the real-time criterion is not fulfilled, as a ranking must be created according to the optimal references, of which there are N_R! options for arrangement depending on the number of references N_R. An optimal reference is characterized by the best approximation of the environmental conditions relative to the monitoring dataset. This is achieved in the classification method, described below, by calculating centroids and distances from the motor control signals.

In order to select a number of suitable references, the system must be trained for many different environmental and operational conditions in the baseline state of the structure. ²⁴ For this reason, reference datasets are analyzed with the SHM algorithms in the training phase and calculated matrixes are saved for subsequent use in the monitoring phase. The decision whether a structure is intact or damaged is made by comparing the DI from the training and monitoring phase. In the case of a damaged structure, the DI would increase.

The algorithms are explained mathematically in the following subsections. Acceleration signals $\mathit{s}$ are assumed for this purpose. Due to the different sensitivity and orientation of the spatial axes of the triaxial accelerometer, the DI are calculated for each spatial axis separately. Moreover, the gravitational acceleration influences $\mathit{s}$ differently for each spatial direction.

Stochastic subspace-based fault detection

This subsection is based on work from Ibrahim et al., ²⁰ Strumpen et al., ²⁶ Basseville et al., ²⁷ and Balmes et al. ²⁸ In order to normalize the acceleration data ${\mathit{s}}^{\prime }$ component-by-component, the first step is to carry out the transformation into a standard normal distribution:

$$\mathit{s}=\frac{{\mathit{s}}^{\prime }-{\mu }_{{\mathit{s}}^{\prime }}}{{\sigma }_{{\mathit{s}}^{\prime }}},$$

(1)

where ${\mu }_{{\mathit{s}}^{\prime }}$ is the expected value and ${\sigma }_{{\mathit{s}}^{\prime }}$ the standard deviation of ${\mathit{s}}^{\prime }$ . After that, the acceleration signals are normalized to the expected mean value 0 ms⁻² and the variance 1 m² s⁻⁴.

In the second step, the Hankel matrix has to be calculated. Generally, it has the form

$${\mathit{H}}_0=\left(\begin{array}{cccc}{R}_0\hfill & R_1\hfill & \cdots \hfill & R_{{N_{\mathrm{cc}}}_{-}1}\hfill \\ \multicolumn{1}{c}{R_1}& R_2\hfill & \cdots \hfill & R_{N_{\mathrm{cc}}}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{R_{{N_{\mathrm{cc}}}_{-}1}}& R_{N_{\mathrm{cc}}}\hfill & \cdots \hfill & R_{2{N_{\mathrm{cc}}}_{-}2}\hfill \end{array}\right),$$

(2)

where

$$R_j(s[t],s[t])=\frac{1}{N_{\mathit{s}}-j}\sum _{i=0}^{N_{\mathit{s}}-j-1}s[t_{i+j}]s[t_i]$$

(3)

are the cross-correlations in time domain for $0\le j\le$ N_cc $-1$ . N_cc stands for the number of considered cross-correlations and is set to $600$ . The structure of ${\mathit{H}}_0$ shows that it is symmetrical:

$${\mathit{H}}_0^T={\mathit{H}}_0$$

(4)

with $(2N_{\mathrm{cc}}-1)$ independent elements.

In the third step, the zero matrix has to be computed that is the kernel of the reference Hankel matrix H ₀ because of its annihilation:

$${\mathit{K}}_0·{\mathit{H}}_0^{\mathrm{R}}=0.$$

(5)

For N_R references, the elements of the same row $k$ and column $l$ of the Hankel matrixes are averaged:

$${\overline{H}}_{0,\{k,l\}}^{\mathrm{R}}=\frac{1}{N_{\mathrm{R}}}\sum _{i=0}^{{N_{\mathrm{R}}}_{-}1}{H}_{0,i,\{k,l\}}^{\mathrm{R}}.$$

(6)

Numerically, ${\mathit{K}}_0$ is determined by a singular value decomposition (SVD) algorithm. It converts the Hankel matrix into the form

$${\mathit{H}}_0=\mathit{US}{\mathit{V}}^T,$$

(7)

where $\mathit{S}$ is a diagonal matrix consisting of singular values arranged in descending order, and $\mathit{U}$ and $\mathit{V}$ are orthogonal transformation matrixes with the property ${\mathit{U}}^T\mathit{U}={\mathit{V}}^T\mathit{V}=\mathit{I}$ ( $\mathit{I}$ is the identity matrix). From the corresponding relations

$${\mathit{H}}_0^{\mathrm{R}}\mathit{v}=\mathit{su},$$

(8a)

$${\left({\mathit{H}}_0^{\mathrm{R}}\right)}^T\mathit{u}=\mathit{sv},$$

(8b)

one can derive the following eigenvalue problems:

$${\left({\mathit{H}}_0^{\mathrm{R}}\right)}^T{\mathit{H}}_0^{\mathrm{R}}\mathit{v}={\mathit{s}}^2\mathit{v},$$

(9a)

$${\mathit{H}}_0^{\mathrm{R}}{\left({\mathit{H}}_0^{\mathrm{R}}\right)}^T\mathit{u}={\mathit{s}}^2\mathit{u}.$$

(9b)

From this, it can be concluded that the singular values are the square roots of the eigenvalues ${\mathit{s}}^2$ and the corresponding eigenvectors $\mathit{u}$ and $\mathit{v}$ are column vectors of $\mathit{U}$ and $\mathit{V}$ . Because of symmetry of ${\mathit{H}}_0$ in Equation (4), it also applies ${\left({\mathit{H}}_0^{\mathrm{R}}\right)}^T{\mathit{H}}_0^{\mathrm{R}}={\mathit{H}}_0^{\mathrm{R}}{\left({\mathit{H}}_0^{\mathrm{R}}\right)}^T={\left({\mathit{H}}_0^{\mathrm{R}}\right)}^2$ , which result in

$$\mathit{U}=\mathit{V}.$$

(10)

The zero matrix ${\mathit{K}}_0$ should finally be composed of linearly independent, transposed column vectors of $\mathit{U}$ in order to reduce the dimension. Mathematically, the kernel of $\mathit{U}$ is considered:

$${\mathit{K}}_0=\ker \left(\mathbf{U}\right).$$

(11)

The rank of $\mathit{U}$ therefore determines the dimension of ${\mathit{K}}_0$ . For reduced computing load, the number of considered singular values can be set manually as a percentage of the highest singular value $s_0$ . For the analysis with the SSFD, the percentage value has been set to p = 0.1%.

Numerically, the SVD algorithm is implemented on the MCU computing Jacobian rotations $\mathit{J}(i,j,\theta )$ . ²⁹ Jacobian rotation matrixes are orthogonal matrixes that contain a two-dimensional orthogonal rotation matrix as a submatrix in an identity matrix. The position of this submatrix is characterized by the indices $i,j\in \mathbb{N},i\ne j$ along the main diagonal. $\theta$ denotes the angle of rotation. The aim is to rotate ${\mathit{H}}_0^{\mathrm{R}}$ until the singular values no longer change significantly compared to the previous iteration step. ${\mathit{H}}_0^{\mathrm{R}}$ is transformed into a diagonal form using the so-called “orthogonal, two-sided Jacobian transformation for diagonalization.” This can be solved for small rotation angles with little effort. For matrixes that are not necessarily symmetrical, Hestenes ³⁰ was able to show that the one-sided transformation for orthogonalizing two column vectors $u_i$ and $u_j$ results in the column vectors of $\mathit{U}$ , whereby the column vectors after each iteration step depend on the previous ones. For symmetric matrixes, $\mathit{J}(i,j,\theta )=\mathit{J}(j,i,\theta )$ also applies so that $\left(\begin{array}{c}N\mathrm{\Delta }t\hfill \\ \multicolumn{1}{c}{2}\end{array}\right)$ Jacobian rotations are calculated for $N_{\mathit{J}}$ rows and columns in one iteration step according to combinatorics.

It was found that the calculation of ${\mathit{K}}_0$ takes about 2 min. Therefore, it is advantageous to perform the calculations in the training phase after the desired number of reference data sets have been recorded.

In the fourth step, the residuals are calculated, which are the column vectors of the matrix multiplication of ${\mathit{K}}_0$ and the Hankel matrix of a damaged structure ${\mathit{H}}_0^{\mathrm{D}}$ . In general, this matrix multiplication is not zero:

$$\mathit{\xi }={\mathit{K}}_0·{\mathit{H}}_0^{\mathrm{D}}\ne 0.$$

(12)

In the last step, the root mean square of $\mathit{\xi }$ is calculated over all elements of the number $N_{\mathit{\xi }}$ to obtain the DI:

$$\mathrm{D}{\mathrm{I}}_{\mathrm{SSFD}}=\sqrt{\frac{\sum _{i,j}{\xi }_{i,j}^2}{N_{\xi }}}.$$

(13)

Histogram difference

For the HD, the frequency distributions of the reference a ^R and damaged structure a ^D are subtracted from each other. By definition, one bin in a histogram represents the number of acceleration data points within a jointed value range. The condition that the bins are identical must be fulfilled. In order not to distort the statistics of the acceleration signals, no normalization is applied here.

The first step is to find the global maximum and minimum from the acceleration data. In the second step, the bin width can be calculated after

$$\mathrm{\Delta }s=\frac{max({\mathit{s}}^{\mathrm{D}},s^{\mathrm{R}})-min({\mathit{s}}^{\mathrm{D}},{\mathit{s}}^{\mathrm{R}})}{N_{\mathrm{B}}},$$

(14)

where N_B is the specified number of bins. In the third step, the N_B + 1 limits of the bins can be calculated:

$$b_i=min({\mathit{s}}^{\mathrm{D}},{\mathit{s}}^{\mathrm{R}})+\mathrm{\Delta }s·(i-1).$$

(15)

Here, N_B has been set to $50$ . The bins are now defined so that the acceleration data points can be sorted within the corresponding bin limits and the counts can be incremented. In the last step, the DI is considered as the area difference between the frequency distributions of the reference and the damaged structure:

$$\mathrm{D}{\mathrm{I}}_{\mathrm{HD}}=\sum _{i=0}^{N_{\mathrm{R}}-1}(b_{i+1}-b_i)(a_i^{\mathrm{D}}-\frac{1}{N_{\mathrm{R}}}\sum _{j=0}^{N_{\mathrm{R}}-1}{a}_{i,j}^{\mathrm{R}}).$$

(16)

If several references are considered, their frequency distributions are averaged bin-by-bin. The factor $(b_{i+1,j}-b_{i,j})$ for the area calculation is redundant because the bar lengths are all the same due to the equidistant subdivision

$$\mathrm{\Delta }s=b_{i+1}-b_i.$$

(17)

Combination of damage indicators

DIs are calculated using two algorithms based on the recorded acceleration data along three spatial directions. All three axes were considered to avoid any assumption related to the impact of a damage for a specific axis. In total, six DIs have been calculated. Combining them into one DI helps to detect changes better in the signal and thus to improve the robustness of the SHM system. First of all, the DIs have to be normalized to 1 for each spatial component and algorithm. After that, the spatial components are combined to a norm as follows:

$$\mathrm{D}{\mathrm{I}}_{3\mathrm{D},\mathrm{SSFD}}=\sqrt{\sum _{i\in \{x,y,z\}}{\mathrm{DI}}_{i,\mathrm{SSFD}}^2},$$

(18)

$$\mathrm{D}{\mathrm{I}}_{3\mathrm{D},\mathrm{HD}}=\sqrt{\sum _{i\in \{x,y,z\}}{\mathrm{DI}}_{i,\mathrm{HD}}^2}.$$

(19)

After another normalization to 1, the algorithms are combined to a norm as well:

$$\mathrm{D}{\mathrm{I}}_{\mathrm{res}}=\sqrt{\sum _{i\in \{\mathrm{SSFD},\mathrm{HD}\}}{\mathrm{DI}}_{3\mathrm{D},i}^2}.$$

(20)

Normalization is necessary as long as the scaling and ranges of the individual DIs differ from each other due to different calculation operations.

Classification method

Classification methods are essential to clearly identify the presence of damage and to compensate environmental parameters that influence the vibration behavior. The focus is on the motor control signals, which can change spontaneously due to wind and balancing maneuvers of the drone. Experience has shown that the outdoor temperature change by ±3°C within a few hours, so they were not included in the evaluation.

Since the control technology of the drone automatically ensures that the planned flight movement is maintained, it is quite possible that the servomotors may rotate at different speeds. This behavior is described by the displacement of an introduced centroid from its original position in the two-dimensional plane of the eight servomotors. This is illustrated in the schematic diagram in Figure 5. By definition, the centroids should point from the origin in the direction which has to be balanced. The adjacent servomotors always rotate in the opposite direction to ensure symmetry and to enable rotation around the vertical axis.

Figure 5.

Schematic representation of the drone in two dimensions. The coordinate system was selected so that the y-axis points to the forward movement of the drone. The blue arrows indicate the direction of rotation of the servomotors.

To determine the centroids of a measurement from $N_{\mathit{n}}$ motor control signals $\mathit{n}$ , the axial components $x$ and $y$ , and the diagonal components $d_1$ and $d_2$ have to be computed by the following equations in the first step:

$$x^{\prime }=\frac{1}{N_{\mathit{n}}}\sum _{i=0}^{N_{\mathit{n}}-1}[(n_{i,2}-n_{i,4})+(n_{i,6}-n_{i,8})+(n_{i,7}-n_{i,3})+(n_{i,5}-n_{i,1})],$$

(21)

$$y^{\prime }=\frac{1}{N_{\mathit{n}}}\sum _{i=0}^{N_{\mathit{n}}-1}[(n_{i,6}-n_{i,7})+(n_{i,2}-n_{i,5})+(n_{i,4}-n_{i,1})+(n_{i,8}-n_{i,3})],$$

(22)

$$d_1=\frac{1}{N_{\mathit{n}}}\sum _{i=0}^{N_{\mathit{n}}-1}[(n_{i,7}-n_{i,5})+(n_{i,6}-n_{i,1})+(n_{i,2}-n_{i,3})+(n_{i,4}-n_{i,8})],$$

(23)

$$d_2=\frac{1}{N_{\mathit{n}}}\sum _{i=0}^{N_{\mathit{n}}-1}[(n_{i,6}-n_{i,2})+(n_{i,7}-n_{i,4})+(n_{i,5}-n_{i,8})+(n_{i,1}-n_{i,3})].$$

(24)

The theorem of Pythagoras is used to project the diagonal components along the axial components:

$$d_1=\sqrt{{{\tilde{x}}^{\prime }}^2+{{\tilde{y}}^{\prime }}^2}\Rightarrow {\tilde{x}}^{\prime }=+{\tilde{y}}^{\prime }=\frac{d_1}{\sqrt{2}},$$

(25)

$$d_2=\sqrt{{{\tilde{x}}^{\prime }}^2+{(-{\tilde{y}}^{\prime })}^2}\Rightarrow {\tilde{x}}^{\prime }=-{\tilde{y}}^{\prime }=\frac{d_2}{\sqrt{2}}.$$

(26)

In the second step, centroids can be calculated taking into account all components:

$$\mathit{c}=\left(\begin{array}{c}{c}_{x^{\prime }}\hfill \\ \multicolumn{1}{c}{c_{y^{\prime }}}\end{array}\right)=\left(\begin{array}{c}{x}^{\prime }+{\tilde{x}}^{\prime }\hfill \\ \multicolumn{1}{c}{y^{\prime }+{\tilde{y}}^{\prime }}\end{array}\right)$$

(27)

with the norm

$$c=\sqrt{c_{x^{\prime }}^2+c_{y^{\prime }}^2}.$$

(28)

In the final step, for optimal reference selection, the minimal Euclidean distance from the centroid of the damaged structure to the centroid of N_R considered references is calculated by

$$\begin{array}{c}\multicolumn{1}{c}{d_{\min }=\min (\sqrt{{(c_{x^{\prime }}^{\mathrm{D}}-\frac{1}{N_{\mathrm{R}}}\sum _{i=0}^{N_{\mathrm{R}}}{c}_{x^{\prime },i}^{\mathrm{R}})}^2}}\\ \multicolumn{1}{c}{\overline{+{(c_{y^{\prime }}^{\mathrm{D}}-\frac{1}{N_{\mathrm{R}}}\sum _{i=0}^{N_{\mathrm{R}}}{c}_{y^{\prime },i}^{\mathrm{R}})}^2}).}\end{array}$$

The smaller the distances, the more similar the environmental and operational conditions are between the reference and a damaged state.

These three steps are performed at the beginning of the monitoring phase (see Figure 4) to select the optimal reference(s) with minimal distances of centroids.

Experimental results: onboard propeller fault diagnosis during flight in real time

This section is divided into the experimental results of the measurements during linear flights in summer and during hovering flights in winter. The results of the linear flight are discussed more in detail, as this case is practically more relevant for delivery drones. All linear flights were carried out in a straight line at a constant speed of 4 ms⁻¹ in the forward direction of the delivery drone.

In addition to the DI, raw data and centroids are also plotted. If DIs of a reference are calculated with themselves, a DI of $0$ is expected for the HD due to the subtraction of identical histograms, while a DI greater than $0$ can be assumed for the SSFD due to numerical uncertainties in the calculation of ${\mathit{K}}_0$ . As regards to the SHM parameters $p$ , N_cc, and N_B, it should be noted that their determination is based on a compromise between accuracy and calculation time to operate in real time.

Linear flight

First, an exemplary measurement dataset of a linear flight is shown, for which the plotted raw data are discussed on plausibility. Figure 6 demonstrates the acceleration data of the x-axis, the corresponding histogram and the motor control signals of servomotor 4 for three different structural states. It can be seen that a rotor blade with $5.2$  cm shortening generates higher vibrations on the drone system than a rotor blade with $2.6$  cm shortening or an intact rotor blade. In the histogram, the strong vibrations are reflected in a broader distribution. With a higher severity, it can also be seen that the servomotors rotate faster and the signals are more jagged, which can be explained by the missing value of buoyancy due to the shorter rotor blade.

Figure 6.

(a) Measured raw acceleration data along x-axis, (b) corresponding histogram and (c) retrieved motor control signals of servomotor 4 for three structural states.

The centroids for all linear flights are plotted in Figure 7. The white arrow in the y-direction marks the forward movement of the delivery drone. As the rear servomotors have to rotate faster for this, the centroids can be found at negative values of $c_{y^{\prime }}$ . A point cloud can be seen in which partial overlaps of the centroids of different structural states exist. The overlaps with reference centroids are particularly interesting, as this ensures that an optimum reference can be found. The measurement series with a rotor blade shortened by $5.2$  cm has the largest scattering behavior of the centroids.

Figure 7.

Centroids for the linear flights, calculated from the appropriate motor control signals.

To display the normalized DI separately for each spatial component and each algorithm, three-dimensional plots were created separately for the SSFD and HD in Figure 8. Only the optimal reference was considered for this purpose. The HD enables a more precise separation of clusters of different structural states.

Figure 8.

Representation of normalized DI in three-dimensional space, calculated by the (a) SSFD and (b) HD method for the linear flights.

The number of the most suitable references was varied for the combined DI. Figure 9 demonstrates the normalized DI with 1, 10 and all ( $48$ ) considered references. The same graph also plots the minimum centroid distances, which increase with higher severity and the number of considered references. The most expected trend is shown in Figure 9(a), where only the optimal reference was considered. Nevertheless, all three cases show that the DI increases as expected with higher severity of the damage.

Figure 9.

Trend of combined, normalized DI and mean of minimal distances to (a) 1, (b) 10, and (c) 48 reference centroids during the linear flights.

Hovering flight

The centroids for the hovering flights are shown in Figure 10. Once again, the centroids of the structural state with the highest severity are particularly remarkable. This makes it significantly more difficult to find optimal references. Another aspect that can be observed during the comparison is the narrower range of the axes, although the wind was more present during the measurements in hovering flight. This shows that the control of the servomotors for a linear flight requires significantly larger rotation speeds than the wind to the thrust of the delivery drone.

Figure 10.

Centroids for the hovering flights, calculated from the appropriate motor control signals.

The combined, normalized DI were calculated with 1, 10, and all ( $60$ ) considered references and plotted in Figure 11. For the first case, the trend corresponds to the expectation, whereas for the other two cases, the reference and rotor blade shortening of $2.3$  cm are barely distinguishable from each other. Looking at Figure 10, the centroids also barely overlap that makes it difficult to find an optimal reference. This can also be seen in the minimum centroid distances. The reason for this lies in the strongly varying wind conditions that dominated on the whole day. It is suggested to consider only one reference that is closest to the reference measurement ( $N_R=1$ ).

Figure 11.

Trend of combined, normalized DI and mean of minimal distances to (a) 1, (b) 10, and (c) 60 reference centroids during the hovering flights.

Conclusions

In this article, it was shown that an onboard evaluation of the airworthiness of the Q8000 delivery drone can be carried out on an MCU in real-time. The benchmark is a DI combined from the SSFD and HD to classify the severity of the cut rotor blade up to $5.2$  cm. To compensate different environmental conditions which are primary forced by wind, an algorithm was proposed that calculates centroids from motor control signals of eight servomotors. It could be figured out that taking into account only the optimum reference gave the best results for linear and hovering flights (see Figures 9(a) and 11(a)).

The developed SHM system can be used for any unmanned aerial vehicle (UAV) as long as up-to-date and suitable reference datasets have been recorded beforehand. For drone development, it can be embedded to initiate an automatic emergency landing in the event of damage. Extending the SHM system with machine learning, local methods and more acceleration sensors can increase detection accuracy, although it is important to ensure that the firmware does not become too complex on a MCU given by the limited computational ressources. An interesting approach for further research could also be the reuse of reference datasets. This would allow the SHM system to be better trained for different environmental parameters and the short battery life would no longer be a limiting factor.