A Spatiotemporal Weighted Discriminant and Least Squares Fitting Method for Identifying Migratory Bird Flyways

Abstract

This invention discloses a spatiotemporal weighted discrimination and least squares fitting method for identifying migratory bird flyways. The method includes preprocessing trajectory data, using three-level adaptive interpolation to repair missing segments and calculate confidence levels; fusing motion features and confidence levels, calculating migration tendency scores using differentiable functions and hidden Markov models to accurately identify migration points and merge them into valid segments; spatial aggregation based on the symmetrical distance between segments, and determining complete flyways through endpoint pairing and trajectory overlap checks; fitting the flyway centerline using weighted least squares B-splines, determining boundaries based on residuals, and quantifying the geometric uncertainty of the flyway through Bootstrap resampling to output the flyway containing confidence intervals. This invention effectively handles missing data and noise, achieving high-precision and robust flyway identification and providing reliable statistical evaluation.

Description

Technical Field

[0001] This invention relates to the interdisciplinary fields of mobile biology, satellite tracking data processing, spatiotemporal data science and statistical fitting analysis, specifically a method for identifying migratory bird migration routes using spatiotemporal weighted discrimination and least squares fitting. Background Technology

[0002] With the widespread adoption and miniaturization of microsatellite locators and Global Positioning System (GPS) technologies, an increasing number of research institutions and conservation organizations are using satellite locators for long-term group tracking of migratory birds to obtain information such as their migration routes, stopover points, and habitat usage. Current technologies typically identify habitats based on cluster analysis of trajectory points, classify behavioral states using methods such as speed / turn angle thresholds or Hidden Markov Models, or employ simple spatial kernel density estimation to identify high-use areas. However, existing technologies have several shortcomings in identifying "migration channels" from discrete positioning points, mainly in the following aspects: First, the original trajectory data generally suffers from irregular sampling intervals, short-term or long-term missing data, and differences in positioning accuracy; traditional simple linear interpolation or time resampling methods cannot simultaneously take into account the continuity of dynamics and the biological rationality of behavioral patterns, and usually ignore the impact of interpolation uncertainty on subsequent analysis; second, existing channel identification relies heavily on point density or direction-based clustering methods, which are easily affected by isolated high-density segments, abnormal routes of single individuals, or equipment interruptions, leading to false channels or false mergers; third, existing estimates of the channel centerline and boundaries are mostly empirical descriptions, lacking a robust fitting process based on weighted least squares and statistical uncertainty quantification, making it difficult to provide managers with quantifiable confidence intervals and decision-making basis; finally, many methods rely heavily on manual verification or empirical thresholds in the channel determination process, making it difficult to promote automated and reproducible pipelines. Based on the above shortcomings, there is an urgent need for a method that can automatically and robustly identify true migration channels on large-scale, heterogeneous trajectory data, and simultaneously provide the channel centerline, width boundary and its uncertainty quantification, in order to support ecological protection and management decisions. Summary of the Invention

[0003] To address the shortcomings of existing technologies, this invention provides a method for identifying migratory bird migration routes using spatiotemporal weighted discrimination and least squares fitting, aiming to solve the problems in the background technology.

[0004] To achieve the above objectives, the present invention provides the following technical solution: a method for identifying migratory bird flyways based on spatiotemporal weighted discrimination and least squares fitting, comprising the following steps:

[0005] Step S1: Data preprocessing and three-level adaptive interpolation; The original migratory bird trajectory data is uniformly formatted and projected. For short-term or long-term missing segments in the original migratory bird trajectory data, a three-level adaptive interpolation strategy is used for interpolation reconstruction, and the interpolation confidence is calculated for each interpolation point.

[0006] Step S2, Spatiotemporal Feature Extraction and Migration Point Determination: Calculate the basic spatiotemporal features of each sampling point in the original migratory bird trajectory data after preprocessing and three-level adaptive interpolation, and introduce a comprehensive confidence level. Use a differentiable S-shaped function to map features, define a point-level migration tendency score, and combine it with the HMM posterior probability to obtain a fusion score. Determine migration points according to the threshold and merge them into migration segments, and select effective migration segments.

[0007] Step S3, Initial Judgment, Spatial Aggregation, Robustness Assessment and Allocation of Candidate Channels: Using the selected effective migration segments as units, calculate the symmetrical average minimum vertical distance between effective migration segments and construct a distance matrix. Aggregate them into candidate clusters according to the merging threshold. Perform endpoint pairing, time pairing and trajectory overlap checks on the segments within the candidate clusters. Determine the complete channel based on the number of strong pairings and the proportion of segment directions.

[0008] Step S4: Channel centerline fitting, merging, boundary decision and uncertainty quantification: For the point set of the candidate cluster determined to be a complete channel, the centerline is fitted using weighted least squares B-spline with point weights; the residual from the point to the centerline is calculated to determine the channel half-width and construct the boundary to form the fitted channel; adjacent fitted channels that meet the conditions are automatically merged; the empirical distribution of the centerline and boundary is obtained through Bootstrap resampling, the confidence band and channel confidence score are calculated, and the migration channel containing the centerline, boundary and confidence interval is output.

[0009] Furthermore, in step S1, the original migratory bird trajectory data is uniformly formatted and projected; the original migratory bird trajectory data includes the first... Original location sequence of a single migratory bird , Indicates the first Only one migratory bird Timestamps of each observation point; Indicates the first Only one migratory bird The planar coordinate projection of each observation point , They represent the first Only one migratory bird The horizontal and vertical coordinates of each observation point on the projection plane; Indicates the first Only one migratory bird Height information of each observation point; Indicates the first Only one migratory bird Positioning quality indicators for each observation point Indicates the first The total number of observation points for individual migratory birds.

[0010] Furthermore, in step S1, the irregularly sampled trajectories are uniformly resampled to the target time step. To address short-term or long-term missing segments in the original migratory bird trajectory data, a three-level adaptive interpolation strategy is employed for interpolation reconstruction, and the interpolation confidence level is calculated for each interpolation point. In the three-level adaptive interpolation strategy, the target time... interpolation position The interpolated continuous trajectory is given by a weighted average of neighborhood samples; neighborhood samples refer to samples taken at the target time. , used as the observation point for interpolation.

[0011] Furthermore, the three-level adaptive interpolation strategy includes:

[0012] When the time interval between the two endpoints of a missing segment is in the original migratory bird trajectory data At that time, cubic Hermite spline interpolation with endpoint velocity constraints is used: let the start endpoint and end endpoint of the missing segment be respectively... , The corresponding time is , Calculated by difference between adjacent points , speed , Introducing normalization parameters The interpolation function is constructed using Hermite basis functions;

[0013] When the missing segment in the original migratory bird trajectory data is of medium length, that is , When the target sampling interval is indicated and there are historical similar situation samples within the time window that meet the preset number, neighborhood weighted interpolation is used;

[0014] When the missing segment in the original migratory bird trajectory data is of medium length, that is Furthermore, there is insufficient data on similar migratory bird tracks. At that time, a Monte Carlo candidate trajectory generation strategy was adopted.

[0015] Furthermore, the specific process for calculating the basic spatiotemporal features of each sampling point in the preprocessed and three-level adaptive interpolation of the original migratory bird trajectory data is as follows: the basic spatiotemporal features of each sampling point in the preprocessed and three-level adaptive interpolation of the original migratory bird trajectory data include displacement, instantaneous velocity, acceleration and altitude change rate.

[0016] Calculate displacement: , Indicates time displacement, express Interpolation position at time;

[0017] Calculate instantaneous velocity: , Indicates time Instantaneous velocity;

[0018] Calculate acceleration: , Indicates time The acceleration;

[0019] Calculate the rate of change of height: , Indicates time The rate of change of height, Indicates time The height.

[0020] Furthermore, a comprehensive confidence score is introduced, and a differentiable sigmoid function is used to map features. A point-level migration tendency score is defined and combined with the posterior probability of the Hidden Markov Model (HMM) to obtain a fusion score. Migration points are determined according to a threshold and merged into migration fragments. The specific process for selecting effective migration fragments is as follows:

[0021] Fusion interpolation confidence and localization quality: , Indicates time The overall confidence level; express Empirical normalization results of the time-level precision factor. express Horizontal precision factor at any given time, Indicates to Reference values ​​for empirical normalization; , All represent weighting coefficients. ; This represents the interpolation confidence level at each interpolation point;

[0022] Displacement and velocity are normalized using a differentiable sigmoid function to obtain the displacement. normalized value and instantaneous speed normalized value ;

[0023] definition Moment-level migration tendency score :

[0024] ;

[0025] In the formula, Indicates the weighting factor; Represents the normalization constant; Weighted by positioning accuracy; Speed ​​weights; Weights for directional consistency; For acceleration weights; High weight; Normalized values ​​representing directional consistency; Normalized value representing acceleration; The normalized value representing the height; Indicators representing directional consistency;

[0026] Hidden Markov Models (HMMs) are used to decode the state of the trajectory sequence to obtain the posterior probability of the migration state. , Indicates time state, Indicates migration status. Representative moment In migration mode, This represents the basic spatiotemporal characteristics of the sampling points;

[0027] Calculate the final fusion score:

[0028] ;

[0029] In the formula, Indicates time The final migration integration score; Indicates the fusion coefficient;

[0030] like If the sampling point is positive, then the corresponding sampling point is marked as a migration point; otherwise, it is a non-migration point. The migration points are merged into migration segments and filtered.

[0031] Furthermore, the specific process of step S3 is as follows:

[0032] Step S3.1: Let the set of migration fragments be... , Indicates the first One migration segment; calculate the set of migration segments. Any two migration segments Symmetric average minimum vertical distance ;

[0033] Step S3.2: Based on any two migration segments Symmetric average minimum vertical distance A distance matrix is ​​constructed, and a "single-link clustering" algorithm is used with a set merging threshold. When any two migration segments Symmetric average minimum vertical distance ≤ The corresponding migration fragments are aggregated into the same candidate cluster. Candidate clusters The set of all migration points within the region is , This indicates a segment of migration;

[0034] Step S3.3: For each candidate cluster The system uses a three-level check of "endpoint pairing - time pairing - trajectory overlap" to determine whether it is a complete channel.

[0035] Furthermore, the specific process of step S3.3 is as follows:

[0036] Endpoint pairing check:

[0037] Select candidate clusters Any pair of migration segments , Define migration segments The starting point ,end Calculate the endpoint distance to obtain the migration segment. Origin and Migration Fragments Spatial distance between endpoints and migration fragments The End Point and Migration Fragments Spatial distance between the starting points ;

[0038] like ≤ First preset threshold and If the endpoint pairing is less than or equal to the second preset threshold, then the endpoint pairing is valid.

[0039] Time-paired test:

[0040] Calculate the migration segments separately , Central Time , Set a reasonable upper limit for the return time window. Lower limit ;like Then the time pairing is valid;

[0041] Trajectory overlap check:

[0042] Set buffer radius Calculate migration segments exist Coverage within the buffer and migration fragments exist Coverage within the buffer ;

[0043] based on , calculate and Symmetrical overlap between ;like If the value is greater than or equal to the third preset threshold, then the trajectory overlap test is successful;

[0044] Complete channel determination:

[0045] candidate clusters The number of migration segment pairs that satisfy the three-level test of "endpoint pairing-time pairing-trajectory overlap" is defined as the strong pair number. ;

[0046] If candidate cluster The memory contains at least one pair of segments that satisfy the three-level test of "endpoint pairing - time pairing - trajectory overlap". It is directly determined to be a complete channel;

[0047] like Then count the candidate clusters Bidirectional proportion of inward migration segments: Density clustering is performed on all origin and destination points within the cluster to obtain the bidirectional direction of "origin A - destination B", and the number of segments in the A→B direction is counted. Number of segments in the B→A direction Candidate clusters Total number of segments ;like and and If it is a complete channel, it is considered a complete channel; otherwise, it is considered an incomplete channel and is discarded. The minimum threshold representing the total number of segments. The minimum threshold representing the percentage of directional segments.

[0048] Furthermore, for candidate clusters that are determined to be complete channels... The point set is used, and the centerline is fitted using weighted least squares B-spline with point weights; the residual from the point to the centerline is calculated to determine the channel half-width and construct the boundary. The specific process of forming the fitted channel is as follows:

[0049] Weight setting: Define the candidate clusters that are judged as complete channels. The point set is , express The Middle One point, This represents the number of points in the point set; the weight of each point is represented as: , express At the corresponding time, express The fusion score, express interpolation confidence level;

[0050] B-spline fitting implementation:

[0051] Central axis parameterization: Let the central axis be... , Represents parameterized variables. Indicates the first A cubic B-spline basis function, Indicates and The corresponding control point coefficients, This represents the number of B-spline basis functions;

[0052] Objective function: The goal is to minimize the weighted orthogonal residuals plus smoothing regularization.

[0053] ;

[0054] In the formula, Represent the objective function; Indicates the smoothing regularization parameter; express The third derivative; Indicates the integral;

[0055] The alternating iterative projection-linear update method is used to solve the problem: first, the initial central axis is obtained through weighted linear least squares or principal component analysis. ;No. In each iteration, calculate each point On the current central axis Recent projection parameters on: ,fixed The control point coefficients were then solved using linear weighted least squares. The updated central axis was obtained. Iterate until the objective function converges or the change in the control point coefficients is less than the threshold. ;

[0056] Calculation points Residual to the central axis , express The final projection parameters, express Find the nearest projection point on the central axis; take the residual. of quantiles As a reference for fitting the channel half-width;

[0057] For parameterized variables Segmented calculation of local fractions Extending along the normal direction on both sides of the central axis Generate closed polygonal boundaries to define the spatial range of the fitted channel;

[0058] The fitting channel is defined as a structured entity containing a central axis and boundaries; the boundaries include a half-width datum and a spatial extent.

[0059] Furthermore, adjacent fitted channels that meet the conditions are automatically merged; the empirical distribution of the central axis and boundaries is obtained through Bootstrap resampling, the confidence band and channel confidence score are calculated, and the migration channel containing the central axis, boundaries, and confidence intervals is output. The specific process is as follows:

[0060] Calculate the average mutual projection distance of the centerlines of two adjacent fitted channels and the envelope overlap rate; if the mutual projection distance is less than the channel merging distance threshold and the envelope overlap rate is greater than the fourth preset threshold, then automatically merge the two adjacent fitted channels and refit them to obtain the merged fitted channels.

[0061] Bootstrap resampling is performed on the merged fitted channels to calculate the confidence intervals of the central axis and boundaries, as well as the channel confidence score. Finally, the complete migration channel containing the central axis, boundaries, and confidence intervals is output.

[0062] Compared with existing technologies, the present invention has the following advantages:

[0063] (1) In view of the problems that are common in the original trajectory data, such as irregular sampling intervals, short or long time missing data, and significant differences in positioning accuracy, this invention proposes a three-level adaptive trajectory interpolation strategy, which integrates time decay, behavior similarity and positioning quality weighting mechanism. It can reconstruct biologically reasonable continuous trajectories under different missing lengths and behavior patterns, and provides a quantifiable confidence assessment for each interpolation point, providing a reliable basis for subsequent analysis.

[0064] (2) In view of the problem that traditional behavior state discrimination methods based on simple thresholds or single features are prone to misjudgment and omission when distinguishing between migration and rest, this invention constructs a comprehensive scoring model of migration tendency that integrates instantaneous motion features, time series context and interpolation confidence, so as to realize the automated and high-precision discrimination of migration points and effectively improve the robustness and interpretability of migration segment extraction.

[0065] (3) In view of the shortcomings of existing migration channel identification methods, which rely too much on spatial clustering results and are easily affected by abnormal trajectories or local high-density interference, resulting in pseudo channels or false mergers, this invention proposes channel aggregation and integrity judgment rules based on symmetric spatial distance measurement, endpoint bidirectional pairing test and trajectory overlap analysis to ensure that the identified channels have structural and functional integrity in an ecological sense. Fourth, in view of the problem that the current channel geometric representation lacks statistical robustness and uncertainty quantification support, this invention proposes a channel centerline and boundary band fitting method based on weighted least squares B-spline fitting and bootstrap resampling to realize accurate reconstruction of channel geometry and automatic generation of confidence intervals, output migration channel products with uncertainty evaluation, and support scientific management decision-making. Attached Figure Description

[0066] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0067] like Figure 1 As shown, the present invention provides a technical solution: a method for identifying migratory bird flyways based on spatiotemporal weighted discrimination and least squares fitting, comprising the following steps:

[0068] Step S1: Data preprocessing and three-level adaptive interpolation; The original migratory bird trajectory data is uniformly formatted and projected. For short-term or long-term missing segments in the original migratory bird trajectory data, a three-level adaptive interpolation strategy is used for interpolation reconstruction, and the interpolation confidence is calculated for each interpolation point.

[0069] This involves uniformly formatting and projecting the original migratory bird trajectory data; the original migratory bird trajectory data includes the first... Original location sequence of a single migratory bird , Indicates the first Only one migratory bird Timestamps of each observation point; Indicates the first Only one migratory bird The planar coordinate projection of each observation point , They represent the first Only one migratory bird The horizontal and vertical coordinates of each observation point on the projection plane; Indicates the first Only one migratory bird Height information of each observation point; Indicates the first Only one migratory bird The positioning quality index of each observation point is used to describe the reliability of the observation. Indicates the first The total number of observation points for individual migratory birds.

[0070] This involves uniformly resampling the irregularly sampled trajectories to the target time step. To address short-term or long-term missing segments in the original migratory bird trajectory data, a three-level adaptive interpolation strategy is employed for interpolation reconstruction, and the interpolation confidence level is calculated for each interpolation point. In the three-level adaptive interpolation strategy, the target time... interpolation position The interpolated continuous trajectory is given by a weighted average of neighborhood samples, expressed as:

[0071] ;

[0072] In the formula, Indicates at the target time , No. Standardized weights of neighborhood samples; Indicates the first The location of each neighboring sample; Indicates the target time The set of neighborhood samples; neighborhood samples refer to samples within the target time. These observation points can be used for interpolation. These samples generate continuous trajectories through weighted averaging, providing a basis for subsequent probability boundary analysis and node identification.

[0073] Among them, standardized weights Represented as:

[0074] ;

[0075] In the formula, Indicates at the target time , No. The original weights of each neighboring sample; Indicates at the target time , No. The original weights of the neighboring samples.

[0076] Among them, the original weights Determined by temporal distance decay, behavioral similarity, and positioning quality, it can be expressed as:

[0077] ;

[0078] In the formula, Indicates the first The time corresponding to each observation point; This represents the time decay bandwidth, a constant greater than 0, used to control the degree to which time distance decays the weights. The larger the distance, the weaker the decay effect of time distance on the weight; Indicates the target time The behavioral feature estimate or reference vector is used to compare with the behavioral feature vector of the observation point to measure behavioral similarity; Indicates the observation point The behavioral feature vector contains feature quantities that reflect the behavior of the observed object, such as instantaneous velocity, rotation angle, and acceleration. Indicates the observation point The normalized positioning quality, with values ​​ranging from (0,1], The closer the value is to 1, the higher the positioning quality of the observation point. The bandwidth representing behavioral similarity is a constant greater than 0, used to control the degree to which differences in behavioral features affect the weights. The larger the value, the weaker the impact of behavioral characteristic differences on the weight. This represents a small constant to prevent division by zero.

[0079] The three-level adaptive interpolation strategy includes:

[0080] 1. When the time interval between the two endpoints of a missing segment is present in the original migratory bird trajectory data. When the missing data is short and the motion state is approximately continuous, cubic Hermite spline interpolation with endpoint velocity constraints is used. Specifically:

[0081] Let the start and end endpoints of the missing segment be respectively , The corresponding time is , Calculated by difference between adjacent points , speed , Introducing normalization parameters The interpolation function is constructed using Hermite basis functions:

[0082] ;

[0083] In the formula, The basis function represents the position of the starting endpoint; The velocity basis function represents the initial endpoint; The basis function represents the position of the endpoint; The velocity basis function represents the endpoint.

[0084] Among them, the normalization parameter Represented as:

[0085] .

[0086] This strategy can satisfy the constraint of position and velocity continuity at the endpoints, thus better preserving the short-term dynamic characteristics of migratory birds.

[0087] 2. When the missing segment in the original migratory bird trajectory data is of medium length ( , When the target sampling interval is specified, and there are a sufficient number of historically similar context samples within the time window (historically similar context samples refer to trajectory segments within the time window that are similar to the context (such as time, spatial environment, migration behavior patterns, etc.) of the current missing segment of the migratory bird trajectory), neighborhood-weighted interpolation is used. Specifically:

[0088] In the original weights Based on the formula, behavioral similarity adjustments are introduced:

[0089] ;

[0090] In the formula, Indicates the adjusted original weights ; express and Euclidean distance, Indicates the observation point The local motion feature vector (including velocity, acceleration, and steering angle). Indicates the target time Local motion feature vectors; This represents the similarity bandwidth parameter.

[0091] This strategy assigns higher weights to historical points with more similar behavioral patterns, thereby improving the biological plausibility of the interpolation results.

[0092] 3. When the missing segment in the original migratory bird trajectory data is of medium length ( Furthermore, similar migratory bird track data (highly comparable migratory bird track data observed at other times or on other migratory bird individuals) is insufficient or has extremely long missing segments. When the time of missing original migratory bird trajectory data is much longer than the target sampling interval, a Monte Carlo candidate trajectory generation strategy is adopted, specifically:

[0093] Estimating prior motion parameters, including the probability distribution of velocity, from historical complete migratory bird migration trajectory data. Angular velocity probability distribution Probability distribution of vertical velocity wait, Indicates speed, Indicates angular velocity. This indicates vertical velocity.

[0094] Random generation based on parameter prior Candidate trajectories , Indicates the first 10 candidate trajectories; calculate the confidence score for each candidate trajectory. , is represented as:

[0095] ;

[0096] In the formula, This represents the set of available observation indexes associated with the missing segment; Indicates the first Measured coordinates of one available observation point; Indicates the first Candidate trajectories at time The predicted coordinates; express exist The probability density is given below.

[0097] Specifically, the interpolation trajectory is selected or weighted and synthesized based on the confidence score, and the interpolation confidence score is calculated for each interpolation point. Quantization interpolation uncertainty.

[0098] Step S2, Spatiotemporal Feature Extraction and Migration Point Determination: Calculate the basic spatiotemporal features of each sampling point (position record point in the original migratory bird trajectory data after preprocessing and three-level adaptive interpolation) in the preprocessed and three-level adaptive interpolation migratory bird trajectory data and introduce a comprehensive confidence score. Use a differentiable S-shaped function to map features, define a point-level migration tendency score and combine it with the HMM posterior probability to obtain a fusion score. Determine migration points according to the threshold and merge them into migration segments, and select effective migration segments.

[0099] The specific process for calculating the basic spatiotemporal features of each sampling point in the original migratory bird trajectory data after preprocessing and three-level adaptive interpolation is as follows: Calculate the basic spatiotemporal features of each sampling point in the original migratory bird trajectory data after preprocessing and three-level adaptive interpolation. The basic spatiotemporal features include displacement, instantaneous velocity, acceleration, and rate of change of altitude.

[0100] Calculate displacement: , Indicates time displacement, express Interpolation position at time;

[0101] Calculate instantaneous velocity: , Indicates time Instantaneous velocity;

[0102] Calculate acceleration: , Indicates time The acceleration;

[0103] Calculate the rate of change of height: , Indicates time The rate of change of height, Indicates time The height.

[0104] The process involves introducing a comprehensive confidence level, employing a differentiable sigmoid function to map features, defining a point-level migration tendency score, and combining it with the posterior probability of the Hidden Markov Model (HMM) to obtain a fusion score. Migration points are then determined based on a threshold and merged into migration segments. The specific process for selecting effective migration segments is as follows:

[0105] Fusion interpolation confidence and localization quality: , Indicates time The overall confidence level is used to comprehensively consider the interpolation confidence level and the positioning quality, and is a comprehensive judgment index for the interpolation and positioning quality; express Empirical normalization results of the Horizontal Dilution of Precision (HDOP) at time points. express Horizontal precision factor at any given time, Indicates to Reference values ​​for empirical normalization; , All represent weighting coefficients. .

[0106] Displacement and velocity are normalized using a differentiable sigmoid function (mapping the features to the [0,1] interval):

[0107] ;

[0108] ;

[0109] In the formula, Displacement The normalized value; Indicates instantaneous velocity The normalized value; The reference median value representing the displacement; Scale parameters representing displacement; The reference median value representing speed; A scale parameter representing velocity.

[0110] definition Moment-level migration tendency score :

[0111] ;

[0112] In the formula, Indicates the weighting factor; Represents the normalization constant; Indicators representing directional consistency; The positioning accuracy weight means that the higher the positioning accuracy of a point, the more reliable its contribution to the score. As a velocity weight, it is associated with the instantaneous velocity calculated from the trajectory. The closer the velocity value is to the typical migration velocity range, the greater its positive contribution to the score. As a weight for directional consistency, the importance of the feature of directional concentration is quantified. The more stable the direction, the greater the contribution of this feature to the migration score. As acceleration weights, different behaviors have different motion accelerations, which determine the weight of acceleration features in distinguishing these states; The altitude is weighted and is related to the flight altitude of migratory birds; Normalized values ​​representing directional consistency; Normalized value representing acceleration; The value represents the normalized value of the altitude; each function maps the indicators to a unified scoring scale representing the "migration tendency".

[0113] Hidden Markov Models (HMMs) are used to decode the state of the preprocessed and three-level adaptive interpolation-based original migratory bird trajectory data to obtain the posterior probability of the migration state. , Indicates time state, Indicates migration status. Representative moment In migration mode, This represents the basic spatiotemporal characteristics of the sampling points.

[0114] Calculate the final fusion score:

[0115] ;

[0116] In the formula, Indicates time The final migration integration score; This represents the fusion coefficient.

[0117] like If the sampling point is not identified, it is marked as a migration point; otherwise, it is considered a non-migration point. Migration points are then merged into migration fragments and filtered.

[0118] Merging rule: If the time difference between adjacent migration points is equal to... Merge directly; if the time difference is ≤2 And the discontinuity (the time difference between two adjacent migration points is greater than) But not exceeding 2 Average interpolation confidence level (for the time period) Exceeding the threshold If they are, then they are merged into the same segment.

[0119] Filtering rule: Calculate the duration of the segment and total displacement Remove <2 hours or The final set of migration segments is obtained by capturing segments shorter than 50km (to avoid misjudging short-term hovering and short-distance movement as migration).

[0120] Among them, duration Represented as: , Indicates the number of migration points within a segment;

[0121] Among them, total displacement Represented as: , Indicates the end time of the segment The spatial location of the corresponding sampling point Indicates the start time of the segment The spatial location of the corresponding sampling point.

[0122] Step S3: Initial judgment, spatial aggregation, robustness assessment and allocation of candidate channels: Using the selected effective migration segments as units, calculate the symmetric average minimum vertical distance between effective migration segments and construct a distance matrix. Aggregate them into candidate clusters according to the merging threshold. For candidate clusters The internal segments are checked for endpoint pairing, time pairing and trajectory overlap, and the complete channel is determined based on the number of strong pairs and the proportion of segment directions.

[0123] Step S3.1: Let the set of migration fragments be... Calculate the set of migration fragments Any two migration segments Symmetric average minimum vertical distance , is represented as:

[0124] ;

[0125] In the formula, Indicates the first A segment of the migration flight; , They represent , Migration points; Indicates the migration point to the migration point The minimum vertical distance; Indicates the migration point to the migration point The minimum vertical distance; express Any one of the migration points; express Any migration point in the middle.

[0126] Step S3.2: Based on any two migration segments Symmetric average minimum vertical distance A distance matrix is ​​constructed, and a "single-link clustering" algorithm is used with a set merging threshold. (Generally set to 10–25 km); when any two migration segments Symmetric average minimum vertical distance ≤ The corresponding migration fragments are aggregated into the same candidate cluster. Candidate clusters The set of all migration points within the region is , This indicates a segment of migration.

[0127] Step S3.3: For each candidate cluster The system employs a three-tiered test—endpoint pairing, time pairing, and trajectory overlap—to determine whether a pathway possesses ecological significance and is a complete pathway. Specifically:

[0128] Endpoint pairing check:

[0129] Select candidate clusters Any pair of migration segments , Define migration segments The starting point ,end Calculate the distance between the endpoints:

[0130] ;

[0131] ;

[0132] In the formula, This indicates a migration segment. Origin and Migration Fragments The spatial distance between the endpoints; This indicates a migration segment. The End Point and Migration Fragments The spatial distance between the starting points.

[0133] like ≤5km and If the distance is ≤5km, then the endpoint pairing is initially established.

[0134] Time-paired test:

[0135] Calculate the migration segments separately , Center Time: , , , These represent migration segments. The timestamps of the start and end points , These represent migration segments. Origin and destination timestamps; setting an upper limit for a reasonable return time window. Lower limit ;like If so, then the time pairing is valid.

[0136] Trajectory overlap check:

[0137] Set buffer radius Calculate migration segments exist Coverage within the buffer :

[0138] ;

[0139] In the formula, Indicates the migration point Migration segment The minimum vertical distance; This represents an exponential function.

[0140] calculate and Symmetrical overlap between , This indicates a migration segment. exist Coverage within the buffer; if If the value is ≥0.5, the trajectory overlap test is valid.

[0141] Complete channel determination:

[0142] candidate clusters The number of migration segment pairs that satisfy the three-level test of "endpoint pairing-time pairing-trajectory overlap" is defined as the strong pair number. ;

[0143] If candidate cluster The memory contains at least one pair of segments that satisfy the three-level test of "endpoint pairing - time pairing - trajectory overlap". It is directly determined to be a complete channel;

[0144] like Then count the candidate clusters Bidirectional proportion of inward migration segments: Density clustering (clustering threshold 5km) is performed on all origin and destination points within the cluster to obtain the bidirectional direction of "origin (A) - destination (B)", and the number of segments in the A→B direction is counted. Number of segments in the B→A direction Candidate clusters Total number of segments ;like ( )and ( )and If it is a complete channel, it is considered a complete channel; otherwise, it is considered an incomplete channel and is discarded. The minimum threshold representing the total number of segments. The minimum threshold representing the percentage of directional segments.

[0145] Step S4, channel centerline fitting, merging, boundary decision and uncertainty quantification: For candidate clusters determined to be complete channels The point set is combined with the point weights and a weighted least squares B-spline is used to fit the central axis; the residual from the point to the central axis is calculated to determine the channel half-width and construct the boundary to form the fitted channel; adjacent fitted channels that meet the conditions are automatically merged; the empirical distribution of the central axis and boundary is obtained by Bootstrap resampling, the confidence band and channel confidence score are calculated, and the migration channel containing the central axis, boundary and confidence interval is output.

[0146] Among them, for candidate clusters that are determined to be complete channels The point set is used, and the centerline is fitted using weighted least squares B-spline with point weights; the residual from the point to the centerline is calculated to determine the channel half-width and construct the boundary. The specific process of forming the fitted channel is as follows:

[0147] Weight setting: Define the candidate clusters that are judged as complete channels. The point set is , express The Middle One point, This represents the number of points in the point set; the weight of each point is represented as: , express At the corresponding time, express The fusion score (reflecting the credibility of migration flight data). express The interpolation confidence level (reflecting the reliability of the location); the higher the weight, the greater the influence of the point on the fitting result.

[0148] B-spline fitting implementation:

[0149] 1. Parameterization of the centerline: Let the centerline be... , Represents parameterized variables. Indicates the first A cubic B-spline basis function, Indicates and The corresponding control point coefficients, This represents the number of B-spline basis functions.

[0150] 2. Objective function: The objective is to minimize the weighted orthogonal residuals plus smoothing regularization.

[0151] ;

[0152] In the formula, Represent the objective function; Indicates the smoothing regularization parameter; express The third derivative; This represents the integral.

[0153] 3. Iterative solution using alternating projection-linear update method: First, obtain the initial central axis through weighted linear least squares or principal component analysis. ;No. In each iteration, calculate each point On the current central axis Recent projection parameters on: ,fixed The control point coefficients were then solved using linear weighted least squares. The updated central axis was obtained. Iterate until the objective function converges or the change in the control point coefficients is less than the threshold. .

[0154] 3. Calculation points Residual to the central axis , express The final projection parameters, i.e., for a specific data point After calculation using orthogonal projection, the projection points are obtained. and The specific parameter value that minimizes the distance between them. express Find the nearest projection point on the central axis; take the residual. of quantiles Construct the envelope region as the reference for fitting the channel half-width. , is represented as:

[0155] ;

[0156] In the formula, Represents any point.

[0157] 4. To reflect the change in channel width with position, parameterized variables can be used. Segmented calculation of local fractions Extending along the normal direction on both sides of the central axis This generates a closed polygon boundary, defining the spatial range of the fitted channel.

[0158] 5. Define the fitting channel as a structured entity containing the centerline and boundaries (half-width datum and spatial range).

[0159] The process involves automatically merging adjacent fitted channels that meet the criteria; obtaining the empirical distribution of the central axis and boundaries through Bootstrap resampling; calculating the confidence band and channel confidence score; and outputting the migration channel containing the central axis, boundaries, and confidence intervals.

[0160] 1. Calculate the average mutual projection distance of the centerlines of two adjacent fitted channels. (Project the centerline point of fitted channel 1 onto the centerline of fitted channel 2, calculate the average distance, and vice versa, taking the average of the two) and the envelope overlap rate. (The area of ​​the intersection of the envelopes of the two fitted channels divided by the smaller envelope area); if Distance threshold for channel merging and If the function is not found, then two adjacent fitted channels will be automatically merged and refitted to obtain the merged fitted channel.

[0161] 2. Quantify the uncertainty of the merged fitted channels (Bootstrap resampling), calculate the 95% confidence intervals of the central axis and the boundary, and the channel confidence score. Finally, output the complete migration channel containing the central axis, the boundary, and the confidence interval.

[0162] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for identifying migratory bird flyways using spatiotemporal weighted discrimination and least squares fitting, characterized in that, Includes the following steps: Step S1: Data preprocessing and three-level adaptive interpolation; The original migratory bird trajectory data is uniformly formatted and projected. For short-term or long-term missing segments in the original migratory bird trajectory data, a three-level adaptive interpolation strategy is used for interpolation reconstruction, and the interpolation confidence is calculated for each interpolation point. Step S2, Spatiotemporal Feature Extraction and Migration Point Determination: Calculate the basic spatiotemporal features of each sampling point in the original migratory bird trajectory data after preprocessing and three-level adaptive interpolation, and introduce a comprehensive confidence level. Use a differentiable S-shaped function to map features, define a point-level migration tendency score, and combine it with the HMM posterior probability to obtain a fusion score. Determine migration points according to the threshold and merge them into migration segments, and select effective migration segments. Step S3: Initial Judgment, Spatial Aggregation, Robustness Assessment and Allocation of Candidate Channels: Using the selected effective migration segments as units, calculate the symmetrical average minimum vertical distance between effective migration segments and construct a distance matrix. Aggregate them into candidate clusters according to the merging threshold. Perform endpoint pairing, temporal pairing, and trajectory overlap checks on the segment pairs within the candidate clusters. Determine complete channels based on the number of strong pairings and the segment direction ratio. Then, allocate the candidate clusters... The number of migration segment pairs that satisfy the three-level test of "endpoint pairing-time pairing-trajectory overlap" is defined as the strong pairing number. ; Step S4: Channel centerline fitting, merging, boundary decision and uncertainty quantification: For the point set of candidate clusters determined to be complete channels, the centerline is fitted using weighted least squares B-splines combined with point weights; the residual from the point to the centerline is calculated to determine the channel half-width and construct the boundary, forming the fitted channel; adjacent fitted channels that meet the conditions are automatically merged; the empirical distribution of the centerline and boundary is obtained through Bootstrap resampling, the confidence band and channel confidence score are calculated, and the migration channel containing the centerline, boundary and confidence interval is output; Introducing a comprehensive confidence level, employing a differentiable sigmoid function mapping feature, defining a point-level migration tendency score, and combining it with the posterior probability of the Hidden Markov Model (HMM) to obtain a fusion score, determining migration points based on a threshold and merging them into migration fragments, the specific process for selecting effective migration fragments is as follows: Fusion interpolation confidence and localization quality: , Indicates time The overall confidence level; express Empirical normalization results of the time-level precision factor. express Horizontal precision factor at any given time, Indicates to Reference values ​​for empirical normalization; , All represent weighting coefficients. ; This represents the interpolation confidence level at each interpolation point; Displacement and velocity are normalized using a differentiable sigmoid function to obtain the displacement. normalized value and instantaneous speed normalized value ; definition Moment-level migration tendency score : ； In the formula, Indicates the weighting factor; Represents the normalization constant; Weighted by positioning accuracy; Speed ​​weights; Weights for directional consistency; For acceleration weights; High weight; Normalized values ​​representing directional consistency; Normalized value representing acceleration; The normalized value representing the height; Indicators representing directional consistency; Hidden Markov Models (HMMs) are used to decode the state of the trajectory sequence to obtain the posterior probability of the migration state. , Indicates time state, Indicates migration status. Representative moment In migration mode, This represents the basic spatiotemporal characteristics of the sampling points; Calculate the final fusion score: ； In the formula, Indicates time The final migration integration score; Indicates the fusion coefficient; like If the sampling point is positive, then the corresponding sampling point is marked as a migration point; otherwise, it is a non-migration point. The migration points are merged into migration segments and filtered.

2. The method for identifying migratory bird flyways using spatiotemporal weighted discrimination and least squares fitting according to claim 1, characterized in that: In step S1, the original migratory bird trajectory data is uniformly formatted and projected; the original migratory bird trajectory data includes the first... Original location sequence of a single migratory bird , Indicates the first Only one migratory bird Timestamps of each observation point; Indicates the first Only one migratory bird The planar coordinate projection of each observation point , They represent the first Only one migratory bird The horizontal and vertical coordinates of each observation point on the projection plane; Indicates the first Only one migratory bird Height information of each observation point; Indicates the first Only one migratory bird Positioning quality indicators for each observation point Indicates the first The total number of observation points for individual migratory birds.

3. The method for identifying migratory bird flyways using spatiotemporal weighted discrimination and least squares fitting according to claim 2, characterized in that: In step S1, the irregularly sampled trajectories are uniformly resampled to the target time step. To address short-term or long-term missing segments in the original migratory bird trajectory data, a three-level adaptive interpolation strategy is employed for interpolation reconstruction, and the interpolation confidence level is calculated for each interpolation point. In the three-level adaptive interpolation strategy, the target time... interpolation position The interpolated continuous trajectory is given by a weighted average of neighborhood samples; neighborhood samples refer to samples taken at the target time. , used as the observation point for interpolation.

4. The method for identifying migratory bird flyways using spatiotemporal weighted discrimination and least squares fitting according to claim 3, characterized in that: Three-level adaptive interpolation strategies include: When the time interval between the two endpoints of a missing segment is in the original migratory bird trajectory data At that time, cubic Hermite spline interpolation with endpoint velocity constraints is used: let the start endpoint and end endpoint of the missing segment be respectively... , The corresponding time is , Calculated by difference between adjacent points , speed , Introducing normalization parameters The interpolation function is constructed using Hermite basis functions; When the missing segment in the original migratory bird trajectory data is of medium length, that is , When the target sampling interval is indicated and there are historical similar situation samples within the time window that meet the preset number, neighborhood weighted interpolation is used; When the missing segment in the original migratory bird trajectory data is of medium length, that is Furthermore, there is insufficient data on similar migratory bird tracks. At that time, a Monte Carlo candidate trajectory generation strategy was adopted.

5. The method for identifying migratory bird flyways using spatiotemporal weighted discrimination and least squares fitting according to claim 4, characterized in that: The specific process for calculating the basic spatiotemporal features of each sampling point in the preprocessed and three-level adaptive interpolation data of migratory birds is as follows: Calculate the basic spatiotemporal features of each sampling point in the preprocessed and three-level adaptive interpolation data of migratory birds. The basic spatiotemporal features include displacement, instantaneous velocity, acceleration, and rate of change of altitude. Calculate displacement: , Indicates time displacement, express Interpolation position at time; Calculate instantaneous velocity: , Indicates time Instantaneous velocity; Calculate acceleration: , Indicates time The acceleration; Calculate the rate of change of height: , Indicates time The rate of change of height, Indicates time The height.

6. The method for identifying migratory bird flyways using spatiotemporal weighted discrimination and least squares fitting according to claim 5, characterized in that: The specific process of step S3 is as follows: Step S3.1: Let the set of migration fragments be... , Indicates the first One migration segment; calculate the set of migration segments. Any two migration segments Symmetric average minimum vertical distance ; Step S3.2: Based on any two migration segments Symmetric average minimum vertical distance Construct a distance matrix, employ the "single-link clustering" algorithm, and set a merging threshold. When any two migration segments Symmetric average minimum vertical distance ≤ The corresponding migration fragments are aggregated into the same candidate cluster. Candidate clusters The set of all migration points within the region is , This indicates a segment of migration; Step S3.3: For each candidate cluster The system uses a three-level check of "endpoint pairing - time pairing - trajectory overlap" to determine whether it is a complete channel.

7. The method for identifying migratory bird flyways using spatiotemporal weighted discrimination and least squares fitting according to claim 6, characterized in that: The specific process of step S3.3 is as follows: Endpoint pairing check: Select candidate clusters Any pair of migration segments , Define migration segments The starting point ,end Calculate the endpoint distance to obtain the migration segment. Origin and Migration Fragments Spatial distance between endpoints and migration fragments The End Point and Migration Fragments Spatial distance between the starting points ; like ≤ First preset threshold and If the endpoint pairing is less than or equal to the second preset threshold, then the endpoint pairing is valid. Time-paired test: Calculate the migration segments separately , Central Time , ; Set a reasonable upper limit for the return time window. Lower limit ;like Then the time pairing is valid; Trajectory overlap check: Set buffer radius Calculate migration segments exist Coverage within the buffer and migration fragments exist Coverage within the buffer ; based on , calculate and Symmetrical overlap between ;like If the value is greater than or equal to the third preset threshold, then the trajectory overlap test is successful; Complete channel determination: If candidate cluster The memory contains at least one pair of fragments that satisfy the three-level check of "endpoint pairing - time pairing - trajectory overlap". It is directly determined to be a complete channel; like Then count the candidate clusters Bidirectional proportion of inward migration segments: Density clustering is performed on all origin and destination points within the cluster to obtain the bidirectional direction of "origin A - destination B", and the number of segments in the A→B direction is counted. Number of segments in the B→A direction Candidate clusters Total number of segments ;like and and If so, it is determined to be a complete channel; Otherwise, it will be considered an incomplete channel and will be removed. The minimum threshold representing the total number of segments. The minimum threshold representing the percentage of directional segments.

8. The method for identifying migratory bird flyways using spatiotemporal weighted discrimination and least squares fitting according to claim 7, characterized in that: For candidate clusters that are determined to be complete channels The point set is used, and the centerline is fitted using weighted least squares B-spline with point weights; the residual from the point to the centerline is calculated to determine the channel half-width and construct the boundary. The specific process of forming the fitted channel is as follows: Weight setting: Define the candidate clusters that are judged as complete channels. The point set is , express The Middle One point, This represents the number of points in the point set; the weight of each point is represented as: , express At the corresponding time, express The fusion score, express interpolation confidence level; B-spline fitting implementation: Central axis parameterization: Let the central axis be... , Represents parameterized variables. Indicates the first A cubic B-spline basis function, Indicates and The corresponding control point coefficients, This represents the number of B-spline basis functions; Objective function: The goal is to minimize the weighted orthogonal residuals plus smoothing regularization. ； In the formula, Represent the objective function; Indicates the smoothing regularization parameter; express The third derivative; The alternating iterative projection-linear update method is used to solve the problem: first, the initial central axis is obtained through weighted linear least squares or principal component analysis. ;No. In each iteration, calculate each point On the current central axis Recent projection parameters on: ,fixed The control point coefficients were then solved using linear weighted least squares. The updated central axis was obtained. Iterate until the objective function converges or the change in the control point coefficients is less than the threshold. ; Calculation points Residual to the central axis , express The final projection parameters, express Find the nearest projection point on the central axis; take the residual. of quantiles As a reference for fitting the channel half-width; For parameterized variables Segmented calculation of local fractions Extending along the normal direction on both sides of the central axis Generate closed polygonal boundaries to define the spatial range of the fitted channel; The fitting channel is defined as a structured entity containing a central axis and boundaries; the boundaries include a half-width datum and a spatial extent.

9. The method for identifying migratory bird flyways using spatiotemporal weighted discrimination and least squares fitting according to claim 8, characterized in that: The process of automatically merging adjacent fitted channels that meet the conditions, obtaining the empirical distribution of the central axis and boundaries through Bootstrap resampling, calculating the confidence band and channel confidence score, and outputting the migration channel containing the central axis, boundaries, and confidence interval is as follows: Calculate the average mutual projection distance of the centerlines of two adjacent fitted channels and the envelope overlap rate; if the mutual projection distance is less than the channel merging distance threshold and the envelope overlap rate is greater than the fourth preset threshold, then automatically merge the two adjacent fitted channels and refit them to obtain the merged fitted channels. Bootstrap resampling is performed on the merged fitted channels to calculate the confidence intervals of the central axis and boundaries, as well as the channel confidence score. Finally, the complete migration channel containing the central axis, boundaries, and confidence intervals is output.

Images