A method for quickly identifying geological disaster hidden points based on unmanned aerial vehicle aerial survey

Abstract

The application discloses a kind of geological disaster hidden danger point fast identification method based on unmanned aerial vehicle aerial survey, it is related to hidden danger point fast identification, by collecting airborne camera image sequence, acceleration and angular velocity of inertial measurement unit, propeller speed, flight height, wind speed and direction and body vibration signal, with image exposure time as unified time axis, form time series sample, establish the momentum entrainment priori of rotor downwash near ground boundary layer based on time series sample, calculate downwash dynamic pressure field time series estimation, particle mobilization probability time series and equivalent micro-attitude disturbance spectrum, sub-pixel level optical flow and phase accumulation analysis are carried out to image sequence, time cross-correlation is carried out with equivalent micro-attitude disturbance spectrum, coupling evidence is obtained, the joint determination of trigger intensity and coupling degree is carried out, to trigger flight parameter adaptive fine-tuning strategy, spatial risk analysis is carried out to the image after compensation reconstruction, to generate spatial risk calibration chart, can effectively identify the geological disaster hidden danger point in critical stable state.

Description

A rapid identification method for geological hazard hazard points based on UAV aerial surveying Technical Field

[0001] This invention relates to the field of uniformity detection and control, specifically a method for rapid identification of geological hazard points based on UAV aerial surveying. Background Technology

[0002] When conducting low-altitude, high-resolution aerial surveys using UAVs at potential sites of loose geological hazards such as landslide deposits and debris slopes, the strong downwash airflow generated by the rotor interacts significantly with the ground surface aerodynamically. This interaction manifests not only as a simple aerodynamic lift effect but also involves a complex rotor downwash momentum entrainment process in the near-surface boundary layer. When the UAV passes over the surface of loose deposits, the dynamic pressure of the downwash field often exceeds the critical wind erosion threshold of fine-grained materials on the surface, leading to resuspension of surface particles or the ejection of tiny particles. This transient physical displacement of the surface medium disrupts the temporal consistency of image texture, resulting in non-rigid local undulations in the image. The rebounding airflow and the lifted particles impact the aircraft and gimbal, inducing high-frequency mechanical vibrations. These vibrations, coupled with rotor aerodynamic disturbances, form a specific equivalent micro-attitude disturbance spectrum, causing micro-flutter in the imaging optical axis that is difficult to completely eliminate using conventional mechanical stabilization.

[0003] Existing UAV-based geological hazard identification technologies typically treat the land surface as a static rigid body and the UAV as a passive observer, neglecting the fluid-solid-aerodynamic coupling effect between the rotor downwash and the loose surface medium. When processing data acquired in such scenarios, traditional methods often treat surface fine-particle disturbances induced by the downwash as image noise or registration errors for smoothing, and filter out the aerodynamic-mechanical coupling vibrations of the UAV as simple flight instability. Current technologies lack the ability to quantify the causal relationship between the downwash dynamic pressure field and the probability of particle mobilization, and cannot utilize the high-frequency components of UAV vibration and the time-frequency coupling characteristics of local image undulations to inversely deduce the degree of surface looseness. Therefore, it is difficult to effectively identify geological hazard hazard points in a critically stable state in the early stages when visual features are not obvious.

[0004] To address the aforementioned shortcomings, a technical solution is provided. Summary of the Invention

[0005] To address the technical problems mentioned in the background section, this invention is proposed. Embodiments of this invention provide a method for rapid identification of geological hazard hazard points based on unmanned aerial vehicle (UAV) aerial surveying.

[0006] The objective of this invention can be achieved through the following technical solution: a method for rapid identification of geological hazard hazard points based on UAV aerial surveying, comprising:

[0007] The system collects image sequences from airborne cameras, acceleration and angular velocity from inertial measurement units, propeller speed, flight altitude, wind speed and direction, and aircraft vibration signals, forming time-series samples with image exposure time as the unified time axis.

[0008] Based on time-series samples, a priori on the momentum entrainment of the near-ground boundary layer during rotor downwash is established. Using propeller speed, flight altitude and wind field data, the time-series estimation of the downwash dynamic pressure field, the time series of particle mobilization probability and the equivalent micro-attitude perturbation spectrum are calculated.

[0009] Subpixel-level optical flow and phase accumulation analysis were performed on the image sequence. Blind source separation was performed on the high-frequency components of body vibration and inertial measurement. Time cross-correlation was performed with the equivalent micro-attitude perturbation spectrum to obtain coupling evidence.

[0010] Based on the time series estimation of the downwash dynamic pressure field, the joint determination of the triggering strength and coupling degree of the particle mobilization probability time series, a flight parameter adaptive fine-tuning strategy is triggered.

[0011] In the triggering section, image plane motion compensation and rolling shutter micro-distortion correction are performed on the image, and the compensated and reconstructed image is output.

[0012] Spatial risk analysis is performed on the compensated and reconstructed imagery to generate a spatial risk calibration map.

[0013] Furthermore, the steps for calculating the timing estimate of the sludge dynamic pressure field are as follows:

[0014] A priori on the momentum entrainment of the near-ground boundary layer during rotor downwash is established. An axisymmetric-shear coupled dynamic pressure field parameterization is constructed based on propeller speed, flight altitude, and wind speed and direction. The time series estimation of the downwash dynamic pressure field is solved by energy conservation and momentum flux closure.

[0015] Furthermore, the calculation steps for the equivalent micro-attitude perturbation spectrum are as follows:

[0016] By establishing a station-grid mapping, surface roughness and zero displacement height parameters are extracted. Combined with wind direction, wind speed and surface zoning information, critical wind erosion and shear thresholds for each particle size are calculated. Uncertainty correction is considered to obtain the particle mobilization probability time series.

[0017] The equivalent torque perturbation is obtained by multiplying the timing estimate of the downwash dynamic pressure field with the transfer matrix of the body / gimbal. Combined with the maneuvering state of the inertial measurement unit, the equivalent micro-attitude perturbation spectrum is solved.

[0018] Furthermore, the calculation steps for the particle mobilization probability time series are as follows:

[0019] By combining the partition weight vector and partition label of each pixel, as well as the surface roughness parameter and zero displacement height parameter, and the wind direction and near-ground wind information that change over time, the partition threshold parameters are weighted to obtain the critical wind erosion threshold and equivalent critical shear threshold for each particle size. Based on turbulent gusts, surface parameter uncertainty and partition heterogeneity, a comprehensive uncertainty scale is generated.

[0020] Based on the critical wind erosion threshold, equivalent critical shear threshold and comprehensive uncertainty scale of each particle size, combined with time-varying frictional wind information, the mobilization probability of each particle size is calculated, and weighted according to the local particle size distribution to obtain the time series of particle mobilization probability.

[0021] Furthermore, the steps for obtaining the partition weight vector and partition label of each pixel are as follows:

[0022] Using the main line-of-sight projection point as the station, reference survey lines are laid out along the main wind direction and side, and the near-ground wind information is unified to the same height benchmark. Combined with the surface type zoning and micro-topographic undulation information, robust fitting and consistency tests are performed on the candidate roughness parameters. The surface roughness parameters, zero displacement height parameters and their uncertainty range corresponding to the spatial location are output, and a station-grid mapping index is established.

[0023] Based on the station-grid mapping index, image grid cells are mapped to surface patches. For each cell, the visible area ratio of each patch is calculated. The ratio is corrected using attitude and spectral consistency to form the partition weight vector and partition label of each cell.

[0024] Furthermore, the steps for obtaining the coupling evidence are as follows:

[0025] Subpixel-level optical flow and phase accumulation analysis were performed on the image sequence to extract the multi-scale vector field and corresponding phase energy density sequence of local motion in the image, and to establish a dynamic observable set of image fluctuations.

[0026] Blind source separation is performed on the vibration signal of the machine body and the high-frequency components of the inertial measurement unit to extract the main vibration mode and noise basis of the structure and form a characterization of the machine body's natural vibration response;

[0027] The dynamic observable set of image fluctuations and the characterization of the body's natural vibration response are cross-correlated and coherently tested with the equivalent micro-attitude perturbation spectrum on a unified time axis to obtain a comprehensive coupling coefficient. Based on the comprehensive coupling coefficient, it is determined that there is significant multi-scale time-frequency coupling, and the coupling evidence is output.

[0028] Furthermore, the steps for triggering the adaptive fine-tuning strategy for flight parameters are as follows:

[0029] By combining frame-level comprehensive triggering probability, coupling evidence, and segment-level triggering results, the triggering intensity and coupling degree of each surface segment are jointly determined to identify high-risk segments.

[0030] The lateral micro-misalignment and micro-elevation commands are calculated based on high-risk sections, and the upper limits of pitch and roll change rates are limited to reduce downwash dynamic pressure and micro-attitude disturbances.

[0031] Furthermore, the steps for obtaining the segment-level triggering result are as follows:

[0032] The comprehensive coupling coefficient is used as the weight of the posterior probability distribution of particle-induced fluctuations in each frame to form a weighted coherence index.

[0033] Using the frame-by-frame particle mobilization probability map, the posterior probability distribution of particle-induced fluctuations in each frame, and the weighted coherence index as ternary evidence, Bayesian updates are used to obtain the frame-level comprehensive triggering probability, and time smoothing is performed within a sliding window to generate segment-level triggering results.

[0034] Furthermore, the steps for obtaining the posterior probability distribution of particle-induced fluctuations in each frame are as follows:

[0035] Based on the time series estimation of the downwash dynamic pressure field, the time series of particle mobilization probability is resampled and spatially interpolated to generate a frame-by-frame particle mobilization probability map, and image-to-ground registration is completed with image coordinates.

[0036] The dynamic observable set of image fluctuations is mapped into a fine-grained disturbance intensity distribution map through line-of-sight projection and ground-shadow geometry. Using the frame-by-frame particle mobilization probability map as a prior, an observation likelihood function is constructed and a posterior inference is performed to obtain the posterior probability distribution of particle-induced fluctuations in each frame.

[0037] Furthermore, the steps for outputting the compensated and reconstructed image are as follows:

[0038] Reacquire flight images of the trigger zone for flight control adjustment, and perform adaptive deflickering and brightness normalization on the acquired images to obtain stable images;

[0039] Using the equivalent micro-attitude perturbation spectrum and the main vibration modes of the structure as priors, the equivalent micro-attitude time series is solved and the image plane motion compensation field is generated.

[0040] Based on the image plane motion compensation field, image plane compensation and rolling shutter micro-distortion correction are performed on the stabilized image, and the compensated and reconstructed image is output for texture analysis.

[0041] Furthermore, the steps for generating the spatial risk calibration map are as follows:

[0042] After compensation and reconstruction, the image is reconstructed using multi-scale structural tensors to generate a texture robust field.

[0043] By combining the immersion pressure field and the particle mobilization probability, a hierarchical conditional random field segmentation is performed on the texture robust field to obtain suspected candidate regions.

[0044] Morphological connectivity detection and stability threshold determination are performed based on suspected candidate regions to generate a spatial risk calibration map.

[0045] Compared with the prior art, the beneficial effects of the present invention are:

[0046] This invention acquires image sequences from airborne cameras, acceleration and angular velocity data from inertial measurement units, propeller speed, flight altitude, wind speed and direction, and aircraft vibration signals. Using image exposure time as a unified time axis, time-series samples are formed. Based on these time-series samples, a priori arithmetic of near-ground boundary layer momentum entrainment during rotor downwash is established. Using propeller speed, flight altitude, and wind field data, time-series estimations of downwash dynamic pressure field, particle mobilization probability time series, and equivalent micro-attitude perturbation spectrum are calculated. Subpixel-level optical flow and phase accumulation analysis are performed on the image sequences. Blind source separation is performed on the high-frequency components of aircraft vibration and inertial measurement, and time cross-correlation is performed with the equivalent micro-attitude perturbation spectrum to obtain coupling evidence. This establishes a physical mapping relationship between the priori arithmetic of near-ground boundary layer momentum entrainment during rotor downwash and the probability of surface particle mobilization, transforming the UAV from a traditional passive image acquisition platform into an active aerodynamic excitation detection source. By calculating the multi-scale time-frequency coupling evidence between the dynamic observable set of image fluctuations and the characterization of the body's natural vibration response, this invention can capture the microscopic non-rigid deformation and fine particle resuspension characteristics of loose deposits under the action of downwash dynamic pressure under non-contact conditions. This mechanism enables the system to identify hidden geological hazard points that have indistinct texture features on static images but are in a critically stable mechanical state, thus realizing the dynamic perception of the looseness of the surface medium.

[0047] This invention employs a joint determination of trigger strength and coupling degree based on the temporal estimation of the downwash dynamic pressure field, coupled evidence, and the time series trigger strength and coupling degree of particle mobilization probability. This triggers an adaptive fine-tuning strategy for flight parameters, performing image plane motion compensation and rolling shutter micro-distortion correction on the image within the triggering segment. The resulting compensated and reconstructed image is then analyzed for spatial risk to generate a spatial risk calibration map. Blind source separation technology decouples the main vibration mode of the airframe structure from aerodynamic environmental noise, inverting the equivalent micro-attitude perturbation spectrum. Based on this spectrum, an image plane motion compensation field for rolling shutter micro-distortion is established, effectively restoring the geometric fidelity of the perturbed image. By introducing a texture robust field and constructing a hierarchical conditional random field model, the temporal estimation of the downwash dynamic pressure field and the posterior probability of particle-induced fluctuations are updated using Bayesian fusion. This method not only suppresses false triggers caused by environmental wind shear or illumination flicker but also quantifies the causal relationship between aerodynamic energy and surface response through a weighted coherence index, thus outputting a spatial risk calibration map with clear physical meaning. Attached Figure Description

[0048] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. The following drawings are not drawn to scale according to the actual size, but are intended to show the main idea of ​​the present invention.

[0049] Figure 1 is a flowchart of the method of the present invention;

[0050] Figure 2 is a flowchart of steps S2021-S2024 of the present invention;

[0051] Figure 3 is a diagram of the cross-shaped survey line layout centered on the survey station according to the present invention;

[0052] Figure 4 is a schematic diagram illustrating the conversion principle of wind speed at different altitudes according to the present invention. Detailed Implementation

[0053] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are also within the scope of protection of the present invention.

[0054] As shown in Figure 1, a method for rapid identification of geological hazard hazard points based on UAV aerial surveying includes:

[0055] Step 1: Collect image sequences from airborne cameras, acceleration and angular velocity from inertial measurement units, propeller speed, flight altitude, wind speed and direction, and aircraft vibration signals, and form time-series samples with image exposure time as a unified time axis.

[0056] Step 2: Based on time-series samples, establish a priori momentum entrainment in the near-ground boundary layer of the rotor downwash. Using propeller speed, flight altitude and wind field data, calculate the time-series estimation of the downwash dynamic pressure field, the time series of particle mobilization probability and the equivalent micro-attitude perturbation spectrum.

[0057] Step S201: Establish the prior of rotor downwash near-ground boundary layer momentum entrainment, construct an axisymmetric-shear coupled dynamic pressure field parameterization based on propeller speed, flight altitude and wind speed and direction, and solve the downwash dynamic pressure field time series estimation through energy conservation and momentum flux closure.

[0058] Establish a dynamic coordinate system centered on the UAV. Specifically, define the origin of the coordinate system as the vertical projection point of the UAV's center of mass onto the ground. The current heading of the UAV is defined as the 0-degree reference direction in the polar coordinate system, thereby constructing a ground polar coordinate system that translates with the UAV. ,in Radial distance, The coordinate system is the azimuth angle. A one-to-one mapping relationship is established between this coordinate system and the airframe inertial coordinate system through the real-time pose of the UAV. Pose refers to position and attitude. The dynamic coordinate system is used to accurately describe the interaction between the rotor downwash and the ground surface. In this coordinate system, an axisymmetric-shear coupled parameterized expression of the rotor downwash dynamic pressure field is established. The prior definition presupposes that the flow field structure is composed of a superposition of an axisymmetric reference field under ideal windless conditions and a shear disturbance field caused by ambient wind. The axial velocity of the axisymmetric reference field... radial Distribution Parameterization can be exemplified using a Gaussian function: The parameters here This represents the maximum speed in the downwash core region and is positively correlated with rotor thrust. It is an exponential function with parameters Represents the radial attenuation scale of the downwash airflow, related to flight altitude and propeller diameter, axisymmetric reference field. The shear disturbance field is represented by an offset vector that reflects the blowing effect of the ambient wind on the axisymmetric flow field. The magnitude of this offset vector is proportional to the wind speed, and its direction is related to the wind direction, causing the actual downwash center to shift from the origin. Offset to point This causes the isobars in the flow field to exhibit a non-concentric shape, with an exemplified wind speed being... The wind direction angle is In Cartesian coordinates, the shear perturbation field Couple these two fields to obtain an arbitrary point. The total velocity vector at that point, and according to the fluid dynamics formula The dynamic pressure is calculated, where air density, The total velocity after coupling is obtained by vector addition of the axisymmetric reference field and the shear perturbation field; the unknown parameters in the above parameterized expression are determined in real time by constructing and solving the closed equations of energy conservation and momentum flux, such as... and This is achieved by constructing and solving a closed set of equations relating energy conservation and momentum flux. The specific method is as follows: Establish constraint equations based on energy conservation: This law refers to the mechanical power output of the rotor system. It should be equal to the increase in the kinetic energy of the entire downwash flow field, and the left side of the equation is: It can be estimated based on propeller speed, pitch, and motor model. For example, the mechanical power output by the motor is equal to its output torque. Multiplied by the angular velocity of the propeller ,Right now The right side of the equation is the kinetic energy density of the parameterized velocity field ( The result obtained by spatial integration is a result containing unknown parameters. and The mathematical expression, by setting both sides equal, yields the first expression about... and The constraint equations are established based on momentum flux closure: this law states that the thrust generated by the rotor to maintain hovering or flight... It should be equal to the vertical momentum flux of the downwash gas passing through the near-surface control surface; the left side of the equation represents the thrust. Its value can be directly obtained from the flight control system. The right side of the equation is the vertical momentum flux density of the parameterized velocity field. The result obtained by performing an area integral is another result containing unknown parameters. and The mathematical expression, by setting both sides equal, yields the second expression about... and The constraint equations are then solved. By simultaneously solving these two nonlinear equations, a unique set of parameters can be obtained at each time step using real-time inputs such as propeller speed, flight altitude, wind speed, and wind direction. ,in, Characterizing time The maximum velocity in the downwash core region, which is positively correlated with rotor thrust. Characterizing time The radial attenuation scale of the downwash airflow, which is related to flight altitude and propeller diameter, is a time-varying parameter obtained by inverse solving of real-time flight data. Substituting the solved parameters into the parameterized expression, a two-dimensional dynamic pressure distribution map on the ground polar coordinate grid at that moment can be generated. Furthermore, by stitching together the dynamic pressure distribution maps with continuous time steps, a time series estimate of the downwash dynamic pressure field is formed. This time series estimate refers to the dynamic pressure field sequence that evolves over time, reflecting in real time the dynamic evolution process of the distribution, intensity, and morphology of ground dynamic pressure under the combined effects of UAV control input and changes in the external environment.

[0059] As shown in Figure 2, step S202: by establishing a station-grid mapping, extract the surface roughness and zero displacement height parameters, combine wind direction, wind speed and surface zoning information, calculate the critical wind erosion and shear thresholds for each particle size, and consider uncertainty correction to obtain the particle mobilization probability time series.

[0060] Step S2021: Using the main line-of-sight projection point as the station, set up reference survey lines along the main wind direction and lateral direction respectively, unify the near-ground wind information to the same height benchmark, combine the surface type zoning and micro-topographic undulation information, perform robust fitting and consistency test on the candidate roughness parameters, output the surface roughness parameters, zero displacement height parameters and their uncertainty range corresponding to the spatial location, and establish the station-grid mapping index.

[0061] Figure 3 shows the layout of the cross-shaped survey lines centered on the survey station;

[0062] Referring to Figure 3, it illustrates how the system establishes a spatial reference for local environmental parameter analysis. The black dot in the center of the figure represents the station, which is the vertical projection of the UAV's main line of sight onto the ground, serving as the origin of the analysis. Using this as the center, the system lays out a main survey line based on the current prevailing wind direction (indicated by the arrow) and its opposite direction, and a lateral survey line perpendicular to the prevailing wind direction, thus forming a cross-shaped sampling path. This layout aims to sample key environmental information such as surface roughness and micro-topographic undulations along specific directions, so as to subsequently establish a mapping index between the station and the geographic raster, and provide a standardized spatial reference for inverting surface aerodynamic parameters.

[0063] Figure 4 shows the conversion principle diagram for wind speeds at different altitudes.

[0064] Referring to Figure 4, the physical process of normalizing wind speed to altitude using a logarithmic wind profile model is described. In the coordinate system, the horizontal axis represents wind speed, and the vertical axis represents altitude. The curve shows that wind speed increases logarithmically with altitude. The figure demonstrates how to calculate the normalized wind speed at a standard reference altitude (e.g., 10 meters) from the measured wind speed of a UAV at low altitude (e.g., 2 meters) using a conversion formula that includes parameters such as surface roughness and zero displacement altitude. This step eliminates data bias caused by different measurement altitudes, ensuring that the wind field data used in subsequent calculations of the critical wind erosion threshold and particle mobilization probability have a unified and accurate physical benchmark.

[0065] As shown in Figures 3 and 4, the main line-of-sight projection point is used as the station. The main line of sight refers to the central axis of the beam emitted by the remote sensing equipment, such as lidar or microwave radar. The projection point is the unique geographical location where this central axis intersects with the surface represented by the pre-acquired high-precision digital elevation model (DEM). This projection point is defined as a temporary station, serving as the reference origin for local environmental parameter analysis. Reference survey lines are laid out along the main wind direction and the lateral direction. Here, the main wind direction refers to the direction of the prevailing wind vector near the ground at the station's location, obtained from external meteorological data sources such as weather forecast models and measured data from nearby meteorological stations. The lateral direction is the direction perpendicular to the main wind direction. The specific method for laying out the reference survey lines is as follows: with the station as the center, a first virtual straight line is laid out along the main wind direction and its opposite direction, and a second virtual straight line is laid out along the lateral direction. These two orthogonal survey lines constitute the path for sampling surrounding environmental information starting from the station. To unify near-ground wind information to a common altitude reference, it's crucial to address the issue that wind speed data from different sources may correspond to different measurement altitudes. For example, a meteorological model might output wind speed at 10 meters above ground, while a drone's onboard sensor might measure wind speed at 2 meters. Since wind speed varies logarithmically or exponentially with altitude, altitude normalization is necessary to ensure the accuracy of subsequent calculations. For instance, the logarithmic law wind profile formula from atmospheric boundary layer theory can be used. All wind speed information is converted to wind speed values ​​at a standard reference height, for example, 10 meters. It is the wind speed at altitude z. It is the friction speed, an intermediate physical quantity. It is the von Kármán constant, approximately equal to 0.4. It is the height above the ground. It is the zero-displacement height parameter, that is, the zero-plane displacement height, which characterizes the vertical lift of the zero point of the wind speed profile. It is a surface roughness parameter, also known as aerodynamic roughness, which characterizes the resistance of surface roughness elements to airflow. It is a natural logarithmic function with the natural constant e as the base. Combining surface type zoning and micro-topographic relief information, robust fitting and consistency checks are performed on the candidate roughness parameters. Here, surface type zoning refers to dividing the area traversed by the survey line into different surface cover categories, such as water bodies, grasslands, forests, and urban built-up areas, based on satellite remote sensing imagery, land use databases, etc. Each category corresponds to an empirical range of candidate surface roughness parameters in the prior knowledge base. Micro-topographic relief information is extracted from a high-precision DEM, reflecting the actual surface undulations and obstacle distribution. The execution process is as follows: a series of sampling points are collected along the reference survey line at certain intervals to obtain the surface type and elevation of each sampling point, with an interval of, for example, 1 meter; based on the surface type of the sampling points, initial surface roughness parameters are selected from the prior knowledge base. and zero displacement height parameter As candidate parameters, using actual wind speed observation profiles along the survey line, a robust fitting algorithm, such as the Random Sample Consensus Algorithm (RANSAC) or the Huber loss function minimization algorithm, is employed to fit the logarithmic law wind profile formula, thus solving for the optimal solution. and The robust algorithm aims to eliminate interference from outlier data points caused by local eddies or measurement noise, resulting in a more robust fitting result; the fitted value... and The values ​​undergo a consistency test, determining whether they fall within the physically reasonable range corresponding to the surface type, and the goodness of fit is evaluated. If the parameters exceed the reasonable range or the fit is poor, they are marked as low-confidence results. The goodness of fit is expressed as the coefficient of determination R². The output includes the surface roughness parameters, zero-displacement height parameters, and their uncertainty ranges corresponding to the spatial location, and establishes a station-grid mapping index. After completing the above steps, a set of core surface aerodynamic parameters for the station location is output, including: the best estimated surface roughness parameters obtained after robust fitting and testing. Zero displacement height parameters The uncertainty range, calculated from the covariance matrix of the fitting algorithm, is based on the statistical results of the robust fitting algorithm and is typically given as a 95% confidence interval. This range provides a direct indication of the reliability of the parameter inversion results, including the surface roughness parameter. Zero displacement height parameter characterizes the Earth's surface's ability to drag airflow. The effective surface uplift height caused by dense obstacles is characterized to quantify the reliability of the estimation results. A station-raster mapping index is established to facilitate rapid retrieval of these local parameters in the global geographic information system. Specifically, the entire operational area is divided into a uniform geographic raster network, such as a 10m x 10m grid. A data lookup table or hash table is created, where the key is the station's unique identifier or geographic coordinates, and the value is a data structure containing a list of raster cells covered by the surrounding area analyzed by the station, as well as the calculated values ​​for these rasters. , And its uncertainty. In this way, when it is necessary to query the surface parameters of any raster, the relevant station calculation results can be quickly located through the index.

[0066] Step S2022: Based on the station-grid mapping index, map the image grid cells to the surface patches. For each cell, calculate the proportion of the visible area of ​​each patch. Correct the proportion using attitude and spectral consistency to form the partition weight vector and partition label of each cell.

[0067] Based on the two-dimensional pixel coordinate system of UAV aerial survey imagery, an image grid coordinate system is established, with the image width set as... Height is Divide the image into For any pixel unit Its center point has coordinates of [coordinates missing] on the image plane. i and j are pixel index variables, representing the horizontal and vertical directions in the image coordinate system. Using the station-raster mapping index, the projected position of the pixel in the surface coordinate system is determined through the ground-image geometric relationship. 'g' represents the surface coordinate variable. The land-image geometric relationship is calculated using a combined method of exterior and interior orientation elements, which will not be elaborated upon here. Based on the location of the projection point, one or more surface patch regions corresponding to the pixel are determined. A surface patch refers to a surface unit with relatively consistent roughness characteristics and material properties within a surface type zoning, such as bare rock areas, gravel areas, vegetated areas, and loose deposit areas. Each patch has a clearly defined boundary polygon in the surface coordinate system. Surface patch data can be generated through expert visual interpretation and manual delineation of high-resolution satellite maps, or it can be converted from existing land use / land cover thematic maps. These acquisition methods are conventional techniques for those skilled in the art. Calculating pixels. The visible area percentage of each surface patch under imaging geometry conditions. Specifically: based on the UAV attitude angles and camera intrinsic parameters, a geometric inverse calculation algorithm based on the attitude matrix and camera intrinsic parameter matrix is ​​used to calculate the set of intersection points between the pixel's line of sight and the ground. The calculation method is existing technology and will not be elaborated here. The UAV attitude angles, such as pitch angle, roll angle, and yaw angle, are used to determine the projected area of ​​the pixel's line of sight within each patch area based on the set of ground intersection points. k represents the surface patch index variable, which represents the total coverage area. As a normalization benchmark, the visible area percentage of each patch is obtained. Attitude and spectral consistency corrections are performed on the visible area percentage. Attitude consistency correction is based on the angle between the camera's optical axis and the ground normal. ,pass Factor correction for local area distortion caused by tilted imaging, the corrected visible area percentage is ,in It can be calculated from the attitude matrix components; spectral consistency correction is used to eliminate the influence of differences in reflectance of different surface materials on the weights, by extracting pixels. spectral reflectance vector It includes red, green, blue, and near-infrared channels, and incorporates typical spectral templates of surface patches. Calculate cosine similarity This similarity coefficient is used to perform spectral correction on the attitude-corrected weights, and the corrected patch weights are: Normalize the weights of all patches so that... The partition weight results for each pixel are then organized into a partition weight vector. ,in This represents the number of patches mapped to this pixel. The primary partition label of the pixel is determined based on the patch number corresponding to the maximum weight. ,in Indicates searching for Take the maximum value The final output for each pixel is two types of results: one is a partition weight vector containing the weights of each patch, used for subsequent threshold weighting calculation; the other is the main partition label, used for regional classification association. The above method completes the multi-source fusion mapping from aerial survey image pixels to surface partitions, providing reliable basic data for the partition weighting calculation of the critical wind erosion threshold for particles in the subsequent step S2023. The entire processing flow can be automatically executed frame-by-frame and pixel-by-pixel in the image processing module, ensuring the spatiotemporal consistency and repeatability of the mapping results under different attitudes and surface conditions.

[0068] Step S2023: Combine the partition weight vector and partition label of each pixel, as well as the surface roughness parameter and zero displacement height parameter, with the wind direction and near-ground wind information that change over time, and weight the partition threshold parameters to obtain the critical wind erosion threshold and equivalent critical shear threshold for each particle size. Based on turbulent gusts, surface parameter uncertainty and partition heterogeneity, generate a comprehensive uncertainty scale.

[0069] Obtain wind direction in time series near-ground wind speed Data, with temporal resolution synchronized with drone image sampling, at every moment The reference shear velocity was calculated based on the logarithmic wind speed profile model. :

[0070] ;

[0071] in Height above the ground At ground height Time series data of near-surface wind speed collected synchronously at the location; For the first The zero-displacement height parameter of the patch-like structure characterizes the vertical lift of the zero point of the wind speed profile; For the first The surface roughness parameters of the patch-like structures characterize the impedance effect of surface roughness elements on airflow. These parameters are calculated independently for each zone to obtain the instantaneous shear velocity at the zone level. For each type of surface patch material, the theoretical critical wind erosion threshold under undisturbed conditions is determined based on the Bagnold wind erosion initiation formula or the Shao particle initiation empirical model. This falls under existing technology in the field of wind erosion mechanics and will not be elaborated upon here. Considering the actual influence of surface roughness and zero displacement height, the theoretical threshold is corrected to an effective threshold:

[0072] ;

[0073] in This is the wind erosion correction factor. The average particle size of the plaque particles. For the main windward direction of this patch, the first part of the correction term describes the impedance effect of surface roughness, and the second part describes the directional correction caused by wind direction deviation. This is based on the weights of each patch within the pixel. For patch-level thresholds By performing a weighted average, the equivalent critical shear rate at the image pixel scale is obtained. : Then, by grouping and statistically analyzing the particles by size, the critical wind erosion threshold for each particle size class was obtained. superscript This indicates the particle size category. To assess the impact of surface and meteorological fluctuations on the threshold results, the obtained... An uncertainty quantification model is constructed, which simultaneously considers three types of error sources: (1) instantaneous fluctuations in wind speed caused by turbulent gusts, using the standard deviation of wind speed. and turbulence intensity Characterization, (2) The measurement uncertainty of surface parameters, including the average near-ground wind speed measured within a given time period; estimation error , Represents surface roughness parameters The estimated standard deviation of the covariance matrix of the robust fit is the square root of the diagonal element of the estimated standard deviation of the corresponding parameter; (3) Spatial heterogeneity index of the partition The underlined value indicates a neighborhood averaging operation, with a range of values ​​of 100. A larger value indicates a more pronounced difference in local land surface types. For pixels The primary surface zoning category identifier, and its weighted weight. The patch number corresponding to the largest one is... To quantify the spatial heterogeneity between adjacent pixels, each primary partition label is... Mapped to feature vectors This vector consists of the surface type encoding of its corresponding patch, and cos is the cosine similarity calculation. (Integrated uncertainty scale) Defined as: ;

[0074] in The empirical weighting coefficient is used to identify the threshold confidence level of each pixel in space. Through the above processing, the distribution maps of critical wind erosion thresholds and equivalent critical shear thresholds for each particle size are obtained, along with the comprehensive uncertainty scale corresponding to each pixel. The results can be used for wind erosion dynamic simulation and wind-blown sand transport risk assessment in subsequent step S2024, providing a repeatable quantitative basis for the surface wind erosion monitoring model under UAV imagery. The entire calculation process can be automatically completed by the data processing module, ensuring that the logical relationship between time series, zoning characteristics, and parameter weighting is consistent and the physical meaning is clear.

[0075] Step S2024: Based on the critical wind erosion threshold, equivalent critical shear threshold and comprehensive uncertainty scale of each particle size, and combined with time-varying frictional wind information, calculate the particle mobilization probability of each particle size, and weight it according to the local particle size distribution to obtain the particle mobilization probability time series.

[0076] In a unified geographic grid In the coordinate system, for each pixel Read its different moments friction speed The sequence, and the corresponding equivalent critical shearing threshold. Comparisons were made for each particle size class. ,when Exceeding its corresponding critical wind erosion threshold At that moment, the particle size is considered to be potentially mobilized by the airflow. To quantify the uncertainty of mobilization, a probability function is used. Perform calculations, where It is either a standard normal distribution function or a cumulative distribution function determined by measured statistics, reflecting the relative exceedance probability between friction speed fluctuations and a threshold. The value range is [0, 1], representing the particle size grade. The probability of a particle being moved by wind at a specified location and time. To obtain the actual overall wind erosion mobilization probability of the earth's surface, the mobilization probability of each particle size class is calculated. Based on local measured or remote sensing inversion particle size distribution By performing a weighted summation, the comprehensive particle mobilization probability at the pixel scale is obtained. To characterize the statistical reliability of the calculation results, the uncertainty scale is used based on the probability distribution. Extract confidence intervals to form 95% confidence upper and lower bound probability fields. ,in These represent the lower and upper limits of the overall particle mobilization probability at a 95% confidence level, respectively, used to define the probability fluctuation range under uncertainty conditions. The probability of each pixel is ordered along the time dimension. The sequence is composed of a time series plot of particle mobilization probability, and the probability is calculated at a set threshold. (like The occurrence frequency under certain conditions is used to generate a wind erosion occurrence frequency raster map based on threshold determination. This is achieved by traversing all time points. With spatial grid The calculation results yield a spatiotemporal distribution dataset of particle mobilization probability with clear temporal and spatial resolution, along with corresponding confidence intervals, providing a unified input basis for subsequent wind erosion flux estimation, geomorphological evolution simulation, and remote sensing inversion verification.

[0077] Steps S2021 to S2024 aim to establish a physical constraint relationship between particle stress and wind erosion response based on surface wind field and underlying surface characteristics, thereby achieving dynamic quantification of the mobilization probability of surface particles. Specifically, step S2021 standardizes complex surface conditions by establishing a station-grid mapping and performing robust inversion of roughness and zero displacement height, providing accurate boundary parameters for subsequent local wind field calculations. Step S2022 maps image pixels to surface zones, using visible area ratio and attitude-spectral consistency correction to construct a zone weighting system linking pixels to the surface, enabling a one-to-one mapping between image information and surface physical properties. Step S2023, based on this, combines time-varying wind direction and speed to obtain the critical wind erosion and shear thresholds for each particle size group, and quantifies the uncertainties caused by turbulence, surface parameters, and spatial heterogeneity. Finally, step S2024 calculates the mobilization probability of particles of different sizes based on the dynamic comparison of thresholds and frictional wind, and generates a time series weighted by particle size distribution. This series of steps achieves the layer-by-layer transfer and fusion of three-dimensional elements—wind field, surface, and particle size—into time-series mobilization probability, enabling the system to accurately characterize the initiation pattern of surface particles under the influence of downwash airflow. Through this process, the scheme not only improves the physical accuracy of wind erosion threshold estimation and aerodynamic response calculation, but also ensures that the particle mobilization probability results are spatially matched with imagery and temporally synchronized with the wind field, thereby enhancing the scientific rigor, traceability, and dynamic judgment capability of hazard point identification.

[0078] Step S203: Multiply the timing estimate of the downwash dynamic pressure field with the transfer matrix of the body / gimbal to obtain the equivalent torque disturbance. Combine the maneuvering state of the inertial measurement unit to solve the equivalent micro-attitude disturbance spectrum.

[0079] ground polar coordinate system The time series estimation results of the hysteresis field obtained below The direction cosine matrix is ​​projected onto the aircraft's aerodynamic coordinate system. This matrix is ​​obtained from the UAV's current attitude angles using existing Euler angle transformation coefficients. The principle of this coordinate transformation is a well-known technique in the field of aircraft aerodynamic analysis and will not be elaborated here. The origin is taken as the UAV's center of mass. The longitudinal axis of the fuselage and the direction in which the nose points are defined as... The axis, the transverse span direction is defined as The axis, defined vertically downwards, is... The axes, along with the other two, form a right-handed rectangular coordinate system. Within this system, the time-series estimation of the projected downwash dynamic pressure field is spatially superimposed according to the geometric model of the UAV's lower surface or wing surface. This involves weighting and accumulating the dynamic pressure estimates at each surface micro-element with their corresponding normal direction and position relative to the center of mass. By comprehensively considering the direction, arm length, and temporal variation of the aerodynamic forces in each region, the torque at the center of mass is obtained. This process essentially transforms the time-varying pressure field distribution into a total torque response acting on the body's center of mass, used to characterize the perturbation effect of transient aerodynamic force uneven distribution caused by sway on the UAV's attitude. Based on the structural connection characteristics of the body-gimbal, a linear transfer matrix is ​​established for the aerodynamic perturbation torque transmitted through the structure to the camera's optical axis. The matrix is ​​established based on existing linearization theories of aircraft structural dynamics. The frequency domain transfer function is derived by linearizing the finite element model of the airframe-gimbal system within a small perturbation range. Experimental modal parameters include the stiffness matrix, damping matrix, and mass matrix. Specific methods for obtaining these parameters can be found in existing standard literature, such as *Flight Dynamics Principles, Cook, 2013*. This is a mature existing technology and will not be elaborated further here. Based on the principle of mechanical equilibrium, the equivalent torque perturbation is calculated. This product form has a clear physical basis and a theoretical foundation in linear systems: under the assumption of small perturbations, the transmission of torque input to output can be approximately expressed as a first-order linear system response relationship. The attitude dynamics equations based on the linearized perturbation model using Euler's dynamics equations are established in machine-system coordinates.

[0080] ;in, The equivalent inertia tensor of the body-gimbal system can be obtained using existing technologies, such as calculating it using conventional moment of inertia calculation methods based on known three-dimensional mass and structural parameters of the entire machine, or measuring it using existing experimental identification methods. It is the angular velocity vector. For gimbal control torque input, Let represent the angular acceleration vector. After linearization using the small perturbation assumption, this equation can be approximated as: This characterizes the linear mapping relationship between the equivalent disturbance torque and angular acceleration. The frequency domain spectrum of the attitude angular disturbance is obtained by performing a fast Fourier transform or power spectrum estimation on the attitude dynamics equations. ,in, Let be the angular response function of the system. The equivalent moment perturbation spectrum is derived from the attitude dynamics equations. The frequency domain transfer characteristic is obtained through Fourier transform derivation and integration. f is the frequency. The imaginary unit is used to construct complex variables in the frequency domain transfer function. The equivalent moment perturbation spectrum is derived from the aforementioned aerodynamic perturbations. Through linear transfer matrix Obtain the equivalent moment The power spectrum is estimated or obtained by performing power spectrum estimation or Fourier transform on the torque time series. Combining the real-time angular velocity and angular acceleration data output by the inertial measurement unit, the obtained... Based on the time window characteristics of the IMU output, a sliding window or wavelet packet decomposition method is used to correct the time-frequency characteristics, ensuring that the disturbance energy distribution in different frequency bands is consistent with the actual aerodynamic response. The core objective is to eliminate energy distortion caused by variations in the sampling interval or signal non-stationarity. The corrected attitude power spectrum is then subjected to energy normalization, i.e., the average spectral energy during the system's steady-state period is used as a reference benchmark to determine the reference energy. The scaling factor is the multiplication of the frequency domain spectrum of the attitude angle perturbation by a scaling factor. This approach conserves total energy while distributing energy across frequency bands according to probability density, preserving the overall shape of the spectrum and standardizing its intensity. This ensures a consistent total energy while balancing the distribution between frequency bands, resulting in an equivalent micro-attitude perturbation spectrum that matches the actual dynamic response of the system. It is used to quantify the impact of motion disturbances on camera-stabilized imaging, and to provide input parameters for attitude compensation and stabilization control algorithms.

[0081] Step 3: Perform subpixel-level optical flow and phase accumulation analysis on the image sequence, perform blind source separation on the high-frequency components of body vibration and inertial measurement, and perform time cross-correlation with the equivalent micro-attitude perturbation spectrum to obtain coupling evidence;

[0082] Step S301: Perform subpixel-level optical flow and phase accumulation analysis on the image sequence to extract the multi-scale vector field and corresponding phase energy density sequence of local motion in the image, and establish a dynamic observable set of image fluctuations;

[0083] A continuous image sequence encompassing the steady-state hovering or flight of a UAV is acquired. This image sequence is obtained at a fixed frame rate by an optical imaging device mounted on the UAV body or gimbal. The image sequence is then subjected to inter-frame registration and geometric correction in chronological order to eliminate lens distortion and global translation effects. On the corrected image sequence, using pixel intensity distribution as input, a sub-pixel-level optical flow calculation method is employed to determine the local motion vector between adjacent frames. Specifically, brightness or gradient changes are calculated within the neighborhood of each pixel. The initial pixel displacement is obtained by solving the optical flow constraint equation. Sub-pixel interpolation or fitting methods are used to refine the displacement results. Fitting methods include, for example, bilinear interpolation or phase-correlation interpolation, to achieve optical flow field calculations with an accuracy of less than one pixel. After acquiring the optical flow vector field between each frame, phase accumulation analysis is performed on the phase change of the same pixel location across multiple consecutive frames. This phase accumulation analysis, based on Fourier transform or complex wavelet transform, represents the brightness change of the time series as a phase signal, and accumulates it frame by frame to obtain the phase drift amount for each time period, thus reflecting the local fluctuation characteristics of the image under the influence of minute structural disturbances or vibrations. To capture the relative motion characteristics at different spatial scales, the optical flow and phase information are processed hierarchically through multi-scale filtering or pyramid decomposition to obtain a multi-scale vector field set ranging from large-scale overall displacement to small-scale local texture jitter. Pyramid decomposition methods include, for example, Gaussian pyramids or Laplace pyramids. The energy density value of the phase change is calculated at each scale to characterize the dynamic energy distribution of the image corresponding to different spatial frequencies. The phase energy density can be characterized by calculating the square of the phase change rate with time or its spectral power. The multi-scale vector field data and the phase energy density data are combined under a unified coordinate and time index to form a dynamic observable set of image fluctuations. This observable set reflects the minute motion direction, displacement amplitude, and corresponding phase energy characteristics of each region in the image, and can be used as input parameters for subsequent steps to identify the coupling effect between body attitude perturbations and aerodynamic forces.

[0084] Step S302: Perform blind source separation on the body vibration signal and the high-frequency components of the inertial measurement unit, extract the main vibration mode and noise basis of the structure, and form a characterization of the body's natural vibration response;

[0085] Piezoelectric vibration sensors installed at the drone's arms, gimbal connections, and core stress points simultaneously collect mechanical vibration data from the drone's surface. This data, combined with raw triaxial acceleration and angular velocity data output from the inertial measurement unit, allows for the setting of a sampling frequency. The frequency should be no lower than 500Hz to cover the natural frequency of the body structure and the high-frequency aerodynamic disturbance band; the acquired multi-channel raw signals are preprocessed by removing DC components, linear drift correction, and bandpass filtering to filter out slowly varying rigid body motion components with frequencies below 1Hz and high-frequency noise higher than 5 times the first-order bending frequency of the structure. The time-series signals from each sensor channel are constructed into an observation signal matrix. ,in For time variables, This is the transpose of the matrix. Indicates the first The time-domain observation data components acquired by each sensor channel; based on the assumption of a linear instantaneous mixing model, the observation signal matrix... Characterized as , in for An unknown set of statistically independent source signals. Indicates the first Each of the statistically independent potential physical source signal components to be separated for An unknown constant mixing matrix of dimension 3 represents the path attenuation and coupling coefficient of each vibration source transmitted to the sensor. The term is additive white Gaussian noise; the FastICA algorithm or the SOBI algorithm based on second-order statistics in independent component analysis is used to analyze it. Blind source separation is performed by iteratively optimizing the time delay correlation matrix by maximizing its non-Gaussianity or diagonalizing it, thus solving for the separation matrix. That is, the mixed matrix The generalized inverse estimator is used to calculate the estimated value of the source signal. For each isolated independent source component Fast Fourier Transform was performed to calculate its power spectral density, and classification and identification were performed based on spectral characteristics: the components with concentrated spectral energy, significant peaks, and peak frequencies matching the natural frequencies obtained from the finite element modal analysis of the UAV airframe structure were defined as the dominant vibration modes of the structure. This mode characterizes the inherent elastic response of the airframe under rotor rotation excitation. Specifically, the spectral energy concentration is obtained by calculating the ratio of local band energy near the main peak frequency to the total energy of the entire frequency band. A preset bandwidth window centered on the main peak frequency is set, and the energy proportion within this window is calculated. If the energy proportion exceeds a preset concentration threshold, it indicates that the energy of this component is highly concentrated at a specific frequency point, and is judged as spectral energy concentration. The peak factor of the spectrum is used to measure this, that is, the ratio of the maximum amplitude in the spectrum to the average amplitude of the entire frequency band. If the ratio is greater than a preset significance threshold, it indicates that the frequency component has a prominent peak characteristic in the background noise, and is judged as having a significant peak. The identified main peak frequency is compared with the pre-acquired set of natural frequencies of the UAV body structure. The set of natural frequencies of the UAV body structure is obtained through finite element simulation or modal experiments, such as the first-order rotor overpass frequency and the arm bending frequency. The relative deviation between the main peak frequency and the frequency with the closest value in the set of natural frequencies is calculated. If the relative deviation is within a preset matching tolerance range, such as within 5%, the two are judged to match, and it is confirmed as the main mode of structural vibration. The residual component with relatively uniform spectral energy distribution over a wide bandwidth, without significant spikes, and whose statistical characteristics approximate a Gaussian distribution, is defined as the noise basis. The substrate includes random airflow disturbances and sensor noise floor. Evaluation is performed by calculating the spectral flatness, which is the ratio of the geometric mean to the arithmetic mean of the spectral amplitude. If this ratio is close to 1, it indicates that the spectral energy distribution is relatively flat across frequency bands, with no obvious spikes and a relatively uniform distribution. In other words, if the calculated peak factor is less than a preset significance threshold, it indicates that there are no prominent narrowband components or significant peaks in the signal. The kurtosis index of the independent source component time-domain signal is used for measurement. Kurtosis is a statistical measure reflecting the steepness of the probability density distribution curve. If the calculated kurtosis value is close to the standard value of a normal distribution (close to 3), then the statistical characteristics of this component conform to a Gaussian distribution, thus confirming it as a noise base. The extracted structural vibration principal modes are then... With noise base Construct a characterization of the body's natural vibration response by aligning and encapsulating the structure along the time axis. This serves as a benchmark for distinguishing between external aerodynamic environmental excitation (such as wash flow rebound) and the mechanical vibration of the machine itself in subsequent steps.

[0086] Step S303: Cross-correlation and coherence tests are performed on the dynamic observable set of screen fluctuations and the characterization of the body's natural vibration response on a unified time axis with the equivalent micro-attitude perturbation spectrum to obtain the comprehensive coupling coefficient. Based on the comprehensive coupling coefficient, it is determined that there is significant multi-scale time-frequency coupling, and coupling evidence is output.

[0087] The dynamic observable set of screen fluctuations and the characterization of the body's natural vibration response are synchronized and registered on a unified time axis. Using the timestamp of the inertial measurement unit as a global reference, signals at different sampling rates are resampled through spline interpolation to ensure that the dynamic observable set of screen fluctuations and the body vibration signal are synchronized. Each value corresponds to a unique value at any given time. Zero-mean and standardization are applied to the synchronized signal to eliminate amplitude differences and dimensional effects. To determine the time delay relationship between the image fluctuation response and the dominant vibration mode, values ​​are calculated at each representative spatial scale. Next, calculate the time-domain cross-correlation function. , Representing each scale as well as Multi-scale vector field of position, This represents the time delay or time shift, and the peak position is determined. As a phase delay estimate at this scale This represents the mathematical expectation operator. The result is... Phase compensation used in subsequent frequency domain analysis, Multiply by the phase correction factor in the frequency domain , This represents the imaginary unit to ensure that the image response at different scales is phase-aligned with the vibration signal within the dominant frequency energy range. The normalized signal for each scale... With the revised Perform a Fast Fourier Transform (FFT) to calculate the power spectrum. , and cross-power spectrum The scale-dependent coherence function is obtained. Its value range is [0, 1], representing the frequency. With scale The degree of linear coupling between the dynamics of the image and the vibration of the organism is shown below. To comprehensively measure the multi-scale coupling strength, the coherence spectra of all scales are sequenced according to their phase energy density. Average energy in the frequency domain The weights are used as weights to superimpose the data, resulting in a full-scale integrated coherence spectrum. The full-scale integrated coherence spectrum With equivalent micro-attitude perturbation spectrum In the effective frequency band The weighted integral is performed to obtain the overall coupling coefficient. ,in and These are the critical frequency points for the rise and fall of attitude disturbance energy, determined within a range of ±20dB from the actual spectral peak. The integration interval covers the main coupling frequency band after compensation to ensure that the calculated energy truly reflects the physical relationship between the image, structure, and attitude. To enhance the reliability of the quantification results, [further details are needed]. and Monte Carlo resampling calculation sample set Find Confidence interval And define the confidence level ,in When performing Monte Carlo statistical tests, the representative... The comprehensive coupling coefficient value corresponding to the random substitute data generated by the secondary resampling, Calculate background coupling energy And define the signal-to-noise ratio enhancement factor. ,by replace Repeat the same cross-correlation, coherence spectrum, and weighted integration processes to obtain the comprehensive coupling coefficient under the background state. .like and If the value is located at the upper limit of the confidence interval, then significant multi-scale time-frequency coupling is determined to exist, and coupling evidence is output. With confidence level The results characterize the energy transfer and response synchronization relationship between dynamic fluctuations in the image, structural vibrations, and attitude disturbances, providing quantitative input for subsequent triggering and discrimination steps. This indicates the number of Monte Carlo resampling attempts and the number of samples. and These represent the lower and upper limits of the confidence interval, respectively. This represents the significance level, corresponding to the residual probability part of the confidence score. Let... , This indicates the corresponding confidence level symbol. This refers to the statistical reliability index.

[0088] Step 4: Based on the time series estimation of the downwash dynamic pressure field, the joint determination of the triggering strength and coupling degree of the particle mobilization probability time series, a flight parameter adaptive fine-tuning strategy is triggered.

[0089] Step S401: Based on the time series estimation of the downwash dynamic pressure field, perform temporal resampling and spatial interpolation on the particle mobilization probability time series to generate a frame-by-frame particle mobilization probability map, and complete image-to-ground registration with the image coordinates.

[0090] Based on the time series results of the downwash dynamic pressure field An instantaneous force equilibrium model was established for surface particles: ,in For fluid density, It is the acceleration due to gravity. The equivalent particle size is... For relative density, The coefficient of friction, This represents the force balance value. A positive force balance value indicates that the fluid dynamic pressure is sufficient to overcome the static resistance of the particles. To express the continuity of this process, the particle mobilization probability is formed using a logistic function. ,in These are control parameters used to adjust the steepness of the probability curve; the resulting values ​​range from... To synchronize with the image timeline, the image sampling time series is resampled three times using piecewise Hermite spline time resampling to obtain synchronized particle mobilization probabilities, ensuring a one-to-one correspondence between the probability series and image frames in the time dimension. Polar coordinates are then used. The projection position of the pixel in the Earth's surface coordinate system is determined by the Earth-image geometry relationship. In the geographic grid The discrete, synchronized particle mobilization probability space is made continuous using bilinear or kriging interpolation, resulting in a planar grid particle mobilization probability matrix. , forming a frame-by-frame probability sequence To achieve image-to-ground registration, a projection mapping is established between the geographic coordinate grid and the image pixel coordinates. Perspective transformation is performed using the camera geometry parameters of the UAV imaging system to achieve spatial alignment between the ground coordinates and the image coordinates. This projection and registration method is a mature existing technology and will not be elaborated upon in this application. Through this registration process, a particle mobilization probability distribution map corresponding to the image frame can be obtained for subsequent display and analysis.

[0091] Step S402: The set of observable dynamic fluctuations in the image is mapped into a fine-grained disturbance intensity distribution map through line-of-sight projection and ground-shadow geometric relationship. Using the frame-by-frame particle mobilization probability map as the prior, the observation likelihood function is constructed and the posterior inference is performed to obtain the posterior probability distribution of particle-induced fluctuations in each frame.

[0092] After transforming the multi-scale vector field and phase energy density sequence to surface coordinates through line-of-sight projection and ground-shadow mapping, the obtained motion and energy response information at each scale is preprocessed for normalization. Specifically, the vector modulus corresponding to each pixel in the multi-scale vector field is extracted to characterize the local motion intensity, and normalized to the [0, 1] range according to the maximum amplitude within each frame. For the phase energy density sequence, the corresponding phase energy response intensity is calculated and normalized within the same interval to eliminate scale and dimensional differences. After normalization, in each frame, the multi-scale results at the same geographic unit are fused according to the line-of-sight overlap relationship within the surface coordinate system to comprehensively reflect the instantaneous response characteristics of the region under airflow excitation, obtaining the fine-grained disturbance intensity distribution map at that moment. This intensity is used to characterize the real-time dynamic disturbance degree of surface fine-grained material caused by the rotor airflow and serves as input data for subsequent frame-by-frame probability inversion. The frame-by-frame particle mobilization probability map is used as the prior probability input, denoted as... The fine-grained disturbance intensity distribution map is used as the input observation data, denoted as... Based on the statistical relationship between disturbance intensity and particle fluctuation response, an observation likelihood function is constructed to describe the conditional probability relationship of observing a certain disturbance intensity under a given mobilization state, defined as follows: ,in This represents the undulating state of particles at this surface unit. This indicates that the particles were induced to fluctuate. This indicates that no fluctuation occurred. The observed likelihood can be expressed as... and The probability density functions for the two scenarios are used to characterize the response distribution of disturbance intensity to different fluctuation states. For each surface unit, the posterior probability is calculated using the Bayesian inference formula within the same frame, i.e.:

[0093] ;

[0094] In the formula, Indicates the observed disturbance strength as The posterior probability that the surface unit particles are induced to undulate; This indicates the intensity of the disturbance observed under conditions of particle fluctuation. The probability of; This indicates the intensity of the disturbance observed under conditions where the particles are not fluctuating. The probability of; This represents the prior probability that the surface unit particle was in a mobilizable state before the current frame. By normalizing the prior probability and the observed likelihood, the sum of the probabilities of all states corresponding to the denominator is guaranteed to be 1. The result after calculation is... This is the posterior probability distribution of particle-induced fluctuations for each frame. This probability distribution reflects the actual fluctuation response probability of particles within an instantaneous frame under the influence of rotor airflow disturbance. A higher posterior probability value indicates that particles in that surface area are more likely to be in a disturbed fluctuation state. Calculated from consecutive frames... The sequence can be used to construct a time series probability plot to characterize the dynamic changes in surface particles caused by rotor disturbance.

[0095] Step S403: Use the comprehensive coupling coefficient as the weight of the posterior probability distribution of particle-induced fluctuations in each frame to form a weighted coherence index.

[0096] Combined coupling coefficient Introduced as a weighting factor, it is used to characterize the structure-aerodynamic energy coupling strength of local surface regions under multi-scale frequency domain characteristics. (Comprehensive coupling coefficient) Derived from the fusion calculation of multi-scale coherent features, a larger value indicates a stronger ground-airflow coupling effect in that region. As a weight, it is compared with the posterior probability of particle-induced fluctuations in the current frame. Weighted fusion is performed to obtain the weighted coherence index of the frame. In the formula, Indicates spatial location Location, Time Frame The weighted coherence value is used to comprehensively reflect the energy-response coupling degree of the surface unit under the influence of airflow disturbance; The comprehensive coupling coefficient extracted from the multi-scale coherence spectrum represents the energy interaction capability between the local surface structure and the aerodynamic field. The particle-induced fluctuation probability, inferred posteriorly at the same time, represents the true dynamic activity probability of surface fine particles at that location. This weighted calculation moderately enhances the contribution of particle activation probability in high-coupling regions while suppressing the response in low-coupling regions, thus forming a spatially coherent comprehensive index that simultaneously reflects both aerodynamic intensity and surface response activity. This is applied across the entire frame sequence. By storing and visualizing the data in chronological order, a weighted coherence distribution sequence of consecutive frames can be obtained. This sequence can be used to characterize the dynamic consistency of rotor-surface interaction and provide a quantitative basis for subsequent disturbance risk assessment and surface inversion model optimization.

[0097] Step S404: Using the frame-by-frame particle mobilization probability map, the posterior probability distribution of particle-induced fluctuations in each frame, and the weighted coherence index as ternary evidence, Bayesian update is used to obtain the frame-level comprehensive trigger probability, and time smoothing is performed within the sliding window to generate segment-level trigger results.

[0098] Using the frame-by-frame particle mobilization probability map, the posterior probability distribution of particle-induced fluctuations in each frame, and the weighted coherence index as ternary evidence inputs, a Bayesian update model for frame-level particle triggering probability is constructed. Characterizes the inherent potential for particles to be mobilized under the initial influence of airflow. Describe the real-time response probability of particle-induced fluctuations after combining observed perturbations with prior mobilization estimates. This reflects the strength of the structure-aerodynamic coupling and the aerodynamic coherence characteristics of the region. Based on these three independent but related probabilistic pieces of evidence, a Bayesian update strategy is employed to fuse prior and observational information at the frame-level to form a comprehensive probabilistic estimate of surface particle triggering. Specifically, for each spatial location... and time frame Define the frame-level integrated trigger probability In the formula, These are normalization coefficients used to ensure that the overall trigger probability at different spatial locations is within a certain range. Within the specified interval and satisfying the overall probability constraint, the above calculation achieves the coordinated fusion of mobilization prior, fluctuation posterior, and coherence weights, enabling the comprehensive triggering probability to simultaneously reflect the dynamic conditions, instantaneous response state, and local coupling weighting effects of particle initiation by airflow. To suppress single-frame noise interference and extract a stable triggering evolution process, a sliding window is set in the time dimension. Within a window (e.g., 5 to 10 frames) Time smoothing is performed. The time smoothing calculation can use a weighted moving average method to obtain the segment-level triggering results. ,in, It is a frame sequence number index in the time direction. This represents the number of frames contained within the sliding window, obtained after smoothing. This reflects the sustained probability distribution of surface particle-triggered activity within a specified time period under the influence of airflow disturbance. It can significantly reduce the uncertainty caused by short-term fluctuations in a single frame and highlight the continuity and stability of the disturbance process. Through the above steps, the segment-level comprehensive triggering results of surface particles during the entire rotor-surface interaction process can be obtained, providing a reproducible quantitative basis for subsequent surface disturbance area identification, risk classification, and dynamic model calibration.

[0099] The design purpose of steps S401 to S404 is to achieve comprehensive perception, fusion judgment, and spatiotemporal stability identification of particle triggering processes in the dynamic environment of UAV rotor-ground interaction. Step S401 combines the downwash dynamic pressure field with time-synchronized frictional wind information to calculate a frame-by-frame particle mobilization probability map, realizing the transformation of aerodynamic conditions into probability space. Step S402, based on the multi-scale perturbation observation results of the image, uses Bayesian inference to fuse the mobilization prior with the actual fluctuation response of the image to obtain the posterior probability distribution of particle-induced fluctuations, thereby establishing a statistical relationship between airflow perturbation and ground response in the image-ground correspondence. Step S403 further introduces the structure-aerodynamic comprehensive coupling coefficient as a weight to form a weighted coherence index, which is used to quantify the coupling strength between aerodynamic energy and particle response in a local area. Step S404 uses the particle mobilization probability, posterior probability, and coherence index as complementary evidence for Bayesian update and smooths it within a time sliding window to obtain a segment-level stable comprehensive triggering probability. Through this hierarchical fusion mechanism, the system can simultaneously consider airflow disturbance dynamics, surface response characteristics, and structural coupling effects, improving the accuracy and confidence level of identifying real particle initiation events. This process enables the solution to stably extract high-risk triggering areas even under complex conditions such as wind field fluctuations, noise interference, and multi-scale heterogeneous surfaces, achieving automatic and reliable determination of potential geological hazard points.

[0100] Step S405: Based on the frame-level comprehensive trigger probability, coupling evidence, and segment-level trigger results, jointly determine the trigger strength and coupling degree of each surface segment to obtain high-risk segments.

[0101] The observation area is divided into multiple segment units based on the surface coordinate system. Specifically, it is divided along the X-axis according to the spatial resolution Δx. A vertical dividing line, along the Y-axis, is used to divide the area according to a resolution Δy. A horizontal dividing line forms × Each surface segment unit is used as the determination object, and the frame-level comprehensive trigger probability is calculated separately. Coupling evidence and the segment-level triggering results obtained by smoothing through a sliding window. As the input for decision-making, the average value of the three indicators within each segment is calculated to obtain the average frame-level trigger probability of that segment. Average coupling strength and average segment-level trigger probability , in Represents a surface segment unit. Three types of judgment thresholds are set. , and The threshold is determined based on the 95% confidence upper limit of the statistical distribution of historical experimental samples or simulated data, and is used to distinguish between normal and significant disturbances. When the segment unit satisfies... , and When all three conditions are met Representing the specific first line, number The surface segment units of the column, It is the first line, number Average coupling strength within column segments, It is the first line, number The column section is Average frame-level trigger probability at time step It is the first line, number The column section is The average segment-level trigger probability at any given time indicates that the region has high instantaneous trigger intensity, significant air-to-ground coupling strength, and continuous triggering characteristics within the current time frame. Based on this, the segment is determined to be a high-risk segment and marked for output. If any condition is not met, it is determined to be a low-risk segment. Through the above joint determination, high-risk segment identification results with temporal and spatial resolution can be formed throughout the monitoring area, which can be used for subsequent model verification, disturbance source localization, and risk prevention and control strategy formulation.

[0102] Step S406: Calculate lateral micro-misalignment and micro-elevation commands based on high-risk sections, and limit the upper limits of pitch and roll change rates to reduce downwash dynamic pressure and micro-attitude disturbances.

[0103] Using the identified high-risk sections as control targets, lateral micro-misalignment commands are calculated based on the UAV's current attitude data and aerodynamic response model. With micro-increase command ,in Used to adjust for small lateral positional deviations of the aircraft. Used to correct vertical altitude position; by making minute adjustments to both, the relative aerodynamic coupling distance between the blade disk and the ground is changed, thereby reducing the peak downwash dynamic pressure above high-risk sections. While calculating this command, the pitch rate of change is also considered. With roll rate of change Set upper limit threshold and This is to prevent excessively rapid attitude response from causing new micro-disturbances or gain oscillations. The threshold is determined based on the aircraft's dynamic characteristic curve and the control system bandwidth, ensuring that control commands are executed within a safe linear range. This achieves simultaneous suppression of downwash pressure and attitude disturbances, and calculates lateral micro-misalignment commands based on the UAV's current attitude data and aerodynamic response model. With micro-increase command This method is used to fine-tune the aircraft's position to achieve aerodynamic coupling adjustment. This calculation method belongs to existing UAV attitude-displacement closed-loop control technology and can be implemented through proportional-integral-derivative control or model predictive control algorithms. The relevant implementation methods have been disclosed in existing literature and commercial flight control systems, and will not be elaborated upon here.

[0104] Step 5: Perform image plane motion compensation and rolling shutter micro-distortion correction on the image in the trigger section, and output the compensated and reconstructed image;

[0105] Step S501: Reacquire flight images of the trigger section for flight control adjustment, and perform adaptive deflickering and brightness normalization on the acquired images to obtain stable images;

[0106] After completing the flight control adjustment of the trigger section, the flight image of the trigger section is re-acquired with a re-stabilized flight attitude. Following step S2022, the image width is set based on the two-dimensional image pixel coordinate system. Height is Divide the image into For any pixel unit, The triggered segment images will be arranged in time frame order. Read and build frame sequence ,in Represents pixel coordinates, For the first The frame acquisition time. For each frame, calculate its global average brightness. This characterizes changes in illumination and flicker intensity using a sliding time window of fixed duration. For each unit, time-domain smoothing and amplitude normalization are performed on the brightness sequence of adjacent frames. Specifically, a weighted moving average is applied to the brightness values ​​of each frame within the time window, and amplitude normalization is performed using the maximum and minimum brightness values ​​of the window. This ensures that the output brightness sequence remains consistent within the dynamic range, thus obtaining the reference brightness benchmark for that window. This is used to compensate for short-term brightness fluctuations caused by minor camera attitude disturbances or changes in illumination, as well as short-term brightness fluctuations caused by minor airborne attitude disturbances, rotor flash reflections, or changes in the angle of sunlight incidence. In the flicker removal process, the frame sequence... High-frequency scintillation components Removal is achieved through adaptive filtering, the filter parameters of which are self-adjusted based on the local brightness change rate: when a brightness change rate is detected... Exceeding the set threshold When flickering occurs, the temporal filter weight is automatically increased to suppress flicker; when it falls below a threshold, detailed textures are preserved, thus balancing brightness stability and texture fidelity. This adaptive weight adjustment method is existing technology and will not be elaborated here. After flicker removal, brightness normalization is performed on each frame of the image. This process uses... Based on this, a normalization transformation is performed on the grayscale value of each pixel. This process unifies the brightness range of each frame of the image to the same baseline level, eliminating brightness deviations caused by ambient light differences or sensor response drift. A second spatial domain local equalization is performed on the normalized image sequence. A contrast-limited adaptive histogram equalization algorithm is used in this process to balance shadows and high-reflectivity areas, maintaining local contrast while ensuring overall image brightness continuity. Spatial domain local equalization is an existing technique and will not be elaborated upon here. After processing, a stable image sequence is obtained. The image exhibits stable brightness changes over time, without short-term flicker, and uniform spatial brightness distribution. This stabilized image serves as input data for subsequent step S502, which solves for the equivalent micro-pose timing and generates the image plane motion compensation field, ensuring photometric consistency and registration accuracy during image plane motion compensation and rolling shutter distortion correction. The brightness normalization and flicker suppression algorithms used in this step are conventional methods for image temporal stabilization, and those skilled in the art can directly implement them based on the parameter definitions and operation sequence.

[0107] Step S502: Using the equivalent micro-attitude perturbation spectrum and the main vibration modes of the structure as priors, solve the equivalent micro-attitude time series and generate the image plane motion compensation field;

[0108] After completing the flight control adjustment and updating the equivalent micro-attitude disturbance spectrum in the trigger zone, With structural vibration principal modes The frequency domain disturbance distribution was re-analyzed for the updated steady-state flight conditions. Based on... The time series of each attitude angle perturbation component is reconstructed from the amplitude spectrum and phase spectrum. ,in This is the inverse Fourier transform. Let be the phase function of the attitude perturbation. For the corresponding frame image or sampling time, The frequency is used to characterize the frequency distribution of the disturbance. Image plane coordinates are established based on the imaging geometry model. The response relationship between attitude perturbations is assumed to be... Let be the spatial vector of the image point, then its instantaneous displacement can be expressed as: ,in Let Jacobian matrix be the image plane matrix determined by focal length and field of view. and These are the lateral and longitudinal physical distances of the pixel center relative to the optical axis, typically expressed in pixel size or sensor length units. The dominant modes of structural vibration... Applying this to the response, we obtain the combined displacement field under structure-attitude coupling:

[0109] The first term reflects the equivalent micro-attitude perturbation spectrum. The energy response to a frequency- and time-domain consistent angular perturbation is given by a given term, and the second term represents the spatial deformation mapping of the structural modes. This method yields the results... This is the image plane motion compensation field, which can be used to dynamically correct the pixel positions of multi-frame sequences and realize image stabilization processing under attitude-structure coupling compensation.

[0110] Step S503: Perform image plane compensation and rolling shutter micro-distortion correction on the stabilized image based on the image plane motion compensation field, and output the compensated and reconstructed image for texture analysis.

[0111] Based on image plane motion compensation field Subsequently, the raw image sequence acquired during the stable flight phase was analyzed. Pixel-level compensation is performed frame by frame. This process uses the image plane coordinate system as the reference coordinate system to establish the compensated and stabilized image plane coordinates. This spatial position correction operation compensates for image point displacement caused by attitude disturbances and structural vibrations in each frame, thereby stabilizing the image plane. To eliminate the temporal micro-distortion introduced by the line exposure time difference in rolling shutter imaging mode, assuming the sensor exposes line by line from top to bottom, the exposure center time of a single line can be expressed as... ,in For the current row index, This is the start time of the exposure for that frame. Add line readout delay between adjacent scan lines. Combined with the image plane motion compensation field. Obtained continuous attitude response Calculate the average angular offset within each exposure interval. ,in This represents the single-line exposure duration. From this, the space-time correction mapping for rolling shutter distortion can be obtained: ,in The camera imaging Jacobian matrix describes the linear response of image point displacement to attitude perturbations. After the above two compensation steps, pixel grayscale is resampled using bilinear or higher-order spline interpolation methods to obtain a distortion-corrected and stabilized image. The image's attitude jitter and rolling shutter timing errors in the spatial domain are eliminated, and the image's geometric position and radiation characteristics are restored to the true imaging state corresponding to the flight path.

[0112] The output of this step As a stabilized image reconstruction sequence, it can be directly input into subsequent steps for operations such as texture extraction, brightness consistency analysis, and surface feature recognition, achieving high-precision remote sensing image geometric calibration and texture feature preservation.

[0113] Step 6: Perform spatial risk analysis on the compensated and reconstructed imagery to generate a spatial risk calibration map.

[0114] Step S601: Perform multi-scale structural tensor reconstruction on the compensated and reconstructed image to generate a texture robust field;

[0115] In the original pixel grid, for the image Calculate its in Direction (lateral direction) and The first-order gray-level gradient components in the longitudinal direction are denoted as follows: and The gradient calculation is implemented using a Gaussian derivative filter, that is:

[0116] ;

[0117] ; among which symbols This represents the convolution operation. The variance is The two-dimensional Gaussian smoothing kernel function, This is a smoothing scale parameter used to control the range of local gradient response and reduce the impact of noise. It is based on the gradient components at each pixel. Constructing local structure tensors Its definition is:

[0118] In the formula, For smoothing scale The Gaussian weighting function is used to perform a weighted average of the second moments of the gradient within the local neighborhood, thereby obtaining a smooth estimate of the texture direction information. The parameters... For the outer layer smoothing scale parameter, and This ensures robustness and directional continuity at texture boundaries. To achieve multi-scale feature fusion, in... Different smoothing scales Parallel computing ,in The multi-scale synthetic structure tensor is obtained by weighted superposition: ;in For the weights at each scale, satisfying This is used to balance the contribution of detail-scale and large-scale texture information. Perform eigenvalue decomposition, and let its eigenvalues ​​be respectively... and (in The corresponding unit eigenvector is and .in Reflecting the salience of the direction of local texture, This reflects the gradient energy intensity. Based on this, the texture robustness function is defined as follows:

[0119] ;in To prevent extremely small positive constants with a denominator of zero, and to maintain numerical stability; The range of values ​​is A larger value indicates a more stable texture orientation and a more pronounced local structure. Texture intensity and orientation information are combined to form a texture robust field. ;in The output is a multi-scale texture robust field dataset. Indicates robustness strength. and This represents the principal direction and orthogonal direction of the local texture of the pixel. Through this multi-scale structural tensor reconstruction process, the local directional structure and texture energy distribution of the image can be stably described at different scale levels. This texture robust field can be used as input features for subsequent texture feature calculation, land cover classification, anomaly detection, and image registration, achieving robust adaptation to illumination changes, noise disturbances, and imaging scale differences.

[0120] Step S602: Combining the immersion pressure field and particle mobilization probability, a hierarchical conditional random field segmentation is performed on the texture robust field to obtain suspected candidate regions;

[0121] After completing the flight control adjustments for the triggering section, the timing estimation results of the downwash dynamic pressure field in the aircraft aerodynamic coordinate system will be used. Combined with frame-level trigger probability Spatial and temporal alignment and normalization are performed. For any observation time... Position on the image raster Construct joint constraint vectors from corresponding pixels ,in For scale-invariant texture robustness, This refers to the instantaneous downwash dynamic pressure at that pixel. This represents the frame-level overall trigger probability. A hierarchical conditional random field model is established with all pixels as nodes. ,in Indicates the first Line number Column pixel nodes For its connection with neighboring nodes, This represents the index of a neighboring pixel node. Node potential function. according to Defined as a prior cost for determining whether a pixel belongs to an anomalous or stable state. in Here, is the dynamic pressure weighting coefficient, and log is the logarithmic function. This represents the globally averaged downwash dynamic pressure field, used to measure the degree of dynamic pressure deviation. (Boundary potential function) Based on the definition of texture similarity, it is used to maintain spatial continuity: ,in To smooth out the weights, For texture standard deviation, The indicator function is used. The total energy function is expressed as follows: ,in 0 represents a stable background, and 1 represents a suspected anomalous region. Optimization is achieved using Graph-Cut or Loopy-Belief Propagation algorithms. And in time dimension Slide upwards smoothly to obtain the optimal label. Candidate anomaly regions are obtained. , indicating at time The spatial region below simultaneously exhibits characteristics of dynamic pressure anomaly, high trigger probability, and texture mutation.

[0122] Step S603: Based on the suspected candidate regions, perform morphological connectivity detection and stability threshold determination to generate a spatial risk calibration map.

[0123] Based on suspected candidate regions To identify potential geological hazards, the area underwent verification: morphological connectivity detection was implemented, specifically using the eight-neighbor connectivity detection method to identify spatially connected pixel sets, with each set defined as a candidate block. And calculate its area. Used to characterize the spatial scale of anomaly regions; morphological factors are calculated to quantitatively describe the morphology of the mass. ,in Let be the perimeter of the outer boundary of the block, when A value close to 1 indicates a compact structure and regular shape. A lower value indicates that the region is elongated, fragmented, and unstable;

[0124] The average dynamic pressure deviation of the calculation block ,in Candidate block The number of pixels contained within. At the same time The baseline mean dynamic pressure over the entire observation area serves as a reference value to reflect the global stable field. This indicates the average deviation of the aerodynamic disturbance within the candidate block from the overall field; Average trigger probability at the frame level and morphological factors Common and preset threshold , and When comparing, and and When the candidate block is determined to be a potentially unstable element, the block that satisfies the condition is calculated. centroid coordinates As a source of geological hazard hazard points, a spatial risk mapping map is generated. This provides spatial constraints for subsequent risk assessment and geomorphological analysis, and presets thresholds. , and The critical judgment levels for characterizing dynamic pressure, trigger probability, and texture perturbation were determined through statistical distribution analysis and empirical calibration of historical sample data.

[0125] It should be understood that although the steps in the flowcharts of the various embodiments of the present invention are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the various embodiments may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least a portion of the sub-steps or stages of other steps.

[0126] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments described above. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), Rambus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0127] The foregoing description is illustrative of the invention and should not be construed as limiting it. Although several exemplary embodiments of the invention have been described, those skilled in the art will readily understand that many modifications can be made to the exemplary embodiments without departing from the novel teachings and advantages of the invention. Therefore, all such modifications are intended to be included within the scope of the invention as defined in the claims. It should be understood that the foregoing description is illustrative of the invention and should not be construed as limiting it to the specific embodiments disclosed, and modifications to the disclosed embodiments and other embodiments are intended to be included within the scope of the appended claims. The invention is defined by the claims and their equivalents.

Claims

1. A method for rapid identification of geological hazard hazard points based on UAV aerial surveying, characterized in that, include: Image sequences from airborne cameras, acceleration and angular velocity data from the inertial measurement unit, propeller speed, flight altitude, wind speed and direction, and airframe vibration signals were collected. Time-series samples were formed using image exposure time as a unified time axis. Based on these time-series samples, a priori on the momentum entrainment of the near-ground boundary layer during rotor downwash was established. Using propeller speed, flight altitude, and wind field data, time-series estimates of the downwash dynamic pressure field, particle mobilization probability time series, and equivalent micro-attitude perturbation spectrum were calculated. Subpixel-level optical flow and phase accumulation analysis were performed on the image sequences. Blind source separation was performed on the high-frequency components of airframe vibration and inertial measurement, and time cross-correlation was performed with the equivalent micro-attitude perturbation spectrum to obtain coupling evidence. Based on the time-series estimation of the downwash dynamic pressure field, the joint determination of coupling evidence and the triggering strength and coupling degree of the particle mobilization probability time series, an adaptive fine-tuning strategy for flight parameters is triggered. Image plane motion compensation and rolling shutter micro-distortion correction are performed on the image in the triggering segment, and the compensated and reconstructed image is output. Spatial risk analysis is performed on the compensated and reconstructed image to generate a spatial risk calibration map. The steps for obtaining coupling evidence are as follows: sub-pixel-level optical flow and phase accumulation analysis are performed on the image sequence to extract the multi-scale vector field of local motion and the corresponding phase energy density sequence, establishing a dynamic observable set of image fluctuations; blind source separation is performed on the high-frequency components of the body vibration signal and the inertial measurement unit to extract the main mode of structural vibration and the noise basis, forming a characterization of the body's self-vibration response; the dynamic observable set of image fluctuations and the characterization of the body's self-vibration response are cross-correlated and coherently tested with the equivalent micro-attitude perturbation spectrum on a unified time axis to obtain a comprehensive coupling coefficient. Based on the comprehensive coupling coefficient, it is determined that there is significant multi-scale time-frequency coupling, and coupling evidence is output.

2. The method for rapid identification of geological hazard hazard points based on UAV aerial surveying according to claim 1, characterized in that, The steps for calculating the time series estimation of the downwash dynamic pressure field are as follows: establish a priori momentum entrainment in the near-ground boundary layer of the rotor downwash, construct an axisymmetric-shear coupled dynamic pressure field parameterization based on propeller speed, flight altitude, and wind speed and direction, and solve the time series estimation of the downwash dynamic pressure field through energy conservation and momentum flux closure.

3. The method for rapid identification of geological hazard hazard points based on UAV aerial surveying according to claim 1, characterized in that, The calculation steps of the equivalent micro-attitude perturbation spectrum are as follows: by establishing a station-grid mapping, the surface roughness and zero displacement height parameters are extracted. Combined with wind direction, wind speed and surface zoning information, the critical wind erosion and shear thresholds of each particle size are calculated, and uncertainty correction is considered to obtain the particle mobilization probability time series. The downwash dynamic pressure field time series estimate is multiplied with the body / gimbal transfer matrix to obtain the equivalent torque perturbation. Combined with the maneuvering state of the inertial measurement unit, the equivalent micro-attitude perturbation spectrum is solved.

4. The method for rapid identification of geological hazard hazard points based on UAV aerial surveying according to claim 3, characterized in that, The calculation steps for the particle mobilization probability time series are as follows: The partition weight vector and partition label of each pixel, along with surface roughness parameters, zero displacement height parameters, and time-varying wind direction and near-surface wind information, are combined to weight the partition threshold parameters to obtain the critical wind erosion threshold and equivalent critical shear threshold for each particle size. A comprehensive uncertainty scale is generated based on turbulent gusts, surface parameter uncertainty, and partition heterogeneity. Based on the critical wind erosion threshold, equivalent critical shear threshold, and comprehensive uncertainty scale for each particle size, combined with time-varying frictional wind information, the particle mobilization probability for each particle size is calculated and weighted according to the local particle size distribution to obtain the particle mobilization probability time series.

5. The method for rapid identification of geological hazard hazard points based on UAV aerial surveying according to claim 4, characterized in that, The steps for obtaining the partition weight vector and partition label of each pixel are as follows: Taking the main line-of-sight projection point as the station, reference survey lines are laid out along the main wind direction and lateral direction respectively, and the near-ground wind information is unified to the same height reference. Combined with the surface type partition and micro-topographic undulation information, robust fitting and consistency checks are performed on the candidate roughness parameters, and the surface roughness parameters, zero displacement height parameters and their uncertainty range corresponding to the spatial location are output. A station-raster mapping index is established. Based on the station-raster mapping index, the image grid pixels are mapped to surface patches. For each pixel, the visible area ratio of each patch is calculated, and the ratio is corrected using attitude and spectral consistency to form the partition weight vector and partition label of each pixel.

6. The method for rapid identification of geological hazard hazard points based on UAV aerial surveying according to claim 1, characterized in that, The steps of the adaptive fine-tuning strategy for triggering flight parameters are as follows: using frame-level integrated trigger probability, coupling evidence, and segment-level trigger results, the trigger intensity and coupling degree of each surface segment are jointly determined to obtain high-risk segments; The lateral micro-misalignment and micro-elevation commands are calculated based on high-risk sections, and the upper limits of pitch and roll change rates are limited to reduce downwash dynamic pressure and micro-attitude disturbances.

7. The method for rapid identification of geological hazard hazard points based on UAV aerial surveying according to claim 6, characterized in that, The steps for obtaining the segment-level triggering result are as follows: the comprehensive coupling coefficient is used as the weight of the posterior probability distribution of particle-induced fluctuations in each frame to form a weighted coherence index; the frame-by-frame particle mobilization probability map, the posterior probability distribution of particle-induced fluctuations in each frame and the weighted coherence index are used as ternary evidence, and Bayesian update is used to obtain the frame-level comprehensive triggering probability, and the segment-level triggering result is generated by time smoothing within the sliding window.

8. The method for rapid identification of geological hazard hazard points based on UAV aerial surveying according to claim 7, characterized in that, The steps for obtaining the posterior probability distribution of particle-induced fluctuations in each frame are as follows: Based on the time series estimation of the downwash dynamic pressure field, the particle mobilization probability time series is resampled and spatially interpolated to generate a frame-by-frame particle mobilization probability map, and image-to-ground registration is completed with the image coordinates; the dynamic observable set of image fluctuations is mapped into a fine-grained disturbance intensity distribution map through line-of-sight projection and ground-to-shadow geometric relationship; using the frame-by-frame particle mobilization probability map as the prior, the observation likelihood function is constructed and posterior inference is performed to obtain the posterior probability distribution of particle-induced fluctuations in each frame.

9. A method for rapid identification of geological hazard hazard points based on UAV aerial surveying according to claim 1, characterized in that, The steps for outputting the compensated and reconstructed image are as follows: re-acquire flight control adjustment images of the trigger section, perform adaptive descintillation and brightness normalization on the acquired images to obtain a stable image; using the equivalent micro-attitude perturbation spectrum and the main mode of structural vibration as priors, solve the equivalent micro-attitude time series and generate an image plane motion compensation field; perform image plane compensation and rolling shutter micro-distortion correction on the stable image based on the image plane motion compensation field, and output the compensated and reconstructed image for texture analysis.

10. A method for rapid identification of geological hazard hazard points based on UAV aerial surveying according to claim 1, characterized in that, The steps for generating the spatial risk calibration map are as follows: after compensation and reconstruction, the image is reconstructed using a multi-scale structural tensor to generate a texture robust field; combined with the downwash dynamic pressure field and particle mobilization probability, the texture robust field is segmented using a hierarchical conditional random field to obtain suspected candidate regions; based on the suspected candidate regions, morphological connected component detection and stability threshold determination are performed to generate the spatial risk calibration map.

Images