An unmanned aerial vehicle state estimation method based on two-stage search

Abstract

This invention discloses a UAV state estimation method based on a two-stage search, relating to the field of low-altitude sensing and control technology. The method includes: constructing a three-dimensional reference coordinate system for a substation area, a set of no-entry flight zones, and a constraint grid table; generating an environmental constraint data packet; configuring and transmitting OTFS sensing frames based on the environmental constraint data packet, acquiring echo sequences, performing time synchronization, frequency offset correction, static clutter suppression, and frequency domain pre-whitening processing to generate a weighted observation sequence and a gridded environmental cost table; performing a coarse search using two-dimensional orthogonal matching pursuit under environmental constraints to obtain a unique candidate state group; performing a continuous domain fine search based on the unique candidate state group to complete position and velocity correction and generate current cycle state data; and generating initial state data for the next cycle based on the current cycle state data. This invention balances search efficiency and estimation accuracy through a two-stage search; and improves the accuracy of UAV state estimation through precise time delay and precise Doppler fitting.

Description

Technical Field

[0001] This invention relates to the field of low-altitude sensing and control technology, and in particular to a method for estimating the state of a UAV based on a two-stage search. Background Technology

[0002] With the development of low-altitude sensing and control and integrated sensing technologies, state estimation methods for UAV targets are gradually being applied to scenarios such as substations, park boundaries, and protection of key facilities. Conventional solutions typically rely on station area geometry information, echo signal acquisition, and time-delay-Doppler parameter calculation to estimate the target's distance, velocity, and motion state. Combined with pilot configuration, echo synchronization, clutter suppression, and parameter search, these processes generate dynamic perception results for UAV targets to support applications such as low-altitude target monitoring, trajectory tracking, and boundary protection.

[0003] However, the above-mentioned conventional methods still have two limitations in complex station area environments: on the one hand, propagation path obstruction, multipath reflection and frequency band interference often act simultaneously on echo observations, making it difficult for conventional search processes to take into account the differences in propagation reliability and spectral interference, thus affecting the stability of state estimation; on the other hand, single-stage parameter refinement methods are difficult to balance between search range control and estimation accuracy, which is not conducive to achieving a fine characterization of UAV state parameters while ensuring processing efficiency. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides a UAV state estimation method based on two-stage search to solve the problems of insufficient utilization of propagation environment constraints and insufficient coordination between search efficiency and estimation accuracy in the existing technology.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution: This invention provides a UAV state estimation method based on a two-stage search, comprising: constructing a three-dimensional reference coordinate system for a substation area, a set of no-entry flight zones, and a constraint grid table; acquiring the background electromagnetic spectrum; calculating the path obstruction rate, multipath attenuation factor, and power spectral density of interference signals; generating an environmental constraint data packet; configuring and transmitting OTFS sensing frames based on the environmental constraint data packet; acquiring echo sequences; sequentially performing time synchronization, frequency offset correction, static clutter suppression, and frequency domain pre-whitening processing; generating a weighted observation sequence and a grid environmental cost table; performing a coarse search of environmental constraint two-dimensional orthogonal matching pursuit using the weighted observation sequence and the grid environmental cost table to obtain a unique candidate state group; constructing a unique fine search interval based on the unique candidate state group; performing a continuous domain fine search to obtain fine time delay, fine Doppler, and fine scattering coefficients; and combining the constraint grid table to complete position correction and velocity correction, generating current cycle state data; calculating the beam correction angle and risk trigger marker based on the current cycle state data, and generating initial state data for the next cycle.

[0007] As a preferred embodiment of the UAV state estimation method based on two-stage search described in this invention, the generation of environmental constraint data packets includes: uniformly converting the substation electronic fence boundary, primary equipment layout, equipment height, and sensor base station installation location data into a unified three-dimensional reference coordinate system for the station area, constructing the equipment spatial outline, no-entry flight zone, and a three-dimensional constraint grid with unique attribute markers; establishing all candidate propagation paths that can enter the grid center, using the phase center of the sensor base station as the propagation starting point, sampling and querying grid attributes along the paths, and calculating the path occlusion rate and multipath attenuation factor; collecting background electromagnetic signals within a fixed pre-sampling window before the sensing period, performing frequency domain transformation, obtaining the interference signal power spectral density table, and finally generating the environmental constraint data packet.

[0008] As a preferred embodiment of the UAV state estimation method based on two-stage search described in this invention, the generation of weighted observation sequences and gridded environmental cost tables includes: calculating the joint optimal value of the path frequency band for each candidate propagation path at each discrete frequency point based on environmental constraint data packets; sorting the joint optimal value of the path frequency band from largest to smallest, and using the top-ranked discrete frequency points as pilot carrying frequency bands; selecting a fixed middle time slot as the pilot carrying time slot within each sensing period to form a pilot resource grid; filling constant modulus reference symbols into the pilot resource grid, filling communication service symbols into the communication data resource grid, filling zero-value symbols into the protection resource grid, and generating and transmitting OTFS sensing frames.

[0009] As a preferred embodiment of the UAV state estimation method based on two-stage search described in this invention, the generation of the weighted observation sequence and gridded environment cost table further includes: acquiring the original echo sequence, using the current period's transmission reference sequence to complete time synchronization and frequency offset correction; combining a static background template to eliminate stable background echo components generated by fixed reflectors, and performing differentiated weighting processing on the dynamic echo spectrum according to the interference signal power spectral density table to obtain the weighted observation sequence; mapping candidate propagation paths to time delay grids and Doppler grids according to the subcarrier spacing and time slot length of the current period's OTFS sensing frame; and summarizing the comprehensive path optimization values ​​falling into each time delay grid and Doppler grid to construct the gridded environment cost.

[0010] As a preferred embodiment of the UAV state estimation method based on two-stage search described in this invention, the step of obtaining a unique candidate state group includes: constructing a time delay and Doppler discrete search grid, and generating a two-dimensional discrete atom set based on the launch reference sequence; using the weighted observation sequence as the initial residual vector, calculating the environmental constraint matching score that combines grid environmental cost for all time delays and Doppler grid points, and selecting the grid index with the highest score to write into the support index set.

[0011] As a preferred embodiment of the UAV state estimation method based on two-stage search described in this invention, the step of obtaining a unique candidate state group further includes: extracting the discrete atom column vectors corresponding to the support index set from the two-dimensional discrete atom set to form a support sub-dictionary; inverting the coarse scattering coefficient vector using the least squares method; reconstructing the coarse echo and updating the residual vector, residual energy, and residual improvement ratio; stopping the coarse search based on the residual convergence threshold, the minimum improvement threshold, and the maximum allowed number of iterations; after the coarse search terminates, calculating the dominant candidate priority value based on the grid environmental cost and the coarse scattering coefficient vector; selecting the grid point with the largest dominant candidate priority value as the dominant grid point, and generating a unique candidate state group.

[0012] As a preferred embodiment of the UAV state estimation method based on two-stage search described in this invention, the acquisition of fine time delay, fine Doppler, and fine scattering coefficients includes: constructing a unique fine search time delay interval and Doppler interval based on the current period coarse time delay and coarse Doppler in the unique candidate state group; constructing a corresponding continuous domain echo response for any continuous time delay and continuous Doppler combination within the unique fine search interval; calculating the fine scattering coefficients based on the continuous domain echo response; calculating the continuous domain fine search cost through the fine scattering coefficients; and acquiring the fine time delay and fine Doppler by adopting a continuous domain single-peak fine search order of time delay first and Doppler second.

[0013] As a preferred embodiment of the UAV state estimation method based on two-stage search described in this invention, the generation of current cycle state data includes: calculating the predicted position of the current cycle according to the continuous position change relationship between adjacent sensing cycles; refitting the fine scattering coefficient based on fine time delay and fine Doppler to calculate the fine scattering coefficient energy and scattering stability of the current cycle; reading the constraint grid table, and extracting all accessible grids within a limited spatial neighborhood around the predicted position of the current cycle as the center, comparing the theoretical radial distance and the fine radial distance of the current cycle to obtain the optimal candidate grid, and taking the center of the optimal candidate grid as the corrected spatial position of the current cycle; calculating the state confidence value using the fine scattering coefficient energy and scattering stability of the current cycle, and generating the current cycle state data using the corrected spatial position and state confidence value of the current cycle.

[0014] As a preferred embodiment of the UAV state estimation method based on two-stage search described in this invention, the step of generating initial state data for the next cycle includes: constructing the line-of-sight direction for the current cycle based on the base station location and the spatial location corrected for the current cycle, and obtaining the spatial velocity corrected for the current cycle; calculating the intrusion trend index and trajectory continuity index based on the spatial velocity corrected for the current cycle; and calculating the boundary proximity based on the spatial location corrected for the current cycle.

[0015] As a preferred embodiment of the UAV state estimation method based on two-stage search described in this invention, the generation of initial state data for the next cycle further includes: calculating a comprehensive risk index based on intrusion trend indicators, trajectory continuity indicators, and boundary proximity, and generating a risk trigger marker for the current cycle; adjusting the coarse search range and fine search window width for the next cycle according to the risk trigger marker for the current cycle to form search window parameters for the next cycle; and generating initial state data for the next cycle using the corrected spatial position, corrected spatial velocity, state confidence value, and search window parameters for the current cycle.

[0016] The beneficial effects of this invention are as follows: By calculating the joint optimal value of the path frequency band, propagation reliability and frequency band interference intensity are considered simultaneously in the propagation path selection, thus optimizing the configuration of the pilot frequency band and the communication data frequency band; by accurately fitting the fine time delay and fine Doppler, high-precision estimation of time delay and Doppler is achieved, which greatly improves the accuracy of UAV state estimation and ensures precise tracking and control; by using a two-stage search method, candidate state intervals are quickly screened in the coarse search stage, ensuring computational efficiency while accurately narrowing the fine search range, thereby greatly improving estimation accuracy and search speed. Attached Figure Description

[0017] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0018] Figure 1 This is a flowchart of a UAV state estimation method based on two-stage search.

[0019] Figure 2 A flowchart for generating environmental constraint data packages.

[0020] Figure 3 A flowchart for generating weighted observation sequences and gridded environmental cost tables.

[0021] Figure 4 This is a flowchart of the coarse search and fine search. Detailed Implementation

[0022] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0023] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0024] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0025] Reference Figures 1-4 This is one embodiment of the present invention, which provides a UAV state estimation method based on two-stage search, including the following steps: S1. Construct a three-dimensional reference coordinate system for the substation area, a set of no-entry flight zones, and a constraint grid table. Collect the background electromagnetic spectrum, calculate the path obstruction rate, multipath attenuation factor, and power spectral density of interference signals, and generate an environmental constraint data package.

[0026] Data on the substation's electronic fence boundary, primary equipment layout, equipment height, and the installation location of the sensing base station are read. A unified three-dimensional reference coordinate system for the substation area is established, using the phase center of the sensing base station as the spatial reference origin. Under this unified three-dimensional reference coordinate system, the original plane coordinates and elevation data of the electronic fence boundary, main transformer, circuit breaker, disconnector, frame, busbar support, and lightning protection equipment are translated and oriented with the phase center of the sensing base station as the unified origin and the north and vertical directions as the unified reference directions. This results in the acquisition of the substation boundary contour and equipment spatial contour.

[0027] After completing the spatial mapping of the station area, an outer envelope is constructed for each type of primary equipment, and the outer surface of the outer envelope is expanded outward according to a safety buffer distance (e.g., 3 meters) to obtain the corresponding safety restricted areas for each equipment. All safety restricted areas are spatially merged to generate a station-wide restricted flight zone. Within the bounded area of ​​the station area's outer boundary, the entire station area space is divided into three-dimensional small volume regions of uniform size layer by layer, column by column, and row by row, according to the method of dividing at equal intervals along the length, width, and height directions. The center point of each three-dimensional small volume region is used as the grid center point to generate a three-dimensional constraint grid. For each grid center point, it is sequentially determined whether it is located within the electronic fence, whether it falls within the restricted flight zone, or whether it is located within the near-boundary buffer area outside the restricted flight zone, and a unique attribute label is assigned accordingly.

[0028] Using the phase center of the sensing base station as the propagation starting point, candidate propagation paths from the base station to the grid center are established for all accessible grid center points. Equal-interval sampling is performed along each candidate propagation path, and the attribute marker of the grid where each sampling point is located is queried segment by segment. Path segments passing through inaccessible grids are considered strong occlusion segments, and path segments passing through near-boundary grids are considered weak occlusion segments. To ensure that the occlusion amount simultaneously reflects both the physical obstruction effect of the equipment and the propagation instability effect caused by boundary proximity, the path occlusion rate is calculated, expressed as: ; in, Indicates the first The path occlusion rate of each candidate propagation path. Indicates the first The length of the strongly occluded segment that passes through an inaccessible grid in each candidate propagation path. Indicates the first The length of the weakly occluded segment that crosses the near-boundary grid in a candidate propagation path. This represents the weak occlusion reduction factor. Indicates the first The total length of the candidate propagation paths.

[0029] It should be noted that, The value is obtained by selecting near-boundary grid crossing samples and entity occlusion crossing samples, calculating the average echo amplitude attenuation of each sample, and then taking the ratio of the two values. The value range is [0,1].

[0030] Furthermore, a multipath attenuation factor is calculated for each candidate propagation path. Specifically, a path neighborhood channel with a fixed radius (e.g., 2 meters) is established with each path as the center line. The occupancy degree of the equipment envelope within the path neighborhood channel is statistically analyzed, and the multipath attenuation is constructed in conjunction with the total path length. Considering that multipath risk in substation scenarios is mainly caused by path growth and increased density of surrounding equipment, the expression for multipath attenuation is as follows: ; in, Indicates the first Multipath attenuation factor for candidate propagation paths Indicates the reference length constant. Indicates the first The occupancy rate of devices within the neighborhood channels of each candidate propagation path.

[0031] It should be noted that the reference length constant is obtained by statistically analyzing the total path length of all candidate propagation paths in the historical sensing cycle and taking the median as the reference length constant.

[0032] The control base station enters a silent reception state and continuously acquires background electromagnetic signals within a fixed pre-sampling window (e.g., 20 milliseconds) before the start of the sensing period (e.g., 100 milliseconds), without transmitting sensing pilots. The acquired background signals are divided into frames, and frequency domain transformation is performed on each frame. The energy at the same frequency point is averaged to obtain the background interference intensity distribution at each frequency point, forming an interference signal power spectral density table. After obtaining the path occlusion rate, multipath attenuation factor, and interference signal power spectral density, the three-dimensional reference coordinate information of the station area, the no-entry flight zone information, the constraint grid table, the candidate propagation path table, the path occlusion rate table, the path multipath attenuation factor table, and the interference signal power spectral density table are uniformly encapsulated to generate an environmental constraint data packet.

[0033] It should be noted that the sensing period refers to the time interval between the start times of two adjacent synesthetic frame transmissions; the fixed presampling window refers to a silent reception time segment intercepted before the transmission of the synesthetic frame in each sensing period.

[0034] S2. Configure and transmit OTFS sensing frames according to the environmental constraint data package, collect echo sequences, and sequentially perform time synchronization, frequency offset correction, static clutter suppression, and frequency domain pre-whitening processing to generate weighted observation sequences and gridded environmental cost tables.

[0035] To ensure that pilot resource allocation considers both propagation reliability and band interference intensity, a joint optimal value for path bands is defined based on environmentally constrained data packets. The expression is as follows: ; in, Indicates the first Candidate propagation paths and the first The joint optimal value of the path frequency band corresponding to each discrete frequency point Indicates the first Interference signal power spectral density at discrete frequency points This represents the interference normalization constant.

[0036] It should be noted that, The interference normalization constant is obtained by taking the average power spectral density of all discrete frequency interference signals.

[0037] The path frequency band joint optimization values ​​of the same discrete frequency point on all candidate propagation paths are aggregated at the frequency point level to obtain the comprehensive optimization value of the frequency point corresponding to each discrete frequency point. The comprehensive optimization values ​​of the frequency points corresponding to each discrete frequency point are sorted from largest to smallest, and several discrete frequency points with the highest ranking are selected as pilot carrier frequency bands. In each sensing period, a fixed middle time slot is selected as the pilot carrier time slot. Specifically, one OTFS sensing frame is transmitted in each sensing period. In a single OTFS sensing frame, a fixed number (e.g., 16 bits) of subcarriers with the highest path frequency band joint optimization values ​​are selected along the frequency direction as pilot frequency bands, and a fixed number (e.g., 4) of consecutive middle time slots are selected along the time direction as pilot time slots to form a pilot resource grid. In the remaining resource grids, the parts near the beginning and end of the frame are used as protection resource grids, and the remaining middle parts are used as communication data resource grids.

[0038] It should be noted that frequency-level aggregation takes the average of the joint optimal values ​​of the path frequency bands corresponding to all candidate propagation paths at the discrete frequency point as the comprehensive optimal value of the frequency point corresponding to the current discrete frequency point.

[0039] After the pilot resource grid is determined, constant modulus reference symbols are filled into the pilot resource grid, communication service symbols are filled into the communication data resource grid, and zero-value symbols are filled into the protection resource grid. The time-frequency resource grid is transformed into a time-domain transmission sequence according to the OTFS frame mapping method, and the transmission is completed based on the start time of the current sensing period. During transmission, the transmission reference sequence of this period is collected. After the OTFS sensing frame is transmitted, the base station immediately switches to the receiving state and continuously collects echo signals in the receiving time window within the current sensing period to obtain the original echo sequence. The original echo sequence is subjected to analog-to-digital conversion and digital down-conversion to obtain the discrete original receiving sequence.

[0040] It should be noted that the transmission reference sequence refers to the known time-domain signal sequence generated and actually transmitted by the sensing base station based on the OTFS sensing frame resource mapping result within the current sensing period, and simultaneously stored in the local cache.

[0041] The transmit reference sequence in the current cycle is read. Using the transmit reference sequence in the current cycle as the correlation template, the original echo sequence is slid point by point according to the sampling points. At each sliding position, the cross-correlation value between the echo segment with the same length as the transmit reference sequence at that position and the transmit reference sequence is calculated. The cross-correlation amplitudes corresponding to all sliding positions are arranged into a correlation curve in time order. The sampling position corresponding to the point with the largest amplitude in the correlation curve is determined as the correlation peak position, and the correlation peak position is used as the echo start position. The original echo sequence is truncated and aligned to complete time synchronization. The frequency offset is estimated by using the phase change of the echo corresponding to the pilot resource grid, and the frequency offset is compensated for on the time-synchronized echo sequence to obtain the time-synchronized and frequency offset corrected echo sequence.

[0042] By continuously collecting echo sequences for multiple sensing cycles when the UAV does not enter the monitoring airspace, and averaging the echo amplitude and phase at the corresponding sampling positions to obtain a static background template, the current cycle-corrected echo sequence and the static background template are aligned with the same length. Then, the current cycle-corrected echo sequence and the static background template are subtracted one by one at the same sampling position to eliminate the stable background echo components generated by fixed reflectors, and obtain the dynamic echo sequence after removing the fixed reflection components.

[0043] A frequency domain transformation is performed on the dynamic echo sequence to obtain the dynamic echo spectrum. Based on the power spectral density table of the interference signal in the environmental constraint data packet, suppression weights of different intensities are applied to different frequency points. Specifically, at each discrete frequency point, the frequency weight of that frequency point is set as the ratio of the interference normalization constant to the interference intensity at that frequency point. The dynamic echo spectrum is then multiplied by the corresponding frequency weight at each frequency point to suppress the echo components at high interference frequencies and relatively preserve the effective target components at low interference frequencies. The processed frequency domain result is then inversely transformed back to the time domain to obtain the weighted observation sequence.

[0044] Using the subcarrier spacing and time slot length of the current OTFS sensing frame, the propagation length corresponding to each candidate propagation path is mapped to a time delay grid, and the relative radial motion interval corresponding to the candidate propagation path is mapped to a Doppler grid. For each time delay and Doppler grid, the joint optimal value of the path frequency band of all candidate propagation paths falling into that grid is summarized, and the environmental cost of that grid is constructed, expressed as: ; in, Representing time delay grid points With Doppler grid points The corresponding grid environment cost, Indicates mapping to grid points The set of candidate propagation path indexes, Indicates the first The comprehensive path optimization value of the candidate propagation paths.

[0045] It should be noted that, Through the first The joint optimal values ​​of the path frequency bands corresponding to each candidate propagation path at each discrete frequency point are compared, and the maximum value is taken.

[0046] S3. Perform a coarse search of two-dimensional orthogonal matching pursuit with environmental constraints using a weighted observation sequence and a gridded environmental cost table to obtain a unique candidate state group.

[0047] Construct a discrete search grid for time delay and Doppler. Specifically, divide the time delay direction into discrete grid points consistent with the OTFS discrete time delay resolution, and divide the Doppler direction into discrete grid points consistent with the OTFS discrete Doppler resolution, forming a two-dimensional search grid. Combine each time delay grid point with each Doppler grid point in pairs to obtain all candidate time delay and Doppler grid point indices.

[0048] First, a discrete shift is applied to the corresponding time delay grid point of the transmitted reference sequence, and then a complex exponential modulation is applied to the corresponding Doppler grid point. The resulting echo response is expanded into a column vector with the same length as the weighted observation sequence, which serves as the discrete atom of that grid point. The discrete atoms corresponding to all grid points are then concatenated column by column to generate a two-dimensional discrete atom set.

[0049] The weighted observation sequence of the current period is directly used as the initial residual vector. The support index set is initialized to an empty set, and the number of iterations is initialized to zero. The environmental constraint matching score is calculated for each time delay and Doppler grid point. The expression is: ; in, Indicates the first Each sensing cycle at the time delay grid point and Doppler grid points Environmental constraint matching score Indicates the first Each sensing cycle at the time delay grid point and Doppler grid points The conjugate transpose of the discrete atomic vector at that location. Indicates the first The residual vector corresponding to the current iteration of each sensing cycle.

[0050] After completing the environmental constraint matching score for all time delays and Doppler grid points, the grid point index with the highest environmental constraint matching score is selected and written into the support index set. Once the support index set contains at least one grid point index, the discrete atom column vector corresponding to the support index set is extracted from the two-dimensional discrete atom set to form a support sub-dictionary. The current coarse scattering coefficient vector is inverted using the least squares method, with the current period weighted observation sequence as the fitted object. By minimizing the sum of squared residuals between the linear combination result of the weighted observation sequence and the support sub-dictionary, the coarse scattering coefficient vector corresponding to each selected grid point is obtained, and the support sub-dictionary is then compiled. The product of the coarse scattering coefficient vector and the coarse echo reconstruction is used as the coarse echo reconstruction for the current iteration. After the coarse echo reconstruction is completed, the difference between the current weighted observation sequence and the coarse echo reconstruction result is used as the updated residual vector. The updated residual energy and residual improvement ratio are further calculated. The current period weighted observation sequence is subtracted from the coarse echo reconstruction result of the current iteration point by point to obtain the updated residual vector. The squares of the amplitudes of each sampling point of the residual vector are summed to obtain the updated residual energy. The ratio of the decrease in residual energy in this iteration relative to the residual energy in the previous iteration to the initial residual energy is used as the residual improvement ratio.

[0051] It should be noted that the data obtained by summing the squares of the amplitudes of each sampling point in the weighted observation sequence is used as the initial residual energy.

[0052] The coarse search stops when the residual energy drops below the residual convergence threshold, or when the residual improvement brought by the current iteration is less than the minimum improvement threshold; if neither condition is met, the search continues until the maximum allowed number of iterations is reached.

[0053] It should be noted that the residual convergence threshold is obtained by statistically analyzing the residual energy of multiple targetless perception cycles before the UAV enters the monitoring airspace, and taking the maximum value as the residual convergence threshold, with a value range of [0,1]. The minimum improvement threshold is obtained by statistically analyzing the residual improvement ratio of two adjacent iterations in the historical target perception cycle, and taking the minimum residual improvement ratio in the continuous stable convergence phase, with a value range of [0,1].

[0054] After the coarse search terminates, in order to select a unique dominant grid point from the final support index set, the priority value of the dominant candidate is calculated, expressed as: ; in, Indicates the first In the first sensing cycle The dominant candidate priority value for each selected grid point This represents the final coarse scattering coefficient corresponding to the selected lattice point. This represents the grid environment cost corresponding to the selected grid point. This represents the final residual energy at the termination of the coarse search. This represents the initial residual energy.

[0055] After obtaining the dominant candidate priority values ​​of all selected grid points, the selected grid point with the largest dominant candidate priority value is selected as the dominant grid point of the current period, and the discrete time delay, discrete Doppler and coarse scattering coefficients corresponding to the dominant grid point are output in a unified manner to generate a unique candidate state group.

[0056] It should be noted that the discrete time delay and discrete Doppler corresponding to the dominant grid point are obtained by reading the time delay grid point index and Doppler grid point index of the dominant grid point in the time delay and Doppler discrete search grid, respectively; the coarse scattering coefficient corresponding to the dominant grid point is obtained by extracting the coefficient component corresponding to the dominant grid point index from the final coarse scattering coefficient vector.

[0057] S4. Construct a unique fine search interval based on the unique candidate state group, perform a continuous domain fine search, obtain the fine time delay, fine Doppler and fine scattering coefficients, and combine with the constraint grid table to complete the position correction and velocity correction, and generate the current period state data.

[0058] The system reads the unique candidate state group output from the previous step, which includes the current cycle coarse delay, coarse Doppler, and coarse scattering coefficient. Simultaneously, it reads the previous cycle's state result from the historical state cache, which includes the spatial position, spatial velocity, state confidence value, and beam pointing information saved in the previous sensing cycle. Based on the previous cycle's spatial position and spatial velocity, and following the continuous positional change relationship between adjacent sensing cycles, the system calculates the current cycle's predicted position using the previous cycle's spatial position in the historical state cache as the starting point and the product of the previous cycle's spatial velocity and the current sensing cycle's duration as the position increment. The system then directly uses the previous cycle's spatial velocity as the current cycle's predicted velocity. Using the sensing base station coordinates as a reference, the system calculates the predicted radial distance and predicted radial velocity, which are then converted into predicted delay and predicted Doppler.

[0059] After obtaining the current period's coarse time delay, coarse Doppler, prediction time delay, and prediction Doppler, two time delay intervals are constructed around the coarse time delay and prediction time delay, respectively, and two Doppler intervals are constructed around the coarse Doppler and prediction Doppler, respectively. Then, the intersection of the two types of intervals is calculated to obtain the unique fine search time delay interval and the unique fine search Doppler interval. Specifically, the width of the time delay interval constructed around the coarse time delay is twice the size of one discrete time delay step, and the width of the time delay interval constructed around the prediction time delay is twice the size of two discrete time delay steps; the width of the Doppler interval constructed around the coarse Doppler is twice the size of one discrete Doppler step, and the width of the Doppler interval constructed around the prediction Doppler is twice the size of two discrete Doppler steps. After this processing, if the coarse estimation center and the prediction center are close, the intersection interval will naturally shrink; if there is a certain deviation between the two, the intersection interval will still retain sufficient search space.

[0060] After establishing the unique fine search interval, the local and current-period weighted observation sequences of the transmitted reference sequence are read. For any combination of continuous time delay and continuous Doppler within the unique fine search interval, the corresponding continuous domain echo response is constructed. Based on the fitting relationship between the continuous domain echo response and the current-period weighted observation sequence, the fine scattering coefficients under this combination are calculated. The fine scattering coefficients are obtained by performing complex amplitude phase fitting on a single continuous domain response and the current-period weighted observation sequence. To ensure that the fine search considers not only the residual fitting effect but also environmental reliability and dynamic continuity, the continuous domain fine search cost is defined as follows: ; in, Indicates the first Each sensing cycle has a continuous time delay and continuous Doppler The cost of a fine search in a continuous domain. Indicates the first A weighted observation sequence for each sensing cycle, Indicates continuous delay and continuous Doppler The fine scattering coefficients obtained by fitting under the given conditions Indicates continuous delay and continuous Doppler The corresponding continuous domain response vector, Indicates the value of continuous environmental replacement. Indicates the dynamic deviation. Indicates the environmental consistency enhancement coefficient. Indicates the dynamic continuity enhancement coefficient. Indicates the environmental compressibility factor. This represents the dynamic compression factor.

[0061] It should be noted that, It is a vector formed by generating echo response signals according to the corresponding propagation path at each combination of time delay and Doppler within the time delay and Doppler search interval, and splicing all echo response signals according to the time domain sequence; By determining continuous delay and continuous Doppler The four adjacent time delay and Doppler discrete grid points that fall into the grid are obtained by bilinear interpolation and normalization of the grid environmental cost corresponding to the four time delay and Doppler discrete grid points according to the relative distance in the time delay direction and the Doppler direction. The sum of squares is obtained by subtracting the current continuous delay and the current continuous Doppler from the predicted delay and the predicted Doppler, respectively, and then normalizing them according to the half-width of the corresponding delay search interval and the half-width of the Doppler search interval. This was achieved by keeping other parameters constant in historical sensing data, gradually adjusting the multiplier before the environmental cost term, and selecting the multiplier that minimizes the fine search error. This was achieved by keeping other parameters constant in historical sensing data, gradually adjusting the multiplier before the dynamic deviation term, and selecting the multiplier that minimizes the trajectory continuity error. and The values ​​range from 0.5 to 2; The environmental cost value index was obtained by gradually adjusting the decay slope of the term in historical sensing data and selecting the slope that maximizes the distinction between high-confidence and low-confidence environmental regions. The dynamic deviation index was obtained by gradually adjusting the decay slope of the term in historical sensing data and selecting the slope that maximizes the distinction between the predicted center neighborhood and the deviation neighborhood. and The values ​​range from 1 to 5.

[0062] After establishing the continuous domain fine search cost, a continuous domain single-peak fine search sequence of first time delay and then Doppler is adopted. Specifically, the Doppler is first fixed near the coarse Doppler, and the time delay search interval is gradually narrowed within the unique fine search time delay interval to find the time delay value that minimizes the continuous domain fine search cost, thus obtaining the fine time delay. The time delay is then fixed as the fine time delay, and the Doppler search interval is gradually narrowed within the unique fine search Doppler interval to find the Doppler value that minimizes the continuous domain fine search cost, thus obtaining the fine Doppler. The fine scattering coefficients are then refitted onto the continuous domain response corresponding to the fine time delay and fine Doppler to obtain the fine scattering coefficients for the current period. The fine time delay is then... The radial distance is converted to the current period's radial distance, and the fine Doppler is converted to the current period's radial velocity. The constrained grid table is read, and all accessible grids are extracted within a limited spatial neighborhood around the predicted position of the current period. For each accessible grid, the theoretical radial distance of the candidate grid center relative to the sensing base station is calculated, and the theoretical radial distance is compared with the current period's fine radial distance to calculate the spatial offset between the candidate grid center and the current period's predicted position. The grid with the smallest radial distance difference and the smallest spatial offset is selected as the optimal candidate grid, and the center of the optimal candidate grid is taken as the corrected spatial position for the current period.

[0063] To unify radial consistency, spatial continuity, and boundary safety into a single position correction value, a position correction cost is defined as follows: ; in, Indicates the first The first sensing cycle The cost of position correction for each local candidate raster. Indicates the first The theoretical radial distance of a local candidate grid relative to the sensing base station. Indicates the first The precise radial distance per sensing cycle, Indicates the distance normalization scale. Indicates the first The square of the Euclidean distance between the center of each local candidate grid and the current periodic prediction position. This represents the square of the local search radius. Indicates the first Attribute tags for local candidate rasters, Represents the predicted neighborhood constraint coefficient. This represents the raster attribute penalty coefficient.

[0064] It should be noted that, The distance normalization scale is obtained by statistically analyzing the deviation distribution between the theoretical radial distance and the precise radial distance of local candidate grids within the historical sensing period, and the median value of the deviation distribution is taken as the distance normalization scale. The value is obtained by gradually adjusting the multiplier before the spatial offset term during the historical sensing cycle and selecting the multiplier that minimizes the position correction error. The value range is from 0.5 to 2. The multiplier before the penalty term for near-boundary grids is gradually adjusted during the historical sensing cycle, and the multiplier that makes the correction result fall into the grids with the highest proportion is selected. The value range is from 0.5 to 3. By judging the first The spatial region of the center point of each grid is obtained. When the grid is located within the restricted flight zone, the grid attribute is marked as 0. When it is located within the buffer zone of the outer boundary of the restricted flight zone, the grid attribute is marked as 1. When it is located in other passable areas, the grid attribute is marked as 2. The local search radius is obtained by statistically analyzing the Euclidean distance distribution between the predicted position of the current period and the final corrected spatial position in the historical sensing period. The 90th percentile of the Euclidean distance distribution is taken as the local search radius. If the quantile is too low, some normal prediction deviations may fall outside the local search area, resulting in the omission of candidate grids. If the quantile is too high, a small number of abnormal deviations may be included, resulting in an excessively wide local search area, which weakens the local constraint of position correction.

[0065] The current period's line-of-sight direction is constructed based on the base station location and the spatially corrected location for the current period. Specifically, the coordinate difference vector between the location of the sensing base station and the spatially corrected location for the current period is used as the current period's line-of-sight direction vector. After normalizing the coordinate difference vector according to the magnitude, the unit direction vector of the current period's line-of-sight direction is obtained. The predicted velocity of the previous period is decomposed into a velocity component along the line-of-sight direction and a velocity component orthogonal to the line-of-sight direction. The line-of-sight direction component in the predicted velocity is replaced with the current period's fine radial velocity, while keeping the velocity component orthogonal to the line-of-sight direction unchanged, to obtain the spatially corrected velocity for the current period.

[0066] Based on the relative relationship between the fine scattering coefficient energy and the final residual energy, the scattering stability of the current period is calculated; the state confidence value of the current period is calculated, expressed as: ; in, This represents the state confidence value for the current period.

[0067] The current period corrected spatial position, current period corrected spatial velocity, current period fine radial distance, current period fine radial velocity, current period fine scattering coefficient, current period scattering stability, and current period state confidence value are all unified as the current period state data.

[0068] It should be noted that the fine scattering coefficient energy is obtained by taking the square of the modulus of the fine scattering coefficient for the current period; the scattering stability for the current period is obtained by calculating the ratio of the fine scattering coefficient energy to the sum of the fine scattering coefficient energy and the reconstruction residual energy per unit length.

[0069] S5. Calculate the beam correction angle and risk trigger marker based on the current cycle state data, and generate the initial state data for the next cycle.

[0070] Using the current period-corrected spatial position as the boundary assessment input, the current period-corrected spatial velocity as the intrusion trend assessment input, and the current period state confidence value and current period scattering stability as control confidence constraint inputs, a nearest boundary search is performed around the grid corresponding to the current period-corrected spatial position. This yields the minimum distance from the current position to the no-entry zone boundary and the minimum distance to the safety buffer boundary, respectively. To increase the risk response when the target approaches the sensitive boundary, the current period boundary proximity is defined, expressed as: ; in, Indicates the first The proximity of the boundary of each perception cycle This represents the minimum distance from the current periodically corrected spatial position to the nearest no-entry flight zone boundary. This represents the minimum distance from the current periodically corrected spatial location to the nearest safe buffer boundary. This represents the normalization constant for the distance to the forbidden boundary. This represents the normalized constant of the safety buffer boundary distance.

[0071] It should be noted that the safety buffer boundary is the outer boundary formed by extending the outer surface of the restricted flight zone outward along the outer normal direction by a safety buffer distance; It is obtained by the average spacing between the centers of adjacent forbidden boundary grids in the statistical constraint grid table; It is obtained by statistically analyzing the average spacing between the centers of adjacent boundary grids on the safety buffer boundary.

[0072] The intrusion direction is constructed by pointing the current period-corrected spatial position to the nearest no-entry zone boundary point, and the normalized projection of the current period-corrected spatial velocity on the intrusion direction is calculated to obtain the current period intrusion trend index. The current period-corrected spatial velocity is compared with the previous period-corrected spatial velocity component by component, and the trajectory continuity index is calculated based on the velocity change amplitude.

[0073] It should be noted that the trajectory continuity index is obtained by calculating the velocity difference between the space velocity after the current cycle correction and the space velocity after the previous cycle correction in each direction, summing the squares of the velocity differences in each direction to obtain the velocity change, and then normalizing and compressing the velocity change to obtain the trajectory continuity index.

[0074] Based on the coordinate difference between the coordinates of the sensing base station and the spatial position after the current cycle correction, calculate the target horizontal angle and target elevation angle for the current cycle; subtract the target horizontal angle for the current cycle from the beam horizontal pointing angle after the previous cycle to obtain the beam horizontal correction angle for the current cycle; subtract the target elevation angle for the current cycle from the beam elevation pointing angle after the previous cycle to obtain the beam elevation correction angle for the current cycle.

[0075] Define the comprehensive risk index as follows: ; in, Indicates the first A comprehensive risk index for each perception cycle. Indicates the first Intrusion trend indicators for each perception cycle Indicates the first State confidence value for each sensing cycle. Indicates the first Scattering stability per sensing cycle, Indicates the first Trajectory continuity index for each sensing cycle.

[0076] When the comprehensive risk index is lower than the first risk threshold, the risk trigger flag is set to 0; when the comprehensive risk index is not lower than the first risk threshold and is lower than the second risk threshold, the risk trigger flag is set to 1; when the comprehensive risk index is not lower than the second risk threshold, the risk trigger flag is set to 2.

[0077] It should be noted that the first risk threshold and the second risk threshold are obtained by statistically analyzing the comprehensive risk index of the historical sensing period, and taking the average risk index at the boundary between normal tracking samples and early warning samples as the first risk threshold, and taking the average risk index at the boundary between early warning samples and high-risk samples as the second risk threshold. The values ​​of the first risk threshold and the second risk threshold are both greater than 0 and less than 1, and the first risk threshold is less than the second risk threshold.

[0078] The beam horizontal correction angle for the current cycle is superimposed on the beam horizontal pointing angle after the previous cycle, and the beam elevation correction angle for the current cycle is superimposed on the beam elevation pointing angle after the previous cycle, thus obtaining the beam horizontal pointing angle and beam elevation pointing angle after the current cycle. Based on the risk trigger flag for the current cycle, the coarse search range and fine search window width for the next cycle are adjusted. Specifically, when the risk trigger flag is 0, the half-width of the coarse search range for the next cycle is maintained at 1 times the discrete time delay step size and discrete Doppler step size of the current cycle, and the half-width of the fine search window for the next cycle is maintained at the current cycle's coarse search range... The coarse search range half-width is increased by 0.5 times; when the risk trigger flag is 1, the coarse search range half-width of the next cycle is increased by a fixed multiple (e.g., 2 times) of the current cycle's discrete delay step size and discrete Doppler step size, and the fine search window half-width of the next cycle is increased by a fixed multiple (e.g., 0.75 times) of the current cycle's coarse search range half-width; when the risk trigger flag is 2, the coarse search range half-width of the next cycle is increased by a fixed multiple (e.g., 3 times) of the current cycle's discrete delay step size and discrete Doppler step size, and the fine search window half-width of the next cycle is increased by a fixed multiple (e.g., 1 time) of the current cycle's coarse search range half-width.

[0079] The corrected spatial position, corrected spatial velocity, state confidence value, scattering stability, horizontal beam pointing angle, and elevation beam pointing angle after the current cycle are all written into the historical state cache to form the initial state data for the next cycle.

[0080] In summary, this invention, through path and frequency band joint optimization value calculation, simultaneously considers propagation reliability and frequency band interference intensity in propagation path selection, optimizing the configuration of pilot and communication data frequency bands; through precise fitting of fine time delay and fine Doppler, it achieves high-precision estimation of time delay and Doppler, significantly improving the accuracy of UAV state estimation and ensuring precise tracking and control; through a two-stage search method, it achieves rapid screening of candidate state intervals in the coarse search stage, ensuring computational efficiency while precisely narrowing the fine search range, thereby significantly improving estimation accuracy and search speed.

[0081] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for estimating the state of a UAV based on a two-stage search, characterized in that, include: Construct a three-dimensional reference coordinate system for the substation area, a set of no-entry flight zones, and a constraint grid table; collect the background electromagnetic spectrum; calculate the path obstruction rate, multipath attenuation factor, and power spectral density of interference signals; and generate environmental constraint data packets. Configure and transmit OTFS sensing frames according to environmental constraint data packets, collect echo sequences, and sequentially perform time synchronization, frequency offset correction, static clutter suppression, and frequency domain pre-whitening processing to generate weighted observation sequences and gridded environmental cost tables. A coarse search for environmental constraints is performed using a weighted observation sequence and a gridded environmental cost table to obtain a unique candidate state group; A unique fine search interval is constructed based on a unique candidate state group. A continuous domain fine search is performed to obtain fine time delay, fine Doppler and fine scattering coefficients. The position correction and velocity correction are completed by combining the constraint grid table to generate the current period state data. The beam correction angle and risk trigger marker are calculated based on the current cycle state data to generate the initial state data for the next cycle.

2. The UAV state estimation method based on two-stage search as described in claim 1, characterized in that, The generated environment constraint data package includes: The data on the substation electronic fence boundary, primary equipment layout, equipment height, and sensor base station installation location are uniformly converted into a unified three-dimensional reference coordinate system for the station area to construct the equipment spatial outline, no-entry flight zone, and three-dimensional constraint grid with unique attribute markers. Using the phase center of the sensing base station as the propagation starting point, establish all candidate propagation paths that can enter the grid center, sample and query grid attributes along the path, and calculate the path occlusion rate and multipath attenuation factor. Background electromagnetic signals are collected within a fixed pre-sampling window before the sensing period, frequency domain transformation is performed, the power spectral density table of interference signals is obtained, and finally, environmental constraint data packets are generated.

3. The UAV state estimation method based on two-stage search as described in claim 2 or 1, characterized in that, The generation of the weighted observation sequence and gridded environment cost table includes: Calculate the joint optimal value of the path frequency band for each candidate propagation path at each discrete frequency point based on environmentally constrained data packets; The path frequency band joint optimization values ​​are sorted from largest to smallest, and the top discrete frequency points are used as pilot carry frequency bands. Within each sensing cycle, a fixed middle time slot is selected as the pilot carry time slot to form a pilot resource grid. Fill in constant modulus reference symbols on the pilot resource grid, fill in communication service symbols on the communication data resource grid, fill in zero-value symbols on the protection resource grid, and generate and transmit OTFS inductive frames.

4. The UAV state estimation method based on two-stage search as described in claim 3, characterized in that, The generation of the weighted observation sequence and gridded environment cost table also includes: The original echo sequence is acquired, and time synchronization and frequency offset correction are completed using the transmission reference sequence of this cycle. By combining a static background template to eliminate the stable background echo components generated by fixed reflectors, and by performing differential weighting processing on the dynamic echo spectrum based on the power spectral density table of the interference signal, a weighted observation sequence is obtained. Based on the subcarrier spacing and slot length of the current OTFS sensing frame, the candidate propagation paths are mapped to delay grid points and Doppler grid points; The comprehensive path optimization values ​​falling into each time delay grid point and Doppler grid point are summarized to construct the grid point environmental cost.

5. The UAV state estimation method based on two-stage search as described in claim 4 or 1, characterized in that, The process of obtaining a unique candidate state group includes: Construct a time-delay and Doppler discrete search grid, and generate a two-dimensional discrete atom set based on the emission reference sequence; Using the weighted observation sequence as the initial residual vector, an environmental constraint matching score is calculated for all time delays and Doppler grid points, combining the grid point environmental cost. The grid point index with the highest score is then written into the support index set.

6. The UAV state estimation method based on two-stage search as described in claim 5, characterized in that, The process of obtaining a unique candidate state group also includes: The discrete atom column vectors corresponding to the support index set are extracted from the two-dimensional discrete atom set to form a support sub-dictionary. The coarse scattering coefficient vector is inverted using the least squares method to reconstruct the coarse echo and update the residual vector, residual energy and residual improvement ratio. The coarse search is stopped based on the residual convergence threshold, the minimum improvement threshold, and the maximum allowed number of iterations. After the coarse search terminates, the dominant candidate priority value is calculated based on the grid environment cost and the coarse scattering coefficient vector; Select the grid point with the highest priority value as the dominant grid point to generate a unique candidate state group.

7. The UAV state estimation method based on two-stage search as described in claim 1, characterized in that, The acquisition of fine time delay, fine Doppler, and fine scattering coefficients includes: Based on the current period coarse time delay and coarse Doppler in the unique candidate state group, construct a unique fine search time delay interval and Doppler interval; For any combination of continuous time delay and continuous Doppler within the unique fine search interval, construct the corresponding continuous domain echo response, and calculate the fine scattering coefficient based on the continuous domain echo response; The cost of fine search in the continuous domain is calculated by using fine scattering coefficients, and fine time delay and fine Doppler are obtained by adopting a continuous domain single-peak fine search order of time delay first and Doppler later.

8. The UAV state estimation method based on two-stage search as described in claim 7, characterized in that, The generation of current cycle status data includes: Calculate the predicted position for the current period based on the continuous positional change relationship between adjacent sensing periods; Based on the fine time delay and fine Doppler refitting of the fine scattering coefficient, the fine scattering coefficient energy and scattering stability of the current period are calculated. Read the constrained raster table, and extract all accessible rasters within the limited spatial neighborhood around the current period prediction position as the center. Filter the optimal candidate rasters by comparing the theoretical radial distance and the current period precise radial distance, and use the center of the optimal candidate rasters as the corrected spatial position for the current period. The state confidence value is calculated using the fine scattering coefficient energy and scattering stability of the current period, and the state data for the current period is generated using the corrected spatial location and state confidence value for the current period.

9. The UAV state estimation method based on two-stage search as described in claim 8, characterized in that, The generation of initial state data for the next cycle includes: Construct the current period line of sight direction based on the base station location and the current period-corrected spatial location, and obtain the current period-corrected spatial velocity; Intrusion trend indicators and trajectory continuity indicators are calculated using the spatial velocity after current cycle correction. The boundary proximity is calculated based on the spatial position after the current period correction.

10. The UAV state estimation method based on two-stage search as described in claim 9, characterized in that, The generation of initial state data for the next cycle also includes: A comprehensive risk index is calculated based on intrusion trend indicators, trajectory continuity indicators, and boundary proximity, and a risk trigger marker for the current period is generated. Adjust the coarse search range and fine search window width for the next period based on the current period's risk trigger flag to form the search window parameters for the next period; The initial state data for the next cycle is generated by using the corrected spatial position, corrected spatial velocity, state confidence value, and search window parameters for the next cycle.

Images