Unmanned aerial vehicle positioning method based on Bayesian network

Abstract

The invention discloses an unmanned aerial vehicle positioning method based on a Bayesian network, and the method comprises the steps: collecting unmanned aerial vehicle data, carrying out the preprocessing, obtaining sound wave data, motion data, infrared image data and environment data, carrying out the spectrum analysis based on the sound wave data, obtaining sound spectrum time sequence data, and carrying out the positioning of the unmanned aerial vehicle. Constructing a probability density thermodynamic diagram and extracting a first range through a sparse variational Gaussian process according to time, frequency and energy three-dimensional information of the sound spectrum time sequence data and environment data, and obtaining first positioning data through a Bayesian network according to the first range and motion data, and performing adaptive wavelet down-sampling on the infrared image data and the sound spectrum time sequence data in the lateral direction to obtain lateral positioning data, and performing deviation correction on the first positioning data by using the lateral positioning data to obtain second positioning. According to the method, through space-time synchronization, adaptive processing and probabilistic reasoning of sound spectrum time sequence data, the positioning precision and stability in a complex environment are improved.

Description

Technical Field

[0001] The present invention relates to the field of UAV positioning, and particularly to a positioning method for UAVs based on Bayesian networks. Background Art

[0002] UAV positioning technology has wide application requirements in fields such as environmental monitoring, logistics distribution, and disaster relief. Traditional positioning methods mainly rely on the Global Navigation Satellite System (GPS), etc. However, in complex electromagnetic environments, urban canyons, or indoor scenarios, problems such as signal occlusion and multipath effects exist, resulting in positioning failures. Therefore, multi-sensor fusion technologies, such as Inertial Measurement Units (IMUs), visual sensors, and acoustic sensors, have gradually become research hotspots, and the positioning robustness is improved by fusing multi-modal data.

[0003] However, the existing multi-sensor fusion positioning methods have the following technical bottlenecks. The lack of spatio-temporal synchronization accuracy leads to different sampling frequencies and time bases for sensors such as acoustic waves and infrared, and spatio-temporal asynchrony causes feature alignment errors, affecting the fusion accuracy. In a strong noise background, there is a lack of an adaptive mechanism for effective signal detection and frequency band noise reduction, which easily introduces pseudo-feature interference in positioning and cannot effectively track and position external sound source information. Traditional Bayesian networks do not fully utilize the correction effect of environmental data on sound speed and sensor noise, and lack directional fusion of the UAV's motion state, resulting in the accumulation of positioning deviations in complex motion scenarios. Therefore, there is an urgent need for a positioning method that can fully consider acoustic wave data and infrared vision data to improve the positioning robustness under multi-modal data. Summary of the Invention

[0004] The object of the present invention is to provide a positioning method for UAVs based on Bayesian networks.

[0005] To achieve the above object, the present invention is implemented according to the following technical solutions: The first aspect of the present invention provides a positioning method for UAVs based on Bayesian networks, including: Collect UAV data for preprocessing to obtain acoustic wave data, motion data, infrared image data, and environmental data, where the environmental data includes temperature, humidity, and air pressure; Extract outliers based on the acoustic wave data, determine effective signals according to the outliers, perform time series slicing on the effective signals, and perform spectral analysis on the slices to obtain acoustic spectrum time series data; Construct a probability density heat map through a sparse variational Gaussian process according to the three-dimensional information of time, frequency, and energy of the acoustic spectrum time series data and the environmental data, and extract a first range based on the probability density heat map; Determine the UAV's deviation state according to the first range and the motion data, obtain the main direction and motion state data according to the deviation state, and input the acoustic spectrum time series data and motion state data of the main direction into the Bayesian network to obtain the first positioning data; Determine the lateral direction based on the motion state data, perform adaptive wavelet downsampling on the infrared image data and the acoustic spectrum time series data in the lateral direction, use the downsampled data to input into the Bayesian network to obtain the lateral positioning data, and use the lateral positioning data to correct the bias of the first positioning data to obtain the second positioning.

[0006] As a further method, the preprocessing method includes: Collect the acoustic wave data, motion data, infrared image data, and environmental data of the drone, synchronize the timestamps through pulses, dynamically compensate the IMU zero bias of the motion data through a temperature sensor; perform dynamic range compression on the infrared image data to distinguish direction features; calibrate the temperature, humidity, and air pressure values of the environmental data through redundant sensor data.

[0007] As a further method, the method for obtaining the effective signal includes: Perform time-domain analysis based on the acoustic wave data to obtain time series data, calculate the mean and standard deviation of the sampling points through a sliding window according to the time series data, set the anomaly threshold to 3 times the standard deviation, mark the sampling points that deviate from the mean by more than the anomaly threshold as outliers, identify continuous anomaly regions based on the outliers using a density clustering algorithm, extrapolate the boundary of the anomaly region by 10% along the time axis as the transition region, and perform linear interpolation smoothing on the signal intensity within the transition region; Calculate the short-time energy of the acoustic wave signal based on the smoothed time series data through a time-domain energy detection algorithm, and use the sum of the mean short-time energy of the silent source period and plus or minus 2 times the standard deviation as the background noise energy. Adaptively determine the energy threshold through the maximum inter-class variance method according to the background noise energy, and use the period higher than the energy threshold as the effective signal segment.

[0008] As a further method, the method for obtaining the acoustic spectrum time series data includes: Correct the speed of sound based on the temperature, humidity, and air pressure of the environmental data, obtain the relative speed according to the absolute value of the difference between the corrected speed of sound and the speed of the drone in the motion data within the corresponding time, use the ratio of the speed of sound to the relative speed as the scaling factor and multiply by 0.05 seconds as the dynamic length. If the dynamic length is less than or equal to the fixed value of 0.2 seconds, use the fixed value as the slice length. If the dynamic length is greater than the fixed value of 0.2 seconds, use the dynamic length as the slice length, and perform time series slicing on the effective signal according to the slice length; Decompose the signal into different frequency bands through wavelet packet transform for the slices, extract the signal-to-noise ratio based on different frequency bands. If the signal-to-noise ratio is less than 5 decibels, expand the threshold weight to 1.5. If the signal-to-noise ratio is greater than or equal to 5 decibels and less than 15 decibels, keep the threshold weight as 1. If the signal-to-noise ratio is greater than or equal to 15 decibels, reduce the threshold weight to 0.8; Noise reduction is performed on different frequency bands according to the threshold weight using an adaptive threshold formula, and the adaptive threshold formula is: Where is the noise reduction threshold of the wavelet coefficient, is the static noise benchmark, is the air temperature, is the slice length, is the threshold weight, is the short-time energy of the current slice, is the background noise energy, which is refreshed every 10 seconds, is the time correlation degree, which is the Pearson correlation coefficient of 10 consecutive sampling points; The spectrogram is calculated by short-time Fourier transform for the noise-reduced slice. The Hanning window is used as the window function, and the window function overlap rate is set to 75% to obtain the acoustic spectrum time series data with three-dimensional distribution of time, frequency, and energy.

[0009] As a further method, the method for obtaining the first range includes: Extract three-dimensional features of time, frequency, and energy based on the acoustic spectrum time series data, calculate the Doppler frequency shift offset according to the corrected sound speed and the corresponding environmental data, calculate the change gradient of the short-time energy of the acoustic spectrum with respect to the timestamp in the time dimension, take the frequency band corresponding to the energy peak in the acoustic spectrum time series data as the main frequency band, calculate the spatio-temporal distribution entropy of the energy proportion of the main frequency band based on the short-time energy, and generate a spatio-temporal spectrum joint feature matrix based on the change gradient, offset, and spatio-temporal distribution entropy; According to the spatio-temporal spectrum joint feature matrix, the attitude angle of the motion data, and the IMU data, use the azimuth angle solution model to map the acoustic features to the geographical grid to obtain the spatial coordinates, and map the spatial coordinates, timestamp, and corresponding eigenvalue of the grid to the observation input of the sparse variational Gaussian process; Based on the observation input, solve the posterior distribution through the variational inference formula, and the variational inference formula is: Where is the observed data is the log marginal likelihood of is the time-frequency feature vector of the th data point, is the slice index, is the total number of slices, is the latent function value of the th data point, that is, the true sound source state, is the latent function is the Gaussian distribution of is the induced variable, which is initialized by the sparse grid method, is the induced variable The variational distribution of is the induced variable The Gaussian prior distribution of, and the covariance matrix of the prior distribution is determined by the estimated result of the noise standard deviation of the acoustic spectrum time series data. is the ambient temperature of the th slice. is the corrected sound speed of the th slice. is the preset regularization coefficient; Regarding the grid with the posterior distribution probability greater than or equal to 0.7 as the continuous region, perform morphological closing operation on the continuous region through a 3×3 pixel structural element to obtain the probability density heat map, and extract the three-dimensional minimum circumscribed bounding box as the first range.

[0010] As a further method, the method for obtaining the main direction and the motion state data according to the bias state includes: Extract the geometric center coordinates of the three-dimensional minimum circumscribed bounding box based on the first range as the estimated sound source direction, extract the synchronized motion data according to the time stamp corresponding to the center coordinates, obtain the current coordinates, yaw angle, pitch angle, and roll angle according to the motion data, and obtain the relative azimuth deviation according to the horizontal azimuth angle and pitch angle between the current coordinates and the center coordinates; Divide the front, left, right, up, and down as the bias states according to the preset thresholds of the relative azimuth deviation and the pitch angle. When the bias state is the front, the main direction is the current yaw angle of the UAV. When the bias state is the left or right, the main direction is the sum of the yaw angle and the relative azimuth deviation; Extract the linear velocity component in the body coordinate system according to the motion data and convert it to the northeast sky coordinate system through the yaw angle, pitch angle, and roll angle. Perform vector decomposition on the converted linear velocity component along the azimuth angle of the main direction and its perpendicular direction to obtain the forward velocity aligned with the estimated sound source direction and the lateral velocity perpendicular to the sound source direction. Take the main direction azimuth angle, forward velocity, lateral velocity, and bias state as the motion state data.

[0011] As a further method, the method for obtaining the first positioning data includes: Perform interpolation alignment at a fixed time interval based on the motion state data of the main direction and the acoustic spectrum time series data. Perform first-order central difference calculation according to the forward velocity and lateral velocity of the motion state data to obtain the velocity change rate. Calculate the absolute value of the azimuth change amount and the standard deviation of the azimuth change amount between adjacent moments according to the main direction azimuth angle sequence of the motion state data; Take the velocity change rate, the absolute value of the azimuth change amount, and the standard deviation of the azimuth change amount as the motion characteristics, take the main frequency band energy, energy distribution entropy, and Doppler frequency shift of the acoustic spectrum time series data as the acoustic spectrum characteristics, and perform normalization processing on the motion characteristics and the acoustic spectrum characteristics. Taking the UAV position and speed as hidden variables, and the motion features and spectrogram features as observation variables, obtaining the conditional probability relationship between the spectrogram features and the UAV position according to the spherical diffusion model, and obtaining the conditional probability relationship between the motion state and the speed according to the UAV dynamics model; Establishing directed edges from hidden variables to observation variables according to the conditional probability relationship, obtaining a directed acyclic graph, constructing a Bayesian network based on the directed acyclic graph, obtaining a particle set of the prior distribution using a probability density heat map, performing particle filtering through the Bayesian network according to the particle set, updating the particle positions according to the motion features, calculating the particle weights based on the spectrogram features, calculating the posterior probability distribution of the UAV position through Monte Carlo integration according to the particle positions and particle weights, and taking the coordinate value with the maximum posterior probability as the first positioning data.

[0012] As a further method, the method for obtaining the lateral positioning data includes: Determining the lateral direction range based on the deviation state of the motion state data, extracting the region of interest of the infrared image data according to the lateral direction range, performing two-layer Daubechies-4 wavelet downsampling on the region of interest, performing wavelet decomposition on the spectrogram time series data in the lateral direction to obtain frequency bands, performing 4-fold downsampling on the frequency bands with a signal-to-noise ratio less than 5 dB, and maintaining the original resolution for the frequency bands greater than or equal to 5 dB; Extracting the histogram of oriented gradient feature vectors from the downsampled infrared image, clustering the histogram of oriented gradient feature vectors, and taking the connected region with a gray value greater than or equal to 1.5 times the average gray value and an area greater than or equal to 10 pixels as an effective heat source to obtain the number of heat sources; Extracting the Mel cepstral coefficient feature vectors from the downsampled spectrogram time series data, and solving for the optimal time offset through the spatio-temporal synchronization constraint equation according to the histogram of oriented gradient feature vectors and Mel cepstral coefficient feature vectors. The spatio-temporal synchronization constraint formula is: where is the optimal time offset, is the histogram of oriented gradient feature vector at time is the speed of sound value under standard atmospheric pressure, is the environmental temperature, is the standard temperature of 273.15 K, is the air density obtained from the pressure data, is the standard density, is the Mel cepstral coefficient feature vector at time is the norm subscript; Normalize, splice the histogram of oriented gradients feature vectors and mel-frequency cepstral coefficients feature vectors after spatio-temporal synchronization, and reduce the dimension to 64 dimensions through principal components to obtain the lateral input features. Take the lateral position offset and the target orientation angle as latent variables, and the lateral input features as observation variables. Input the lateral input features into the Bayesian network, use the probability density heat map to obtain a particle set of the prior distribution, perform particle filtering through the Bayesian network according to the particle set, update the particle positions according to the motion features, calculate the particle weights based on the acoustic spectrum features, calculate the posterior probability distribution of the UAV position through Monte Carlo integration according to the particle positions and particle weights, and take the coordinate value with the maximum posterior probability as the lateral positioning information. Take the number of heat sources and the lateral positioning information as the lateral positioning data.

[0013] As a further method, the method for obtaining the second positioning includes: Based on the lateral positioning data, take the maximum probability value of the posterior probability distribution output by the Bayesian network as the data confidence. If the data confidence is greater than or equal to 0.6 or the number of heat sources in the lateral positioning data is less than or equal to 2, then trigger the bias correction process and use the lateral positioning data to correct the first positioning data; If the bias state extracted in the correction process is above or below, skip the horizontal correction and directly use the height coordinate of the first positioning data as the second positioning height value; If the bias state is left or right, then enter the horizontal correction process, calculate the horizontal correction vector according to the difference between the first positioning data coordinates and the lateral positioning data coordinates, and project the horizontal correction vector onto the vertical direction of the main direction azimuth angle to obtain the lateral offset component; Calculate the spatial distribution entropy based on the integral two-dimensional projection in the vertical axis direction of the main direction and the lateral direction regions in the probability density heat map to obtain the correction weight, and obtain the second positioning by adding the vector dot product of the correction weight and the lateral offset component to the horizontal coordinate of the first positioning data.

[0014] The second aspect of the present invention provides a positioning system for a UAV using a Bayesian network, including: Data acquisition module: used to collect UAV data for preprocessing to obtain acoustic wave data, motion data, infrared image data and environmental data, and the environmental data includes temperature, humidity, and air pressure; Acoustic spectrum data processing module: used to extract outliers based on the acoustic wave data, determine the effective signal according to the outliers, perform time series slicing on the effective signal, and perform spectrum analysis on the slices to obtain acoustic spectrum time series data; First range acquisition module: used to construct a probability density heat map through a sparse variational Gaussian process according to the three-dimensional information of time, frequency, and energy of the acoustic spectrum time series data and the environmental data, and extract the first range based on the probability density heat map; The first positioning data module: It is used to determine the deviation state of the UAV according to the first range and motion data, obtain the main direction and motion state data according to the deviation state, input the spectrogram time series data and motion state data of the main direction into the Bayesian network, and obtain the first positioning data; The second positioning acquisition module: It is used to determine the side direction based on the deviation state of the motion state data, perform adaptive wavelet downsampling on the infrared image data and spectrogram time series data of the side direction, input the downsampled data into the Bayesian network to obtain the lateral positioning data, and use the lateral positioning data to correct the deviation of the first positioning data to obtain the second positioning.

[0015] Compared with the prior art, the embodiments of the present invention have at least the following advantages or beneficial effects: The present invention realizes the effective alignment of infrared HOG features and spectrogram MFCC features through the combination of spatio-temporal synchronization constraints and the correction of the sound speed by the ambient temperature, solves the problem of multi-sensor time asynchrony, improves the time correlation of the fused features, downsamples the spectrogram frequency band, retains the features with high signal-to-noise ratio, and at the same time, through wavelet packet transform and dynamic adjustment of threshold weights, in the environment where the noise power spectral density increases, the effective signal retention rate is improved, and at the same time, the positioning calculation overhead in the side direction is reduced, the positioning real-time performance is improved. Through the positioning of external sound sources, tracking and searching can also be carried out. A three-dimensional probability density heat map is constructed through sparse variational Gaussian processes. Especially in the environment where indoor information such as GNSS is denied, the effective sampling area of the initial particle set of the Bayesian network is reduced and the convergence speed of particle filtering is improved. Through temperature, humidity, and air pressure calibration of the IMU zero bias and sound speed correction, the positioning accumulation error can be effectively reduced in the environment of temperature change and air pressure fluctuation. Description of the Drawings

[0016] Figure 1 It is a step flow chart of a positioning method for a UAV with a Bayesian network in an embodiment of the present invention. Detailed Embodiments

[0017] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0018] Refer to Figure 1 As shown, the present invention provides a positioning method for a UAV with a Bayesian network, including: Collect UAV data for preprocessing to obtain acoustic wave data, motion data, infrared image data, and environmental data. The environmental data includes temperature, humidity, and air pressure; In the actual evaluation, the drone collects acoustic data, motion data, and environmental data through various devices such as a microphone, an MPU6050 IMU sensor, and an infrared camera. The acoustic data includes various audio information such as engine noise signals and environmental noise as acoustic data, with a duration of 10 s and 441,000 sampling points. The motion data includes acceleration, angular velocity, etc. output by the IMU. The infrared camera collects image data with a resolution of 640×480 and a frame rate of 30 fps. The environmental data includes 25°C and 50% RH output by the temperature and humidity sensor, and 101.3 kPa air pressure data output by the air pressure sensor. The above data is time-stamped and calibrated through a hardware synchronization pulse of 1PPS, and the synchronization error is less than or equal to 1 μs. The zero bias of the IMU is compensated by looking up a table using the temperature sensor data. The acceleration zero bias is reduced from ±50 μg to ±10 μg. The contrast of the infrared image is enhanced by logarithmic transformation to achieve dynamic range compression, and the direction features in the image are extracted by the Sobel operator.

[0019] Based on the acoustic data, outliers are extracted, valid signals are determined according to the outliers, the valid signals are sliced in time series, and the sliced data is subjected to spectral analysis to obtain acoustic spectrum time series data; Among them, during the flight of the drone, various specific noise signals are often encountered, which may come from the environment or be the internal noise of the drone. By detecting the energy outliers of the noise background, most natural Gaussian noises can be filtered out, avoiding the interference of engine vibration noise and environmental sudden noise on positioning, improving the signal-to-noise ratio of the valid signal segment, and also judging the distance of the noise source by the level of energy to avoid damage to the drone itself.

[0020] In the actual evaluation, time-domain analysis is performed based on the acoustic data to obtain time series data. According to the time series data, the mean value of the sampling points is calculated as 0.1 V and the standard deviation is 0.02 V through a sliding window with a width of 512 points. The abnormal threshold is set to 3 times the standard deviation. The sampling points that deviate from the mean value by more than the abnormal threshold of 0.06 V are marked as outliers. Based on the outliers, a density clustering algorithm is used to identify continuous abnormal regions. The boundaries of the abnormal regions are extrapolated by 10% along the time axis as the transition region. For example, for points 1000 - 2000, the boundary is extrapolated by 10% to 800 - 2200. The signal intensity in the transition region is linearly interpolated and smoothed. Based on the smoothed time series data, the short-time energy of the acoustic signal is calculated by the time-domain energy detection algorithm. The mean value of the short-time energy in the silent source period and the sum of plus and minus 2 times the standard deviation are used as the background noise energy , and the energy threshold is adaptively determined by the maximum inter-class variance method according to the background noise energy. The time period higher than the energy threshold is used as the valid signal segment. In the actual evaluation, the sound speed is corrected based on the temperature, humidity, and air pressure of the environmental data, where the absolute temperature , the corrected sound speed is , , the relative speed is obtained according to the absolute value of the difference between the corrected sound speed and the speed of the drone in the motion data within the corresponding time , use the ratio of the sound speed to the relative speed as the scaling factor, multiply it by 0.05 seconds as the dynamic length of 0.0508s. Since the dynamic length is less than or equal to the fixed value of 0.2 seconds, the fixed value is used as the slice length. According to the slice length, the effective signal is sliced in time series. The sliced signal is decomposed into 16 frequency bands through wavelet packet transform. For example, the SNR of frequency band 1 is 4dB, which is less than 5dB, so its threshold weight is 1.5. The SNR of frequency band 2 is 10dB, which is greater than 5dB and less than 15dB, so the threshold weight is 1. According to the threshold weight, different frequency bands are denoised through the adaptive threshold formula. For example, for frequency band 1, its static noise reference is 0.01, the air temperature is 25°C, the slice length is 0.2s, and the short-time energy is , the time correlation is 0.8, the calculated denoising threshold is 0.019. The denoised slice is subjected to short-time Fourier transform, set as a Hanning window, with an overlap rate of 75%, to generate a time-frequency spectrogram, whose time resolution is 11.6ms and frequency resolution is 43Hz, and the acoustic spectrum time series data with three-dimensional distribution of time, frequency, and energy is obtained.

[0021] According to the three-dimensional information of time, frequency, and energy of the acoustic spectrum time series data and the environmental data, a probability density heat map is constructed through sparse variational Gaussian process, and the first range is extracted based on the probability density heat map; Among them, the core is to realize the probabilistic expression of the sound source area through Bayesian non-parametric modeling, provide a spatial prior constraint for subsequent positioning, map the three-dimensional acoustic spectrum features of time-frequency-energy to the geographical grid, and compress the high-dimensional data into a three-dimensional space probability distribution through sparse variational Gaussian process. This can effectively reduce the computational complexity and adapt to the real-time processing ability of the drone, so that the heat map can reflect the acoustic wave propagation characteristics under different altitudes and climates. This operation essentially constructs a spatial prior model of data-driven-physical constraint-probability inference, solves the ambiguity problem of sound source localization in complex environments, and provides a reliable spatial reference for multi-modal fusion positioning.

[0022] In the actual evaluation, based on the acoustic spectrum time series data, the three-dimensional features of time, frequency, and energy are extracted, and the Doppler frequency shift offset is calculated according to the corrected sound speed and the corresponding environmental data as , calculate the change gradient of the short-time energy of the sound spectrum with respect to the time stamp according to the time stamp in the time dimension, take the frequency band 2 corresponding to the energy peak in the sound spectrum time series data as the main frequency band, calculate the spatio-temporal distribution entropy of the energy proportion of the main frequency band based on the short-time energy to be 0.7, and generate a spatio-temporal spectrum joint feature matrix based on the change gradient, offset, and spatio-temporal distribution entropy; In the actual evaluation, according to the spatio-temporal spectrum joint feature matrix, the attitude angle of the motion data, and the azimuth angle solution model of the IMU data, map the acoustic features to the geographical grid, where the geographical grid resolution is 1m×1m×1m, with a total of 1000 grid points, obtain the spatial coordinates, map the spatial coordinates, time stamp, and corresponding eigenvalue of the grid to the observation input of the sparse variational Gaussian process, initialize 100 inducing points by the sparse grid method, and use the estimation result of the square of the noise standard deviation of the sound spectrum time series data to determine the diagonal elements of the prior covariance matrix. For example, the variance of frequency band 1 is , solve the posterior distribution through the variational inference formula based on the observation input, take the grid with the posterior distribution probability greater than or equal to 0.7 as the continuous region, perform morphological closing operation on the continuous region with a 3×3 pixel structuring element to obtain the probability density heat map, and extract the three-dimensional minimum bounding box (100m, 200m, 10m), with a side length of 5m as the first range.

[0023] Determine the deviation state of the UAV according to the first range and the motion data, obtain the main direction and motion state data according to the deviation state, input the sound spectrum time series data and motion state data of the main direction into the Bayesian network, and obtain the first positioning data; In the actual evaluation, based on the first range, extract the geometric center coordinates (100m, 200m, 10m) of the three-dimensional minimum bounding box as the estimated sound source direction, extract the synchronized motion data according to the time stamp corresponding to the center coordinates, obtain the current coordinates (0, 0, 0), yaw angle, pitch angle, and roll angle according to the motion data, obtain the relative azimuth deviation of -26.6° according to the horizontal azimuth angle of 63.4° and the pitch angle between the current coordinates and the center coordinates, divide the azimuth angle ±25° and pitch angle ±15° as the front according to the preset thresholds of the relative azimuth deviation and pitch angle, the azimuth angle from -180° to -25° as the left side, the azimuth angle from +25° to +180° as the right side, the pitch angle from +15° to +90° as the upper side, and the pitch angle from -90° to -15° as the lower side as the deviation state. The current deviation state is the left side. Extract the linear velocity components in the body coordinate system according to the motion data and convert them to the northeast celestial coordinate system through the yaw angle, pitch angle, and roll angle. Perform vector decomposition on the converted linear velocity components along the azimuth angle of the main direction and its perpendicular direction to obtain the forward velocity of 4.47m / s aligned with the estimated sound source direction and the lateral velocity of 2.24m / s perpendicular to the sound source direction. Take the azimuth angle of the main direction, forward velocity, lateral velocity, and deviation state as the motion state data.

[0024] In the actual evaluation, interpolation alignment at a fixed time interval is performed based on the motion state data in the main direction and the spectrogram time series data. First-order central difference calculation is carried out according to the forward velocity and lateral velocity of the motion state data to obtain a velocity change rate of 0.1 m / s². According to the main direction azimuth angle sequence of the motion state data, the absolute value of the azimuth angle change between adjacent moments is 5°, and the standard deviation of the azimuth angle change is 2°. In the actual evaluation, the velocity change rate, the absolute value of the azimuth angle change, and the standard deviation of the azimuth angle change are used as motion characteristics, and the main frequency band energy , the energy distribution entropy of 0.6, and the Doppler frequency shift of 2.17 Hz are used as spectrogram characteristics. The motion characteristics and spectrogram characteristics are normalized. The UAV position and velocity are used as latent variables, and the motion characteristics and spectrogram characteristics are used as observed variables. According to the spherical diffusion model, the conditional probability relationship between the spectrogram characteristics and the UAV position is obtained. According to the UAV dynamics model, the conditional probability relationship between the motion state and the velocity is obtained. Directed edges from the latent variables to the observed variables are established according to the conditional probability relationship to obtain a directed acyclic graph. A Bayesian network is constructed based on the directed acyclic graph. 1000 particle sets of the prior distribution are obtained using the probability density heat map. Particle filtering is performed through the Bayesian network according to the particle set. The particle positions are updated according to the motion characteristics, and the particle weights are calculated based on the spectrogram characteristics. The posterior probability distribution of the UAV position is calculated through Monte Carlo integration according to the particle positions and particle weights. After 50 iterations, the coordinate with the maximum posterior probability is (98 m, 195 m, 12 m), which is used as the first positioning data.

[0025] The side direction is determined based on the deviation state of the motion state data. Adaptive wavelet downsampling is performed on the infrared image data and spectrogram time series data in the side direction. The downsampled data is input into the Bayesian network to obtain the side direction positioning data. The first positioning data is corrected for deviation using the side direction positioning data to obtain the second positioning.

[0026] Among them, since the drone has a rigid body coordinate system with itself as the coordinate origin in addition to the globally positioned coordinates during movement, which is used to represent the orientation of the nose and the azimuth of the body equipment, these azimuth information will also affect the sensing ability of the body sensors. Therefore, it is necessary to judge the true position of the drone through data of multiple azimuths and orientations, delimit the side direction range according to the deviation state of the drone itself, and only process the regional data related to the current movement direction, so as to greatly reduce the amount of infrared ROI data and the processing amount of the sound spectrum frequency band, and significantly improve the real-time performance. For example, when the drone flies sideways, only the corresponding side sensors are activated to avoid noise interference in the invalid area, and downsampling is performed on the high-frequency noise area of the sound wave with low signal-to-noise ratio in the side direction data, so as to effectively suppress the noise pollution degree of the sound spectrum MFCC feature, provide a pure input for the Bayesian network, and the side direction positioning data can correct the positioning deviation caused by occlusion in the main direction. For example, when the sound wave in the main direction is interfered by building reflection, the heat source position detected by the side infrared can provide an independent positioning reference, and the positioning reliability in complex maneuvering scenarios is improved through entropy value weighting.

[0027] In actual evaluation, the lateral direction range (45° - 135°) is determined based on the motion state data. The region of interest (ROI) of 320×480 pixels in the left half of the infrared image data is extracted according to the lateral direction range. The ROI is downsampled to 80×120 pixels through two - layer Daubechies - 4 wavelet downsampling. The spectrogram time - series data in the lateral direction is wavelet - decomposed to obtain frequency bands. The frequency band 1 with a signal - to - noise ratio less than 5 dB is downsampled by a factor of 4, and the 50 Hz resolution is retained. The HOG feature vector is extracted from the downsampled infrared image. For example, the HOG feature is a vector with 9 directions, 8×8 cell units, and 128 dimensions. The HOG feature vector is clustered, and the connected region with a gray - scale value greater than or equal to 1.5 times the average gray - scale and an area greater than or equal to 10 pixels is regarded as an effective heat source, and the number of heat sources is obtained as 1. The MFCC feature vector is extracted from the downsampled spectrogram time - series data. For example, the MFCC feature is a vector with 13 dimensions. The optimal time offset is solved through the spatio - temporal synchronization constraint equation based on the histogram of oriented gradients (HOG) feature vector and the Mel - frequency cepstral coefficients (MFCC) feature vector. The spatio - temporally synchronized HOG feature vector and MFCC feature vector are normalized, concatenated, and reduced to 64 dimensions through principal component analysis to obtain the lateral direction input feature. The lateral position offset and the target orientation angle are used as latent variables, and the lateral direction input feature is used as the observed variable. The lateral direction input feature is input into the Bayesian network, and a particle set of the prior distribution is obtained using the probability density heat map. Particle filtering is performed through the Bayesian network based on the particle set. The particle positions are updated according to the motion characteristics, and the particle weights are calculated based on the spectrogram characteristics. The posterior probability distribution of the UAV position is calculated through Monte Carlo integration based on the particle positions and particle weights. The coordinate values with the maximum posterior probability (95m, 205m, 11m) are taken as the lateral positioning information. The number of heat sources and the lateral positioning information are used as the lateral positioning data.

[0028] In the actual evaluation, based on the lateral positioning data, the maximum probability value of the posterior probability distribution output by the Bayesian network is used as the data confidence. Since the data confidence of 0.7 is greater than 0.6, the bias correction process is triggered. Using the lateral positioning data to correct the first positioning data, and the bias state is the left side, then enter the horizontal correction process. Calculate the horizontal correction vector as (95 - 98, 205 - 195) = (-3m, 10m) according to the difference between the coordinates of the first positioning data and the lateral positioning data. Project the horizontal correction vector onto the vertical direction of the main direction azimuth angle = (95 - 98, 205 - 195) = (-3m, 10m) to obtain the lateral offset component (-3×sin63.4° + 10×cos63.4°) = (-2.68 + 4.47) = 1.79m. Calculate the spatial distribution entropy based on the integral two-dimensional projection in the vertical axis direction of the main direction and the side direction regions in the probability density heat map to obtain the correction weight 0.6 / (0.4 + 0.6) = 0.6. According to the vector dot product of the correction weight and the lateral offset component plus the horizontal coordinate of the first positioning data, obtain the second positioning coordinate: (98 + 0.6×(-2.68), 195 + 0.6×4.47, 12) ≈ (96.4m, 197.7m, 12m).

[0029] In this embodiment, the method of the preprocessing includes: Collect the acoustic wave data, motion data, infrared image data, and environmental data of the unmanned aerial vehicle, perform timestamp synchronization through pulses, dynamically compensate the IMU zero bias of the motion data through a temperature sensor; perform dynamic range compression on the infrared image data to distinguish direction features; calibrate the temperature, humidity, and air pressure values of the environmental data through redundant sensor data.

[0030] In this embodiment, the method of obtaining the effective signal includes: Perform time-domain analysis on the acoustic wave data to obtain time-series data. Calculate the mean and standard deviation of the sampling points through a sliding window according to the time-series data. Set the anomaly threshold as 3 times the standard deviation. Mark the sampling points that deviate from the mean by more than the anomaly threshold as outliers. Identify continuous anomaly regions based on the outliers using a density clustering algorithm. Extrapolate the boundary of the anomaly region by 10% along the time axis as the transition region, and perform linear interpolation smoothing on the signal intensity within the transition region; Calculate the short-time energy of the acoustic wave signal based on the smoothed time-series data through a time-domain energy detection algorithm. Take the mean of the short-time energy in the silent source period plus and minus 2 times the standard deviation as the background noise energy, Adaptively determine the energy threshold according to the background noise energy through the maximum inter-class variance method, and take the period higher than the energy threshold as the effective signal segment.

[0031] In this embodiment, the method of obtaining the acoustic spectrum time-series data includes: Based on environmental data, correct the sound speed according to temperature, humidity and air pressure. Obtain the relative speed based on the absolute value of the difference between the corrected sound speed and the UAV speed in the motion data within the corresponding time. Use the ratio of the sound speed to the relative speed as the scaling factor, multiply it by 0.05 seconds to obtain the dynamic length. If the dynamic length is less than or equal to the fixed value of 0.2 seconds, use the fixed value as the slice length. If the dynamic length is greater than the fixed value of 0.2 seconds, use the dynamic length as the slice length. Perform time series slicing on the effective signal according to the slice length; Decompose the sliced signal into different frequency bands through wavelet packet transform, extract the signal-to-noise ratio based on different frequency bands. If the signal-to-noise ratio is less than 5 dB, expand the threshold weight to 1.5. If the signal-to-noise ratio is greater than or equal to 5 dB and less than 15 dB, keep the threshold weight as 1. If the signal-to-noise ratio is greater than or equal to 15 dB, reduce the threshold weight to 0.8; Perform noise reduction on different frequency bands according to the threshold weight through the adaptive threshold formula. The adaptive threshold formula is: where is the noise reduction threshold of the wavelet coefficient, is the static noise benchmark, is the air temperature, is the slice length, is the threshold weight, is the short-time energy of the current slice, is the background noise energy, which is refreshed every 10 seconds, is the time correlation degree, which is the Pearson correlation coefficient of 10 consecutive sampling points; Calculate the time-frequency spectrum diagram of the noise-reduced slice through short-time Fourier transform, use the Hanning window as the window function, set the window function overlap rate to 75%, and obtain the acoustic spectrum time series data with three-dimensional distribution of time, frequency and energy.

[0032] In this embodiment, the method for obtaining the first range includes: Extract three-dimensional features of time, frequency and energy from the acoustic spectrum time series data, calculate the Doppler frequency shift offset according to the corrected sound speed and the corresponding environmental data, calculate the change gradient of the short-time energy of the acoustic spectrum with respect to the time stamp in the time dimension, use the frequency band corresponding to the energy peak in the acoustic spectrum time series data as the main frequency band, calculate the spatio-temporal distribution entropy of the energy proportion of the main frequency band based on the short-time energy, and generate a spatio-temporal spectrum joint feature matrix based on the change gradient, offset and spatio-temporal distribution entropy; According to the spatio-temporal spectrum joint feature matrix, the attitude angle of the motion data, and the IMU data, use the azimuth angle resolution model to map the acoustic features to the geographical grid, obtain the spatial coordinates, and map the spatial coordinates, time stamp and corresponding eigenvalue of the grid as the observation input of the sparse variational Gaussian process; Solving the posterior distribution through the variational inference formula based on the observed input, and the variational inference formula is: where is the observed data of the log marginal likelihood, is the th time-frequency feature vector of the data point, is the slice index, is the total number of slices, is the th potential function value of the data point, that is, the true sound source state, is the potential function of the Gaussian distribution, is the induced variable, obtained by initializing through the sparse grid method, is the induced variable of the variational distribution, is the induced variable of the Gaussian prior distribution, and the covariance matrix of the prior distribution is determined by the estimated result of the noise standard deviation of the acoustic spectrum time series data, is the th environmental temperature of the slice, is the th corrected sound speed of the slice, is the preset regularization coefficient; Taking the grid with the posterior distribution probability greater than or equal to 0.7 as the continuous region, performing morphological closing operation on the continuous region through a 3×3 pixel structural element, obtaining the probability density heat map, and extracting the three-dimensional minimum circumscribed bounding box as the first range.

[0033] In this embodiment, the method for obtaining the main direction and the motion state data according to the deviation state includes: Extracting the geometric center coordinates of the three-dimensional minimum circumscribed bounding box based on the first range as the estimated sound source direction, extracting the synchronized motion data according to the time stamp corresponding to the center coordinates, obtaining the current coordinates, yaw angle, pitch angle, and roll angle according to the motion data, and obtaining the relative azimuth deviation according to the horizontal azimuth angle and pitch angle between the current coordinates and the center coordinates; Dividing the front, left, right, above, and below as the deviation states according to the preset thresholds of the relative azimuth deviation and the pitch angle. When the deviation state is the front, the main direction is the current yaw angle of the UAV. When the deviation state is the left or right, the main direction is the sum of the yaw angle and the relative azimuth deviation; Extract the linear velocity component in the body coordinate system based on the motion data and convert it to the northeast celestial coordinate system through the yaw angle, pitch angle, and roll angle. Perform vector decomposition on the converted linear velocity component along the azimuth angle of the main direction and its perpendicular direction to obtain the forward velocity aligned with the estimated sound source direction and the lateral velocity perpendicular to the sound source direction. Take the azimuth angle of the main direction, forward velocity, lateral velocity, and deviation state as the motion state data.

[0034] In this embodiment, the method for obtaining the first positioning data includes: Perform interpolation alignment at a fixed time interval based on the motion state data of the main direction and the spectrogram time series data. Perform first-order central difference calculation based on the forward velocity and lateral velocity of the motion state data to obtain the velocity change rate. Calculate the absolute value of the azimuth angle change and the standard deviation of the azimuth angle change between adjacent moments based on the azimuth angle sequence of the motion state data; Take the velocity change rate, absolute value of the azimuth angle change, and standard deviation of the azimuth angle change as motion features, take the main frequency band energy, energy distribution entropy, and Doppler frequency shift of the spectrogram time series data as spectrogram features, and perform normalization processing on the motion features and spectrogram features; Take the UAV position and velocity as latent variables, take the motion features and spectrogram features as observed variables, obtain the conditional probability relationship between the spectrogram features and the UAV position according to the spherical diffusion model, and obtain the conditional probability relationship between the motion state and the velocity according to the UAV dynamics model; Establish a directed edge from the latent variable to the observed variable according to the conditional probability relationship to obtain a directed acyclic graph. Construct a Bayesian network based on the directed acyclic graph. Use the probability density heat map to obtain a particle set of the prior distribution. Perform particle filtering through the Bayesian network according to the particle set. Update the particle position according to the motion features, calculate the particle weights based on the spectrogram features, calculate the posterior probability distribution of the UAV position through Monte Carlo integration according to the particle position and particle weights, and take the coordinate value with the maximum posterior probability as the first positioning data.

[0035] In this embodiment, the method for obtaining the lateral positioning data includes: Determine the lateral direction range based on the deviation state of the motion state data. Extract the region of interest of the infrared image data according to the lateral direction range. Perform two-layer Daubechies-4 wavelet downsampling on the region of interest. Perform wavelet decomposition on the spectrogram time series data in the lateral direction to obtain frequency bands. Perform 4-fold downsampling on the frequency bands with a signal-to-noise ratio less than 5 dB, and keep the original resolution for the frequency bands greater than or equal to 5 dB; Extract the histogram of oriented gradients feature vector from the downsampled infrared image, cluster the histogram of oriented gradients feature vector, and take the connected region with a gray value greater than or equal to 1.5 times the average gray value and an area greater than or equal to 10 pixels as the effective heat source to obtain the number of heat sources; Extract the Mel cepstral coefficient feature vectors based on the downsampled spectrogram time series data. Solve for the optimal time offset according to the histogram of oriented gradients (HOG) feature vectors and the Mel cepstral coefficient feature vectors through the spatio-temporal synchronization constraint equation. The spatio-temporal synchronization constraint formula is as follows: where is the optimal time offset, is the HOG feature vector at time is the speed of sound value under standard atmospheric pressure, is the ambient temperature, is the standard temperature of 273.15 K, is the air density obtained from the pressure data, is the standard density, is the Mel cepstral coefficient feature vector at time is the norm subscript; Normalize, splice the spatio-temporally synchronized HOG feature vectors and Mel cepstral coefficient feature vectors, and reduce the dimension to 64 through principal component analysis to obtain the side direction input features. Use the lateral position offset and the target orientation angle as latent variables, and the side direction input features as observed variables. Input the side direction input features into the Bayesian network, use the probability density heat map to obtain the particle set of the prior distribution, perform particle filtering according to the particle set through the Bayesian network, update the particle positions according to the motion characteristics, calculate the particle weights based on the spectral features, calculate the posterior probability distribution of the UAV position through Monte Carlo integration according to the particle positions and particle weights, and take the coordinate value with the maximum posterior probability as the lateral positioning information. Use the number of heat sources and the lateral positioning information as the lateral positioning data.

[0036] In this embodiment, the method for obtaining the second positioning includes: Based on the lateral positioning data, use the maximum probability value of the posterior probability distribution output by the Bayesian network as the data confidence. If the data confidence is greater than or equal to 0.6 or the number of heat sources in the lateral positioning data is less than or equal to 2, trigger the bias correction process and use the lateral positioning data to correct the first positioning data; If the bias state extracted in the correction process is above or below, skip the horizontal correction and directly use the height coordinate of the first positioning data as the second positioning height value; If the bias state is left or right, enter the horizontal correction process, calculate the horizontal correction vector according to the difference between the first positioning data coordinates and the lateral positioning data coordinates, and project the horizontal correction vector onto the vertical direction of the main direction azimuth angle to obtain the lateral offset component; The spatial distribution entropy is calculated based on the integral two-dimensional projection in the vertical axis direction of the main direction and the side direction regions in the probability density heat map to obtain a corrected weight. The second positioning is obtained by adding the vector dot product of the corrected weight and the lateral offset component to the horizontal coordinate of the first positioning data. The second aspect of the present invention also provides a positioning system for an unmanned aerial vehicle with a Bayesian network, including: Data acquisition module: used to collect unmanned aerial vehicle data for preprocessing to obtain acoustic wave data, motion data, infrared image data, and environmental data, where the environmental data includes temperature, humidity, and air pressure; Acoustic spectrum data processing module: used to extract outliers based on the acoustic wave data, determine effective signals according to the outliers, perform time series slicing on the effective signals, and perform spectrum analysis on the slices to obtain acoustic spectrum time series data; First range acquisition module: used to construct a probability density heat map through a sparse variational Gaussian process based on the three-dimensional information of time, frequency, and energy of the acoustic spectrum time series data and the environmental data, and extract the first range based on the probability density heat map; First positioning data module: used to determine the deviation state of the unmanned aerial vehicle according to the first range and the motion data, obtain the main direction and motion state data according to the deviation state, input the acoustic spectrum time series data and the motion state data of the main direction into the Bayesian network, and obtain the first positioning data; Second positioning acquisition module: used to determine the side direction according to the deviation state of the motion state data, perform adaptive wavelet downsampling on the infrared image data and the acoustic spectrum time series data of the side direction, input the downsampled data into the Bayesian network to obtain lateral positioning data, and use the lateral positioning data to correct the deviation of the first positioning data to obtain the second positioning.

[0037] The above content is only an example and explanation of the structure of the present invention. Those skilled in the art of this technology can make various modifications or supplements to the described specific embodiments or use similar methods to replace them. As long as they do not deviate from the structure of the invention or exceed the scope defined by this claim book, they should fall within the protection scope of the present invention.

Claims

1. A positioning method for an unmanned aerial vehicle of a Bayesian network, characterized in that, The method includes the following steps: Collect drone data for preprocessing to obtain acoustic wave data, motion data, infrared image data, and environmental data. The environmental data includes temperature, humidity, and air pressure; Extract outliers based on the acoustic wave data, determine effective signals according to the outliers, perform time series slicing on the effective signals, and perform spectral analysis on the slices to obtain acoustic spectrum time series data; Construct a probability density heat map through a sparse variational Gaussian process based on the three-dimensional information of time, frequency, and energy of the acoustic spectrum time series data and the environmental data, and extract a first range based on the probability density heat map; Determine the deviation state of the drone according to the first range and the motion data, obtain the main direction and motion state data according to the deviation state, input the acoustic spectrum time series data and motion state data of the main direction into a Bayesian network, and obtain first positioning data; Determine the side direction based on the deviation state of the motion state data, perform adaptive wavelet downsampling on the infrared image data and acoustic spectrum time series data of the side direction, input the downsampled data into a Bayesian network to obtain side positioning data, and use the side positioning data to correct the deviation of the first positioning data to obtain a second positioning.

2. The positioning method of the UAV with a Bayesian network according to claim 1, characterized in that, The method for the preprocessing includes: Collect the acoustic wave data, motion data, infrared image data, and environmental data of the drone, perform timestamp synchronization through pulses, dynamically compensate the IMU zero bias of the motion data through a temperature sensor; perform dynamic range compression on the infrared image data to distinguish direction features; calibrate the temperature, humidity, and air pressure values of the environmental data through redundant sensor data.

3. The positioning method of the UAV with a Bayesian network according to claim 1, characterized in that The method for obtaining the effective signal includes: Perform time domain analysis on the acoustic wave data to obtain time series data, calculate the mean and standard deviation of the sampling points through a sliding window according to the time series data, set the anomaly threshold to 3 times the standard deviation, mark the sampling points that deviate from the mean by more than the anomaly threshold as outliers, identify continuous anomaly regions based on the outliers using a density clustering algorithm, extrapolate the boundaries of the anomaly regions by 10% along the time axis as a transition zone, and perform linear interpolation smoothing on the signal intensity within the transition zone; Calculate the short-time energy of the acoustic wave signal based on the smoothed time series data through a time domain energy detection algorithm, and use the sum of the short-time energy mean of the silent source period and plus or minus 2 times the standard deviation as the background noise energy; Adaptively determine the energy threshold according to the background noise energy through the maximum inter-class variance method, and use the period higher than the energy threshold as the effective signal segment.

4. A positioning method for an unmanned aerial vehicle of a Bayesian network according to claim 1, characterized in that, The method for obtaining the acoustic spectrum time series data includes: Correct the speed of sound based on the temperature, humidity, and air pressure of the environmental data, obtain the relative speed according to the absolute value of the difference between the corrected speed of sound and the speed of the drone in the motion data within the corresponding time, use the ratio of the speed of sound to the relative speed as a scaling factor multiplied by 0.05 seconds as the dynamic length. If the dynamic length is less than or equal to the fixed value of 0.2 seconds, use the fixed value as the slice length. If the dynamic length is greater than the fixed value of 0.2 seconds, use the dynamic length as the slice length, and perform time series slicing on the effective signal according to the slice length; The signal of the slice is decomposed into different frequency bands through wavelet packet transform. The signal-to-noise ratio is extracted based on different frequency bands. When the signal-to-noise ratio is less than 5 decibels, the threshold weight is expanded to 1.

5. When the signal-to-noise ratio is greater than or equal to 5 decibels and less than 15 decibels, the threshold weight is maintained at 1. When the signal-to-noise ratio is greater than or equal to 15 decibels, the threshold weight is reduced to 0.8; Noise reduction is performed on different frequency bands according to the threshold weight through the adaptive threshold formula. The adaptive threshold formula is: ； where is the noise reduction threshold of the wavelet coefficient, is the static noise benchmark, is the air temperature, is the slice length, is the threshold weight, is the short-time energy of the current slice, is the background noise energy, refreshed every 10 seconds, is the time correlation degree, which is the Pearson correlation coefficient of 10 consecutive sampling points; The time-frequency spectrogram is calculated for the denoised slice through short-time Fourier transform. The Hanning window is used as the window function, and the window function overlap rate is set to 75% to obtain the acoustic spectrum time series data with three-dimensional distribution of time, frequency, and energy.

5. A positioning method for an unmanned aerial vehicle of a Bayesian network according to claim 1, characterized in that, The method for obtaining the first range includes: Extract the three-dimensional features of time, frequency, and energy based on the acoustic spectrum time series data. Calculate the Doppler frequency shift offset according to the corrected sound speed and the corresponding environmental data. Calculate the change gradient of the short-time energy of the acoustic spectrum with respect to the time stamp in the time dimension. The frequency band corresponding to the energy peak in the acoustic spectrum time series data is used as the main frequency band. Calculate the spatio-temporal distribution entropy of the energy proportion of the main frequency band based on the short-time energy. Generate a spatio-temporal spectrum joint feature matrix based on the change gradient, offset, and spatio-temporal distribution entropy; Map the acoustic features to the geographical grid according to the spatio-temporal spectrum joint feature matrix, the attitude angle of the motion data, and the IMU data using the azimuth resolution model to obtain the spatial coordinates. Map the spatial coordinates, time stamp, and corresponding eigenvalue of the grid to the observation input of the sparse variational Gaussian process; Solve the posterior distribution based on the observation input through the variational inference formula. The variational inference formula is: ； where is the log marginal likelihood of the observed data , is the time-frequency feature vector of the -th data point is the slice index is the total number of slices is the value of the potential function of the -th data point, i.e., the true sound source state is the potential function with a Gaussian distribution are the induced variables, obtained by initializing through the sparse grid method is the induced variable with a variational distribution is the induced variable with a Gaussian prior distribution, and the covariance matrix of the prior distribution is determined by the estimation result of the noise standard deviation of the acoustic spectrum time series data is the environmental temperature of the -th slice is the corrected sound speed of the -th slice is the preset regularization coefficient; The grid with the posterior distribution probability greater than or equal to 0.7 is used as the continuous region. Perform morphological closing operation on the continuous region through a 3×3 pixel structural element to obtain the probability density heat map, and extract the three-dimensional minimum bounding box as the first range.

6. A positioning method for an unmanned aerial vehicle of a Bayesian network according to claim 1, characterized in that, The method for obtaining the main direction and motion state data according to the deviation state includes: Extract the geometric center coordinates of the three-dimensional minimum bounding box based on the first range as the estimated sound source direction. Extract the synchronized motion data according to the time stamp corresponding to the center coordinates. Obtain the current coordinates, yaw angle, pitch angle, and roll angle according to the motion data. Obtain the relative azimuth deviation according to the horizontal azimuth angle and pitch angle between the current coordinates and the center coordinates; Divide the front, left, right, up, and down as the deviation states according to the preset thresholds of the relative azimuth deviation and pitch angle. When the deviation state is the front, the main direction is the current yaw angle of the unmanned aerial vehicle. When the deviation state is the left or right, the main direction is the sum of the yaw angle and the relative azimuth deviation; Extract the linear velocity components in the body coordinate system according to the motion data and convert them to the northeast sky coordinate system through the yaw angle, pitch angle, and roll angle. Perform vector decomposition on the converted linear velocity components along the azimuth angle of the main direction and its perpendicular direction to obtain the forward velocity aligned with the estimated sound source direction and the lateral velocity perpendicular to the sound source direction. The azimuth angle of the main direction, forward velocity, lateral velocity, and deviation state are used as the motion state data.

7. A positioning method for an unmanned aerial vehicle of a Bayesian network according to claim 1, characterized in that The method for obtaining the first positioning data includes: Interpolate and align the motion state data and spectrogram time series data based on the main direction at fixed time intervals. Calculate the first-order central difference of the forward velocity and lateral velocity of the motion state data to obtain the velocity change rate. Calculate the absolute value of the azimuth change and the standard deviation of the azimuth change between adjacent moments based on the main direction azimuth angle sequence of the motion state data; Take the velocity change rate, the absolute value of the azimuth change, and the standard deviation of the azimuth change as motion features, and take the main frequency band energy, energy distribution entropy, and Doppler frequency shift of the spectrogram time series data as spectrogram features. Normalize the motion features and spectrogram features; Take the UAV position and velocity as latent variables, and take the motion features and spectrogram features as observed variables. Obtain the conditional probability relationship between the spectrogram features and the UAV position according to the spherical diffusion model, and obtain the conditional probability relationship between the motion state and the velocity according to the UAV dynamics model; Establish a directed edge from the latent variable to the observed variable according to the conditional probability relationship to obtain a directed acyclic graph. Construct a Bayesian network based on the directed acyclic graph. Use the probability density heat map to obtain a particle set of the prior distribution. Perform particle filtering through the Bayesian network according to the particle set. Update the particle position according to the motion features, calculate the particle weights based on the spectrogram features, calculate the posterior probability distribution of the UAV position through Monte Carlo integration according to the particle position and particle weights, and take the coordinate value with the maximum posterior probability as the first positioning data.

8. A positioning method for an unmanned aerial vehicle of a Bayesian network according to claim 1, characterized in that, The method for obtaining the lateral positioning data includes: Determine the lateral direction range based on the deviation state of the motion state data. Extract the region of interest of the infrared image data according to the lateral direction range. Perform two-layer Daubechies-4 wavelet downsampling on the region of interest. Decompose the spectrogram time series data in the lateral direction to obtain frequency bands. Downsample the frequency bands with a signal-to-noise ratio less than 5 dB by a factor of 4, and keep the original resolution for the frequency bands greater than or equal to 5 dB; Extract the histogram of oriented gradients feature vector from the downsampled infrared image. Cluster the histogram of oriented gradients feature vector. Take the connected region with a gray value greater than or equal to 1.5 times the average gray value and an area greater than or equal to 10 pixels as the effective heat source, and obtain the number of heat sources; Extract the mel-frequency cepstral coefficients feature vector from the downsampled spectrogram time series data. Solve the optimal time offset according to the histogram of oriented gradients feature vector and the mel-frequency cepstral coefficients feature vector through the spatio-temporal synchronization constraint equation. The spatio-temporal synchronization constraint formula is: ； wherein is the optimal time offset; is the histogram of oriented gradients feature vector at time the speed of sound value under standard atmospheric pressure; is the ambient temperature; is the standard temperature 273.15K; is the air density obtained from the pressure data; is the standard density; is the mel-frequency cepstral coefficients feature vector at time is the norm subscript; Normalize, splice the histogram of oriented gradients feature vectors and mel-frequency cepstral coefficients feature vectors after spatio-temporal synchronization, and reduce the dimension to 64 dimensions through principal component analysis to obtain the side-direction input features. Use the lateral position offset and the target orientation angle as latent variables, and the side-direction input features as observation variables. Input the side-direction input features into the Bayesian network, use the probability density heat map to obtain a particle set of the prior distribution, perform particle filtering through the Bayesian network according to the particle set, update the particle positions according to the motion features, calculate the particle weights based on the acoustic spectrum features, calculate the posterior probability distribution of the UAV position through Monte Carlo integration according to the particle positions and particle weights, and take the coordinate value with the maximum posterior probability as the lateral positioning information. Use the number of heat sources and the lateral positioning information as the lateral positioning data.

9. A positioning method for an unmanned aerial vehicle of a Bayesian network according to claim 1, characterized in that, The method for obtaining the second positioning includes: Based on the lateral positioning data, use the maximum probability value of the posterior probability distribution output by the Bayesian network as the data confidence. If the data confidence is greater than or equal to 0.6 or the number of heat sources in the lateral positioning data is less than or equal to 2, then trigger the bias correction process and use the lateral positioning data to correct the first positioning data; If the bias state extracted in the correction process is above or below, skip the horizontal correction and directly use the height coordinate of the first positioning data as the second positioning height value; If the bias state is left or right, then enter the horizontal correction process. Calculate the horizontal correction vector according to the difference between the coordinates of the first positioning data and the coordinates of the lateral positioning data, and project the horizontal correction vector onto the vertical direction of the main direction azimuth angle to obtain the lateral offset component; Calculate the spatial distribution entropy based on the integral two-dimensional projection in the vertical axis direction of the main direction and side-direction regions in the probability density heat map to obtain the correction weight. Obtain the second positioning by adding the vector dot product of the correction weight and the lateral offset component to the horizontal coordinate of the first positioning data.

10. A positioning system for an unmanned aerial vehicle of a Bayesian network, which is used to execute the positioning method of the unmanned aerial vehicle of a Bayesian network according to any one of claims 1 to 9, characterized in that, The system includes: Data acquisition module: used to collect UAV data for preprocessing to obtain acoustic wave data, motion data, infrared image data, and environmental data. The environmental data includes temperature, humidity, and air pressure; Acoustic spectrum data processing module: used to extract outliers based on the acoustic wave data, determine the effective signal according to the outliers, perform time series slicing on the effective signal, and perform spectral analysis on the slices to obtain acoustic spectrum time series data; First range acquisition module: used to construct a probability density heat map through a sparse variational Gaussian process according to the time, frequency, and energy three-dimensional information of the acoustic spectrum time series data and the environmental data, and extract the first range based on the probability density heat map; First positioning data module: used to determine the UAV bias state according to the first range and motion data, obtain the main direction and motion state data according to the bias state, input the acoustic spectrum time series data and motion state data of the main direction into the Bayesian network, and obtain the first positioning data; Second positioning acquisition module: used to determine the side direction according to the bias state of the motion state data, perform adaptive wavelet downsampling on the infrared image data and acoustic spectrum time series data of the side direction, input the downsampled data into the Bayesian network to obtain the lateral positioning data, and use the lateral positioning data to perform bias correction on the first positioning data to obtain the second positioning.

Images