extfieldsofstudy:
Medicine
Computer Science

text:
Radar and Visual Odometry Integrated System Aided Navigation for UAVS in GNSS Denied Environment

Drones are becoming increasingly significant for vast applications, such as firefighting, and rescue. While flying in challenging environments, reliable Global Navigation Satellite System (GNSS) measurements cannot be guaranteed all the time, and the Inertial Navigation System (INS) navigation solution will deteriorate dramatically. Although different aiding sensors, such as cameras, are proposed to reduce the effect of these drift errors, the positioning accuracy by using these techniques is still affected by some challenges, such as the lack of the observed features, inconsistent matches, illumination, and environmental conditions. This paper presents an integrated navigation system for Unmanned Aerial Vehicles (UAVs) in GNSS denied environments based on a Radar Odometry (RO) and an enhanced Visual Odometry (VO) to handle such challenges since the radar is immune against these issues. The estimated forward velocities of a vehicle from both the RO and the enhanced VO are fused with the Inertial Measurement Unit (IMU), barometer, and magnetometer measurements via an Extended Kalman Filter (EKF) to enhance the navigation accuracy during GNSS signal outages. The RO and VO are integrated into one integrated system to help overcome their limitations, since the RO measurements are affected while flying over non-flat terrain. Therefore, the integration of the VO is important in such scenarios. The experimental results demonstrate the proposed system’s ability to significantly enhance the 3D positioning accuracy during the GNSS signal outage.

Introduction
Over the past 10 years, Unmanned Aerial Vehicles (UAVs have become an essential tool for wide civil and military applications, such as first aid delivery, firefighting, rescue, downtown traffic monitoring, disaster management, border mentoring, and reconnaissance. The deployment of these UAVs in such aspects can help minimize human exposure to hazards and risks in dangerous circumstances. Most of the UAVs are commonly dependent on their onboard Global Navigation Satellite System (GNSS)/Inertial Navigation System (INS) integrated navigation systems to estimate where they are flying while performing their missions. The main challenge for this system configuration while flying in such cluttered environments is that the satellite signals might be disturbed in various scenarios, such as signals blockage, attenuation, multipath effects, jamming, spoofing in hostile areas, and flying in urban areas or natural canyons. In such harsh scenarios, the navigation solution is obtained only by the INS prior to the GNSS signal retrieval. Hence, the navigation solution will

Frequency Modulated Continuous Wave Radar Odometry
The proposed RO estimates the vehicle forward velocity and the height above the ground from its Range Doppler Map (RDM). This RO consists of data acquisition, target detection, velocity, and height above ground extraction.

Data Acquisition
This RO acquires the ranges and velocity measurements for the reflected ground scatterer based on an internal radar signal processing. During the flight, the attached micro radar is continuously transmitting a frequency modulated sawtooth chirps toward the ground objects : where is the sweep rate, is the frequency sweep time, and 0 is the initial transmitted frequency. This sweep rate can be obtained as follows:

Frequency Modulated Continuous Wave Radar Odometry
The proposed RO estimates the vehicle forward velocity and the height above the ground from its Range Doppler Map (RDM). This RO consists of data acquisition, target detection, velocity, and height above ground extraction.

Data Acquisition
This RO acquires the ranges and velocity measurements for the reflected ground scatterer based on an internal radar signal processing. During the flight, the attached micro radar is continuously transmitting a frequency modulated sawtooth chirps toward the ground objects f RF TX : where K f is the sweep rate, T is the frequency sweep time, and f 0 is the initial transmitted frequency. This sweep rate can be obtained as follows: where BW represents the transmitted chirp signal bandwidth. A roundtrip propagation time delay ∆t and a small frequency shift between the transmitted and received signals take place due to the range propagation effect. The time propagation delay for the reflected signals from each scatterer i is: where c is the speed of light, and r i is the range between the micro radar and each object which is located within the radar beamwidth. The originally transmitted signal is then mixed with the received signal in a low-pass filter to estimate the video signal x(t) that has a Beat frequency f b , which can be written as: The phase change of this video signal is then utilized to obtain the Doppler frequency f doppler . In each chirp, this video signal is sampled with a 264-ns sampling rate. A two-dimensional Fast Fourier Transform (FFT) is then applied on the sampled signal [28] to form the radar RDM for each of the three receiving antennae. A mean RDM is generated by averaging the RDM from all antennae. This mean RDM has 256 × 256 pixels with a 32-Bit amplitude value for each pixel, which is used to obtain the radar range and velocity measurements for different scatterers. Figure 2 illustrates an RDM image where the X-axis represents the velocity measurements and the Y-axis provides the range measurements. The 32-Bit amplitude value at each pixel represents the received echo signals strength from different ground scatterers.
where represents the transmitted chirp signal bandwidth. A roundtrip propagation time delay ∆ and a small frequency shift between the transmitted and received signals take place due to the range propagation effect. The time propagation delay for the reflected signals from each scatterer is: where is the speed of light, and is the range between the micro radar and each object which is located within the radar beamwidth. The originally transmitted signal is then mixed with the received signal in a low-pass filter to estimate the video signal ( ) that has a Beat frequency , which can be written as: The phase change of this video signal is then utilized to obtain the Doppler frequency . In each chirp, this video signal is sampled with a 264-ns sampling rate. A two-dimensional Fast Fourier Transform (FFT) is then applied on the sampled signal [28] to form the radar RDM for each of the three receiving antennae. A mean RDM is generated by averaging the RDM from all antennae. This mean RDM has 256 × 256 pixels with a 32-Bit amplitude value for each pixel, which is used to obtain the radar range and velocity measurements for different scatterers. Figure 2 illustrates an RDM image where the X-axis represents the velocity measurements and the Y-axis provides the range measurements. The 32-Bit amplitude value at each pixel represents the received echo signals strength from different ground scatterers. This RDM is acquired from the micro radar through an Ethernet cable after the radar signal processing inside the radar is performed. This generated map is then utilized for the target detection, velocity estimation, and height extraction purposes.

Target Detection and Data Extraction
The target detection process of radar is commonly achieved by utilizing the Constant False Alarm Rate (CFAR) algorithm. This algorithm adaptively estimates the threshold power level by locally comparing the power for the cell-under-test (CUT) against its neighborhood cells (background) This RDM is acquired from the micro radar through an Ethernet cable after the radar signal processing inside the radar is performed. This generated map is then utilized for the target detection, velocity estimation, and height extraction purposes.

. Target Detection and Data Extraction
The target detection process of radar is commonly achieved by utilizing the Constant False Alarm Rate (CFAR) algorithm. This algorithm adaptively estimates the threshold power level by locally comparing the power for the cell-under-test (CUT) against its neighborhood cells (background) [29]. The target detection is achieved when the power level of the CUT exceeds the estimated threshold level. Many ground radar stations and airborne radars are relying on the CFAR for the aircraft detection purpose. Although the CFAR is an effective detection approach for these applications, it is not effective for the proposed system since the ground scatterers have approximately the same power level on the RDM. The CFAR-detected targets for the RDM image of Figure 2 is shown in Figure 3 where CFAR could detect a part of the main reflected arc of the ground objects while it misses the rest of it. On the other hand, the CFAR can detect a false target (noise), which has a prominent power level from its background.  [29]. The target detection is achieved when the power level of the CUT exceeds the estimated threshold level. Many ground radar stations and airborne radars are relying on the CFAR for the aircraft detection purpose. Although the CFAR is an effective detection approach for these applications, it is not effective for the proposed system since the ground scatterers have approximately the same power level on the RDM. The CFAR-detected targets for the RDM image of Figure 2 is shown in Figure 3 where CFAR could detect a part of the main reflected arc of the ground objects while it misses the rest of it. On the other hand, the CFAR can detect a false target (noise), which has a prominent power level from its background. The first issue occurs when the CUT is picked inside a patched area of real ground scatterers. In this case, the CUT has approximately the same power level as its surrounding background and the CFAR could not detect all the targets since the CUT has the same power level as its neighborhood, which makes it not distinctive.
The second issue results from random noises, which have relatively high-power levels compared to its local neighborhood cells. Unlike the CFAR, which depends on the local neighborhood cells for target detection, an efficient target detection technique is proposed based on the Gaussian kernel and local maxima to avoid problems associated with the CFAR. The Gaussian kernel is convolved with the RDM as a filtering step to depress the isolated noise responses, which may affect the RO estimated height and velocity. A local maxima-based approach is then applied as a global threshold to obtain the strongest five candidates on the whole image which are considered to be the detected targets, as shown in Figure 4. The target detection process takes approximately 1.3 ms, which is convenient for real-time applications.
The radial velocity and the range for each detected target are then acquired directly from the Xand Y-axes of this image, respectively. These radial velocities are projected by the radar tilt angle, which is 60° with respect to the flight path direction in addition to the vehicle pitch angle. These projected components are then averaged to obtain the resultant forward velocity of the vehicle in the body frame. On the other hand, the ranges are projected and compensated by the same angles with respect to the vertical direction and then averaged to obtain the height above the ground level for the vehicle. The first issue occurs when the CUT is picked inside a patched area of real ground scatterers. In this case, the CUT has approximately the same power level as its surrounding background and the CFAR could not detect all the targets since the CUT has the same power level as its neighborhood, which makes it not distinctive.
The second issue results from random noises, which have relatively high-power levels compared to its local neighborhood cells. Unlike the CFAR, which depends on the local neighborhood cells for target detection, an efficient target detection technique is proposed based on the Gaussian kernel and local maxima to avoid problems associated with the CFAR. The Gaussian kernel is convolved with the RDM as a filtering step to depress the isolated noise responses, which may affect the RO estimated height and velocity. A local maxima-based approach is then applied as a global threshold to obtain the strongest five candidates on the whole image which are considered to be the detected targets, as shown in Figure 4. The target detection process takes approximately 1.3 ms, which is convenient for real-time applications.

Enhanced Moncular Visual Odometry
A GoPro HERO4 camera is attached to the quadcopter with a resolution of 1080 × 1920, and a video measurement rate of 30 frames per second. This camera is mounted on the UAV to have a downward facing orientation. The proposed VO consists of two main steps, which are the monocular VO based on the optical flow, and regression tree-based approach for the velocity drift error compensation.

Monocular Visual Odometry
The proposed VO extracts the optical flow vectors in the X and Y directions by detecting the features of interest from the video frames using the Speeded Up Robust Features (SURF) detector [30]. These detected features are then matched between successive frames. The main reason for utilizing the SURF algorithm is that it has a low computational load, which makes it more convenient for real-time operation. An M-estimator Sample Consensus (MSAC) algorithm is employed for outlier (incorrect match) rejection purpose. Suppose a point = [ , , ] in the space is projected by a pinhole camera model to an image plane at point = [ , , ], which is described as: where is the camera's focal length. Since the camera is mounted perpendicularly to a vehicle body, the coordinate is approximately equal to the distance between ground and camera's projection origin. The estimated ground distance using radar measurement is shown in Figure 5.
By having this ground distance, the displacement in the image plane (∆ , ∆ ) can be converted to a real-world displacement (∆ , ∆ ) as following: The displacement in the image plane is obtained after removing the outliers from the matched points between two successive frames. As the computed displacement (∆ , ∆ ) is usually in pixels, it is required to convert it into real-world units (e.g., meters), which is written as: The radial velocity and the range for each detected target are then acquired directly from the X-and Y-axes of this image, respectively. These radial velocities are projected by the radar tilt angle, which is 60 • with respect to the flight path direction in addition to the vehicle pitch angle. These projected components are then averaged to obtain the resultant forward velocity of the vehicle in the body frame. On the other hand, the ranges are projected and compensated by the same angles with respect to the vertical direction and then averaged to obtain the height above the ground level for the vehicle.

Enhanced Moncular Visual Odometry
A GoPro HERO4 camera is attached to the quadcopter with a resolution of 1080 × 1920, and a video measurement rate of 30 frames per second. This camera is mounted on the UAV to have a downward facing orientation. The proposed VO consists of two main steps, which are the monocular VO based on the optical flow, and regression tree-based approach for the velocity drift error compensation.

Monocular Visual Odometry
The proposed VO extracts the optical flow vectors in the X and Y directions by detecting the features of interest from the video frames using the Speeded Up Robust Features (SURF) detector [30]. These detected features are then matched between successive frames. The main reason for utilizing the SURF algorithm is that it has a low computational load, which makes it more convenient for real-time operation. An M-estimator Sample Consensus (MSAC) algorithm is employed for outlier (incorrect match) rejection purpose. Suppose a point P = [X, Y, Z] in the space is projected by a pinhole camera model to an image plane at point p = x pix , y pix , f cam , which is described as: where f cam is the camera's focal length. Since the camera is mounted perpendicularly to a vehicle body, the coordinate Z is approximately equal to the distance between ground and camera's projection origin. The estimated ground distance using radar measurement is shown in Figure 5.
the X and Y directions. These compensated vectors are then utilized to estimate the vehicle forward velocity V VO , which is described as: where ∆t is the time between two consecutive frames.

Velocity Compensation
During this velocity compensation phase, a regression tree-based approach is employed to enhance the accuracy of the estimated vehicle velocity from the monocular VO and compensate the employed approximations in the previous velocity estimation step. This training process is achieved during the availability of the GNSS signals. The GNSS signals are assumed to be available during the first 50 and 140 s of the first and second flights, respectively. The training process takes place during the first 40 s (from 1 to 40 s) for the first flight and the remaining 10 s (40 to 50 s) of the GNSS signals availability are utilized to obtain the weights for the predicted velocity from the regression trees. These weights are then utilized for the averaging process between the predicted velocity from the regression trees and the estimated vehicle forward velocity from the monocular VO. During the second flight, the whole collected data during the first flight and the first 125 s (from 1 to 125 s of the GNSS signal availability) of the second flight are utilized for the training purpose. The remaining 15 s (125 to 140 s of the GNSS signal availability) are utilized to obtain the weights for the predicted velocity from the regression trees.
The video frames are divided into 3 × 3 cells, and the optical flow vectors in the X and Y directions are then averaged inside each cell to offer a fixed number of vectors to the regression trees during this stage, as shown in Figure 6. Theses averaged optical flow vectors along each cell, the estimated forward velocity from the monocular VO, roll, pitch, and RO height are then utilized as inputs for the regression trees during the training session while the ground truth forward velocity, which is obtained from the GNSS/INS/magnetometer/barometer integration, is utilized as an output during this training phase. Figure 7 illustrates the proposed system architecture during the training session.
When the GNSS signals are unavailable, the trained regression trees start to estimate the vehicle forward velocity in an attempt to compensate the monocular VO drift errors caused by different reasons, such as featureless area, inconsistent matches, and the change of the camera calibration By having this ground distance, the displacement in the image plane (∆u, ∆v) can be converted to a real-world displacement (∆X, ∆Y) as following: The displacement in the image plane is obtained after removing the outliers from the matched points between two successive frames. As the computed displacement (∆u, ∆v) is usually in pixels, it is required to convert it into real-world units (e.g., meters), which is written as: where s is the pixel size. The optical vectors (u, v) are obtained by multiplying these image displacements by the camera's measurements rate. The gyro's measurements ω x , ω y are then utilized to compensate the vehicle rotational motion effect from the estimated optical flow vectors in the X and Y directions. These compensated vectors are then utilized to estimate the vehicle forward velocity V VO , which is described as: where ∆t is the time between two consecutive frames.

Velocity Compensation
During this velocity compensation phase, a regression tree-based approach is employed to enhance the accuracy of the estimated vehicle velocity from the monocular VO and compensate the employed approximations in the previous velocity estimation step. This training process is achieved during the availability of the GNSS signals. The GNSS signals are assumed to be available during the first 50 and 140 s of the first and second flights, respectively. The training process takes place during the first 40 s (from 1 to 40 s) for the first flight and the remaining 10 s (40 to 50 s) of the GNSS signals availability are utilized to obtain the weights for the predicted velocity from the regression trees. These weights are then utilized for the averaging process between the predicted velocity from the regression trees and the estimated vehicle forward velocity from the monocular VO. During the second flight, the whole collected data during the first flight and the first 125 s (from 1 to 125 s of the GNSS signal availability) of the second flight are utilized for the training purpose. The remaining The video frames are divided into 3 × 3 cells, and the optical flow vectors in the X and Y directions are then averaged inside each cell to offer a fixed number of vectors to the regression trees during this stage, as shown in Figure 6. Theses averaged optical flow vectors along each cell, the estimated forward velocity from the monocular VO, roll, pitch, and RO height are then utilized as inputs for the regression trees during the training session while the ground truth forward velocity, which is obtained from the GNSS/INS/magnetometer/barometer integration, is utilized as an output during this training phase. Figure 7 illustrates the proposed system architecture during the training session.
Sensors 2018, 18, 2776 9 of 27 parameters due to the strong vehicle vibration effect during the flight. A weighted average between the predicted velocity from the regression trees and the estimated velocity from the monocular VO is utilized as a measurements update for the EKF. These weights are obtained by computing the Root Mean Square (RMS) errors for the predicted velocity from the regression trees and the estimated velocity from the monocular VO with respect to the reference forward velocity from 40 to 50 s for the first flight and from 125 to 140 s for the second flight. The estimated forward velocity from the RO is also utilized as a measurement update during this GNSS outage period. Figure 8 illustrates the proposed system architecture during the prediction session. The regression tree-based approach is employed since it has the ability to handle the situation of missing inputs (optical flow vectors) in some image parts due to the lack of the observed features or inconsistent matches caused by repeated patterns.   parameters due to the strong vehicle vibration effect during the flight. A weighted average between the predicted velocity from the regression trees and the estimated velocity from the monocular VO is utilized as a measurements update for the EKF. These weights are obtained by computing the Root Mean Square (RMS) errors for the predicted velocity from the regression trees and the estimated velocity from the monocular VO with respect to the reference forward velocity from 40 to 50 s for the first flight and from 125 to 140 s for the second flight. The estimated forward velocity from the RO is also utilized as a measurement update during this GNSS outage period. Figure 8 illustrates the proposed system architecture during the prediction session. The regression tree-based approach is employed since it has the ability to handle the situation of missing inputs (optical flow vectors) in some image parts due to the lack of the observed features or inconsistent matches caused by repeated patterns.   When the GNSS signals are unavailable, the trained regression trees start to estimate the vehicle forward velocity in an attempt to compensate the monocular VO drift errors caused by different reasons, such as featureless area, inconsistent matches, and the change of the camera calibration parameters due to the strong vehicle vibration effect during the flight. A weighted average between the predicted velocity from the regression trees and the estimated velocity from the monocular VO is utilized as a measurements update for the EKF. These weights are obtained by computing the Root Mean Square (RMS) errors for the predicted velocity from the regression trees and the estimated velocity from the monocular VO with respect to the reference forward velocity from 40 to 50 s for the first flight and from 125 to 140 s for the second flight. The estimated forward velocity from the RO is also utilized as a measurement update during this GNSS outage period. Figure 8 illustrates the proposed system architecture during the prediction session. The regression tree-based approach is employed since it has the ability to handle the situation of missing inputs (optical flow vectors) in some image parts due to the lack of the observed features or inconsistent matches caused by repeated patterns.

Data Fusion
This section describes the data fusion between INS, magnetometer, barometer, enhanced VO, RO, and GNSS measurements during its availability in an EKF. The navigation states, which are the position, velocity, and attitudes in the navigation frame (n-frame) are derived from the IMU raw measurements through a mechanization process. The EKF error states vector include 21 states as follows: where r n , v n , n are the (position, velocity and attitude) error vector of INS mechanization, respectively, and b and d are the accelerometers bias and gyros drift, respectively. Finally, and are the accelerometers and gyros scale factor, respectively. The whole sensors are fused in loosely coupled fashion through two main steps inside the EKF, the prediction phase and the measurements update phase.

Prediction Model
The system model, which describes how different INS error states evolve with time, is obtained by linearly perturbing the mechanization equations and can be represented [31,32] as follows:

Data Fusion
This section describes the data fusion between INS, magnetometer, barometer, enhanced VO, RO, and GNSS measurements during its availability in an EKF. The navigation states, which are the position, velocity, and attitudes in the navigation frame (n-frame) are derived from the IMU raw measurements through a mechanization process. The EKF error states vector include 21 states as follows: where δr n , δv n , ε n are the (position, velocity and attitude) error vector of INS mechanization, respectively, and b and d are the accelerometers bias and gyros drift, respectively. Finally, s a and s g are the accelerometers and gyros scale factor, respectively. The whole sensors are fused in loosely coupled fashion through two main steps inside the EKF, the prediction phase and the measurements update phase.

Prediction Model
The system model, which describes how different INS error states evolve with time, is obtained by linearly perturbing the mechanization equations and can be represented [31,32] as follows: where Φ k,k+1 is the state transition matrix,x + k is the error states from the previous epoch,x − k+1 is the predicted error states, and w k is the system noise. The inertial sensor stochastic errors are modeled as a first-order Gauss-Markov process. The prediction of state-error covariance matrix P − k+1 at a certain epoch is: where Q k is the covariance matrix of the system noise, P + k is the error states covariance matrix from the previous epoch, and G k is the noise coefficient matrix.

Observation Model
The differences between the GNSS positioning measurements P where e P is the vector of GNSS position measurement errors and can be expressed as: where e n , e e , e d are the positioning errors in the North, East and Down, respectively.
where R is the equatorial earth radius, e is the eccentricity of the earth ellipsoid, and ϕ is the latitude. The positioning measurements which are utilized for updating the EKF are: (17) where the parameter l b GNSS × is the skew-symmetric form of the lever-arm between the GNSS antenna and the IMU in the body frame, C l b is the rotation matrix between the body and local level frame, δP ned IMU are the positions error states, and ∅ is the attitude errors which can be expressed in a skew-symmetric form E ned (or cross product ∅×) form [33] as: where n , e and d are the attitude errors in the North, East and Down directions, respectively. From Equation (17), the design matrix can be expressed as: Finally, the innovation sequence between the measurement updates z GNSS and the estimated measurementsẑ GNSS is obtained as: The heading measurements ψ mag is acquired from the magnetometer raw measurements (3D magnetic field components) as: where M x M y M z are the unit-vector measurements of magnetic north in the body frame, Φ is the roll angle, θ is the pitch angle and δ mag is the magnetic declination, which is the bearing the difference between the true north and magnetic north, and e ψ mag is the magnetometer heading measurement error. The estimated Direction Cosine Matrix (DCM)Ĉ l b , which is derived from the mechanization process during the EKF prediction stage, is utilized to compute the heading angleψ. The heading measurement, which is used for updating the EKF, is the difference between the estimated headingψ and the measured heading ψ mag and can be written as: The error equation can be obtained by perturbing Equation (22) as: Hence, the design matrix can be expressed as: The innovation sequence between the measurement update z mag and the estimated measurementŝ z mag is obtained as: The barometer is utilized to measure the height of the vehicle h baro as: where e h baro is the barometer height measurement error. The design matrix can be expressed as: The innovation sequence between the measurement update z baro and the estimated measurementŝ z baro is obtained as: The estimated velocity in the flight path direction, which is obtained from the RO, can be written as: where V s RO is the velocity acquired from the RO represented in the sensor frame and e V is the velocity measurement noise. The difference between the predicted forward velocity from the INS in the sensor frame and the extracted forward velocity from the RO is utilized as a measurement update during the GNSS signal outage. The attitude errors and the angular rate errors must be incorporated in the derived velocity in the mechanization process. The EKF measurement updates can be expressed as: whereV s RO is the computed RO forward velocity from the INS represented in the sensor frame, V ned IMU × is the skew-symmetric form of the velocity of the vehicle at the center of the IMU represented in the navigation frame, l b RO × is the skew-symmetric form of the lever arm between the RO and the IMU represented in the body frame, C b l is the rotation matrix between the navigation and body frame. δV ned IMU is the velocity error of the vehicle at the center of the IMU and δω b ib is the angular rate measurement error. From Equation (30), the design matrix can be expressed as: The innovation sequence between the measurement updates z V RO and the estimated measurementsẑ V RO is obtained as: The enhanced monocular VO velocity update is derived in the same manner as the RO velocity update. After the innovation sequence δz k and design matrix H k computation for each sensor, the Kalman gain K k , the updated statesx + k , and its updated covariance matrix P + k can be obtained as:

Experimental Results
The results shown in the upcoming subsections are obtained from two-real flights with 3DR Solo quadcopter equipped with a micro-FMCW radar and a GoPro Hero4 camera.

Hardware Setup
A GoPro Hero4 camera with a fisheye lens is attached to the UAV to get a High Definition (HD) video measurement with 30 frames per second. Two real flights are performed in different places.
During the first flight, the camera's field of view is adjusted to be wide with a resolution of 1080 × 1920, while the field of view is adjusted to be medium with the same resolution during the second flight. A K-MD2 radar module is attached to the quadcopter belly through a wooden frame. This radar has a 24-GHz transmitted frequency with one transmitter and three receiver antennae. The radar range is measured up to 300 m with a 1-m range resolution and velocity measurements up to 40 m/s with a 0.6-m/s resolution accuracy. The UAV has a compact size as the dimensions of the radar are 120 × 72 × 15 mm, the radar weight is 165 g which makes it suitable for small UAVs, and the power consumption of this radar is 7.2 watt at 12 V. Figure 9 illustrates the utilized radar.
The results shown in the upcoming subsections are obtained from two-real flights with 3DR Solo quadcopter equipped with a micro-FMCW radar and a GoPro Hero4 camera.

Hardware Setup
A GoPro Hero4 camera with a fisheye lens is attached to the UAV to get a High Definition (HD) video measurement with 30 frames per second. Two real flights are performed in different places. During the first flight, the camera's field of view is adjusted to be wide with a resolution of 1080 × 1920, while the field of view is adjusted to be medium with the same resolution during the second flight. A K-MD2 radar module is attached to the quadcopter belly through a wooden frame. This radar has a 24-GHz transmitted frequency with one transmitter and three receiver antennae. The radar range is measured up to 300 m with a 1-m range resolution and velocity measurements up to 40 m/s with a 0.6-m/s resolution accuracy. The UAV has a compact size as the dimensions of the radar are 120 × 72 × 15 mm, the radar weight is 165 g which makes it suitable for small UAVs, and the power consumption of this radar is 7.2 watt at 12 V. Figure 9 illustrates the utilized radar.

Figure 9. Utilized Frequency Modulated Continuous Wave (FMCW) radar (RF-beam).
This radar is attached to a UDOO X86 single board computer. This computer is based on Quad Core 64-bit new-generation X86 processors made by Intel ® . The attached camera and radar to the SOLO quadcopter are shown in Figure 10. A 3DR Solo Quadcopter is utilized during the flights 'missions. This UAV has a Pixhawk-2 autopilot with an InvenSense MPU-6000 MEMS IMU, an MS5611 barometer and a U-blox GPS. Figure 11 illustrates a block diagram of the hardware configuration.
The experiments were conducted in two different places, with different trajectories, and the radar was pitched by 60° from the quadcopter body. The first flight was performed over a football playground while the second flight was performed over multiple objects with different altitudes, such as a house, trees, grass, cars, and hangars. Figure 12 shows aerial images, which were obtained from the first and second flights. This radar is attached to a UDOO X86 single board computer. This computer is based on Quad Core 64-bit new-generation X86 processors made by Intel ® . The attached camera and radar to the SOLO quadcopter are shown in Figure 10. A 3DR Solo Quadcopter is utilized during the flights "missions". This UAV has a Pixhawk-2 autopilot with an InvenSense MPU-6000 MEMS IMU, an MS5611 barometer and a U-blox GPS. Figure 11 illustrates a block diagram of the hardware configuration.  The experiments were conducted in two different places, with different trajectories, and the radar was pitched by 60 • from the quadcopter body. The first flight was performed over a football playground while the second flight was performed over multiple objects with different altitudes, such as a house, trees, grass, cars, and hangars. Figure 12 shows aerial images, which were obtained from the first and second flights.

First Experiment
The first flight trajectory is composed of 10 waypoints during the total flight time of 393 s with a maximum speed of 5 m/s, as shown in Figure 13.  In Figure 14, the estimated velocity from the typical closed form monocular VO, the enhanced monocular VO, and the RO, with RMS error values of 1.29, 1.02, and 0.49 m/s, respectively, are compared to the UAV reference forward velocity in the body frame.
Three outage scenarios were carried, with different outage periods, ranging from 30 s to 113 s. The first outage period was performed for 30 s. Figures 15 and 16 show comparisons between the estimated 2D flight trajectory outage segments from the GNSS/INS integration (ground truth segment) and an enhanced monocular VO aided navigation system during the first flight for 30 and 113 s of GNSS signal outage, respectively. Figures 17 and 18 show comparisons between the estimated 2D flight trajectory outage segments from the GNSS/INS integration (ground truth segment) and the proposed integrated system aided navigation system during the first flight for 30 and 113 s of GNSS signal outage, respectively. Figure 19 demonstrates the ability of the proposed system (RO/enhanced monocular VO/magnetometer/barometer) to mitigate the INS drift errors when the GNSS signals get lost, and to enhance the 3D Root Mean Square Error (RMSE) positioning accuracy to be 3.2 m in 113 s. Figure 20 shows a comparison between the estimated 2D flight trajectory outage segments from the GNSS/INS integration (ground truth segment) and INS in a standalone mode during the first flight for 30 s of GNSS signal outage. Figure 21 illustrates the navigation errors for the INS in the standalone mode in the North and East directions during 30 s of GNSS signal outage. In Figure 14, the estimated velocity from the typical closed form monocular VO, the enhanced monocular VO, and the RO, with RMS error values of 1.29, 1.02, and 0.49 m/s, respectively, are compared to the UAV reference forward velocity in the body frame.  Three outage scenarios were carried, with different outage periods, ranging from 30 s to 113 s. The first outage period was performed for 30 s. Figures 15 and 16 show comparisons between the estimated 2D flight trajectory outage segments from the GNSS/INS integration (ground truth segment) and an enhanced monocular VO aided navigation system during the first flight for 30 and 113 s of GNSS signal outage, respectively. Figures 17 and 18 show comparisons between the estimated 2D flight trajectory outage segments from the GNSS/INS integration (ground truth segment) and the proposed integrated system aided navigation system during the first flight for 30 and 113 s of GNSS signal outage, respectively.                     Table 1 provides comparisons of the RMS errors for the position states, which are obtained from the INS in a standalone mode, and enhanced monocular VO, and the proposed integrated system aided navigation during the GNSS outage periods. The results demonstrate the ability of the proposed integrated system to reduce the 3D positioning errors to 2.06% for 30 s, and 0.13% during 113 s of the INS drift errors in the standalone mode during the GNSS signals outages period.
The results also demonstrate the proposed integrated system's capability of reducing the 3D positioning RMS errors to 78.04%, 71.65%, and 65.71% of the enhanced monocular VO. Figure 22 shows a comparison between the RMS errors of the 3D positioning for the proposed integrated system aided navigation during three GNSS outage periods ranged from 30 s to 113 s.   Table 1 provides comparisons of the RMS errors for the position states, which are obtained from the INS in a standalone mode, and enhanced monocular VO, and the proposed integrated system aided navigation during the GNSS outage periods. The results demonstrate the ability of the proposed integrated system to reduce the 3D positioning errors to 2.06% for 30 s, and 0.13% during 113 s of the INS drift errors in the standalone mode during the GNSS signals outages period. The results also demonstrate the proposed integrated system's capability of reducing the 3D positioning RMS errors to 78.04%, 71.65%, and 65.71% of the enhanced monocular VO. Figure 22 shows a comparison between the RMS errors of the 3D positioning for the proposed integrated system aided navigation during three GNSS outage periods ranged from 30 s to 113 s.

Second Experiment
The second flight trajectory is composed of 18 waypoints of total flight time 393 s, with a maximum speed of 5 m/s, as shown in Figure 23. In Figure 24, the estimated velocity from the typical closed form monocular VO, the enhanced monocular VO, and the RO, with RMS error values of 0.61, 0.53, and 0.75 m/s, respectively, are compared to the UAV reference forward velocity in the body frame.
The accuracy for the estimated velocity from the RO was slightly less than those of the monocular VO, and the enhanced monocular VO, because the flight was performed over multiple objects with different altitudes, ranges, and angles.

Second Experiment
The second flight trajectory is composed of 18 waypoints of total flight time 393 s, with a maximum speed of 5 m/s, as shown in Figure 23.

Second Experiment
The second flight trajectory is composed of 18 waypoints of total flight time 393 s, with a maximum speed of 5 m/s, as shown in Figure 23. In Figure 24, the estimated velocity from the typical closed form monocular VO, the enhanced monocular VO, and the RO, with RMS error values of 0.61, 0.53, and 0.75 m/s, respectively, are compared to the UAV reference forward velocity in the body frame.
The accuracy for the estimated velocity from the RO was slightly less than those of the monocular VO, and the enhanced monocular VO, because the flight was performed over multiple objects with different altitudes, ranges, and angles. In Figure 24, the estimated velocity from the typical closed form monocular VO, the enhanced monocular VO, and the RO, with RMS error values of 0.61, 0.53, and 0.75 m/s, respectively, are compared to the UAV reference forward velocity in the body frame.   The accuracy for the estimated velocity from the RO was slightly less than those of the monocular VO, and the enhanced monocular VO, because the flight was performed over multiple objects with different altitudes, ranges, and angles.
Four outage scenarios were carried, with different outage periods, ranging from 60 s to 240 s. The first outage period was performed for 60 s. Figures 25 and 26 show comparisons between the estimated 2D flight trajectory outage segments from the GNSS/INS integration (ground truth segment) and enhanced monocular VO aided navigation during the second flight for 60 s and 240 s of GNSS signal outage, respectively.                       Table 2 provides comparisons of the RMS errors for the position states, which are obtained from the INS in a standalone mode, and the enhanced monocular VO, and the proposed integrated system aided navigation during the GNSS outages periods. The results demonstrate the ability of the proposed integrated system to reduce the 3D positioning errors to 3.09 % for 60 s, and 0.1% during 240 s of the INS drift errors in a standalone mode during the period of GNSS signal outages.     Table 2 provides comparisons of the RMS errors for the position states, which are obtained from the INS in a standalone mode, and the enhanced monocular VO, and the proposed integrated system aided navigation during the GNSS outages periods. The results demonstrate the ability of the proposed integrated system to reduce the 3D positioning errors to 3.09 % for 60 s, and 0.1% during 240 s of the INS drift errors in a standalone mode during the period of GNSS signal outages.   Table 2 provides comparisons of the RMS errors for the position states, which are obtained from the INS in a standalone mode, and the enhanced monocular VO, and the proposed integrated system aided navigation during the GNSS outages periods. The results demonstrate the ability of the proposed integrated system to reduce the 3D positioning errors to 3.09 % for 60 s, and 0.1% during 240 s of the INS drift errors in a standalone mode during the period of GNSS signal outages. The results demonstrate the proposed integrated system aided navigation is capable of reducing the positioning RMSEs to 92.53%, 99.41%, and 92.38% of the enhanced monocular VO during 60, 120, and 180 s of GNSS signal outages, respectively, while it has almost the same accuracy as the enhanced monocular VO during 240 s. Larger data than that for the first flight is used to train the enhanced monocular VO, which helps with an accurate estimate for the vehicle forward velocity. Figure 32 shows a comparison of the 3D positioning RMSEs of the proposed integrated system aided navigation during four GNSS outage periods ranged from 60 s to 240 s. The results demonstrate the proposed integrated system aided navigation is capable of reducing the positioning RMSEs to 92.53%, 99.41%, and 92.38% of the enhanced monocular VO during 60, 120, and 180 s of GNSS signal outages, respectively, while it has almost the same accuracy as the enhanced monocular VO during 240 s. Larger data than that for the first flight is used to train the enhanced monocular VO, which helps with an accurate estimate for the vehicle forward velocity. Figure 32 shows a comparison of the 3D positioning RMSEs of the proposed integrated system aided navigation during four GNSS outage periods ranged from 60 s to 240 s.

Conclusions
An integrated navigation system for UAVs in GNSS denied environments is proposed based on micro FMCW radar and a single camera. The vehicle's forward velocity is estimated from both the RO and the enhanced VO to enhance the INS navigation accuracy during GNSS signal outages. In addition, the estimated height from the RO is utilized to resolve the monocular VO's scale ambiguity. An efficient target detection approach is proposed to detect the ground objects based on a Gaussian kernel and local maxima algorithms. These detected targets are then utilized to estimate the forward velocity and the height above ground level for the vehicle. An optical flow and regression tree-based approaches are utilized to implement the enhanced monocular VO. The optical flow is employed for

Conclusions
An integrated navigation system for UAVs in GNSS denied environments is proposed based on micro FMCW radar and a single camera. The vehicle's forward velocity is estimated from both the RO and the enhanced VO to enhance the INS navigation accuracy during GNSS signal outages.
In addition, the estimated height from the RO is utilized to resolve the monocular VO's scale ambiguity. An efficient target detection approach is proposed to detect the ground objects based on a Gaussian kernel and local maxima algorithms. These detected targets are then utilized to estimate the forward velocity and the height above ground level for the vehicle. An optical flow and regression tree-based approaches are utilized to implement the enhanced monocular VO. The optical flow is employed for the forward velocity estimation purpose, while its associated drift errors are compensated based on a trained regression tree model. Theses estimated forward velocities from the RO and enhanced VO are then fused with the IMU, barometer, and magnetometer measurements via an EKF. The experimental results demonstrate the proposed system's ability to enhance the average 3D positioning errors for the first flight by 98.65%, and to by 99.04% for the second flight with respect to the INS on the stand alone mode. In addition, it improves the average 3D positioning accuracy by 28.21% for the first flight, and 4.72% for the second flight with respect to the enhanced monocular VO.
Unlike other proposed RO methods which utilize a large unmanned aircraft or large radar or even simulate the flight missions, the proposed system utilizes a small SOLO quadcopter and a lightweight micro radar during a real flight. Such small quadcopters are typically utilized in many missions, such as search, rescue, and disaster management. The proposed algorithm has been evaluated in a generic and typical maneuvering scenario. It also avoids the various assumptions imposed by many other researches, such as straight flights, constant velocities, leveled flights, and flying over flat terrain. The RO and VO are integrated into one integrated system to help overcome their limitations. The proposed RO provides a more accurate forward velocity estimation than the enhanced monocular VO while flying over a flat terrain and it is slightly worse than the enhanced monocular VO while flying over non-flat terrain. On the other hand, the radar is immune against the environmental changes and can operate in featureless areas.


extfieldsofstudy:
Computer Science
Medicine
Mathematics

s2fieldsofstudy:
Computer Science
Mathematics

text:
Causal Inference on Discrete Data via Estimating Distance Correlations

In this paper, we deal with the problem of inferring causal directions when the data is on discrete domain. By considering the distribution of the cause $P(X)$ and the conditional distribution mapping cause to effect $P(Y|X)$ as independent random variables, we propose to infer the causal direction via comparing the distance correlation between $P(X)$ and $P(Y|X)$ with the distance correlation between $P(Y)$ and $P(X|Y)$. We infer"$X$ causes $Y$"if the dependence coefficient between $P(X)$ and $P(Y|X)$ is smaller. Experiments are performed to show the performance of the proposed method.

Introduction
Inferring the causal direction between two variables from observational data becomes a hot research topic. Additive noise models (ANMs) (Zhang and Hyvärinen, 2009a;Shimizu et al., 2006;Peters et al., 2011;Hoyer and Hyttinen, 2009;Zhang and Hyvärinen, 2009b;Shimizu et al., 2011;Hyvärinen and Smith, 2013;Hyvärinen et al., 2010;Hyttinen et al., 2012;Mooij et al., 2011) are preliminary trials to solve this problem. They assume that the effect is governed by the cause and an additive noise, and the causal inference is done by finding the direction that admits such a model. Recently, under another view of exploiting the asymmetry between cause and effect, the linear trace method (LTr) (Zscheischler et al., 2011; and information geometric causal inference (IGCI)  are proposed. Suppose X is the cause and Y is the effect. Based on the fact that the generating of P (X) is independent with that of P (Y |X) Lemeire and Janzing, 2013;Schölkopf et al., 2012), LTr suggests that the trace condition is fulfilled in the causal direction while violated in the anti-causal direction, and IGCI shows that the density of the cause and the log slope of the function transforming cause to effect are uncorrelated while the density of the effect and the log slope of the inverse of the function are positively correlated. By accessing these so-called cause-effect asymmetries, one can determine the causal direction. Then a kernel method using the framework of IGCI to deal with high dimensional variables is developed (Chen et al., 2014), and nonlinear extensions of trace method is presented (Chen et al., 2013).
In some situations, the variables of interest are on discrete domains, and researchers have adopted additive noise models to discrete data for causal inference (Peters et al., 2010(Peters et al., , 2011. Given observations of the variable pair, they do regressions in both directions, and test the independence between the residuals and the regressors. The direction that admits an additive noise model is inferred as the causal direction. However, ANMs may not be suitable for modeling discrete variables in some situations. For example, it is not natural to adopt ANMs to modeling categorical variables. Methods with wider applicability would be valuable for causal inference on discrete data. Motivated by the postulate that the generating of P (X) is independent with that of P (Y |X), we suggest that the (P (x), P (Y |x)) is an observation of a pair of variables that are independent with each other (here P (x) = P (X = x), referring to the probability at one specified point). To infer the causal direction, we calculate the dependence coefficients between P (X) and P (Y |X), and P (Y ) and P (X|Y ). Then the direction that induces smaller correlation is inferred as the causal direction. Without a functional causal model assumption, our method is with wider applicability than traditional ANMs. Various experiments are conducted to demonstrate the performance of the proposed method. This paper is organized as follows. Section 2 defines the problem. Section 3 presents the causal inference principle. Section 4 gives the detailed causal inference method. Section 5 shows the experiments, and section 6 concludes the whole paper.

Problem Description
Suppose we have a set of observed variables X and Y with support domain X and Y respectively. X is the cause and Y is the effect but we do not have prior knowledge about the causal direction. We assume that they are on discrete domain (for continuous variable we can perform discretization) and there are no latent confounders. We want to identify the causal direction ("X causes Y" or "Y causes X"). For clarity, we list the symbols that may appear in the following sections in table 1. Since we constrain variables in finite range, P (X, Y ) and P (Y |X) can be written as matrices, and P (Y |x) is a vector in the matrix P (Y |X). We use these representations in the rest of the paper.

Causal Inference Principle
In this section, we will present the principle that we are using for causal inference on discrete data. We start with the basic idea.

Basic Idea
The basic idea is to consider the (P (x), P (Y |x)) as a realization of a variable pair, and the two variables (one is one dimensional and the other is high dimensional) are independent with each other. See figure 1 for an example. Figure 1 shows an example of P (X) and P (Y |X). Suppose |X | = M and |Y| = L (here M = 7 and L = 8). P (X) is a vector in R M and P (Y |X) is a matrix in R L×M . The highlighted (red bars) grids are a pair (P (x), P (Y |x)). Consider P (X) and P (Y |X) as two independent random variables. The generating of the P (X) and P (Y |X) is done by drawing realizations (P (x), P (Y |x)) at each possible value of x (shifting the red bars from right to left). We have |X | realizations. We formalize this in postulate 1.
Postulate 1. P (X) and P (Y |X) are both random variables taking realizations at different x. P (X) is independent with P (Y |X).
Once we have this postulate in mind, one could seek for some properties induced by it for causal discovery. We will discuss this in the next section.

Distance Correlation
If we want to characterize the dependence between P (X) and P (Y |X), one measurement is the correlation. However, P (Y |X) is a high dimensional random vector. Adopting traditional dependence coefficients like Pearson correlations would cause certain estimation bias when sample size is not large. Moreover, it would be useful if the independence between variables corresponds to 0 correlation. This is not true if we use traditional correlations. Here we propose to use distance correlation (Székely et al., 2007) as the dependence measurement.
Distance correlation is a measurement of dependence between two random variables (one dimensional or high dimensional). Suppose we have two random variables (α, β), with characteristic functions f α (t) and f β (s) respectively. Their joint characteristic function is f α,β (t, s). The distance covariance is defined as below.
Definition 1. The distance covariance C(α, β) between two random variables (α, β) is Here · refers the weighted L 2 norm, and similarly we can define distance variance C 2 (α, α) (Székely et al., 2007). Then the distance correlation is defined as and This dependence measurement is a distance metric between the joint and multiply of the marginal characteristic functions. There are other methods to measure the dependence, like mutual information and kernel independence measurements (Gretton et al., 2005). However, mutual information is hard to estimate given finite sample size. Kernel methods involve a few parameters (kernel functions, kernel widths) which is not easy to choose. So we use this metric in our paper. Then we discuss how to estimate the distance correlation empirically from data (Székely et al., 2007). Suppose we have n observations of two random variables (α, β) as For variable α and β, we can construct and then we construct matrices A and B, with its entries as Then we can estimate the empirical distance covariance as follows (Székely et al., 2007).
We can estimate the empirical distance correlation using the empirical distance covariance. The distance correlation has a property that D(α, β) = 0 implies independence between α and β. We show that this helps to identify the causal direction in the next section.

Inferring Causal Directions
In this section we discuss how to infer the causal directions. Suppose we have the joint distribution of the variable pair as P (X, Y ). We are able to factorize it in two directions and get (P (X), P (Y |X)) and (P (Y ), P (X|Y )). Each of them is a random variable pair. We define the dependence measurements of them as below.
Definition 4. The dependence measurements are defined as If postulate 1 is accepted, then in the causal direction, the distance correlation between P (X) and P (Y |X) reaches the lower bound as Since the distance correlation is nonnegative, in the anti-causal direction we have and now we get the causal inference principle as Intuitively speaking, in the causal direction we get smaller dependence coefficient between the marginal distribution and the conditional distribution than that in the anticausal direction. One thing worth attention is that the domain size should be reasonably large to generate reliable statistics. We will give the detailed causal inference method in the next section.

Causal Inference Method
In this section we give a causal inference method which identifies the causal direction via estimating the distance correlations. If X causes Y, the D X→Y should be smaller than D Y →X . However, estimating the coefficients from samples may induce random errors. So we introduce a threshold . They are significantly different if their difference is larger than , and we can decide the causal direction. Otherwise we stay undecided. We detail the inference method in table 2. From table 2, one could see that our method identifies the cause and the effect by factorizing the joint distribution in two directions and comparing the dependence coefficients (D X→Y and D Y →X ). The one with smaller distance correlation (between the marginal distribution and the conditional distribution) is inferred as the causal direction. We name it causal inference via estimating distance correlations (abbreviated as DC). Next we analyze the computational cost of our method. Suppose the sample size is n. The time for constructing the matrix recording the joint distribution is O(n), and the times for calculating D X→Y and D Y →X are O(|X | 2 ) and O(|Y| 2 ) respectively. So the total time is O(n + |X | 2 + |Y| 2 ). One can see that this method is of low computational complexity (linear with respect to sample size). This would be verified in experiments. Input: sample of the discrete variables X and Y , threshold 1. Construct the vector recording the distribution P (X) and the matrix recording the conditional distribution P (Y |X). Calculate D X→Y = D(P (X), P (Y |X)). 2. Construct the vector recording the distribution P (Y ) and the matrix recording the conditional distribution P (X|Y ). Calculate D Y →X = D(P (Y ), P (X|Y )). 3. Decide the causal direction: Else, output "No decision made".

Experiments
In this section we test the performance of our method (DC). The compared algorithm is the discrete regression (DR) (Peters et al., 2011). We perform experiments under various settings. Section 5.1 shows the performance of DC on identifying ANMs with different N . Section 5.2 presents the performance of DC when the distribution of the cause and the conditional distribution mapping cause to effect are randomly generated. Section 5.3 tests the efficiency of the algorithms. Section 5.4 discusses the choice of the threshold parameter . Section 5.5 shows the performance of DC at different decision rates. In section 5.6, we apply DC to real world cause-effect pairs (with discretization) to show its capability in solving practical problems.

Additive Noise Models
We first evaluate the accuracies of DC on identifying ANMs with different ranges of noise term N . The model is written as The probability distributions of the cause X are chosen by: (1) randomly generate a vector (length |X |) with each entry being an integer between [1, |X |/4].
(2) Normalize it to unit sum. We generate the probability distributions of the noise N using the same way. In each trial, the algorithms are forced to make a decision. For each noise setting, we randomly generate 500 functions. Thus we have 500 additive noise models. Y for each model could be different due to the randomness of the mappings. Then we sample 200, 300, 500, 1000, 2000, 4000 points for each model, and apply DC and DR to the samples. The plots showing the accuracies of the algorithms are given in figure 2. For example, when sample size is 4000, the accuracy of DR is 0.82 while that of DC is 0.78. We observe that ANMs with small |N | could sometimes yield small distance correlations in both directions. In these situations, the decision made by DC is close to a random guess. DC performs better than DR when |N | is larger. The accuracies of DC become around 0.9 when sample size is large. But DC does not correctly identify all models. This is because the difference between the estimated distance correlations is sometimes small due to estimation random errors, and this may make the decision wrong.

Models with Randomly Generated P (X) and P (Y |X)
We now test the algorithms on the models with P (X) and P (Y |X) being randomly generated. To be specific, we generate P (X) using the method in section 5.1. Then we generate |N |/4 distributions on Y as a reference set. For x ∈ X , we generate P (Y |x) by randomly taking one of the distributions in the reference set. We choose the domain size to be 1. (|X |, |Y|) = (12, 12). For each setting, we generate 500 models. For each model, we sample 200, 300, 500, 1000, 2000, 4000 points, and apply DC and DR to them. The performance is showed in figure 3.
One can see that DR has unsatisfactory performance in these scenarios. This is because DR often makes a random guess since the models do not satisfy the ANM in either direction. DC has satisfactory performance when sample size goes large. This shows that DC works in these scenarios while DR does not.

On Efficiency
This section investigates the efficiency of the algorithms. We use the same experimental setting as that in section 5.2 (setting 2 and 4). For each sample size, we run them 100 times and record the total running time (seconds) of the algorithms. The records are showed in figure 4. Figure 4 tells that DC has a higher efficiency than DR. For example, when (|X |, |Y|) = (15, 15) and sample size is 1000, DR uses around 400 seconds to finish the experiments while DC uses only 2 seconds. This is because DR searches the whole domain iteratively to find a function that yields minimum dependence between residuals and regressors, which could be time-consuming in practice.

On Parameter
In the sections above, DC is forced to make a decision at each trial ( = 0). In this section, we examine the influence of on the performance of DC. This may help to set the values of in practice. We use the same experimental setting as that in section 5.2. The domain size is choose as (|X |, |Y|) = (15, 15) and (|X |, |Y|) = (20, 20). The sample size is fixed to be 4000. We choose the parameter to be: 1. = 0.01.
For each setting, we generate 500 models and apply DC to them. The proportion of correctly identified models, proportion of wrongly identified models and proportion of non-identified models are showed in figure 5. The proportion of non-identified models becomes large when the threshold is 0.1. This is because DC becomes conservative in this situation. The proportion of correctly identified models and that of wrongly identified models decrease as the threshold goes larger. However, the accuracy (the number of correctly identified models divided by the number of correctly identified models plus the number of wrongly identified models) increases. This is reasonable since the decisions made by DC under a higher threshold are more reliable. Based on the plotted results, we observe that the accuracy and decision date of DC are acceptable when = 0.05. We suggest that 0.05 is a reasonable choice for the parameter.

On Decision Rates
In our algorithm, the threshold parameter controls the decision rates of DC. In other words, if we increase the (from 0 to 1), the decision rates decrease (from 100% to 0%). In this section, we study the influence of the parameter on the performance of DC. We follow the experimental setting in section 5.2, and we choose (|X |, |Y|) = (15, 15). We fix the sample size to be 5000, and we vary the decision rates (by changing the from 0 to 1). The plots showing the percentage of correct decisions versus the decision rate are available in figure 6. Obviously, the percentage of correct decisions decreases as the decision rate increases. For example, the percentage of correct decisions is 100% when the decision rate is less than 20%. But it becomes 77% when the decision rate is 100%. This is acceptable since the decision would be more reliable if the algorithm makes decisions based on higher .

On Real World Data
We apply DC to real world benchmark cause-effect pairs 1 . The dataset contains records of 88 cause-effect pairs. We exclude the pairs numbered 17,44,45,52,53,54,55,68,71,75. They are either multivariate causal pairs or pairs that cannot be fitted into memory using our algorithm. We apply DC, DR (Peters et al., 2011), IGCI , LiNGAM (Shimizu et al., 2006), ANM , PNL (Zhang and Hyvärinen, 2009b), CURE (Sgouritsa et al., 2015). The last 5 methods can be directly applied to the data. DC and DR could be applied to discretized data. To make the variables discrete, we process them using the following method. For a variable 1 Available at https://webdav.tuebingen.mpg.de/cause-effect/ X, if the maximum absolute value of all observations is less than 1, we process it by round(20 * X); else we process it by round(X) 2 . For each pair, we generate 50 replicates using resampling techniques, and then apply the algorithms to the causal pairs. The boxplots showing the accuracies are in figure 7. Determining causal directions on real world data is challenging since the causal mechanisms are often complex and the data records could be noisy (Pearl, 2000;Spirtes et al., 2000). However, figure 7 shows that DC has satisfactory performance on this task. The average accuracy of DC is around 72%, which is highest among all algorithms. ANM based methods (DR, LiNGAM, ANM) do not have a good performance. This may be because the assumptions of additive noise models restrict their applicability.

Conclusions
In this paper, we deal with the causal inference problem on discrete data. We consider the distribution of the cause and the conditional distribution mapping cause to effect as independent random variables, and propose to discover the causal direction via estimating and comparing the distance correlations. Encouraging experimental results are reported. This shows that inferring the causal direction using the independence postulate is a promising research direction. In future we will try to extend this method to deal with high dimensional data.
2 For pair 65, 66, 67 which contain stock returns, we process the variables by round(100 * X)