Real-Time Terrain-Following of an Autonomous Quadrotor by Multi-Sensor Fusion and Control

: For the application of the autonomous guidance of a quadrotor from conﬁned undulant ground, terrain-following is the major issue for ﬂying at a low altitude. This study has modiﬁed the open-source autopilot based on the integration of a multi-sensor receiver (a Global Navigation Satellite System (GNSS)), a Lidar-lite (a laser-range-ﬁnder device), a barometer and a low-cost inertial navigation system (INS)). These automatically control the position, attitude and height (a constant clearance above the ground) to allow terrain-following and avoid obstacles based on multi-sensors that maintain a constant height above ﬂat ground or with obstacles. The INS/Lidar-lite integration is applied for the attitude and the height stabilization, respectively. The height control is made by the combination of an extended Kalman ﬁlter (EKF) estimator and a cascade proportional-integral -derivative (PID) controller that is designed appropriately for the noise characteristics of low accuracy sensors. The proposed terrain-following is tested by both simulations and real-world experiments. The results indicate that the quadrotor can continuously navigate and avoid obstacles at a real-time response of reliable height control with the adjustment time of the cascade PID controller improving over 50% than that of the PID controller.


Introduction
In recent years, advances in constructions, electronics, sensors and control algorithms have fueled a growth in the development of unmanned aerial vehicles (UAVs) [1]. The automatic applications of UAVs rely on the continuous control of the position, attitude, speed and other high level commands [2,3]. Specifically, the attitude and height (the clearance from the ground) controllers allow terrain-following of an autonomous quadrotor, which is useful for topography prediction, navigation and obstacle avoidance [4][5][6]. An integrated navigation framework is an autonomous system used for determining the position, attitude and velocity of UAVs using multi-sensors [7,8], i.e., a global positioning system (GPS), an inertial measurement unit ((IMU) e.g., an accelerometer, a gyroscope and a magnetometer) [9] and radar or distance sensors [10][11][12][13]. To achieve a more reliable travel in real environments, fruitful research efforts have focused on the height control problem [14].
The height control of common UAVs applies a proportional-integral-derivative (PID) control structure for the altitude that handles high impact and changed mass [15][16][17]. These controllers compare the behavior of the UAVs with a model and adapt the output of the controller to minimize the error between the model and the system. The design, implementation, simulation and evaluation of this controller will be presented in this paper divided into algorithms based on a classic control theory and a modern control theory. Control algorithms based on a modern control theory include a dynamic inverse control [18], a linear-quadratic regulator (LQR) control [19], an intelligent control (such as a neural network control [20,21]) and an adaptive control [22,23]. Ngo used the dynamic inversion based adaptive control method to control the height of the UAVs in the climbing stage and verified that the method had good effects through simulation [24]. However, the long-term flight of UAVs is a time-varying nonlinear system [25,26] whose external environment, speed, load and other parameters change continuously. For dynamic inverse control strategies, the control effect is often not satisfactory because there is no effective method to obtain the precise parameters of the model. Due to the requirement of the rapidity of response, the intelligent control algorithm, which needs complex calculations, is less used in practical engineering applications though it has a strong environmental adaptability [16]. Thus, in engineering practice, methods based on a classic control (such as the PID control) are mostly used because of the uncertainty of system parameters and the requirement of real-time control. The PID control is an effective and reliable flight control strategy that has been studied for some time and has a mature design strategy. For the UAVs' terrain-following system, the problem to be solved is that the traditional single-stage PID control cannot meet the requirements in terms of response speed and stability [27]. In order to improve the quality of control, Idres used the cascade PID control method to control the UAVs and verified the efficiency of the cascade PID controlled UAVs in tracking slopes, circles and sinusoidal trajectories through a simulation [28]. Ren analyzed the system based on the classical control theory and proved that the cascade PID was better than the traditional single-stage PID control from an amplitude margin and a phase angle margin [29]. A cascade PID is one of the effective methods used to improve the quality of traditional PID control and it has been widely used in the field of process control [30][31][32]. For the cascade PID control of terrain-following UAVs, how to set the control feedback state quantity and how to adjust its parameters are important issues. The Ziegler-Nichols method [33], the critical proportionality method [34] and the attenuation curve method [35] are commonly used in engineering for tuning. Choosing a suitable method for PID tuning can also help improve the quality of control [36,37].
The main goal of this study is to improve the autonomous terrain-following and navigation performance of a quadrotor system. We address the detailed mathematical modeling of the quadrotor's kinematics and dynamics of the quadrotor. The autonomous controller comprises of an integrated control for the quadrotor's position, attitude and height control with the multi-sensor based control loops for terrain-following depicted in Figure 1. The control methodologies are conducted in MATLAB Simulink and the experiments on a quadrotor test-bed where interference such as wind conditions and sensor noise are incorporated in both the simulation and real-world experiments. A quadcopter measures its altitude relatively to the ground by using a range-finder. However, when the surface of the ground is fluctuating, a traditional PID control algorithm cannot meet the requirements of accurate terrain-following. Work in this innovation is the adaptive control function using the integration of a Global Navigation Satellite System (GNSS) (GPS) and a BeiDou Navigation Satellite System (BDS) [38], an inertial navigation system (INS) [39], a barometer [40] and a Lidar-lite sensor [41], especially the combination of an extended Kalman filter (EKF) and the cascade PID controller. In this paper, a dynamic model of the quadcopter is built and then linearized for simplification. The cascade PID is modeled and simulated in Simulink based on a tuning method according to the transfer function. The rest of this paper is organized as follows. In II, dynamic models of the quadrotor are described. In III, we design the integration of the position and height control based on the combination of an extended Kalman filter and a cascade PID. IV verifies the reliability of the proposed models by MATLAB Simulink and real-world quadrotor flying avoiding a box above the ground. Conclusions are given in V.

Quadrotor Dynamic Models
The quadrotor changes its position and attitude by modifying the speed of the four motors [42,43]. The movements of back and forth and left and right are produced by pitch and roll, respectively, and the movements of up and down are achieved by the four motors increasing or decreasing the speed at the same time.
In order to study the kinematic law of aircraft, two coordinate systems [44] depicted in Figure 2 are defined as follows: 1. The geographic coordinate system ( E E E E O X Y Z )is used to study the motion state of aircraft relative to the earth and the space position coordinates in which is used to study the rotational motion of aircraft relative to the center of gravity in which The rest of this paper is organized as follows. In II, dynamic models of the quadrotor are described. In III, we design the integration of the position and height control based on the combination of an extended Kalman filter and a cascade PID. IV verifies the reliability of the proposed models by MATLAB Simulink and real-world quadrotor flying avoiding a box above the ground. Conclusions are given in V.

Quadrotor Dynamic Models
The quadrotor changes its position and attitude by modifying the speed of the four motors [42,43]. The movements of back and forth and left and right are produced by pitch and roll, respectively, and the movements of up and down are achieved by the four motors increasing or decreasing the speed at the same time.
In order to study the kinematic law of aircraft, two coordinate systems [44] depicted in Figure 2 are defined as follows: 1.
The geographic coordinate system (O E X E Y E Z E ) is used to study the motion state of aircraft relative to the earth and the space position coordinates in which Y E -X E -Z E point to East-North-Up, respectively.

2.
The is used to study the rotational motion of aircraft relative to the center of gravity in which According to the conversion relationship between the Euler angle and the attitude matrix, the attitude transfer matrix R b E from the geographic coordinate system to the airframe coordinate system can be obtained: cos θ sin ψ sin ψ cos θ − sin θ sin φ sin θ cos ψ − cos φ sin ψ sin φ sin θ sin ψ + cos φ cos ψ sin φ cos θ cos φ sin θ cos ψ + sin φ sin ψ cos φ sin ψ sin θ − sin φ cos ψ cos φ cos θ   (1) The vector of x y z φ θ ψ T represents the position and angle of the quadrotor in the geographic coordinate system and the vector of u v w p q r T represents the linear velocity and angular velocity of the quadrotor in the airframe coordinate system. The conversion relationship of these twelve-dimensional variables in the two coordinate systems is: In Equation (2), v E = . x is the direction cosine matrix from the airframe coordinate system to the geographic coordinate system and T is the conversion matrix between the airframe rotation angular velocity and the attitude angle change rate. The dynamic model of the quadrotor is: According to Newton's law of motion, the total force acting on the quadrotor is: where m is the mass of the quadcopter, ⊗ represents the cross product operation and f b = f x f y f z ∈ R 3 is the total force acting on the quadrotor. According to Euler's equation, the total torque produced by the four motors is: In (5), τ b = τ x τ y τ z ∈ R 3 represents the total torque, I is the inertia matrix of the quadrotor, which can be simplified as follows according to the symmetry of the quadrotor: where I xx , I yy , I zz are called the central principal moments of inertia and can be measured through experiments. Therefore, the dynamic model of the quadrotor in the airframe coordinate system is: pI xx − qrI yy + qrI zz τ y = . qI yy + prI xx − prI zz τ z = . rI zz − pqI xx + pqI yy .

of 19
The above equation is established under the condition that the origin and axis of the airframe coordinate system coincide with the centroid and main axis of the quadrotor.
The total external force received by the quadrotor in the airframe coordinate system is: whereê z is the unit vector of the z-axis of the geographic coordinate system,ê 3 is the unit vector of the z-axis of the airframe coordinate system, g is the acceleration of gravity, f t is the total pulling force generated by motors and f w is the force generated by the wind acting on the quadrotor. Ignoring the rotational torque in the control equation, the total external torque T b in the airframe coordinate system can be expressed as: where τ b = τ x τ y τ z T ∈ R 3 , τ w = τ wx τ wy τ wz T ∈ R 3 are the control torques generated by the difference in motor speed and the movement generated by the wind acting on the quadrotor. Substituting Equations (8) and (9) into Equation (7), the complete dynamic model of the quadrotor in the airframe coordinate system is obtained as: rI zz − pqI xx + pqI yy (10) Considering that the force and torque input by the control system are proportional to the square of the motor speed, the parameters in (10) are as follows: where l is the distance between the center of the motor and the center of the airframe, c T is the tensile coefficient of the motor, and c M is the torque coefficient. All of the above parameters can be measured through experiments. Substituting (11) into (10), the complete dynamic model of the quadrotor in the airframe coordinate system can be further expressed as: The parameter setting of the quadrotor models is listed in Table 1. The state vector is expressed as follows: According to Equations (3) and (10), the dynamic equation of the quadrotor can be rewritten as: In hovering mode, the attitude of the aircraft changes very little. It can be made that cos α = 1, sin α= tanα = 0 and the state Equation (14) can be approximately expressed as: The above equation can be written in the form of a differential equation: In Equation (16), The balance point of the system can be set as follows: When the system is in equilibrium, the input of the system is: The state equation is further written as: Ignoring the influence of wind, it can be written as: That is, the linearized model of the hovering state equation can be written as: where A and B can be expressed as follows:

Quadrotor State Estimations
When the laser-range-finder is used for distance measurement in terrain-following flight, the large vertical drop of the terrain will cause the delay of altitude control and the decrease of the stability of the quadrotor. Therefore, it is necessary to ameliorate the altitude control channel of the quadrotor to improve the stability and real-time performance of the quadrotor in a large maneuvering and large overload terrain-following flight.
Combined with the linearized dynamics model of the quadrotor in the hovering state in Section 2.2, the state equation of the height channel control model can be written as: where A, B, x and u are the same as those defined in section A, C = 0 0 0 0 0 0 0 0 1 0 0 0 . Using a Laplace transform on both sides of Equation (24), the transfer function of the vertical speed control can be obtained as: According to the voltage balance equation of the motor armature circuit, the motor model can be approximated as a first-order inertia link and the transfer function of the motor torque is:  (26), the final transfer function of vertical velocity is: . (27) Adding the integral link to Equation (27) to obtain the transfer function of the UAVs' height control channel gives: Using the parameter measurement method proposed in the literature [45], the inertial time constant T m = 0.153 and the open-loop gain K = 1542 of the 2212KV910 brushless motor were measured by experiment. Substituting T m and K K into Equations (27) and (28), Equation (29) can be obtained:

Height and Vertical Velocity Estimations
As the error of the accelerometer will accumulate with the integration, it is not reliable to obtain the vertical velocity by using the accelerometer integration alone. Therefore, the optimal vertical velocity is estimated by integrating the differential data of the range-finder and the integral data of the accelerometer. The sampling period is set as T = 1 and the block diagram of the discretized fusion speed measurement system is shown in Figure 3.  (27) and (28), (29) can be obtained:

Height and Vertical Velocity Estimations
As the error of the accelerometer will accumulate with the integration, it is not reliable to obtain the vertical velocity by using the accelerometer integration alone. Therefore, the optimal vertical velocity is estimated by integrating the differential data of the rangefinder and the integral data of the accelerometer. The sampling period is set as T = 1 and the block diagram of the discretized fusion speed measurement system is shown in Figure  3. In Figure 3, k a is the speed increment from time is the optimal estimation of vertical velocity. k Z is the observations of multi-sensor increment from time to time k t , which are used as the estimated value of the observation correction of the speed measurement system. Assuming that both the system and the measurement noise are Gaussian white noise with an expectation of zero, the discretized extended Kalman filter equation of the velocity measurement system [46,47] is: where k X is the estimated value before correction at k t , 1 k X  is the best estimate at 1 k t  , k P is the mean square error before correction, k K is the Kalman gain, k H is the observation matrix, 1 k Q  and k R are the noise variance of the system and the measurement process, respectively [48].  In Figure 3, a k is the speed increment from time t k−1 to time t k ,v(k − 1) is the optimal estimation of vertical velocity. Z k is the observations of multi-sensor increment from time t k−1 to time t k , which are used as the estimated value of the observation correction of the speed measurement system.
Assuming that both the system and the measurement noise are Gaussian white noise with an expectation of zero, the discretized extended Kalman filter equation of the velocity measurement system [46,47] is: where X k is the estimated value before correction at t k ,X k−1 is the best estimate at t k−1 , P k is the mean square error before correction, K k is the Kalman gain, H k is the observation matrix, Q k−1 and R k are the noise variance of the system and the measurement process, respectively [48].

Height and Vertical Velocity Based on a PID Controller
The multi-sensor control method is a classic control method mainly designed by the PID control principle. The differential equation of the PID controller is as follows: 1.
Continuous form where u(t) is the output of the control system, k p , k i and k d is the proportional, integral and the differential gain. And e is the error between the set point and the feedback value, t represents the current time, and τ is an integral variable of which the value is from 0 to the current time.

2.
Discrete form The PID control system of the quadrotor is based on a hierarchical control of inner and outer loops [49,50] and it can be further subdivided into four levels including position control, attitude control, control distribution and power distribution as shown in Figure 4.

Height and Vertical Velocity Based on a PID Controller
The multi-sensor control method is a classic control method mainly designed by the PID control principle. The differential equation of the PID controller is as follows: where ( ) u t is the output of the control system, p k , i k and d k is the proportional, integral and the differential gain. And e is the error between the set point and the feedback value, t represents the current time, and  is an integral variable of which the value is from 0 to the current time.

Discrete form
The PID control system of the quadrotor is based on a hierarchical control of inner and outer loops [49,50] and it can be further subdivided into four levels including position control, attitude control, control distribution and power distribution as shown in Figure  4. As the height control of the quadrotor has nothing to do with the attitude of the body, there is an independent feedback input z in the height control system and the height correction is output after the height PID controller. Due to the balance between gravity and lift when the quadrotor is hovering, it is necessary to add gravity compensation to the height controller. When there is no error in the height, the controller can still output a stable value to control the motor rotation to maintain the hovering state. The simulation structure of the height controller is shown in Figure 5. As the height control of the quadrotor has nothing to do with the attitude of the body, there is an independent feedback input z in the height control system and the height correction is output after the height PID controller. Due to the balance between gravity and lift when the quadrotor is hovering, it is necessary to add gravity compensation to the height controller. When there is no error in the height, the controller can still output a stable value to control the motor rotation to maintain the hovering state. The simulation structure of the height controller is shown in Figure 5

Height and Vertical Velocity Based on a Cascade PID Controller
The cascade control system is one of the effective methods to improve the control quality and it has been widely used in process controls. The cascade PID controller is composed of a single loop PID regulator (as the main regulator) and an external set regulator (as a secondary regulator) connected in series to form a dual loop regulation system. The control output of the main regulator is used as the external set regulator and the control output of the external set regulator is sent to the control regulation structure.
In order to improve the quality of height control effectively, a cascade PID based height controller system is designed as follows. The error between the expected height and the height measured by the range-finder is set as the outer loop feedback of the cascade PID controller, the expected vertical speed is obtained after the outer loop PID controller and then the difference between the expected vertical velocity and the actual vertical speed is set as the inner loop feedback input to accomplish the adjustment of the vertical speed, which finally realizes the dual loop control of the height and vertical speed. The control structure and Simulink are shown in Figures 6 and 7.

Height and Vertical Velocity Based on a Cascade PID Controller
The cascade control system is one of the effective methods to improve the control quality and it has been widely used in process controls. The cascade PID controller is composed of a single loop PID regulator (as the main regulator) and an external set regulator (as a secondary regulator) connected in series to form a dual loop regulation system. The control output of the main regulator is used as the external set regulator and the control output of the external set regulator is sent to the control regulation structure.
In order to improve the quality of height control effectively, a cascade PID based height controller system is designed as follows. The error between the expected height and the height measured by the range-finder is set as the outer loop feedback of the cascade PID controller, the expected vertical speed is obtained after the outer loop PID controller and then the difference between the expected vertical velocity and the actual vertical speed is set as the inner loop feedback input to accomplish the adjustment of the vertical speed, which finally realizes the dual loop control of the height and vertical speed. The control structure and Simulink are shown in Figures 6 and 7.

Simulink Parameter Tuning by Attenuation Curve
In order to obtain the optimal control effect, the parameters of the PID regulator need to be tuned. There are three common methods for PID parameter tuning in engineering: the Ziegler-Nichols tuning method, the critical proportionality method and the attenuation curve method. These three methods have their own characteristics: the Ziegler-Nichols tuning method is suitable for systems with a delay link in their open loop transfer function while the critical proportionality method requires the order of the system to be third-order or above and the attenuation curve method adjusts controller parameters by attenuation frequency characteristics [28]. If we suppose that the initial state of the quadrotor is

Simulink Parameter Tuning by Attenuation Curve
In order to obtain the optimal control effect, the parameters of the PID regulator need to be tuned. There are three common methods for PID parameter tuning in engineering: the Ziegler-Nichols tuning method, the critical proportionality method and the attenuation curve method. These three methods have their own characteristics: the Ziegler-Nichols tuning method is suitable for systems with a delay link in their open loop transfer function while the critical proportionality method requires the order of the system to be third-order or above and the attenuation curve method adjusts controller parameters by attenuation frequency characteristics [28]. If we suppose that the initial state of the quadrotor is

Simulink Parameter Tuning by Attenuation Curve
In order to obtain the optimal control effect, the parameters of the PID regulator need to be tuned. There are three common methods for PID parameter tuning in engineering: the Ziegler-Nichols tuning method, the critical proportionality method and the attenuation curve method. These three methods have their own characteristics: the Ziegler-Nichols tuning method is suitable for systems with a delay link in their open loop transfer function while the critical proportionality method requires the order of the system to be third-order or above and the attenuation curve method adjusts controller parameters by attenuation frequency characteristics [28]. If we suppose that the initial state of the quadrotor is x = [ 0 0 0 0 0 0 0 0 0 0 0 1 0] T and Ω i = 4100 4100 4100 4100 T , this means that the drone hovers at a height of 10 m and the simulation duration is set to 15 s. After 2 s, a position control command is given to make the aircraft rise to a height of 20 m and remain hovering. The simulation results are shown in Figure 8. make the aircraft rise to a height of 20 m and remain hovering. The simulation results are shown in Figure 8.  (31) of the height control system, the attenuation curve method is selected to tune the PID controller parameters and build the single-stage/cascade PID controller simulation block in Figure  8.
After setting up the control system structure in Simulink, first set the regulator parameters in the control system to a pure proportional action ( = ∞, = 0) put the system into operation and then gradually adjust the proportionality from large to small until the ratio of the amplitude of the two adjacent periodic curves is 4:1 attenuation.
At this time, the ratio is 4:1, the attenuation ratio is , the rise time is and the time interval between two adjacent peaks is . According to , or , use the empirical formula in Table 2 to calculate each parameter value of the regulator [28]. Table 2. PID regulator parameter tuning based on the attenuation curve method.

Controller
Proportional band The sequence of the PID controller parameter tuning is a proportional link first, then an integral link and then a derivative link. The parameters after tuning are: In addition to the former requirements, the cascade PID controller must follow the inner first and then the outer setting rule. The parameters after tuning are: Inner loop: Outer loop:  (31) of the height control system, the attenuation curve method is selected to tune the PID controller parameters and build the single-stage/cascade PID controller simulation block in Figure 8.
After setting up the control system structure in Simulink, first set the regulator parameters in the control system to a pure proportional action (Ti = ∞, τ = 0) put the system into operation and then gradually adjust the proportionality δ from large to small until the ratio of the amplitude of the two adjacent periodic curves is 4:1 attenuation.
At this time, the ratio is 4:1, the attenuation ratio is δs, the rise time is tr and the time interval between two adjacent peaks is Ts. According to δs, tr or Ts, use the empirical formula in Table 2 to calculate each parameter value of the regulator [28]. Table 2. PID regulator parameter tuning based on the attenuation curve method.

Controller
Proportional Band δ

Integration Time Ti
Derivative Time τ The sequence of the PID controller parameter tuning is a proportional link first, then an integral link and then a derivative link. The parameters after tuning are: In addition to the former requirements, the cascade PID controller must follow the inner first and then the outer setting rule. The parameters after tuning are: Inner loop: Outer loop: Set the expected value of the height to 3 m and then input the tuned parameters into the simulation system. After the system has run, the response curves of the single-stage PID controller and the cascade PID controller are shown in Figure 9. 3.8235 5.1282 0.04875 Set the expected value of the height to 3 m and then input the tuned parameters into the simulation system. After the system has run, the response curves of the single-stage PID controller and the cascade PID controller are shown in Figure 9. It can be analyzed from Figure 9 that after the parameters are set by the attenuation curve method, the time it takes for the cascade PID controller to reach the expected height is reduced by 70% compared with the PID controller. The comparison of specific indicators is shown in Table 3 below.

Carton Box Avoidance Experiment
To simulate the drastic changes of relative height to the ground, we carried out a carton box avoidance experiment to validate the proposed control algorithm. The assembled quadrotor as the experimental platform used a core controller stm32f427 cortex M4 and output inertial data with an mpu9250 nine axis IMU. Lidar-lite V3, a high-precision laser-range-finder, was used in the ranging sensor. Its data update frequency can reach 400 Hz, which can meet the requirements of high dynamic ranging. The cascade PID control algorithm was transplanted to flight control. The laser-range-finder height was taken as the output of the outer loop and the change rate of the measured altitude with time was taken as the input of the inner loop. The control performance test was carried out in an open outdoor field in Southeast University, China. The wheelbase of our quadrotor was 400 mm and the parameters of its brushless motor are shown in Table 4. It can be analyzed from Figure 9 that after the parameters are set by the attenuation curve method, the time it takes for the cascade PID controller to reach the expected height is reduced by 70% compared with the PID controller. The comparison of specific indicators is shown in Table 3 below.

Carton Box Avoidance Experiment
To simulate the drastic changes of relative height to the ground, we carried out a carton box avoidance experiment to validate the proposed control algorithm. The assembled quadrotor as the experimental platform used a core controller stm32f427 cortex M4 and output inertial data with an mpu9250 nine axis IMU. Lidar-lite V3, a high-precision laserrange-finder, was used in the ranging sensor. Its data update frequency can reach 400 Hz, which can meet the requirements of high dynamic ranging. The cascade PID control algorithm was transplanted to flight control. The laser-range-finder height was taken as the output of the outer loop and the change rate of the measured altitude with time was taken as the input of the inner loop. The control performance test was carried out in an open outdoor field in Southeast University, China. The wheelbase of our quadrotor was 400 mm and the parameters of its brushless motor are shown in Table 4. In order to simulate the height mutation of ground obstacles, a 0.48 m high carton was placed on flat ground as shown in Figure 10. Considering the great influence of high winds before take-off, the expected height of the UAV was set as 1.9 m, the horizontal flight speed was set as 0.5 m/s and the UAV was manipulated to fly over the carton repeatedly for 10 times. The controlled variables and experimental scenarios in the experiment are shown in Table 5 and Figure 10, respectively. When the quadrotor was holding its altitude, the proposed controller should have detected the edge of the box and manipulated it to fly over the box.

0.153
In order to simulate the height mutation of ground obstacles, a 0.48 m high carton was placed on flat ground as shown in Figure 10. Considering the great influence of high winds before take-off, the expected height of the UAV was set as 1.9 m, the horizontal flight speed was set as 0.5 m/s and the UAV was manipulated to fly over the carton repeatedly for 10 times. The controlled variables and experimental scenarios in the experiment are shown in Table 5 and Figure 10, respectively. When the quadrotor was holding its altitude, the proposed controller should have detected the edge of the box and manipulated it to fly over the box. Figure 10. Terrain-following experiment of the quadrotor above the ground with the cardboard box as the obstacle or the ground undulation.

Variables
Values/m Height of the obstacle (the box) 0.48 Width of the obstacle (the box) 1.2 Travel speed of the quadrotor (m/s) 0.5 The resulting height data of the multi-sensor based controller can be seen in Figure  11 with the target height of 10 m above flat terrain. The difference in the performance of the height control clearly illustrates a better stability and sensitivity of the multi-sensor control. Therefore, the multi-sensor based control can react quickly especially with a reliable position and the velocity control in the z-axis.  The resulting height data of the multi-sensor based controller can be seen in Figure 11 with the target height of 10 m above flat terrain. The difference in the performance of the height control clearly illustrates a better stability and sensitivity of the multi-sensor control. Therefore, the multi-sensor based control can react quickly especially with a reliable position and the velocity control in the z-axis. For the height control of the quadrotor, we implemented with a PID and a cascade PID controller, respectively. The laser-range-finder measurements for the two controllers and the target height (1.9 m) are shown in Figure 12. For the height control of the quadrotor, we implemented with a PID and a cascade PID controller, respectively. The laser-range-finder measurements for the two controllers and the target height (1.9 m) are shown in Figure 12. For the height control of the quadrotor, we implemented with a PID and a cascade PID controller, respectively. The laser-range-finder measurements for the two controllers and the target height (1.9 m) are shown in Figure 12.  From the comparison of the fixed height flight results in Figure 12a,b, it can be seen that the adjustment time of the cascade PID controller was much less than that of the PID controller in the case of sudden changes of the ground altitude. In addition, the stability of the cascade PID control was better than that of the PID. After considering the statistics of all of the height measurement data, the shortest time from detecting a height change to reaching the desired height was 1.13 s and the longest time was 1.66 s when using the PID controller; when using the cascade PID controller, the shortest time was 0.43 s and the longest time was 0.75 s. The adjustment time of the cascade PID controller improved over 50% than that of the PID controller.

Conclusions
This study investigated the position, attitude, velocity and height estimators based on GPS, accelerometers, gyroscopes and Lidar-lite onboard sensors, which were operated on a quadrotor of the Pixhawk open-source autopilot. For the terrain-following of the quadrotor, we have presented real-time height control, navigation and obstacle avoidance above outdoor undulant ground. In our system, GNSS and INS information were used to verify the position, velocity and attitude estimates with a constant height control for autonomous flights. The navigation and estimations depended on an available multi-sensor From the comparison of the fixed height flight results in Figure 12a,b, it can be seen that the adjustment time of the cascade PID controller was much less than that of the PID controller in the case of sudden changes of the ground altitude. In addition, the stability of the cascade PID control was better than that of the PID. After considering the statistics of all of the height measurement data, the shortest time from detecting a height change to reaching the desired height was 1.13 s and the longest time was 1.66 s when using the PID controller; when using the cascade PID controller, the shortest time was 0.43 s and the longest time was 0.75 s. The adjustment time of the cascade PID controller improved over 50% than that of the PID controller.

Conclusions
This study investigated the position, attitude, velocity and height estimators based on GPS, accelerometers, gyroscopes and Lidar-lite onboard sensors, which were operated on a quadrotor of the Pixhawk open-source autopilot. For the terrain-following of the quadrotor, we have presented real-time height control, navigation and obstacle avoidance above outdoor undulant ground. In our system, GNSS and INS information were used to verify the position, velocity and attitude estimates with a constant height control for autonomous flights. The navigation and estimations depended on an available multi-sensor fusion. Apart from the absolute altitude estimation, the terrain-following required relative distance measurements from the Lidar-lite sensor. The proposed control algorithm was tested extensively in MATLAB Simulink simulations and then implemented on the real-world quadrotor vehicle. The results confirmed that the proposed combination of the multisensor fusion and cascade PID controller could ensure autonomous navigation, obstacle avoidance and height control, achieving better accuracy, response time and sensitivity. Future developments will fuse trajectory and velocity smoothing that enable the balance between control efficiency and terrain-following criteria.  Data Availability Statement: The following are available online at https://github.com/yangyuancsi/ Altitude-control: PID: code for the PID controllers, Altitude estimation: code for the Kalman filter, Video: video for carton box avoidance experiment. The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest:
The authors declare no conflict of interest.