Unmanned aerial vehicle threat assessment and defense decision method and system combined with environmental perception

Abstract

This application provides a method and system for drone threat assessment and defense decision-making that combines environmental perception, relating to the field of low-altitude security defense technology. It decouples environmental elements from real-time environmental perception data sets to obtain an electromagnetic propagation attenuation potential field and an atmospheric turbulence disturbance potential field. Based on geographical obstacle boundary contour information, a three-dimensional airspace mobility constraint boundary is established. Using this as a spatial constraint, the two potential fields are superimposed using multiphysics fields to generate multipath environmental masking effect distribution data for intruding drones. A passive detection feasibility assessment is performed on the preset flight paths of the intruding drones, generating threat detectability level identifiers and detection blind zone distribution boundary coordinates. Based on the threat detectability level identifiers and coordinates, a defense resource scheduling decision network is invoked to output defense resource collaborative deployment instructions. This invention achieves environmental perception-driven dynamic threat assessment and adaptive collaborative deployment of defense resources.

Description

Technical Field

[0001] This application relates to the field of low-altitude security defense technology, and more specifically, to a method and system for drone threat assessment and defense decision-making that combines environmental perception. Background Technology

[0002] With the rapid development of the low-altitude economy, the widespread application of drones in both civilian and military fields has led to increasingly severe threats to airspace security, especially unauthorized intrusion drones posing serious security risks to key areas. Existing drone threat assessment and defense decision-making technologies typically rely on fixed-deployment radar detection systems or electro-optical tracking equipment to detect and track intruding targets. These solutions have significant shortcomings when facing complex and ever-changing real-world airspace environments. Specifically, existing technologies often only consider the intruding drone's flight parameters and radar cross-section characteristics during threat assessment, neglecting the coupling effects of multiple environmental factors such as the electromagnetic spectrum, weather conditions, and geographical obstacles on the detection link and flight path. This results in a significant discrepancy between threat assessment results and actual detection effectiveness. At the defense decision-making level, existing solutions typically respond based on pre-set static defense resource deployment plans, lacking the ability to dynamically adjust the beam parameters of detection nodes and the deployment positions of interceptor units based on real-time environmental perception data. This prevents the full utilization of tactical advantages brought by environmental masking effects and makes it difficult to accurately deploy resources to fill detection blind spots, severely limiting the response efficiency and interception success rate of the defense system in complex environments. Summary of the Invention

[0003] In view of this, the purpose of this application is to provide a method and system for drone threat assessment and defense decision-making that combines environmental perception.

[0004] In conjunction with the first aspect of this application, a method for drone threat assessment and defense decision-making that incorporates environmental perception is provided, applied to a drone threat assessment and defense decision-making system that incorporates environmental perception, the method comprising: Acquire a set of real-time environmental perception data collected by an airspace monitoring sensor network. The set of real-time environmental perception data includes electromagnetic spectrum signal distribution information, meteorological fluid dynamic field distribution information, and geographic obstacle boundary contour information. The real-time environmental perception data set is subjected to environmental element decoupling processing to obtain the electromagnetic propagation attenuation potential field and the atmospheric turbulence disturbance potential field, and a three-dimensional airspace navigability constraint boundary is established based on the geographical obstacle boundary contour information. Using the three-dimensional airspace accessibility constraint boundary as a spatial constraint condition, the electromagnetic propagation attenuation potential field and the atmospheric turbulence disturbance potential field are superimposed with a multi-physics field to generate multipath environmental shielding effect distribution data of intruding UAVs within the three-dimensional airspace accessibility constraint boundary. Based on the multipath environmental masking effect distribution data, a passive detection feasibility assessment is performed on the preset flight path of the intrusion drone, generating a threat detectability level identifier and detection blind zone distribution boundary coordinates for the intrusion drone. Based on the threat detectability level identifier and the boundary coordinates of the detection blind zone distribution, the defense resource scheduling decision network is invoked to output a multi-dimensional defense resource collaborative deployment instruction. The multi-dimensional defense resource collaborative deployment instruction includes the beamforming parameter group of the detection node and the sequence of maneuvering and stationing waypoints of the interception unit.

[0005] In conjunction with the second aspect of this application, a drone threat assessment and defense decision-making system incorporating environmental awareness is provided. The drone threat assessment and defense decision-making system incorporating environmental awareness includes a machine-readable storage medium and a processor. The machine-readable storage medium stores machine-executable instructions. When the processor executes the machine-executable instructions, the drone threat assessment and defense decision-making system incorporating environmental awareness implements the aforementioned drone threat assessment and defense decision-making method incorporating environmental awareness.

[0006] In conjunction with a third aspect of this application, a computer-readable storage medium is provided, wherein computer-executable instructions are stored therein, and when the computer-executable instructions are executed, the aforementioned method for drone threat assessment and defense decision-making combined with environmental perception is implemented.

[0007] Combining any of the above aspects, firstly, by acquiring a real-time environmental perception data set containing electromagnetic spectrum signal distribution information, meteorological fluid dynamic field distribution information, and geographic obstacle boundary contour information collected by an airspace monitoring sensor network, and performing environmental element decoupling processing to obtain the electromagnetic propagation attenuation potential field and atmospheric turbulence disturbance potential field, a three-dimensional airspace mobility constraint boundary is established based on the geographic obstacle boundary contour information. Then, using the three-dimensional airspace mobility constraint boundary as a spatial constraint condition, multi-physics field superposition is performed on the electromagnetic propagation attenuation potential field and the atmospheric turbulence disturbance potential field to generate multi-path environmental masking effect distribution data of the intruding UAV within the constraint boundary. This allows for precise quantification of the comprehensive impact of environmental elements on the detection link and flight path. Based on this, the intruding UAV... The system performs a passive detection feasibility assessment on the aircraft's preset flight path, generating threat detectability level identifiers and detection blind zone distribution boundary coordinates. This allows the threat assessment results to directly reflect the detection effectiveness under real-world environmental conditions. More importantly, based on the threat detectability level identifiers and detection blind zone distribution boundary coordinates, the system invokes the defense resource scheduling decision network to output multi-dimensional defense resource collaborative deployment instructions, including beamforming parameter sets for detection nodes and maneuvering waypoint sequences for interceptor units. This eliminates reliance on static preset schemes for defense resource deployment. Instead, it dynamically optimizes beam pointing and interceptor node positions based on real-time environmental masking effects. This enables collaborative adaptive scheduling of detection and interception resources driven by environmental perception, significantly improving the accuracy of UAV threat assessment, the timeliness of defense decisions, and the interception success rate in complex airspace environments. Attached Figure Description

[0008] Figure 1 This application provides a flowchart illustrating the method for drone threat assessment and defense decision-making that combines environmental perception. Detailed Implementation

[0009] To make the objectives, technical solutions, and advantages of this application clearer, the application will be further described in detail below with reference to the accompanying drawings. The described embodiments should not be regarded as limitations on this application. All other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0010] In the following description, the terms "first, second, third" are used merely to distinguish similar objects and do not represent a specific ordering of objects. It is understood that "first, second, third" may be interchanged in a specific order or sequence where permitted, so that the embodiments of this application described herein can be implemented in an order other than that illustrated or described herein.

[0011] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. The terminology used herein is for the purpose of describing embodiments of this application only and is not intended to limit this application.

[0012] Figure 1 This illustration shows a flowchart of a drone threat assessment and defense decision-making method incorporating environmental awareness, provided in an embodiment of this application. It should be understood that in other embodiments, the order of some steps in this drone threat assessment and defense decision-making method incorporating environmental awareness may be shared based on actual needs, or some steps may be omitted or maintained. The detailed components of this drone threat assessment and defense decision-making method incorporating environmental awareness include: Example Implementation Section: This embodiment provides a method for drone threat assessment and defense decision-making that combines environmental perception. This method is applied to airspace security management scenarios involving threat assessment and defense resource deployment for intruding drones. In this embodiment, an airspace monitoring sensor network is deployed, consisting of multiple detection nodes and sensing devices, capable of collecting airspace environmental data in real time. The data collection activities of this airspace monitoring sensor network have been legally authorized by relevant management departments, and the collected data is used solely for security defense purposes and does not involve personal privacy information. For potentially sensitive data such as personal communication signals, desensitization processing is performed at the collection end using homomorphic encryption and differential privacy injection techniques to ensure that the original sensitive information is irreversible. Data storage employs attribute-based access control policies for permission isolation. The specific implementation of this method will be described in detail below with reference to the accompanying drawings.

[0013] As attached Figure 1 As shown, the method in this embodiment includes the following steps.

[0014] Step S110: Obtain the real-time environmental perception data set collected by the airspace monitoring sensor network. The real-time environmental perception data set includes electromagnetic spectrum signal distribution information, meteorological fluid dynamic field distribution information, and geographical obstacle boundary contour information.

[0015] In this embodiment, the airspace monitoring sensor network includes a spectrum monitoring receiver array, a Doppler weather radar station, a lidar scanning unit, and a three-dimensional geographic information system database. The spectrum monitoring receiver array scans the entire monitored airspace at a preset frequency scanning step size, recording the signal strength, phase, and estimated direction of arrival angle at each frequency point. At data acquisition time t, the electromagnetic spectrum signal distribution information output by the i-th receiving channel in the spectrum monitoring receiver array is a four-dimensional data tensor, denoted as E(x, y, z, f, t), where (x, y, z) represents the three-dimensional coordinates of the spatial sampling point, f represents the frequency value, and the value of E represents the electromagnetic signal power spectral density at frequency f at spatial location (x, y, z), with dimensions in watts per square meter per hertz. The meteorological hydrodynamic field distribution information is collected collaboratively by the Doppler weather radar station and wind speed and direction sensors deployed at different altitudes. The weather radar emits a pulse beam, receives the backscattered echoes from aerosol particles in the atmosphere, measures the Doppler frequency shift of the echoes, and inverts the radial wind speed component. The radial wind speed components from multiple radar stations are vector-synthesized to obtain three-dimensional wind vector field distribution information, denoted as W(x, y, z, t), with component form W(x, y, z, t) = [u(x, y, z, t), v(x, y, z, t), w(x, y, z, t)], where u, v, and w represent the eastward, northward, and skyward wind speed components, respectively, with dimensions in meters per second. Geographic obstacle boundary contour information is acquired through a lidar scanning unit. The lidar emits pulsed laser beams into the surrounding space, measures the time interval from emission to reflection from the obstacle surface back to the receiver, and calculates the three-dimensional coordinates of obstacle surface points by combining the laser beam's pointing angle. The coordinates of all surface points belonging to the same obstacle are aggregated to form a vertex coordinate sequence and triangular facet data for each obstacle, stored as geographic obstacle boundary contour information, denoted as O(x, y, z), where O(x, y, z) is a binary labeling function, taking a value of 1 for the space occupied by the obstacle and 0 otherwise.

[0016] Step S120: Perform environmental element decoupling processing on the real-time environmental perception data set to obtain the electromagnetic propagation attenuation potential field and the atmospheric turbulence disturbance potential field, and establish a three-dimensional airspace navigability constraint boundary based on the geographical obstacle boundary contour information.

[0017] Step S121: Perform spectral domain segmentation analysis on the electromagnetic spectrum signal distribution information, extract the signal intensity spatial gradient distribution and spectrum occupancy time evolution sequence in each spectrum segment, and perform time-frequency co-mapping of the signal intensity spatial gradient distribution and the spectrum occupancy time evolution sequence to generate an electromagnetic propagation attenuation state field.

[0018] The frequency range covered by the electromagnetic spectrum signal distribution information E(x, y, z, f, t) is from fa to fb. According to the preset frequency segmentation scheme, this frequency range is divided into Q spectrum segments, and the frequency range of the qth spectrum segment is [fqL, fqH].

[0019] For each spectral segment, the spatial gradient distribution of signal intensity within that segment is extracted as follows: Within the spectral segment q, the power spectral density value is integrated along the frequency dimension to obtain the spatial power distribution Pq(x, y, z, t), where Pq(x, y, z, t) = ∫E(x, y, z, f, t) df, with the lower limit of integration being fqL and the upper limit being fqH, and the dimension being watts per square meter. Next, the spatial derivatives of the spatial power distribution Pq(x, y, z, t) are calculated in the x, y, and z directions to obtain the spatial gradient distribution of signal intensity Gq(x, y, z, t), where Gq(x, y, z, t) = [∂Pq / ∂x, ∂Pq / ∂y, ∂Pq / ∂z], and the dimension being watts per cubic meter.

[0020] The method for extracting the spectral occupancy time evolution sequence is as follows: within the spectral segment q, the power spectral density value at each spatial sampling position (x, y, z) is continuously sampled in the time dimension. The sampling time window length is set to ΔT, containing K sampling moments, with the sampling moment index k ranging from 1 to K. For the spatial position (x, y, z), the time evolution sequence Sq(x, y, z, k) of the spectral segment q is defined as: at the k-th sampling moment within the time window, the power spectral density value within the frequency range [fqL, fqH] is integrated, Sq(x, y, z, k) = ∫E(x, y, z, f, tk)df, with a lower limit of integration of fqL and an upper limit of integration of fqH. Thus, the spectral occupancy time evolution sequence of the spectral segment q at the spatial position (x, y, z) is obtained as a one-dimensional array of length K.

[0021] The electromagnetic propagation attenuation state field is generated by performing a time-frequency co-mapping of the signal strength spatial gradient distribution Gq(x, y, z, t) with the spectrum occupancy time evolution sequence Sq(x, y, z, k). First, within each spectrum segment q, a time-frequency co-mapping kernel function Mq is constructed. At sampling time tk, Mq(x, y, z, tk) = ||Gq(x, y, z, tk)|| × (Sq(x, y, z, k+1) - Sq(x, y, z, k)) / Δt, where ||Gq|| = √((∂Pq / ∂x) 2 +(∂Pq / ∂y) 2 +(∂Pq / ∂z) 2), where Δt is the time interval between adjacent sampling times. Then, the time-frequency co-mapping kernel functions Mq of all Q spectral segments are spatially weighted and superimposed. The superposition weight coefficient αq is determined according to the frequency interval between the center frequency of each spectral segment and the center frequency of the frequency band to be evaluated. The final generated electromagnetic propagation attenuation state field H1(x,y,z,t) is: H1(x,y,z,t)=∑(q=1 to Q)αq×Mq(x,y,z,t), with the dimension of watts per cubic meter per second. This electromagnetic propagation attenuation state field H1(x,y,z,t) characterizes the distribution of the comprehensive attenuation degree of electromagnetic waves propagating in the spatial domain due to atmospheric absorption, scattering, and spectrum occupancy competition in the four-dimensional space-time.

[0022] Step S122: Perform wind vector field curl divergence decomposition on the meteorological fluid dynamic field distribution information, extract the spatiotemporal distribution of the rotation component and the spatiotemporal distribution of the divergence component of the wind vector field, and input the spatiotemporal distribution of the rotation component and the spatiotemporal distribution of the divergence component into the atmospheric refractive index structure constant field reconstruction step to generate the atmospheric turbulent disturbance potential field.

[0023] The wind vector field curl and divergence decomposition is performed on the meteorological fluid dynamic field distribution information W(x, y, z, t) = [u(x, y, z, t), v(x, y, z, t), w(x, y, z, t)]. The curl decomposition is calculated as follows: curl vector R(x, y, z, t) = ∇ × W = [(∂w / ∂y - ∂v / ∂z), (∂u / ∂z - ∂w / ∂x), (∂v / ∂x - ∂u / ∂y)], with dimensions per second. The divergence decomposition is calculated as follows: divergence scalar D(x, y, z, t) = ∇·W = ∂u / ∂x + ∂v / ∂y + ∂w / ∂z, with dimensions per second. By continuously sampling the curl vector R(x, y, z, t) in the time dimension, the spatiotemporal distribution of the rotation component R(x, y, z, tk) is obtained. By continuously sampling the divergence scalar D(x, y, z, t) in the time dimension, the spatiotemporal distribution of the divergence component D(x, y, z, tk) is obtained.

[0024] The spatiotemporal distributions of the rotational component R and the divergent component D are input into the atmospheric refractive index structure constant field reconstruction step to generate the atmospheric turbulent perturbation potential field H2(x, y, z, t). The atmospheric refractive index structure constant field reconstruction step adopts the following calculation process: First, the turbulent kinetic energy TKE of the wind vector field is calculated, TKE(x, y, z, t) = 0.5 × (u′ 2 +v′ 2 +w′ 2), where u′=u-ū, v′=v-ū, w′=w-ū, and ū is the time-averaged wind speed. The calculation window length for time averaging is set according to the time scale of the atmospheric turbulence inertial sub-region. Next, the temperature fluctuation structure function is calculated. Combining the empirical relationship between atmospheric refractive index and temperature and humidity, the temperature fluctuation structure function is converted into the atmospheric refractive index structure constant Cn. 2 Cn 2 The relationship between (x, y, z, t) and TKE is established using the Monin-Obukhov similarity theory, Cn 2 (x, y, z, t) = β × TKE(x, y, z, t) × L0^(4 / 3), where β is the empirical scaling factor and L0 is the external scale of the turbulence. Then, according to Cn 2 The spatial distribution of atmospheric turbulence is calculated, and the mean square value of the phase perturbation of electromagnetic wave propagation by atmospheric turbulence is determined, thereby generating the atmospheric turbulence perturbation potential field H2(x, y, z, t). H2(x, y, z, t) = γ × Cn 2 (x, y, z, t) × Δz, where γ is a coefficient related to the wavelength of the electromagnetic wave, and Δz is the infinitesimal length of the propagation path. The dimension of H2 is the dimensionless phase variance.

[0025] Step S123: Perform 3D convex hull reconstruction on the geographic obstacle boundary contour information, extract the vertex coordinate sequence and face normal vector set of each obstacle in the geographic obstacle boundary contour information, and stitch the vertex coordinate sequence and face normal vector set into a triangular mesh to generate a 3D spatial traversability constraint boundary.

[0026] Extract the set of coordinates of all spatial locations with a value of 1 from the geographic obstacle boundary contour information O(x, y, z). Perform a 3D convex hull reconstruction algorithm on this set of coordinates to generate a triangular mesh for the convex hull surface of each obstacle. For each obstacle, its convex hull surface consists of M triangular faces, each defined by 3 vertex coordinates and 1 unit outward normal vector. Extract the vertex coordinate sequence V={(x1, y1, z1), (x2, y2, z2), ..., (xN, yN, zN)}, where N is the total number of vertices of the obstacle's convex hull. Extract the set of facet normal vectors N={n1, n2, ..., nM}, where the j-th facet normal vector nj=[nxj, nyj, nzj], with a magnitude of 1, and pointing outward from the obstacle's convex hull. Concatenate the vertex coordinate sequence V and the facet normal vector set N into a triangular mesh to generate the 3D spatial traversability constraint boundary B(x, y, z). B(x, y, z) is a three-dimensional binary function that takes the value 1 inside the space enclosed by the convex hull of the obstacles and takes the value 0 in the traversable space outside the convex hull of the obstacles. The union of the three-dimensional traversability constraint boundaries of all obstacles constitutes the traversability constraint boundary of the entire monitored space.

[0027] Step S124: Using the electromagnetic propagation attenuation potential field as the first physical field input data, the atmospheric turbulence disturbance potential field as the second physical field input data, and the three-dimensional spatial mobility constraint boundary as the spatial filtering condition, the first physical field input data and the second physical field input data are superimposed on a voxel-by-voxel basis to generate a superimposed physical field data set.

[0028] The passable space outside the obstacle defined by the three-dimensional spatial passability constraint boundary B(x, y, z) is taken as the effective voxel space Ω, Ω = {(x, y, z) | B(x, y, z) = 0}. For each voxel position (x, y, z) within the effective voxel space Ω, the electromagnetic attenuation feature value at that position is extracted from the electromagnetic propagation attenuation potential field H1(x, y, z, t), and the atmospheric turbulence disturbance feature value at that position is extracted from the atmospheric turbulence disturbance potential field H2(x, y, z, t). The voxel-by-voxel superposition method is as follows: at each voxel position (x, y, z), a superimposed physical field feature vector F(x, y, z, t) = [H1(x, y, z, t), H2(x, y, z, t)] is constructed. This vector is a two-dimensional array, with the first dimension storing the electromagnetic attenuation feature value and the second dimension storing the atmospheric turbulence disturbance feature value. By traversing all voxels within the effective voxel space Ω, we obtain the superimposed physical field data set FΩ = {F(x, y, z, t) | (x, y, z) ∈ Ω}. The superimposed physical field data set FΩ is a data array of shape [NΩ, 2], where NΩ is the total number of voxels within the effective voxel space Ω, and 2 indicates that each voxel position is associated with two physical field eigenvalues.

[0029] Step S125: Extract the spatial gradient direction field from the superimposed physical field data set to obtain the electromagnetic attenuation gradient direction and atmospheric refraction gradient direction at each voxel position. Spatially register the electromagnetic attenuation gradient direction and the atmospheric refraction gradient direction, analyze the spatial angle relationship and relative change trend between the two, and generate composite situation distribution data reflecting the combined influence of the atmospheric environment on electromagnetic propagation.

[0030] Spatial gradient calculation is performed on the electromagnetic attenuation eigenvalue H1(x,y,z,t) at each voxel position (x,y,z) in the superimposed physical field data set FΩ, yielding the electromagnetic attenuation gradient direction vector g1(x,y,z,t)=[∂H1 / ∂x, ∂H1 / ∂y, ∂H1 / ∂z]. Spatial gradient calculation is also performed on the atmospheric turbulence disturbance eigenvalue H2(x,y,z,t) at the same voxel position, yielding the atmospheric refraction gradient direction vector g2(x,y,z,t)=[∂H2 / ∂x, ∂H2 / ∂y, ∂H2 / ∂z]. Spatial registration is achieved by unifying the two gradient direction vectors to the same spatial coordinate system. In this embodiment, they are naturally located in the same three-dimensional spatial coordinate system, requiring no additional coordinate transformation. Calculate the spatial angle θ(x, y, z, t) between the electromagnetic attenuation gradient direction g1 and the atmospheric refraction gradient direction g2, θ = arccos((g1·g2) / (‖g1‖×‖g2‖)), where g1·g2 = g1x×g2x + g1y×g2y + g1z×g2z, and ‖g1‖ = √(g1x × g2x + g1y × g2y + g1z × g2z). 2 +g1y 2 +g1z 2 ), ||g2||=√(g2x 2 +g2y 2 +g2z 2 Analyze the relative trends of the two, and calculate the ratio of the directional derivatives of g1 and g2, ρ(x, y, z, t) = ||g1|| / (||g2||+ε), where ε is a preset minimum positive number to prevent the denominator from being zero. Generate composite situational distribution data C(x, y, z, t) reflecting the combined influence of the atmospheric environment on electromagnetic propagation, C(x, y, z, t) = [θ(x, y, z, t), ρ(x, y, z, t)]. This data is represented by a two-dimensional vector at each voxel position, with the first dimension being the spatial angle and the second dimension being the ratio of the directional derivatives.

[0031] Step S126: Detect the path intersection of the environmental shading effect direction vector distribution data and the preset flight path of the intruding drone. Take each intersection point of the preset flight path and the environmental shading effect direction vector distribution data as a sampling point, extract the magnitude of the environmental shading effect direction vector at each sampling point, and generate a shading effect intensity sequence along the path.

[0032] In the aforementioned composite situational distribution data C(x, y, z, t), the environmental masking effect direction vector is defined as V(x, y, z, t) = ρ(x, y, z, t) × g1(x, y, z, t). The preset flight path of the intruding UAV consists of a sequence of R waypoints, denoted as Path = {(xp1, yp1, zp1), (xp2, yp2, zp2), ..., (xpR, ypR, zpR)}. A continuous path segment is generated between adjacent waypoints using linear interpolation, and this continuous path segment is discretized into S equally spaced sampling points. For each sampling point s, its spatial coordinates are (xs, ys, zs). The environmental masking effect direction vector Vs at that sampling point is extracted from V(x, y, z, t) using trilinear interpolation. The magnitude of Vs is calculated as ms = ||Vs|| = √(Vsx) 2 +Vsy 2 +Vsz 2 Arrange the magnitude values ​​ms of the S sampling points sequentially from front to back along the path to generate the shading effect intensity sequence L1={m1, m2, ..., mS} along the path.

[0033] Step S130: Using the three-dimensional airspace mobility constraint boundary as a spatial constraint condition, superimpose the electromagnetic propagation attenuation potential field and the atmospheric turbulence disturbance potential field using a multi-physics field to generate multipath environmental masking effect distribution data of the intruding UAV within the three-dimensional airspace mobility constraint boundary.

[0034] Step S131: Using the three-dimensional spatial mobility constraint boundary as a spatial filtering condition, extract the data inside the electromagnetic propagation attenuation state field located inside the three-dimensional spatial mobility constraint boundary from the electromagnetic propagation attenuation state field, and use the data inside the electromagnetic propagation attenuation state field as the first physical field valid data.

[0035] The three-dimensional spatial mobility constraint boundary B(x, y, z) divides the entire monitoring space into two parts: a traversable space Ω and an obstacle-occupied space. Using B(x, y, z) = 0 as the spatial filtering condition, all voxel positions in the electromagnetic propagation attenuation state field H1(x, y, z, t) are traversed, and only the H1 values ​​at voxel positions that satisfy B(x, y, z) = 0 are retained, forming the first physical field effective data H1eff = {H1(x, y, z, t) | (x, y, z) ∈ Ω}.

[0036] Step S132: Using the three-dimensional spatial mobility constraint boundary as a spatial filtering condition, extract the data inside the atmospheric turbulence disturbance potential field located inside the three-dimensional spatial mobility constraint boundary from the atmospheric turbulence disturbance potential field, and use the data inside the atmospheric turbulence disturbance potential field as the second physical field effective data.

[0037] Using B(x, y, z) = 0 as the spatial filtering condition, all voxel positions in the atmospheric turbulent perturbation potential field H2(x, y, z, t) are traversed, and only the H2 values ​​at the voxel positions that satisfy B(x, y, z) = 0 are retained to form the effective data of the second physical field H2eff = {H2(x, y, z, t)|(x, y, z)∈Ω}.

[0038] Step S133: Perform spatial registration and alignment processing on the first physical field effective data and the second physical field effective data. Based on the preset physical field coupling mapping relationship, analyze the degree of interaction between electromagnetic propagation attenuation characteristics and atmospheric turbulence disturbance characteristics at each voxel position, and generate a comprehensive environmental interference situation data set.

[0039] The effective data H1eff of the first physics field and the effective data H2eff of the second physics field naturally share the same spatial sampling grid. Spatial registration is completed by one-to-one correspondence of voxel indices. A preset physics field coupling mapping relationship defines a coupling function Φ, which maps the electromagnetic propagation attenuation value h1 and the atmospheric turbulence disturbance value h2 to a comprehensive environmental disturbance value J. Φ(h1, h2) = a1×h1 + a2×h2 + a3×h1×h2, where a1, a2, and a3 are preset coupling coefficients. a1 and a2 represent the linear contributions of electromagnetic propagation attenuation and atmospheric turbulence disturbance to the comprehensive environmental disturbance, respectively, and a3 represents the weight of the nonlinear interaction between the two. For each voxel position (x, y, z) within Ω, calculate J(x, y, z, t) = Φ(H1eff(x, y, z, t), H2eff(x, y, z, t)). Traverse all voxel positions to generate a comprehensive environmental disturbance situation data set JΩ={J(x,y,z,t)|(x,y,z)∈Ω}.

[0040] Step S134: Based on the comprehensive environmental interference situation data set, extract the spatial variation characteristics of electromagnetic propagation attenuation characteristics and atmospheric turbulence disturbance characteristics respectively. By analyzing the spatial correlation and relative orientation relationship between the two, generate comprehensive distribution data of environmental shielding effect.

[0041] From the effective data of the first physical field H1eff corresponding to the comprehensive environmental interference situation data set JΩ, the spatial variation characteristics of electromagnetic propagation attenuation are extracted. This characteristic is characterized by the spatial autocorrelation function R1(Δd) of H1eff within Ω, R1(Δd) = E[(H1eff(x, y, z, t) - μ1) × (H1eff(x + Δx, y + Δy, z + Δz, t) - μ1)], where μ1 is the spatial mean of H1eff within Ω, and Δd = √(Δx 2 +Δy 2 +Δz 2The spatial lag distance is denoted as Ψ(x, y, z, t). Spatial variation characteristics of atmospheric turbulence disturbances are extracted from H2eff, characterized by the spatial autocorrelation function R2(Δd) of H2eff within Ω, calculated in the same way as R1. The spatial correlation between H1eff and H2eff is analyzed, and the spatial cross-correlation function R12(Δd) is calculated: R12(Δd) = E[(H1eff(x, y, z, t) - μ1) × (H2eff(x + Δx, y + Δy, z + Δz, t) - μ2)], where μ2 is the spatial mean of H2eff. The relative orientation relationship between the two is analyzed, and the angle between the gradient directions of H1eff and H2eff at each voxel position is calculated to obtain the directional consistency field Ψ(x, y, z). The comprehensive distribution data of the environmental shading effect is generated as Z(x, y, z, t) = [J(x, y, z, t), R12(Δd), Ψ(x, y, z)]. This data contains three components at each voxel location: the comprehensive environmental disturbance value, the spatial cross-correlation function value, and the directional consistency value.

[0042] Step S135: Based on the comprehensive distribution data of the environmental masking effect and the preset flight path of the intruding drone, perform path intersection detection, take each intersection point of the preset flight path and the distribution data of the direction vector of the environmental masking effect as a sampling point, and extract the magnitude and direction angle of the direction vector of the environmental masking effect at each sampling point.

[0043] The preset flight path of the intruding drone is sampled at equal intervals to obtain S sampling points. For each sampling point s, its spatial coordinates are (xs, ys, zs). The environmental masking effect direction vector Zs at that sampling point is extracted from the comprehensive distribution data Z of the environmental masking effect using trilinear interpolation. The magnitude of the environmental masking effect direction vector ms = ||Zs||. The direction angle φs of the environmental masking effect direction vector is determined by the projection direction angles of Zs in the horizontal and vertical planes, φs = [arctan(Zsy / Zsx), arctan(Zsz / √(Zsx)], [arc ... 2 +Zsy 2 ))).

[0044] Step S136: Arrange the magnitude values ​​of the environmental shading effect direction vectors at each sampling point in order along the preset flight path to generate a sequence of shading effect intensity along the path; arrange the direction angles of the environmental shading effect direction vectors at each sampling point in order along the preset flight path to generate a sequence of shading effect direction angles along the path.

[0045] Arrange the modulus values ms of the S sampling points in order from front to back along the path direction, and generate a path-wise shielding effect intensity sequence L2={m1, m2, ..., mS}. Arrange the direction angles φs of the S sampling points in order from front to back along the path direction, and generate a path-wise shielding effect direction angle sequence A={φ1, φ2, ..., φS}, where each φs is a two-dimensional vector comprising a horizontal direction angle and a vertical direction angle.

[0046] Step S137: Based on the path-wise shielding effect intensity sequence and the path-wise shielding effect direction angle sequence, perform two-dimensional marking of shielding effect amplitude and shielding effect direction on each sampling point on the preset flight path, so as to generate multipath environment shielding effect distribution data.

[0047] For each sampling point s on the preset flight path, construct a two-dimensional mark Ds=[ms, φs]. Combine the two-dimensional marks of all S sampling points on the path in the path order, so as to generate multipath environment shielding effect distribution data MP={D1, D2, ..., DS}. MP is a list with a length of S, and each list element is a three-dimensional vector comprising a modulus value and a direction angle.

[0048] Step S140: Perform passive detection feasibility assessment on the preset flight path of the intruding unmanned aerial vehicle according to the multipath environment shielding effect distribution data, so as to generate a threat detectability level identifier and detection blind area distribution boundary coordinates of the intruding unmanned aerial vehicle.

[0049] Step S141: Extract the path-wise shielding effect intensity sequence and the path-wise shielding effect direction angle sequence from the multipath environment shielding effect distribution data, perform point-by-point comparison between the quantized value of the comprehensive interference intensity of each sampling point in the path-wise shielding effect intensity sequence and a preset detection sensitivity threshold, and mark a sampling point with a shielding effect intensity lower than the preset detection sensitivity threshold as a detection blind area candidate point.

[0050] Extract the path-wise shielding effect intensity sequence L2={m1, m2, ..., mS} from the multipath environment shielding effect distribution data MP. The preset detection sensitivity threshold is set to Td, and the threshold is determined according to the minimum detectable signal power of the receiver of the detection node. Perform point-by-point comparison on each element ms in L2, and if ms<Td, mark the sampling point s as a detection blind area candidate point. The indices of all detection blind area candidate points constitute a set C1.

[0051] Step S142: Perform point-by-point included angle calculation on the dominant interference direction feature of each sampling point in the path-wise shielding effect direction angle sequence and the preset detection beam coverage direction range, and mark a sampling point with a shielding effect direction angle beyond the preset detection beam coverage direction range as a direction unreachable candidate point.

[0052] The masking effect direction angle sequence A = {φ1, φ2, ..., φS} along the path is extracted from the multipath environmental masking effect distribution data MP. The preset detection beam coverage direction range is defined by the azimuth range [αmin, αmax] and the elevation range [βmin, βmax]. For each sampling point s, the horizontal direction angle αs and the vertical direction angle βs are extracted from φs. If αs is not within the interval [αmin, αmax], or βs is not within the interval [βmin, βmax], then the sampling point s is marked as a candidate point of unreachable direction. The indices of all candidate points of unreachable direction are used to form a set C2.

[0053] Step S143: Perform spatial intersection analysis on the candidate points of the detection blind zone and the candidate points of the direction unreachable, and extract the sampling points that simultaneously meet the conditions of the candidate points of the detection blind zone and the conditions of the candidate points of the direction unreachable, as the set of effective points of the detection blind zone.

[0054] Find the intersection of sets C1 and C2 to obtain a set of sampling point indices C3 = C1 ∩ C2 that simultaneously satisfy the conditions for candidate points in the detection blind zone and the condition for candidate points that are unreachable in the direction. The sampling point corresponding to each index in C3 is the effective point in the detection blind zone. These sampling points have a signal strength below the detection threshold and are outside the coverage area of ​​the detection beam in the direction, so the detection node cannot effectively detect the intruding UAV. The set of effective points in the detection blind zone is denoted as Pblind.

[0055] Step S144: Perform spatial connectivity analysis on the set of effective points in the detection blind zone, aggregate the effective points in the detection blind zone that are spatially adjacent into a connected domain of the detection blind zone, extract the boundary coordinate sequence of each connected domain of the detection blind zone, and generate the boundary coordinates of the detection blind zone distribution.

[0056] For the sampled points in the effective point set Pblind of the detection blind zone, connected component analysis is performed based on their spatial coordinates. A three-dimensional connected component labeling algorithm is used, setting the connected neighborhood radius to r. Two effective points in the detection blind zone with a spatial distance less than r are considered adjacent points. A breadth-first search is used to traverse all effective points in the detection blind zone, aggregating interconnected points into a single connected component for the detection blind zone. For each connected component, the coordinates of its boundary points are extracted. The criteria for determining a boundary point are that it belongs to the connected component but at least one of its neighboring points does not. The boundary point coordinates are sorted according to spatial adjacency to generate the boundary coordinate sequence of the connected component for that detection blind zone. The boundary coordinate sequences of all connected components in the detection blind zone constitute the detection blind zone distribution boundary coordinates Bblind.

[0057] Step S145: The proportion of the number of valid points in the detection blind zone to the total number of sampling points on the preset flight path is taken as the detection blind zone coverage rate. A threat detectability level identifier is generated according to the preset coverage rate interval into which the detection blind zone coverage rate falls.

[0058] The detection blind zone coverage rate Rb = |Pblind| / S is calculated, where |Pblind| is the number of elements in the effective point set of the detection blind zone, and S is the total number of sampling points on the preset flight path. The correspondence between the preset coverage rate intervals and threat detectability level identifiers is as follows: if 0 ≤ Rb < η1, the threat detectability level identifier is L1, indicating high detectability; if η1 ≤ Rb < η2, the threat detectability level identifier is L2, indicating medium detectability; if η2 ≤ Rb < η3, the threat detectability level identifier is L3, indicating low detectability; if Rb ≥ η3, the threat detectability level identifier is L4, indicating extremely low detectability. η1, η2, and η3 are the preset coverage rate thresholds.

[0059] Step S150: Based on the threat detectability level identifier and the boundary coordinates of the detection blind zone distribution, the defense resource scheduling decision network is invoked to output a multi-dimensional defense resource collaborative deployment instruction. The multi-dimensional defense resource collaborative deployment instruction includes a beamforming parameter group for detection nodes and a sequence of maneuvering and stationing waypoints for interception units.

[0060] Step S151: Input the threat detectability level identifier into the threat level embedding layer of the defense resource scheduling decision network, map the threat detectability level identifier into a threat level embedding vector, input the detection blind zone distribution boundary coordinates into the spatial coordinate encoding layer of the defense resource scheduling decision network, and map the detection blind zone distribution boundary coordinates into a spatial coordinate encoding vector.

[0061] The defense resource scheduling decision network is a deep neural network model. This network includes a threat level embedding layer, a spatial coordinate encoding layer, a feature cross-fusion layer, a detection resource allocation branch, and an interception resource allocation branch. The threat level embedding layer is a trainable embedding matrix E with dimensions Nlevel × Demb, where Nlevel is the total number of threat detectability level identifiers, and Demb is the dimension of the embedding vector. Using the threat detectability level identifier L as an index, the vector in the Lth row of the embedding matrix E is extracted as the threat level embedding vector eL, with dimension Demb. The spatial coordinate encoding layer is a multilayer perceptron structure containing two fully connected hidden layers, each followed by a ReLU activation function. The boundary coordinates of the detection blind zone distribution, Bblind, are flattened into a one-dimensional vector, input to the spatial coordinate encoding layer, and after forward propagation, output a spatial coordinate encoding vector eS, with dimension Denc.

[0062] Step S152: Input the threat level embedding vector and the spatial coordinate encoding vector into the feature cross-fusion layer of the defense resource scheduling decision network, perform feature fusion and correlation analysis on the threat level embedding vector and the spatial coordinate encoding vector, and generate a joint threat representation vector containing spatial distribution characteristics.

[0063] The feature cross-fusion layer receives the threat level embedding vector eL and the spatial coordinate encoding vector eS as input. First, eL and eS are concatenated along the feature dimensions to obtain the concatenated vector econcat = [eL, eS], with dimensions Demb + Denc. Then, the concatenated vector econcat is input into a feature fusion subnetwork containing multiple fully connected layers. The first layer weight matrix W1 of this subnetwork has dimensions (Demb + Denc) ​​× Dh, the second layer weight matrix W2 has dimensions Dh × Dh, and the third layer weight matrix W3 has dimensions Dh × Djoint, where Dh is the hidden layer dimension and Djoint is the dimension of the joint threat representation vector. Each fully connected layer is followed by a batch normalization layer and a ReLU activation function. The forward propagation outputs the joint threat representation vector ejoin, which contains spatial distribution features and has dimensions Djoint.

[0064] Step S153: Input the joint threat characterization vector into the detection resource allocation branch of the defense resource scheduling decision network, and output the detection node beamforming parameter set, which includes the beam main lobe pointing angle, beamwidth range and transmit power gain coefficient of each detection node.

[0065] The detection resource allocation branch is a fully connected output head that takes the joint threat representation vector ejoin as input and outputs the beamforming parameter set Pbeam for each detection node. Assuming Nd detection nodes are deployed, the output layer dimension of the detection resource allocation branch is Nd × Nparam, where Nparam is the number of beamforming parameters for each detection node. The beamforming parameters for the i-th detection node include the beam main lobe pointing azimuth angle αi, the beam main lobe pointing elevation angle βi, the beam horizontal width Δαi, the beam vertical width Δβi, and the transmit power gain coefficient gi. Pbeam = {(α1, β1, Δα1, Δβ1, g1), (α2, β2, Δα2, Δβ2, g2), ..., (αNd, βNd, ΔαNd, ΔβNd, gNd)}.

[0066] Step S154: Input the threat joint characterization vector into the interception resource allocation branch of the defense resource scheduling decision network, and output the interception unit maneuvering and stationing waypoint sequence. The interception unit maneuvering and stationing waypoint sequence includes the three-dimensional coordinates of the stationing point and the stationing duration of each interception unit.

[0067] The intercept resource allocation branch is another fully connected output head, taking the joint threat representation vector ejoin as input and outputting the intercept unit maneuvering and dwelling waypoint sequence Tintercept. Assume there are Ni intercept units deployed, each allocated Nw dwelling waypoints. The output layer dimension of the intercept resource allocation branch is Ni×Nw×4, where 4 represents the four parameters of each waypoint: x-coordinate, y-coordinate, z-coordinate, and dwell time τ. The maneuvering and dwelling waypoint sequence of the j-th intercept unit is Tj={(xj1, yj1, zj1, τj1), (xj2, yj2, zj2, τj2), ..., (xjNw, yjNw, zjNw, τjNw)}. Tintercept={T1, T2, ..., TNi}.

[0068] Step S155: Synchronize and align the beamforming parameter group of the detection node and the maneuvering waypoint sequence of the interception unit in time to generate a multi-dimensional defense resource collaborative deployment command.

[0069] Time-stamping is applied to the parameters of each detector node in the beamforming parameter group Pbeam, setting the start time tbeam_start and the effective period Tbeam for each parameter. An arrival time stamp calculated using preset speed constraints is appended to each waypoint in the maneuvering waypoint sequence Tintercept of the interceptor unit. The effective time windows of the parameters in the detector node beamforming parameter group are aligned with the arrival time windows of each waypoint in the maneuvering waypoint sequence of the interceptor unit, establishing a temporal correspondence based on the same absolute time reference. The time-aligned Pbeam and Tintercept are encapsulated into a multi-dimensional defense resource collaborative deployment instruction Cmd={Pbeam, Tintercept, tsync}, where tsync is the synchronization timestamp.

[0070] The training process of the defense resource scheduling decision network is implemented in steps S210 to S230.

[0071] Step S210: Construct a training sample dataset. Collect historical airspace monitoring data. Following the methods in steps S110 to S145, extract threat detectability level identifiers and detection blind zone distribution boundary coordinates as training input features. Simultaneously, construct training labels based on manually labeled or simulated optimal deployment schemes of defense resources from historical defense missions. The training labels contain the optimal detection node beamforming parameter set and the optimal interceptor unit maneuvering and stationing waypoint sequence corresponding to the input features. The training sample dataset contains U training samples, where the input feature of the u-th training sample is Xu=(Lu, Bu), and the label is Yu=(Pbeam_u, Tintercept_u).

[0072] Step S220: Construct the loss function and configure the training hyperparameters. The loss function L is composed of a weighted sum of the detection resource allocation loss Ldet and the interception resource allocation loss Lint, L = λ1 × Ldet + λ2 × Lint, where λ1 and λ2 are weight coefficients. Ldet uses the mean square error between the predicted beamforming parameter set and the optimal beamforming parameter set. Lint uses the weighted dynamic time warping distance between the predicted waypoint sequence and the optimal waypoint sequence. The Adam optimizer is selected for training, with an initial learning rate of lr0, a batch size of Bs, and Ep training ep. An early stopping strategy is adopted during training, with the training termination condition being that the loss function value on the validation set no longer decreases. Training is terminated if the validation set loss does not decrease in consecutive Pat rounds.

[0073] Step S230: Perform iterative updates using forward and backpropagation. In each training epoch, the input features Xu of the training samples are input into the defense resource scheduling decision network. After forward propagation through the threat level embedding layer, spatial coordinate encoding layer, feature cross-fusion layer, probe resource allocation branch, and interception resource allocation branch, the predicted output Ypred_u is obtained. The gradient of L with respect to Ypred_u is calculated, and the gradient is backpropagated layer by layer using the backpropagation algorithm to update the weight parameters and bias parameters of each layer in the network. The above process is repeated until the training termination condition is met, and the network weight parameters are saved to obtain the trained defense resource scheduling decision network.

[0074] Step S610: Obtain a set of real-time feedback data during the execution of the defense mission. The set of real-time feedback data includes flight trajectory offset data of the intrusion drone and signal reception power change data of the detection node.

[0075] During the execution of the defense mission, the actual flight trajectory of the intruding drone is continuously tracked through an airspace monitoring sensor network. The actual spatial coordinates of the intruding drone at each moment are recorded and compared with the preset flight path Path to obtain the flight trajectory offset data of the intruding drone ΔPath(t)={(Δx1, Δy1, Δz1, t1), (Δx2, Δy2, Δz2, t2), …}. At the same time, the actual signal received power time series Precv_i(t) of each detection node is obtained from the receiver of the detection node and compared with historical reference values ​​or theoretical calculation values ​​to obtain the signal received power change data of the detection node ΔPrecv(t)={ΔPrecv_1(t), ΔPrecv_2(t), …, ΔPrecv_Nd(t)}.

[0076] Step S620: Calculate the actual flight path deviation of the intrusion drone based on the flight trajectory offset data of the intrusion drone, input the actual flight path deviation into the multipath environmental masking effect distribution data for rematching, and obtain the updated multipath environmental masking effect distribution data corresponding to the actual flight path of the intrusion drone.

[0077] Calculate the actual flight path deviation of the intruding drone: For each time t, the spatial deviation distance d(t) between the actual flight path and the preset flight path Path is √(Δx(t)). 2 +Δy(t) 2 +Δz(t) 2 The original preset flight path Path is replaced by the sampling point sequence Path_actual(t) of the actual flight path. Following the methods in steps S135 to S137, the shading effect data along the actual flight path is re-extracted from the comprehensive distribution data Z of the environmental shading effect to obtain the updated multipath environmental shading effect distribution data MPupdated.

[0078] Step S630: Calculate the actual detection signal gain fluctuation of each detection node based on the signal received power change data of the detection node, compare the actual detection signal gain fluctuation with the updated multipath environmental shading effect distribution data, and identify the prediction deviation area of ​​the environmental shading effect distribution data.

[0079] Calculate the actual detection signal gain fluctuation of each detection node: For detection node i, define its signal gain fluctuation δgi = (max(Precv_i(t)) - min(Precv_i(t))) / Pref within the time period [t1, t2], where t∈[t1, t2] and Preref is the reference signal received power value. Compare the actual detection signal gain fluctuation δgi of each detection node with the predicted value of the occlusion effect intensity at the corresponding spatial location in the updated multipath environment occlusion effect distribution data MPupdated, and calculate the difference Δi = |δgi - m(xi, yi, zi)|, where m(xi, yi, zi) is the occlusion effect intensity value of MPupdated at the location of detection node i. Mark the spatial region where Δi exceeds the preset deviation threshold Tbias as the prediction deviation region Rbias.

[0080] Step S640: Based on the predicted deviation area, locally correct the multipath environmental masking effect distribution data to obtain corrected multipath environmental masking effect distribution data. Based on the corrected multipath environmental masking effect distribution data, re-evaluate the passive detection feasibility of the actual flight path of the intruding UAV, and generate a corrected threat detectability level identifier and corrected detection blind zone distribution boundary coordinates.

[0081] Within the prediction bias region Rbias, the masking effect intensity value in MPupdated is locally corrected. The correction method is: mcorrected(x, y, z) = m(x, y, z) + κ × Δ(x, y, z), where κ is the correction gain coefficient, and Δ(x, y, z) is the difference between the actual detection signal gain fluctuation and the predicted masking effect intensity value at that spatial location. The corrected multipath environment masking effect distribution data MPcorrected replaces the original m value with the mcorrected value within the Rbias region, while the original m value remains unchanged outside the Rbias region. Based on MPcorrected, the actual flight path of the intruding UAV is reassessed for passive detection feasibility according to step S140, generating the corrected threat detectability level identifier Lcorrected and the corrected detection blind zone distribution boundary coordinates Bcorrected.

[0082] Step S650: Input the corrected threat detectability level identifier and the corrected detection blind zone distribution boundary coordinates back into the defense resource scheduling decision network, and output the corrected multi-dimensional defense resource collaborative deployment instruction.

[0083] Using Lcorrected and Bcorrected as inputs, and following the method in step S150, the forward propagation calculation of the defense resource scheduling decision network is used to output the corrected probe node beamforming parameter group Pcorrected and the corrected interceptor unit maneuvering and resident waypoint sequence Tcorrected. Then, the timing synchronization alignment is performed according to the method in step S155 to generate the corrected multi-dimensional defense resource collaborative deployment instruction Cmdcorrected.

[0084] Step S660: Update the beamforming parameter group of the currently executing detection node and the maneuvering waypoint sequence of the interception unit with the modified multidimensional defense resource collaborative deployment instruction.

[0085] Replace the beamforming parameter set of the currently executing probe node with Pcorrected, and replace the maneuvering waypoint sequence of the currently executing interception unit with Tcorrected.

[0086] Step S710: Perform signal source separation processing on the electromagnetic spectrum signal distribution information to identify the electromagnetic radiation characteristic signals of the intruding drone and the environmental background electromagnetic radiation characteristic signals contained in the electromagnetic spectrum signal distribution information.

[0087] The electromagnetic spectrum signal distribution information E(x, y, z, f, t) is processed for signal source separation. Using independent component analysis (ICA) or a cluster-based blind source separation algorithm, E is decomposed into Ns independent signal source components. Among these Ns independent signal source components, their modulation characteristics and spectral morphology are matched against a pre-defined UAV electromagnetic fingerprint database. Signal components with a matching degree exceeding a pre-defined threshold are identified as the electromagnetic radiation characteristic signal Euav of the intruding UAV, while the remaining signal components are identified as the environmental background electromagnetic radiation characteristic signal Ebg. The mathematical expression for signal source separation is E(x, y, z, f, t) = ∑(s=1 to Ns)cs × As(f, t), where cs is the spatial mixing coefficient, and As(f, t) is the s-th signal source component.

[0088] Step S720: Extract the signal modulation mode feature and signal pulse repetition interval feature of the electromagnetic radiation characteristic signal of the intruding drone, compare the signal modulation mode feature and the signal pulse repetition interval feature with the preset drone electromagnetic fingerprint database, and determine the drone type identifier of the intruding drone.

[0089] Modulation pattern identification is performed on the electromagnetic radiation characteristic signal Euav of intrusive drones. An automatic modulation identification algorithm based on high-order cumulants and cyclostationary analysis is used to extract second-order, fourth-order, and sixth-order cumulant features of the signal. These cumulant features are input into a trained support vector machine classifier, which outputs the signal modulation pattern feature Mtype. Pulse repetition interval extraction is performed on Euav, and the autocorrelation function of the signal envelope is calculated. The time interval between the peaks of the autocorrelation function is extracted as the signal pulse repetition interval feature Tpri. The extracted Mtype and Tpri are compared with the electromagnetic feature templates of various types of drones in a pre-set drone electromagnetic fingerprint database. The comparison method is to calculate the Euclidean distance between the feature vectors, and the drone type identifier corresponding to the template with the smallest distance is taken as the drone type identifier IDuav of the intrusive drone.

[0090] Step S730: Based on the spatial propagation attenuation distribution of the electromagnetic radiation characteristic signal of the intruding drone in the electromagnetic propagation attenuation state field, calculate the signal coverage boundary of the electromagnetic radiation characteristic signal of the intruding drone within the three-dimensional airspace traversability constraint boundary, and generate the signal radiation coverage range of the intruding drone.

[0091] Taking the current spatial position (xuav, yuav, zuav) of the intruding drone as the signal radiation source point, the spatial propagation attenuation distribution centered on the source point is extracted from the electromagnetic propagation attenuation situation field H1(x, y, z, t). Using the ray tracing method, rays are emitted from the source point in all directions in space, and the electromagnetic propagation attenuation value is integrated along each ray path. When the integrated attenuation value reaches the preset signal-to-noise ratio detection threshold, the end point of the ray is the signal coverage boundary point. The passable airspace defined by the three-dimensional airspace trafficability constraint boundary B(x, y, z)=0 is taken as the effective spatial range for signal propagation. The set of coverage boundary points of rays in all directions within the effective spatial range constitutes the signal radiation coverage area Ruav of the intruding drone.

[0092] Step S740: Perform spatial correlation analysis on the signal radiation coverage area of the intruding drone and the multipath environment occlusion effect distribution data, and determine a first signal attenuation area affected by the environment occlusion effect and a second signal fidelity area not affected by the environment occlusion effect within the signal radiation coverage area of the intruding drone.

[0093] Perform spatial superposition analysis on the signal radiation coverage area Ruav of the intruding drone and the multipath environment occlusion effect distribution data MP. Within the spatial area covered by Ruav, the occlusion effect intensity value m at each position in MP is extracted. An occlusion effect intensity threshold Tm is set, the area where m≥Tm is marked as the first signal attenuation area Ratten, and the area where m<Tm is marked as the second signal fidelity area Rfidel.

[0094] Step S750: Adjust the pointing angle of the main lobe of the beam of each detection node in the detection node beamforming parameter group according to the spatial distribution of the first signal attenuation area and the second signal fidelity area, so that the coverage of the detection beam is biased toward the second signal fidelity area.

[0095] For each detection node, the initial pointing angle of its beam main lobe is (αi0, βi0). The spatial center position (μx_fidel, μy_fidel, μz_fidel) of the second signal fidelity area Rfidel is calculated. The pointing angle (αi, βi) of the beam main lobe of the detection node i is adjusted so that the direction of the beam main lobe points to the direction of (μx_fidel, μy_fidel, μz_fidel). The adjustment mode is: αi=arctan((μy_fidel-yi) / (μx_fidel-xi)), βi=arctan((μz_fidel-zi) / √((μx_fidel-xi) 2 +(μy_fidel-yi) 2 )), where (xi, yi, zi) is the spatial position coordinate of the detection node i.

[0096] Step S760: Update the multi-dimensional defense resource collaborative deployment command based on the adjusted probe node beamforming parameter group.

[0097] Replace the original Pbeam in the multidimensional defense resource collaborative deployment instruction with the adjusted probe node beamforming parameter group Pbeam_adjusted to generate the updated multidimensional defense resource collaborative deployment instruction.

[0098] Step S810: Perform time evolution trend analysis on the atmospheric turbulence disturbance potential field, extract the spatial structure change characteristics of the atmospheric turbulence disturbance potential field in continuous time segments, and generate the time evolution sequence of the atmospheric turbulence disturbance potential field.

[0099] The atmospheric turbulent perturbation potential field H2(x, y, z, t) is continuously sampled in the time dimension, with a sampling time interval of Δt2, for L consecutive time segments. The atmospheric turbulent perturbation potential field of the l-th time segment is H2(x, y, z, tl), where l = 1, 2, ..., L. The atmospheric turbulent perturbation potential field of adjacent time segments is differentiated to obtain the spatial structure change characteristic ΔH2(x, y, z, tl) = H2(x, y, z, tl+1) - H2(x, y, z, tl). The temporal evolution sequence of the atmospheric turbulent perturbation potential field is H2seq = {H2(x, y, z, t1), ΔH2(x, y, z, t1), H2(x, y, z, t2), ΔH2(x, y, z, t2), ..., H2(x, y, z, tL)}.

[0100] Step S820: Calculate the moving direction parameters and moving velocity parameters of the peak turbulence intensity region in the atmospheric turbulence disturbance potential field based on the time evolution sequence of the atmospheric turbulence disturbance potential field, and generate migration trend data of the peak turbulence region.

[0101] In each time segment tl of the temporal evolution sequence of the atmospheric turbulent perturbation potential field, the peak turbulence intensity region is extracted. This region is a connected region where the value of H2(x, y, z, tl) exceeds a preset peak threshold Tpeak. The geometric center coordinates (cx_peak(tl), cy_peak(tl), cz_peak(tl)) of this peak region are calculated. The direction parameter is determined by the displacement vector of the geometric center of the peak region of the adjacent time segment, and the direction vector v(tl) = [cx(tl+1)-cx(tl), cy(tl+1)-cy(tl), cz(tl+1)-cz(tl)]. The velocity parameter is the magnitude of the displacement vector divided by the time interval, s(tl) = ||v(tl)|| / Δt2. Generate turbulence peak region migration trend data Ttrend={(cx(tl), cy(tl), cz(tl), v(tl), s(tl))|l=1,2,…,L-1}.

[0102] Step S830: Perform spatiotemporal intersection analysis on the migration trend data of the turbulence peak region and the preset flight path of the intruding UAV to predict the time window and spatial location of the intruding UAV crossing the peak region of the atmospheric turbulence disturbance potential field in the subsequent period.

[0103] The spatiotemporal intersection analysis is performed between the geometric center (cx(tl), cy(tl), cz(tl)) and the movement direction vector v(tl) of the peak region in each time segment tl of the turbulent peak region migration trend data Ttrend and the preset flight path Path of the intruding UAV. Starting from the geometric center of the peak region, the direction of the movement direction vector v(tl) is used as the propulsion direction, and the movement speed s(tl) is used as the propulsion speed to construct the spatiotemporal trajectory of the peak region. The time ts at which the intruding UAV arrives at each sampling point (xs, ys, zs) on the Path is calculated by appending the preset flight speed of the intruding UAV to that sampling point. If there exists a time segment tl such that the spatiotemporal distance between the spatiotemporal trajectory of the peak region and the spatiotemporal coordinates of the intruding UAV arriving at a certain sampling point is less than the preset intersection threshold Dcross, then it is predicted that the intruding UAV will traverse the peak region of the atmospheric turbulent disturbance potential field at that sampling point location and time ts. The time window for traversal is [ts-Δtw, ts+Δtw], and the spatial location is (xs, ys, zs). Δtw is the preset half-width of the time window.

[0104] Step S840: Based on the time window and spatial location of the intruding UAV crossing the peak region of the atmospheric turbulence disturbance potential field, the multipath environmental shading effect distribution data is marked with timeliness, and the part that will be affected by the peak region of the atmospheric turbulence disturbance potential field is marked as the shading effect decay region.

[0105] In the multipath environment shading effect distribution data MP, all sampling points within the neighborhood of the peak region crossing point (xs, ys, zs) are extracted. The spatial radius Rneighbor of this neighborhood is set according to the characteristic scale of the turbulence peak region. The sampling points within the above neighborhood are marked as the shading effect decay region Rdecay, and a predicted crossing time window [ts-Δtw, ts+Δtw] is attached as a timeliness marker.

[0106] Step S850: Based on the labeling results of the shading effect decay region, perform attenuation prediction correction on the shading effect intensity sequence along the preset flight path in the multipath environmental shading effect distribution data to obtain the multipath environmental shading effect distribution data after the prediction period.

[0107] Within the shading effect decay region Rdecay, the shading effect intensity value m in the multipath environmental shading effect distribution data MP is subject to attenuation prediction correction. The attenuation prediction correction uses an exponential decay model: mdecay(t) = m0 × e^(-λ × (t-ts)), where m0 is the shading effect intensity value at the current time, λ is the attenuation coefficient, t is the prediction time, and ts is the start time of the crossing time window. For each sampling point within Rdecay, its shading effect intensity value m at the current time is replaced with the attenuation correction value mdecay(tpred) at the prediction time tpred. After the prediction period, the multipath environmental shading effect distribution data MPpred uses the corrected mdecay value within the Rdecay region, while maintaining the original m value in other regions.

[0108] Step S860: Based on the multipath environmental shading effect distribution data after the prediction period, adjust the boundary coordinates of the detection blind zone distribution in advance to generate the adjusted boundary coordinates of the detection blind zone distribution.

[0109] Based on MPpred's preset flight path for the intruding drone, a new passive detection feasibility assessment is performed according to step S140, generating the detection blind zone distribution boundary coordinates Bpred after the predicted time period. The advanced adjusted detection blind zone distribution boundary coordinates Bpred_adjusted=Bpred.

[0110] Step S870: Input the adjusted boundary coordinates of the detection blind zone distribution into the defense resource scheduling decision network, and output the adjusted multi-dimensional defense resource collaborative deployment command.

[0111] Input Bpred_adjusted into the defense resource scheduling decision network, and according to the method in step S150, output the advanced-adjusted detection node beamforming parameter group Ppred and the interception unit maneuvering and stationing waypoint sequence Tpred, and generate the advanced-adjusted multi-dimensional defense resource collaborative deployment instruction Cmdpred.

[0112] Step S910: Obtain the characteristic data of the interference signal actively emitted by the intrusion drone during flight. The interference signal characteristic data includes the center frequency parameter, bandwidth parameter, and power spectral density distribution of the interference signal.

[0113] Interference signals actively emitted by intruding drones are detected using a spectrum monitoring receiver array in an airspace monitoring sensor network. Time-frequency analysis is performed on the received signals to identify broadband high-power signal components that do not match the characteristics of known communication and radar signals, marking them as interference signals. The center frequency parameter fjam, the bandwidth parameter Bjam, and the power spectral density distribution Pjam(f) of the interference signal are extracted, where f∈[fjam-Bjam / 2, fjam+Bjam / 2].

[0114] Step S920: Based on the interference signal center frequency parameter and interference signal bandwidth parameter in the interference signal feature data, perform interference suppression processing on the corresponding frequency band in the electromagnetic spectrum signal distribution information to obtain the interference-free electromagnetic spectrum signal distribution information.

[0115] In the electromagnetic spectrum signal distribution information E(x, y, z, f, t), the portion of frequency f within the interference signal frequency band [fjam-Bjam / 2, fjam+Bjam / 2] is contaminated by the interference signal. Interference suppression processing is applied to the E value within this frequency band using a frequency domain notch filter method. Within the interference frequency band, E(x, y, z, f, t) is multiplied by a notch mask function Mnotch(f). Mnotch(f) takes the value εnotch within the interference frequency band and 1 outside the interference frequency band, where εnotch is a positive number much less than 1. The de-interference electromagnetic spectrum signal distribution information Eclean(x, y, z, f, t) = E(x, y, z, f, t) × Mnotch(f).

[0116] Step S930: The electromagnetic spectrum signal distribution information after interference removal is re-analyzed in the spectrum domain, the spatial gradient distribution of signal intensity in each spectrum segment after interference removal is extracted, the spatial gradient distribution of signal intensity in the electromagnetic spectrum signal distribution information is replaced by the spatial gradient distribution of signal intensity after interference removal, and the electromagnetic propagation attenuation state field after interference removal is regenerated.

[0117] Replace the original E(x,y,z,f,t) with Eclean(x,y,z,f,t), and re-perform the spectrum domain segmentation analysis, signal strength spatial gradient distribution extraction and time-frequency co-mapping according to the method of step S121 to generate the electromagnetic propagation attenuation state field H1clean(x,y,z,t) after interference removal.

[0118] Step S940: Using the three-dimensional spatial mobility constraint boundary as the spatial constraint condition, the electromagnetic propagation attenuation potential field after interference removal and the atmospheric turbulence disturbance potential field are superimposed again using multi-physics fields to generate multi-path environmental shielding effect distribution data after interference removal.

[0119] H1clean(x, y, z, t) is used as the first physical field effective data, and the atmospheric turbulence disturbance potential field H2(x, y, z, t) is used as the second physical field effective data. The multi-physical fields are superimposed again according to the method in step S130 to generate the multipath environmental shading effect distribution data MPclean after interference removal.

[0120] Step S950: Based on the multipath environmental masking effect distribution data after interference removal, re-evaluate the passive detection feasibility of the preset flight path of the intruding UAV, and generate the threat detectability level identifier and the boundary coordinates of the detection blind zone distribution after interference removal.

[0121] Based on the preset flight path of the intruding drone using MPclean, the feasibility of passive detection is reassessed according to the method in step S140, generating the threat detectability level identifier Lclean after interference removal and the boundary coordinates of the detection blind zone distribution Bclean after interference removal.

[0122] Step S960: Input the threat detectability level identifier after interference removal and the boundary coordinates of the detection blind zone distribution after interference removal into the defense resource scheduling decision network, and output the multi-dimensional defense resource collaborative deployment instruction after interference removal.

[0123] Input Lclean and Bclean into the defense resource scheduling decision network, and according to the method in step S150, output the de-interference detection node beamforming parameter group Pclean and the interception unit maneuvering and stationing waypoint sequence Tclean, and generate the de-interference multi-dimensional defense resource collaborative deployment instruction Cmdclean.

[0124] Step S970: Replace the currently executing multi-dimensional defense resource collaborative deployment instruction with the interference-free multi-dimensional defense resource collaborative deployment instruction.

[0125] Replace the currently executing multi-dimensional defense resource collaborative deployment command with Cmdclean.

[0126] For example, step S1010: perform path segmentation processing on the preset flight path of the intrusion drone, and divide the preset flight path into multiple flight path segments according to spatial turning points.

[0127] The preset flight path of the intruding drone, Path={(xp1, yp1, zp1), (xp2, yp2, zp2), ..., (xpR, ypR, zpR)}, is segmented. The direction vector between any two adjacent waypoints on the Path is calculated. If the angle between adjacent direction vectors exceeds a preset turning angle threshold θturn, the waypoint at the angle is marked as a spatial turning point. Using the spatial turning point as the segmentation point, the Path is divided into multiple flight path segments, Segment={seg1, seg2, ..., segNseg}, where each segment is a continuous subsequence of the Path.

[0128] Step S1020: Extract the corresponding multipath environment shading effect distribution data segment for each flight path segment to obtain the segment shading effect intensity sequence and segment shading effect direction angle sequence for each flight path segment.

[0129] For each flight path segment segr, r=1,2,…,Nseg, extract the occlusion effect intensity and orientation angle values ​​corresponding to the sampling points covered by the segment from the multipath environment occlusion effect distribution data MP. The segment occlusion effect intensity sequence is Lseg_r={m_r1, m_r2,…, m_rNr}, and the segment occlusion effect orientation angle sequence is Aseg_r={φ_r1, φ_r2,…, φ_rNr}, where Nr is the number of sampling points within the segment segr.

[0130] Step S1030: Calculate the average occlusion effect intensity of each flight path segment based on the segment occlusion effect intensity sequence corresponding to each flight path segment, and calculate the occlusion effect direction consistency of each flight path segment based on the segment occlusion effect direction angle sequence corresponding to each flight path segment.

[0131] Calculate the average occlusion effect intensity of segment segr: μr = (1 / Nr) × ∑(i=1 to Nr)m_ri. Calculate the occlusion effect direction consistency σr of segment segr, where σr is the circular square of all directional angles in the occlusion effect direction angle sequence Aseg_r, σr = 1 - ||(1 / Nr) × ∑(i=1 to Nr)[cosφ_ri, sinφ_ri] ||.

[0132] Step S1040: Compare the average shielding effect intensity and shielding effect direction consistency of each flight path segment with the preset segment threat determination intensity threshold and the preset segment threat determination direction threshold, respectively, and mark the flight path segments that exceed the preset segment threat determination intensity threshold and are lower than the preset segment threat determination direction threshold as high shielding risk segments.

[0133] Set the preset segment threat determination intensity threshold as Tμ, and the preset segment threat determination direction threshold as Tσ. For each segment segr, if μr>Tμ and σr<Tσ, this segment is marked as a high occlusion risk segment. The set of all high occlusion risk segments is denoted as RiskSeg.

[0134] Step S1050: Extract the start point coordinates and end point coordinates of the high occlusion risk segments according to the spatial position distribution of the high occlusion risk segments on the preset flight path, and generate the spatial range of the high occlusion risk area.

[0135] For each high occlusion risk segment segr in RiskSeg, extract the spatial coordinates of its start sampling point (xr_start, yr_start, zr_start) and the spatial coordinates of its end sampling point (xr_end, yr_end, zr_end). The start-end coordinate pairs of all high occlusion risk segments constitute the spatial range of the high occlusion risk area RiskArea={(xr_start, yr_start, zr_start, xr_end, yr_end, zr_end)|segr∈RiskSeg}.

[0136] Step S1060: Based on the spatial range of the high occlusion risk area, add a defense resource deployment density weight for the spatial range of the high occlusion risk area in the defense resource scheduling decision network, the defense resource scheduling decision network outputs a detection node beamforming parameter set containing enhanced deployment parameters, and the enhanced deployment parameters include a beam main lobe pointing angle adjustment amount and a transmission power gain increase amount.

[0137] A spatial attention weight mechanism is introduced in the detection resource allocation branch of the defense resource scheduling decision network. For the output of the beam main lobe pointing angle in the detection node beamforming parameter set, a spatial attention weight matrix Wattn is multiplied before the output layer. The weight value corresponding to the spatial range of the high occlusion risk area RiskArea in Wattn is set as an enhancement coefficient γattn higher than 1, and the weight values of other areas remain 1. The enhanced beam main lobe pointing angle αi_boost=αi+Δαi, where Δαi is the adjustment amount pointing to the high occlusion risk area. The transmission power gain coefficient gi_boost=gi+Δgi, and Δgi is the gain increase amount. The defense resource scheduling decision network outputs a detection node beamforming parameter set Pboost containing the enhanced deployment parameters.

[0138] Step S1110: Spatial layering of the electromagnetic propagation attenuation state field is performed. The electromagnetic propagation attenuation state field is divided into multiple horizontal propagation attenuation slices along the vertical height direction. The contour distribution of the electromagnetic propagation attenuation state field values ​​in each horizontal propagation attenuation slice is extracted to generate a family of electromagnetic propagation attenuation contours.

[0139] The electromagnetic propagation attenuation state field H1(x, y, z, t) is spatially segmented along the z-axis. The z-axis coordinate range is divided into Znum slices with a layer height of Δzslice. The height range of the k-th horizontal propagation attenuation slice is [zk, zk+Δzslice], and the electromagnetic propagation attenuation state field slice within this layer is H1slice_k(x, y) = H1(x, y, z, t), where z ∈ [zk, zk+Δzslice]. Contour extraction is performed on each slice H1slice_k(x, y), with the contour level set to Cnum and the contour interval set to ΔC. The c-th contour line is the trajectory line of H1slice_k(x, y) = c × ΔC. All contour lines of all slices constitute the electromagnetic propagation attenuation contour family Contour1.

[0140] Step S1120: Spatially segment the atmospheric turbulence disturbance potential field, divide the atmospheric turbulence disturbance potential field into multiple horizontal turbulence disturbance slices along the vertical height direction, extract the contour distribution of the atmospheric turbulence disturbance potential field values ​​in each horizontal turbulence disturbance slice, and generate a family of atmospheric turbulence disturbance contour lines.

[0141] Using the same spatial layering scheme as in step S1110, the atmospheric turbulence disturbance potential field H2(x, y, z, t) is sliced ​​to obtain the atmospheric turbulence disturbance contour family Contour2.

[0142] Step S1130: Perform spatial overlay analysis on the electromagnetic propagation attenuation contour families and atmospheric turbulence disturbance contour families corresponding to the same altitude layer, identify the spatial proximity relationship and intersection region of the two types of contours, extract the coordinates of the intersection region with significant electromagnetic attenuation and strong turbulence characteristics, and generate a set of intersection regions with physical field characteristics.

[0143] Spatial overlay analysis is performed between the electromagnetic propagation attenuation contour families Contour1_k and atmospheric turbulence disturbance contour families Contour2_k corresponding to altitude level zk. For each contour line in Contour1_k and each contour line in Contour2_k, the shortest spatial distance between them is calculated. If the shortest spatial distance between two contour lines is less than a preset proximity threshold Dprox, the intersection region of these two contour lines is marked as a physical field feature intersection region. The coordinates of the intersection region are taken as the geometric center coordinates of the nearest points of the two contour lines. All physical field feature intersection regions at all altitude levels constitute the physical field feature intersection region set Intersect.

[0144] Step S1140: Quantify the local environmental shielding effect of each intersection region in the set of physical field feature intersection regions. Based on the preset coupling evaluation rules, comprehensively analyze the synergistic effect strength of electromagnetic propagation attenuation characteristics and atmospheric turbulence disturbance characteristics at the intersection region, and generate a local environmental interference evaluation value for each intersection region.

[0145] For each intersection region in the set of physical field feature intersection regions Intersect, the electromagnetic propagation attenuation value h1 and the atmospheric turbulence disturbance value h2 are extracted from its corresponding spatial location. Based on the preset coupling evaluation rules, the local environmental disturbance evaluation value Jlocal of the intersection region is calculated, Jlocal=ω1×h1+ω2×h2+ω3×h1×h2, where ω1, ω2, and ω3 are coupling evaluation weight coefficients.

[0146] Step S1150: Spatial interpolation diffusion of the local environmental interference assessment values ​​of all intersection regions in the set of physical field feature intersection regions within the three-dimensional spatial mobility constraint boundary to generate environmental shading effect coupled field distribution data.

[0147] Using the local environmental disturbance assessment value Jlocal of each intersection region in the set of physical field feature intersection regions Intersect as known sample points, Jlocal is spatially diffused using the Kriging space interpolation method within the traversable space Ω defined by the three-dimensional spatial traversability constraint boundary B(x, y, z) = 0. The interpolation result is the environmental shading effect coupled field distribution data Jcouple(x, y, z), which has a coupled disturbance assessment value at each voxel location within Ω.

[0148] Step S1160: Extract the coordinates of local maxima points in the environmental shading effect coupled field distribution data based on the environmental shading effect coupled field distribution data, and mark the coordinates of the local maxima points as sensitive nodes of the environmental shading effect.

[0149] Local maxima detection is performed on the coupled field distribution data Jcouple(x, y, z) of the environmental shading effect. At each voxel location within Ω, its Jcouple value is compared with the Jcouple values ​​of all neighboring voxels in its 3D neighborhood. If the Jcouple value of the voxel is greater than the Jcouple values ​​of all neighboring voxels, then the voxel location is marked as a local maximum. The coordinates of all local maxima points are then labeled as NodeSens, sensitive nodes for the environmental shading effect.

[0150] Step S1170: Perform spatial connectivity clustering on the environmental shading effect sensitive nodes, merge environmental shading effect sensitive nodes with a spatial distance less than a preset connectivity threshold into environmental shading effect sensitive regions, and extract the region boundary contour and region centroid coordinates of the environmental shading effect sensitive regions.

[0151] For all nodes in the set NodeSens, which is sensitive to environmental shading effects, calculate the spatial Euclidean distance between nodes. Connect nodes whose spatial distance is less than a preset connectivity threshold Dconn to form the same connected component. Use a disjoint-set data structure algorithm to perform connectivity clustering, and each connected component constitutes an environmental shading effect sensitive region SensArea. For each SensArea, extract its region boundary contour, which is the triangular mesh of the convex hull surface of the nodes in that region; calculate the region's centroid coordinates as the arithmetic mean of the coordinates of all nodes.

[0152] Step S1180: Based on the environmental shielding effect sensitive area, perform path risk assessment on the preset flight path of the intruding drone, calculate the proportion of the path length that the preset flight path passes through or is close to the environmental shielding effect sensitive area, and generate a path environmental sensitivity assessment value.

[0153] Determine whether each sampling point on the preset flight path Path of the intruding drone falls within the spatial range of any SensArea sensitive area due to environmental masking effects, or whether the distance to the boundary of the SensArea is less than the preset proximity distance Dclose. Count the number of sampling points Nsens that meet the above conditions, and calculate the path length ratio Rsens = Nsens / S, where S is the total number of sampling points. The path environmental sensitivity assessment value Esens = Rsens.

[0154] Step S1190: Correct the initial assessment result of the threat detectability level identifier based on the path environment sensitivity assessment value, and superimpose the path environment sensitivity assessment value into the generation process of the threat detectability level identifier to obtain the corrected threat detectability level identifier.

[0155] In the original threat detectability level identifier generation process, the path environment sensitivity assessment value Esens is used as an adjustment factor and fused with the detection blind zone coverage rate Rb to obtain the corrected comprehensive risk index Rbmod = Rb + ηsens × Esens, where ηsens is the sensitivity adjustment weight coefficient. Based on the preset coverage range in which Rbmod falls, the corrected threat detectability level identifier Lmod is redefined.

[0156] Step S1210: Perform electromagnetic propagation direction inversion on the electromagnetic propagation attenuation state field, extract the main direction vector of electromagnetic energy propagation at each spatial location in the electromagnetic propagation attenuation state field, and generate an electromagnetic energy propagation direction vector field.

[0157] Electromagnetic propagation direction inversion is performed on the electromagnetic propagation attenuation state field H1(x, y, z, t). At each spatial location (x, y, z), the electromagnetic energy propagation direction is determined by the time-averaged direction of the Poynting vector. In this embodiment, the Poynting vector direction is approximately estimated by the negative gradient direction of H1, i.e., the principal direction vector of electromagnetic energy propagation p(x, y, z) = -∇H1 = -[∂H1 / ∂x, ∂H1 / ∂y, ∂H1 / ∂z]. Normalizing p(x, y, z) yields the unit principal direction vector pu(x, y, z) = p / ‖p‖. pu at all spatial locations constitutes the electromagnetic energy propagation direction vector field Pfield.

[0158] Step S1220: Extract the atmospheric refractive gradient direction from the atmospheric turbulence disturbance potential field, extract the atmospheric refractive index gradient direction vector at each spatial location in the atmospheric turbulence disturbance potential field, and generate an atmospheric refractive gradient direction vector field.

[0159] For the atmospheric turbulent perturbation potential field H2(x, y, z, t), calculate its spatial gradient ∇H2=[∂H2 / ∂x, ∂H2 / ∂y, ∂H2 / ∂z]. The atmospheric refractive index gradient direction vector q(x, y, z)=∇H2 / ‖∇H2‖. q at all spatial locations constitutes the atmospheric refractive index gradient direction vector field Qfield.

[0160] Step S1230: Perform spatial vector relationship analysis on the electromagnetic energy propagation direction vector field and the atmospheric refractive index gradient direction vector field. Analyze the spatial angle and relative orientation between the electromagnetic energy propagation main direction vector and the atmospheric refractive index gradient direction vector at each spatial location to generate a directional coupling situation field that reflects the degree of directional coupling between the two.

[0161] At each spatial location (x, y, z), calculate the spatial angle δ(x, y, z) = arccos(pu·q) between the principal direction vector of electromagnetic energy propagation pu and the direction vector of atmospheric refractive index gradient q. The directional coupled situational field K(x, y, z) = δ(x, y, z) is in radians.

[0162] Step S1240: Perform spatial distribution analysis of coupling intensity on the directional coupling situation field, and extract the region where the directional coupling characteristics significantly exceed the preset threshold as the electromagnetic-atmospheric coupling anomaly region.

[0163] Set a directional coupling threshold δth. In the directional coupling situational field K(x,y,z), the spatial region where δ(x,y,z)>δth is marked as the electromagnetic-atmospheric coupling anomaly region Kabnormal.

[0164] Step S1250: Perform internal structure analysis on the electromagnetic-atmospheric coupling anomaly region, compare the angle between the curl vector direction distribution in the electromagnetic-atmospheric coupling anomaly region and the preset flight path direction of the intruding UAV, and extract sub-regions where the angle between the curl vector direction and the preset flight path direction is less than a preset parallel threshold as path co-directional anomaly regions.

[0165] Within the abnormal electromagnetic-atmospheric coupling region Kabnormal, calculate the curl vector Rcurl = ∇ × Pfield of the electromagnetic energy propagation direction vector field Pfield. Extract the tangent direction vector Tpath of the intruding UAV's preset flight path Path at each spatial location. Calculate the angle ζ between Rcurl and Tpath = arccos((Rcurl·Tpath) / (‖Rcurl‖×‖Tpath‖)). Set a preset parallelism threshold ζth. Mark the region within Kabnormal that satisfies ζ < ζth as the path co-directional abnormal region Rparallel.

[0166] Step S1260: Match the spatial coordinates of the path co-directional anomaly region with the multipath environmental shading effect distribution data point by point, and extract the shading effect intensity distribution data located in the path co-directional anomaly region from the multipath environmental shading effect distribution data.

[0167] Spatial intersection is performed between the spatial coordinate range of the path co-directional anomaly region Rparallel and the multipath environment shading effect distribution data MP. The shading effect intensity value m corresponding to each location within Rparallel is extracted to obtain the shading effect intensity distribution data Mparallel within the path co-directional anomaly region.

[0168] Step S1270: Identify abnormal fluctuation features of the shading effect intensity distribution data in the path co-directional anomaly region, extract the spatial variation gradient sequence of the shading effect intensity along the preset flight path direction in the path co-directional anomaly region, and mark the position where the gradient sign in the spatial variation gradient sequence is flipped as the shading effect singularity position.

[0169] Within the path-co-directional anomaly region Rparallel, the directional derivative of the shading effect intensity value m, dm / ds = ∇m·Tpath, is calculated along the tangent direction Tpath of the preset flight path Path. The directional derivatives at each position within Rparallel are arranged in path order to obtain the spatially varying gradient sequence Gseq. The position where the gradient sign flips in Gseq, i.e., the position where the signs of Gseq(k) and Gseq(k+1) are opposite, is detected and marked as the shading effect singularity position, Posingularity.

[0170] Step S1280: Based on the singularity position of the occlusion effect, perform singularity correction on the boundary coordinates of the detection blind zone distribution, incorporate the singularity position of the occlusion effect into the generation basis of the boundary coordinates of the detection blind zone distribution, redefine the boundary line direction of the boundary coordinates of the detection blind zone distribution, and obtain the boundary coordinates of the detection blind zone distribution after singularity correction.

[0171] During the generation of the original detection blind zone distribution boundary coordinates, when the boundary line of the connected domain of the detection blind zone passes near the singularity position of the occlusion effect (Posingularity), the boundary line is forcibly adjusted to pass through Posingularity. Specifically, within the preset neighborhood radius (rsingularity) of Posingularity, the original boundary line is replaced with a polyline segment connecting the neighboring boundary point to Posingularity and then to the neighboring boundary point on the other side. The redefined boundary coordinates (Bsingularity) of the detection blind zone distribution are the result after singularity correction.

[0172] Step S1310: Analyze the electromagnetic reflection characteristics of the obstacle surface based on the boundary contour information of the geographic obstacle. Determine the electromagnetic reflection coefficient distribution of the obstacle surface based on the surface material properties of each obstacle in the boundary contour information of the geographic obstacle, and generate electromagnetic reflection coefficient distribution data of the obstacle surface.

[0173] The geographic obstacle boundary contour information O(x, y, z) is associated with obstacle surface material attribute data. A preset material-electromagnetic reflection coefficient mapping table stores the electromagnetic reflection coefficient values ​​of common obstacle materials in specific frequency bands. Based on the surface material attribute of each obstacle, the electromagnetic reflection coefficient Γ of the obstacle surface is obtained by looking up the table. Electromagnetic reflection coefficient distribution data Γ(x, y, z) of the obstacle surface is generated, where a reflection coefficient value is defined at each point on the obstacle surface, and the value is 0 at non-obstacle surface locations.

[0174] Step S1320: Analyze the surface roughness characteristics of the geographical obstacle boundary contour information, determine the root mean square parameter of the surface undulation of the obstacle surface based on the surface morphology characteristics of each obstacle in the geographical obstacle boundary contour information, and generate obstacle surface roughness distribution data.

[0175] For each obstacle's surface triangular mesh data, calculate the root mean square parameter (hrms) of the surface undulation for each triangular facet. hrms is calculated by extracting the perpendicular distance from each vertex on the triangular facet and within its neighboring triangular facets to the fitted plane, and then calculating the root mean square value of this distance. Generate obstacle surface roughness distribution data hrms(x, y, z), which defines a roughness height value at each point on the obstacle surface, and sets the value to 0 at locations on non-obstacle surfaces.

[0176] Step S1330: Input the electromagnetic reflection coefficient distribution data and the surface roughness distribution data of the obstacle into the obstacle electromagnetic scattering directionality modeling step to generate electromagnetic scattering pattern distribution data for each obstacle surface. The electromagnetic scattering pattern distribution data includes the scattering main lobe direction and scattering side lobe distribution range of the scattered beam on each obstacle surface.

[0177] The Kirchhoff approximation method is used to model the electromagnetic scattering directionality of obstacles. For each triangular facet on the obstacle surface, the scattering pattern of the incident electromagnetic wave on that facet is determined by the electromagnetic reflection coefficient Γ of the facet, the surface roughness hrms, the incident wavelength λ, and the incident angle. The main lobe direction of scattering is the specular reflection direction, and the specular reflection direction vector rspecular = rincident - 2 × (rincident·n) × n, where rincident is the unit direction vector of the incident wave, and n is the normal vector of the facet. The distribution range of the scattering sidelobes is determined by the surface roughness, and the half-power beamwidth θHPBW of the sidelobes is related to the ratio of hrms to λ, with θHPBW increasing as hrms / λ increases. For all triangular facets on each obstacle surface, scattering pattern distribution data is generated, including the main lobe direction rspecular and the distribution range of the scattering sidelobes θHPBW.

[0178] Step S1340: Based on the electromagnetic energy propagation main direction vector at each spatial location in the electromagnetic propagation attenuation state field, determine the direct propagation path of the detection beam emitted from the detection node location to each sampling point on the preset flight path of the intrusion UAV.

[0179] The detection node is located at (xnode, ynode, znode). For each sampling point (xs, ys, zs) on the preset flight path Path of the intrusion drone, the direct propagation path is a straight line segment Ls from (xnode, ynode, znode) to (xs, ys, zs), with direction vector vs = (xs-xnode, ys-ynode, zs-znode) and length ds = ‖vs‖.

[0180] Step S1350: Spatial intersection of the direct propagation path and the boundary contour information of the geographical obstacle is performed, and the coordinates of the reflection point where the direct propagation path intersects with the surface of the obstacle are extracted. The coordinates of the reflection point are used as the multipath reflection source point.

[0181] For each direct propagation path Ls, spatial intersection detection is performed with the geographic obstacle boundary contour information O(x, y, z). If Ls intersects with any obstacle triangular facet, occlusion exists, and the coordinates of the intersection point (xrefl, yrefl, zrefl) are calculated. This intersection point is the multipath reflection source point. The set of all multipath reflection source points is denoted as ReflPoints.

[0182] Step S1360: At each multipath reflection source point, determine the main lobe direction and side lobe distribution range of the scattering on the surface of the obstacle corresponding to the multipath reflection source point based on the electromagnetic scattering pattern distribution data. Perform directional alignment analysis between the main lobe direction and side lobe distribution range and the preset flight path of the intruding UAV, and mark the multipath reflection source points whose scattering energy points to the preset flight path as effective multipath reflection source points.

[0183] For each multipath reflection source, the main lobe direction (rspecular) and side lobe distribution range (θHPBW) of the triangular facet containing the reflection source are extracted from the electromagnetic scattering pattern distribution data of the corresponding obstacle surface. A direction vector is drawn from the reflection source to each sampling point on the preset flight path. If the angle between this direction vector and rspecular is less than θHPBW, the scattered energy of the reflection source is determined to point towards the preset flight path, and the reflection source is marked as a valid multipath reflection source (ReflPointsEff).

[0184] Step S1370: Input the spatial coordinates of each effective multipath reflection source point and the corresponding scattering pattern distribution data into the multipath effect analysis step, evaluate the contribution of the scattered signal to each sampling point on the preset flight path after multipath reflection propagation and the characteristics of the direct signal of the direct propagation path, calculate the composite field strength distribution under the influence of multipath effect based on the electromagnetic wave propagation model, and generate the equivalent electromagnetic propagation attenuation state field after considering the multipath reflection effect.

[0185] For each effective multipath reflection source, the multipath effect is calculated using an electromagnetic wave propagation model. The direct signal field strength Edirect(xs, ys, zs) of the direct propagation path is calculated using the free-space propagation loss formula. The scattered signal field strength Escatter(xs, ys, zs) of the multipath reflection path is calculated by superimposing the obstacle surface reflection loss and the scattering pattern gain on top of the direct propagation loss. The reflection loss is determined by the reflection coefficient Γ, and the scattering pattern gain is determined by the corresponding gain value in the scattering pattern of the deviation angle between the scattering main lobe direction and the direction from the reflection source to the sampling point. The composite field strength Etotal(xs, ys, zs) is the complex vector superposition of the direct signal field strength and the field strengths of each multipath scattered signal, Etotal = Edirect + ∑Escatter. The equivalent electromagnetic propagation attenuation Lmulti(xs, ys, zs) after considering the multipath reflection effect is calculated by taking the logarithm of the ratio of the composite field strength to the emitted field strength. Replace the electromagnetic propagation attenuation value at the corresponding sampling point in the original H1 with Lmulti to generate the equivalent electromagnetic propagation attenuation state field H1multi.

[0186] Step S1380: Using the three-dimensional spatial mobility constraint boundary as the spatial constraint condition, the equivalent electromagnetic propagation attenuation potential field and the atmospheric turbulence disturbance potential field are superimposed again using a multi-physics field to generate equivalent multipath environmental shielding effect distribution data after considering the multipath reflection effect.

[0187] Using H1multi as the first physical field effective data and the atmospheric turbulence disturbance potential field H2 as the second physical field effective data, the multi-physical fields are superimposed again according to the method in step S130 to generate the equivalent multipath environmental shading effect distribution data MPmulti after considering the multipath reflection effect.

[0188] Step S1390: Based on the equivalent multipath environmental masking effect distribution data, re-evaluate the passive detection feasibility of the preset flight path of the intruding UAV, and generate an updated threat detectability level identifier and updated detection blind zone distribution boundary coordinates.

[0189] Based on the preset flight path of the intruding drone using MPmulti, the passive detection feasibility assessment is re-performed according to the method in step S140, generating an updated threat detectability level identifier Lmulti and updated detection blind zone distribution boundary coordinates Bmulti.

[0190] Step S13100: Based on the updated threat detectability level identifier and the updated detection blind zone distribution boundary coordinates, call the defense resource scheduling decision network to output the updated multi-dimensional defense resource collaborative deployment instructions.

[0191] Input Lmulti and Bmulti into the defense resource scheduling decision network, and according to the method in step S150, output the updated probe node beamforming parameter group Pmulti and the interception unit maneuvering and stationing waypoint sequence Tmulti, and generate the updated multi-dimensional defense resource collaborative deployment instruction Cmdmulti.

[0192] In some embodiments, the drone threat assessment and defense decision-making system incorporating environmental awareness for performing the above-described methods can be any electronic device with data computing, processing, and storage capabilities. This drone threat assessment and defense decision-making system incorporating environmental awareness can be used to implement the drone threat assessment and defense decision-making method incorporating environmental awareness provided in the above embodiments.

[0193] Typically, a drone threat assessment and defense decision-making system incorporating environmental perception includes a processor and memory. The processor may include one or more processing cores, such as a quad-core processor or an octa-core processor. The processor can be implemented using at least one hardware form of DSP (Digital Signal Processing), FPGA (Field Programmable Gate Array), or PLA (Programmable Logic Array). The processor may also include a main processor and coprocessors. The main processor, also known as a CPU (Central Processing Unit), processes data in the wake-up state; the coprocessor is a low-power processor that processes data in the standby state. In some embodiments, the processor may integrate a GPU (Graphics Processing Unit), which is responsible for rendering and drawing the content to be displayed on the screen. In some embodiments, the processor may also include an AI (Artificial Intelligence) processor, which handles computational operations related to machine learning.

[0194] The memory may include one or more computer-readable storage media, which may be non-transitory. The memory may also include high-speed random access memory and non-volatile memory, such as one or more disk storage devices or flash memory devices. In some embodiments, the non-transitory computer-readable storage media in the memory is used to store a computer program configured to be executed by one or more processors to implement the aforementioned method for unmanned aerial vehicle threat assessment and defense decision-making combined with environmental awareness.

[0195] In an illustrative embodiment, a computer-readable storage medium is also provided, wherein a computer program is stored in the storage medium. When executed by a processor of a computer device, the computer program implements the aforementioned method for drone threat assessment and defense decision-making combined with environmental perception. Optionally, the aforementioned computer-readable storage medium may be ROM (Read-Only Memory), RAM (Random Access Memory), CD-ROM (CompactDisc Read-Only Memory), magnetic tape, floppy disk, and optical data storage device, etc.

[0196] This application provides a computer program product, which includes computer-executable instructions or a computer program. When the computer-executable instructions or the computer program are executed by a processor, the processor will execute the UAV threat assessment and defense decision-making method combined with environmental awareness provided in this application.

[0197] This application provides a computer-readable storage medium storing computer-executable instructions or computer programs. When the computer-executable instructions or computer programs are executed by a processor, the processor will execute the UAV threat assessment and defense decision-making method combined with environmental awareness provided in this application.

[0198] In some embodiments, the computer-readable storage medium may be a read-only memory (ROM), random access memory (RAM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), flash memory, magnetic surface memory, optical disk, or CD-ROM, etc.; or it may be a device that includes one or any combination of the above-mentioned memories.

[0199] In some embodiments, computer-executable instructions may take the form of programs, software, software modules, scripts, or code, written in any form of programming language (including compiled or interpreted languages, or declarative or procedural languages), and may be deployed in any form, including as stand-alone programs or as modules, components, subroutines, or other units suitable for use in a computing environment.

[0200] As an example, computer-executable instructions may, but do not necessarily, correspond to files in a file system. They may be stored as part of a file that holds other programs or data, for example, in one or more scripts in a HyperText Markup Language (HTML) document, in a single file dedicated to the program in question, or in multiple co-located files (e.g., files that store one or more modules, subroutines, or code sections).

[0201] As an example, computer-executable instructions can be deployed to execute on a single electronic device, or on multiple electronic devices located at one location, or on multiple electronic devices distributed across multiple locations and interconnected via a communication network.

[0202] Finally, it should be noted that the above-disclosed embodiments are merely preferred embodiments of the present invention and are only used to illustrate the technical solutions of the present invention, not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for drone threat assessment and defense decision-making that combines environmental perception, characterized in that, The method includes: Acquire a set of real-time environmental perception data collected by an airspace monitoring sensor network. The set of real-time environmental perception data includes electromagnetic spectrum signal distribution information, meteorological fluid dynamic field distribution information, and geographic obstacle boundary contour information. The real-time environmental perception data set is subjected to environmental element decoupling processing to obtain the electromagnetic propagation attenuation potential field and the atmospheric turbulence disturbance potential field, and a three-dimensional airspace navigability constraint boundary is established based on the geographical obstacle boundary contour information. Using the three-dimensional airspace accessibility constraint boundary as a spatial constraint condition, the electromagnetic propagation attenuation potential field and the atmospheric turbulence disturbance potential field are superimposed with a multi-physics field to generate multipath environmental shielding effect distribution data of intruding UAVs within the three-dimensional airspace accessibility constraint boundary. Based on the multipath environmental masking effect distribution data, a passive detection feasibility assessment is performed on the preset flight path of the intrusion drone, generating a threat detectability level identifier and detection blind zone distribution boundary coordinates for the intrusion drone. Based on the threat detectability level identifier and the boundary coordinates of the detection blind zone distribution, the defense resource scheduling decision network is invoked to output a multi-dimensional defense resource collaborative deployment instruction. The multi-dimensional defense resource collaborative deployment instruction includes the beamforming parameter group of the detection node and the sequence of maneuvering and stationing waypoints of the interception unit.

2. The method for UAV threat assessment and defense decision-making combining environmental perception as described in claim 1, characterized in that, The process of decoupling environmental elements from the real-time environmental perception data set yields the electromagnetic propagation attenuation potential field and the atmospheric turbulence disturbance potential field. Furthermore, a three-dimensional airspace navigability constraint boundary is established based on the geographical obstacle boundary contour information, including: The electromagnetic spectrum signal distribution information is analyzed by spectral domain segmentation, and the signal intensity spatial gradient distribution and spectrum occupancy time evolution sequence in each spectral segment are extracted. The signal intensity spatial gradient distribution and the spectrum occupancy time evolution sequence are then mapped in a time-frequency manner to generate an electromagnetic propagation attenuation state field. The wind vector field curl divergence decomposition is performed on the meteorological fluid dynamic field distribution information to extract the spatiotemporal distribution of the rotation component and the spatiotemporal distribution of the divergence component of the wind vector field. The spatiotemporal distribution of the rotation component and the spatiotemporal distribution of the divergence component are input into the atmospheric refractive index structure constant field reconstruction step to generate the atmospheric turbulence disturbance potential field. The boundary contour information of the geographic obstacles is reconstructed using a three-dimensional convex hull. The vertex coordinate sequence and the set of face normal vectors of each obstacle in the boundary contour information are extracted. The vertex coordinate sequence and the set of face normal vectors are then stitched together with a triangular mesh to generate a three-dimensional spatial traversability constraint boundary. The electromagnetic propagation attenuation potential field is used as the first physical field input data, the atmospheric turbulence disturbance potential field is used as the second physical field input data, and the three-dimensional spatial mobility constraint boundary is used as the spatial filtering condition. The first physical field input data and the second physical field input data are superimposed on a voxel-by-voxel basis to generate a superimposed physical field data set. Spatial gradient direction field extraction is performed on the superimposed physical field data set to obtain the electromagnetic attenuation gradient direction and atmospheric refraction gradient direction at each voxel position. The electromagnetic attenuation gradient direction and the atmospheric refraction gradient direction are spatially registered, and the spatial angle relationship and relative change trend between the two are analyzed to generate composite situation distribution data that reflects the combined influence of the atmospheric environment on electromagnetic propagation. Based on the distribution data of the direction vector of the environmental masking effect and the preset flight path of the intruding drone, the path intersection point is detected. Each intersection point of the preset flight path and the distribution data of the direction vector of the direction vector of the environmental masking effect is taken as a sampling point. The magnitude of the direction vector of the environmental masking effect at each sampling point is extracted to generate a sequence of masking effect intensity along the path.

3. The method for UAV threat assessment and defense decision-making combining environmental perception as described in claim 1, characterized in that, The process involves using the three-dimensional airspace mobility constraint boundary as a spatial constraint condition, superimposing the electromagnetic propagation attenuation potential field and the atmospheric turbulence disturbance potential field using a multiphysics field, and generating multipath environmental masking effect distribution data for intruding UAVs within the three-dimensional airspace mobility constraint boundary, including: Using the three-dimensional spatial mobility constraint boundary as a spatial filtering condition, the electromagnetic propagation attenuation state field data located inside the three-dimensional spatial mobility constraint boundary is extracted from the electromagnetic propagation attenuation state field, and the electromagnetic propagation attenuation state field data is used as the first physical field effective data. Using the three-dimensional spatial mobility constraint boundary as a spatial filtering condition, the data inside the atmospheric turbulence disturbance potential field located inside the three-dimensional spatial mobility constraint boundary is extracted from the atmospheric turbulence disturbance potential field, and the data inside the atmospheric turbulence disturbance potential field is used as the effective data of the second physical field. Spatial registration and alignment are performed on the first and second physical field effective data. Based on the preset physical field coupling mapping relationship, the interaction between electromagnetic propagation attenuation characteristics and atmospheric turbulence disturbance characteristics at each voxel position is analyzed to generate a comprehensive environmental interference situation data set. Based on the comprehensive environmental interference situation data set, the spatial variation characteristics of electromagnetic propagation attenuation and atmospheric turbulence disturbance characteristics are extracted respectively. By analyzing the spatial correlation and relative orientation relationship between the two, comprehensive distribution data of environmental shielding effect is generated. Based on the comprehensive distribution data of the environmental masking effect and the preset flight path of the intruding drone, path intersection point detection is performed. Each intersection point of the preset flight path and the distribution data of the direction vector of the environmental masking effect is taken as a sampling point, and the magnitude and direction angle of the direction vector of the environmental masking effect at each sampling point are extracted. The magnitude of the environmental shading effect direction vector at each sampling point is arranged in order along the preset flight path to generate a sequence of shading effect intensity along the path. The direction angle of the environmental shading effect direction vector at each sampling point is arranged in order along the preset flight path to generate a sequence of shading effect direction angle along the path. Based on the shading effect intensity sequence and the shading effect direction angle sequence along the path, each sampling point on the preset flight path is marked with a dual-dimensional label of shading effect amplitude and shading effect direction to generate multi-path environmental shading effect distribution data.

4. The method for UAV threat assessment and defense decision-making combining environmental perception as described in claim 1, characterized in that, The step of passively detecting the intrusion drone's preset flight path based on the multi-path environmental masking effect distribution data, and generating the intrusion drone's threat detectability level identifier and detection blind zone distribution boundary coordinates, includes: Extract the shading effect intensity sequence and the shading effect direction angle sequence along the path from the multipath environment shading effect distribution data. Compare the comprehensive interference intensity quantization value of each sampling point in the shading effect intensity sequence along the path with a preset detection sensitivity threshold point by point. Mark the sampling points with shading effect intensity lower than the preset detection sensitivity threshold as candidate points for detection blind zone. The dominant interference direction feature of each sampling point in the occlusion effect direction angle sequence along the path is calculated point by point with the preset detection beam coverage direction range, and sampling points whose occlusion effect direction angle exceeds the preset detection beam coverage direction range are marked as unreachable candidate points. Spatial intersection analysis is performed between the candidate points of the detection blind zone and the candidate points of the direction unreachable, and sampling points that simultaneously meet the conditions of the candidate points of the detection blind zone and the conditions of the candidate points of the direction unreachable are extracted as the set of effective points of the detection blind zone. Spatial connectivity analysis is performed on the set of effective points in the detection blind zone. Spatially adjacent effective points in the detection blind zone are aggregated into connected domains of the detection blind zone. The boundary coordinate sequence of each connected domain of the detection blind zone is extracted to generate the boundary coordinates of the detection blind zone distribution. The proportion of the number of valid points in the detection blind zone to the total number of sampling points on the preset flight path is used as the detection blind zone coverage rate. A threat detectability level identifier is generated based on the preset coverage rate interval into which the detection blind zone coverage rate falls.

5. The method for UAV threat assessment and defense decision-making combining environmental perception according to claim 1, characterized in that, The step of invoking the defense resource scheduling decision network to output multi-dimensional defense resource collaborative deployment instructions based on the threat detectability level identifier and the boundary coordinates of the detection blind zone distribution includes: The threat detectability level identifier is input into the threat level embedding layer of the defense resource scheduling decision network, and the threat detectability level identifier is mapped into a threat level embedding vector. The detection blind zone distribution boundary coordinates are input into the spatial coordinate encoding layer of the defense resource scheduling decision network, and the detection blind zone distribution boundary coordinates are mapped into a spatial coordinate encoding vector. The threat level embedding vector and the spatial coordinate encoding vector are input into the feature cross-fusion layer of the defense resource scheduling decision network. Feature fusion and correlation analysis are performed on the threat level embedding vector and the spatial coordinate encoding vector to generate a joint threat representation vector containing spatial distribution characteristics. The threat joint characterization vector is input into the detection resource allocation branch of the defense resource scheduling decision network, and the detection node beamforming parameter set is output. The detection node beamforming parameter set includes the beam main lobe pointing angle, beamwidth range and transmit power gain coefficient of each detection node. The joint threat representation vector is input into the interception resource allocation branch of the defense resource scheduling decision network, and the interception unit maneuvering and stationing waypoint sequence is output. The interception unit maneuvering and stationing waypoint sequence includes the three-dimensional coordinates of the stationing point and the stationing duration of each interception unit. The beamforming parameter group of the detection node and the maneuvering waypoint sequence of the interception unit are synchronized in time to generate a multi-dimensional defense resource collaborative deployment command.

6. The method for UAV threat assessment and defense decision-making combining environmental perception according to claim 1, characterized in that, The method further includes: Acquire a set of real-time feedback data during the execution of the defense mission, the set of real-time feedback data including flight trajectory deviation data of the intrusion drone and signal reception power change data of the detection node; The actual flight path deviation of the intrusion drone is calculated based on the flight trajectory offset data of the intrusion drone. The actual flight path deviation is then input into the multipath environmental masking effect distribution data for rematching, resulting in updated multipath environmental masking effect distribution data corresponding to the actual flight path of the intrusion drone. The actual detection signal gain fluctuation of each detection node is calculated based on the signal reception power change data of the detection node. The actual detection signal gain fluctuation is compared with the updated multipath environmental shading effect distribution data to identify the prediction deviation area of ​​the environmental shading effect distribution data. Based on the predicted deviation area, the multipath environmental masking effect distribution data is locally corrected to obtain the corrected multipath environmental masking effect distribution data. Based on the corrected multipath environmental masking effect distribution data, the actual flight path of the intruding UAV is re-evaluated for passive detection feasibility, and the corrected threat detectability level identifier and the corrected detection blind zone distribution boundary coordinates are generated. The revised threat detectability level identifier and the revised detection blind zone distribution boundary coordinates are input again into the defense resource scheduling decision network, and the revised multi-dimensional defense resource collaborative deployment instruction is output. The modified multi-dimensional defense resource collaborative deployment command is used to update the beamforming parameter group of the currently executing detection node and the maneuvering waypoint sequence of the interception unit.

7. The method for UAV threat assessment and defense decision-making combining environmental perception as described in claim 1, characterized in that, The method further includes: The electromagnetic spectrum signal distribution information is processed to separate the signal sources and identify the electromagnetic radiation characteristic signals of the intruding drone and the electromagnetic radiation characteristic signals of the environmental background contained in the electromagnetic spectrum signal distribution information. The signal modulation mode features and signal pulse repetition interval features of the electromagnetic radiation characteristic signals of the intruding drone are extracted. The signal modulation mode features and signal pulse repetition interval features are compared with a preset drone electromagnetic fingerprint database to determine the drone type identifier of the intruding drone. Based on the spatial propagation attenuation distribution of the electromagnetic radiation characteristic signal of the intruding drone in the electromagnetic propagation attenuation state field, the signal coverage boundary of the electromagnetic radiation characteristic signal of the intruding drone within the three-dimensional airspace traversability constraint boundary is calculated, and the signal radiation coverage range of the intruding drone is generated. Spatial correlation analysis is performed between the signal radiation coverage area of ​​the intrusion drone and the distribution data of the multipath environmental shielding effect to determine the first signal attenuation area affected by the environmental shielding effect and the second signal fidelity area unaffected by the environmental shielding effect within the signal radiation coverage area of ​​the intrusion drone. Based on the spatial distribution of the first signal attenuation region and the second signal fidelity region, the beam main lobe pointing angle of each detection node in the detection node beamforming parameter group is adjusted so that the detection beam coverage is biased towards the second signal fidelity region. The multi-dimensional defense resource collaborative deployment command is updated based on the adjusted probe node beamforming parameter set.

8. The method for UAV threat assessment and defense decision-making combining environmental perception according to claim 1, characterized in that, The method further includes: The temporal evolution trend of the atmospheric turbulent disturbance potential field is analyzed, and the spatial structure change characteristics of the atmospheric turbulent disturbance potential field in continuous time segments are extracted to generate a temporal evolution sequence of the atmospheric turbulent disturbance potential field. Based on the time evolution sequence of the atmospheric turbulence disturbance potential field, the moving direction parameters and moving velocity parameters of the peak turbulence intensity region in the atmospheric turbulence disturbance potential field are calculated to generate migration trend data of the peak turbulence region. By performing a spatiotemporal intersection analysis on the migration trend data of the turbulence peak region and the preset flight path of the intruding drone, the time window and spatial location of the intruding drone crossing the peak region of the atmospheric turbulence disturbance potential field in the subsequent period can be predicted. Based on the time window and spatial location of the intruding drone crossing the peak region of the atmospheric turbulence disturbance potential field, the distribution data of the multipath environmental shielding effect is marked with timeliness, and the part that will be affected by the peak region of the atmospheric turbulence disturbance potential field is marked as the shielding effect decay region. Based on the marking results of the shading effect decay region, the shading effect intensity sequence along the preset flight path in the multipath environmental shading effect distribution data is modified by attenuation prediction to obtain the multipath environmental shading effect distribution data after the prediction period. Based on the multipath environmental shading effect distribution data after the prediction period, the boundary coordinates of the detection blind zone distribution are adjusted in advance to generate the adjusted boundary coordinates of the detection blind zone distribution. The adjusted boundary coordinates of the detection blind zone are input into the defense resource scheduling decision network, and the adjusted multi-dimensional defense resource collaborative deployment instructions are output.

9. The method for UAV threat assessment and defense decision-making combining environmental perception according to claim 1, characterized in that, The method further includes: The interference signal characteristic data actively emitted by the intrusion drone during flight is obtained, and the interference signal characteristic data includes the center frequency parameter, bandwidth parameter and power spectral density distribution of the interference signal; Based on the center frequency parameter and bandwidth parameter of the interference signal in the interference signal feature data, interference suppression processing is performed on the corresponding frequency band in the electromagnetic spectrum signal distribution information to obtain the interference-free electromagnetic spectrum signal distribution information. The electromagnetic spectrum signal distribution information after interference removal is re-analyzed in the spectrum domain, and the spatial gradient distribution of signal intensity in each spectrum segment after interference removal is extracted. The spatial gradient distribution of signal intensity after interference removal is used to replace the original spatial gradient distribution of signal intensity in the electromagnetic spectrum signal distribution information, and the electromagnetic propagation attenuation state field after interference removal is regenerated. Using the three-dimensional airspace navigability constraint boundary as the spatial constraint condition, the electromagnetic propagation attenuation potential field after interference removal and the atmospheric turbulence disturbance potential field are superimposed again by multi-physics fields to generate multi-path environmental shielding effect distribution data after interference removal. Based on the multipath environmental masking effect distribution data after interference removal, the feasibility of passive detection of the preset flight path of the intruding UAV is reassessed, and the threat detectability level identifier and the boundary coordinates of the detection blind zone distribution after interference removal are generated. The interference-free threat detectability level identifier and the interference-free detection blind zone distribution boundary coordinates are input into the defense resource scheduling decision network, and the interference-free multi-dimensional defense resource collaborative deployment instruction is output. Replace the currently executing multidimensional defense resource collaborative deployment instruction with the interference-free multidimensional defense resource collaborative deployment instruction.

10. A drone threat assessment and defense decision-making system combining environmental perception, characterized in that, The device includes a processor and a computer-readable storage medium storing machine-executable instructions, which, when executed by a computer, implement the UAV threat assessment and defense decision-making method combining environmental perception as described in any one of claims 1-9.

Images