Lake region ship trajectory anomaly identification and early warning method and system based on data fusion

Abstract

The invention provides a lake region ship trajectory anomaly identification and early warning method and system based on data fusion, and relates to the technical field of ship safety, and the method comprises the steps: obtaining multi-source trajectory data, carrying out the space-time alignment processing, dividing the fusion data into time windows, and constructing space-time data blocks; identifying candidate abnormal time periods in the high-dimensional feature space by using space mapping; phase-space reconstruction and causal analysis are carried out on the candidate abnormal time periods, and real anomalies are confirmed; and finally generating graded early warning information according to the types and degrees of the anomalies. According to the invention, accurate identification and timely early warning of ship trajectory abnormity are realized, and the safety level of ships in lake areas is improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of ship safety, and in particular to a lake area ship trajectory anomaly identification and early warning method and system based on data fusion. BACKGROUND

[0002] With the rapid development of inland waterway transportation, the number of lake area ship passages is increasing, and ship safety supervision is facing great challenges. Ship trajectory anomaly identification and early warning technology, as a key component of intelligent water transportation supervision, plays an important role in preventing navigation risks and reducing water accidents. Currently, ship trajectory monitoring mainly relies on automatic identification systems, radars, and video monitoring to obtain ship operation state information. With the development of Internet of Things and big data technology, the fusion processing of multi-source trajectory data provides more comprehensive data support for ship behavior analysis.

[0003] Traditional ship abnormal behavior identification methods mainly judge whether a ship deviates from a normal channel or has abnormal speed through preset rules. With the development of machine learning and deep learning technologies, existing technologies automatically learn normal behavior patterns by analyzing historical navigation data, and then identify abnormal behaviors that deviate from these patterns.

[0004] However, the existing technology still has the problems of difficult data integration, inability to fully utilize the advantages of multi-source data for comprehensive anomaly analysis, neglect of the influence of environmental factors on ship behavior, inability to effectively distinguish between normal behavior adjustment caused by environmental condition changes and real abnormal behavior, easy generation of a large number of false positives, lack of comprehensive evaluation of the persistence and severity of abnormal behavior, inability to provide graded early warning according to different types of abnormal behavior, and difficulty in supporting the targeted emergency response decisions of regulatory departments. SUMMARY

[0005] The embodiments of the present application provide a lake area ship trajectory anomaly identification and early warning method and system based on data fusion, which can at least solve some of the problems in the prior art.

[0006] In a first aspect, the embodiments of the present application provide a lake area ship trajectory anomaly identification and early warning method based on data fusion, comprising:

[0007] Obtain multi-source trajectory data and perform spatio-temporal alignment processing to obtain fusion trajectory data by combining pre-set timestamp association and spatial coordinate mapping rules. Divide the fusion trajectory data into multiple time windows and construct a spatio-temporal data block corresponding to each time window.

[0008] mapping the position sequence in the spatio-temporal data block to a high-dimensional feature space through embedding space mapping, comparing the distribution positions of the current spatio-temporal data block and the spatio-temporal data block corresponding to the historical normal navigation in the topological space and identifying a candidate abnormal time period;

[0009] extracting the motion state sequence and the environment information sequence corresponding to the candidate abnormal time period and performing phase space reconstruction in the time dimension, calculating the mutual information transfer strength and transfer direction between the motion state sequence and the environment information sequence in the reconstructed phase space and constructing a causal influence matrix, performing spectral decomposition on the causal influence matrix, obtaining a dominant eigenvalue and judging the dominant action degree of the environment change on the motion state change, judging the candidate abnormal time period based on the dominant action degree to obtain a confirmed abnormal event set;

[0010] iterating the confirmed abnormal event set in chronological order and identifying sustained abnormalities and sudden abnormalities, generating hierarchical warning information according to the duration of the sustained abnormalities and the deviation degree of the sudden abnormalities and outputting.

[0011] In an optional implementation,

[0012] Obtaining multi-source trajectory data and performing spatio-temporal alignment processing to obtain fused trajectory data in combination with pre-set timestamp association and spatial coordinate mapping rules, dividing the fused trajectory data into multiple time windows and constructing a spatio-temporal data block corresponding to each time window includes:

[0013] Obtaining the position sequence and the motion state sequence corresponding to the ship and the environment information sequence corresponding to the lake area and combining to obtain multi-source trajectory data, extracting the timestamp corresponding to the multi-source trajectory data as an association key value to establish a synchronization index table, grouping and associating data records with timestamp difference within a tolerance range through the synchronization index table to obtain an associated data group;

[0014] Performing coordinate system conversion on the position sequence in the associated data group to unify to a standard coordinate system and generating an aligned position sequence in combination with interpolation correction, performing resampling on the motion state sequence in the associated data group and unifying to a standard sampling frequency to generate an aligned motion sequence, performing resampling on the environment information sequence in the associated data group and unifying to a standard sampling frequency to generate an aligned environment sequence, and merging the aligned position sequence, the aligned motion sequence and the aligned environment sequence according to the timestamp to generate fused trajectory data;

[0015] The fusion trajectory data is divided into a preliminary window sequence according to a fixed step length, a fusion trajectory data segment corresponding to each preliminary window in the preliminary window sequence is extracted, a complexity index of an alignment position sequence in the fusion trajectory data segment is calculated, and a time span of the preliminary window is adjusted to generate an adaptive window sequence according to the complexity index;

[0016] A fusion trajectory data segment corresponding to each adaptive window in the adaptive window sequence is extracted, and a plurality of spatiotemporal data blocks are obtained through difference calculation.

[0017] In an optional implementation,

[0018] The position sequence in the spatiotemporal data block is mapped to a high-dimensional feature space through an embedding space mapping, and the distribution positions of the current spatiotemporal data block and the spatiotemporal data block corresponding to the historical normal navigation in the topological space are compared in the high-dimensional feature space to identify a candidate abnormal time period, including:

[0019] For the position sequence in the spatiotemporal data block, a neighborhood matrix is constructed and a reconstruction weight is calculated through a local linear embedding algorithm, the position sequence is projected to a high-dimensional feature space to obtain a current embedding vector based on the reconstruction weight, and a metric tensor is constructed through a Riemannian geometry algorithm based on the distribution of the current embedding vector in the high-dimensional feature space.

[0020] The spatiotemporal data block corresponding to the historical normal navigation is obtained, and the position sequence in the historical spatiotemporal data block is projected to a high-dimensional feature space to obtain a historical embedding vector set, the topological distance between each historical embedding vector in the historical embedding vector set and the nearest neighbor historical embedding vector is calculated through an adaptive density clustering algorithm based on the metric tensor, and a local density radius is determined in combination with a preset quantile position, and a hypersphere surface point set is constructed based on the historical embedding vector as the center and in combination with the local density radius to obtain a normal distribution boundary point set.

[0021] The topological distance between the current embedding vector and each boundary point in the normal distribution boundary point set is calculated based on the metric tensor, and the minimum value is extracted as a deviation metric value, the abnormal score is solved based on the statistical characteristics of the topological distance between the deviation metric value and the historical embedding vector set, and the corresponding time window is determined as a candidate abnormal time period according to the abnormal score.

[0022] In an optional implementation,

[0023] The topological distance between each historical embedding vector in the historical embedding vector set and the nearest neighbor historical embedding vector is calculated through an adaptive density clustering algorithm based on the metric tensor, and a local density radius is determined in combination with a preset quantile position, including:

[0024] Based on the metric tensor, geodesic path integration between different historical embedding vectors in the set of historical embedding vectors is calculated to obtain a topological distance, for each historical embedding vector, the historical embedding vector corresponding to the minimum value of the topological distance is extracted as the nearest neighbor historical embedding vector, and the topological distance from each historical embedding vector to the nearest neighbor historical embedding vector is recorded to obtain a set of nearest neighbor distances;

[0025] The topological distances in the set of nearest neighbor distances are arranged in ascending order to obtain an ordered distance sequence, and based on the preset quantile position, the corresponding index position in the ordered distance sequence is located and the topological distance at the index position is extracted as the global quantile distance;

[0026] The number of historical embedding vectors contained in the hypersphere with the global quantile distance as the radius in the set of historical embedding vectors is calculated to obtain a neighborhood density value, and the nearest neighbor topological distance corresponding to each historical embedding vector is weighted and fused with the reciprocal of the neighborhood density value to obtain a local density radius.

[0027] In an optional implementation,

[0028] The motion state sequence and the environmental information sequence corresponding to the candidate abnormal period are extracted and reconstructed in the phase space in the time dimension, and the mutual information transfer strength and transfer direction between the motion state sequence and the environmental information sequence in the reconstructed phase space are calculated, and a causal influence matrix is constructed, including:

[0029] Based on the start time point and the end time point of the candidate abnormal period, the motion state sequence and the environmental information sequence within the range from the start time point to the end time point are intercepted from the spatio-temporal data block;

[0030] For each sub-sequence in the motion state sequence and the environmental information sequence, the optimal delay time is calculated by the time delay embedding algorithm and the pre-set mutual information function, the optimal embedding dimension is calculated based on the false nearest neighbor point algorithm, and the phase space is reconstructed based on the optimal delay time and the optimal embedding dimension to obtain a set of motion state phase space vectors and a set of environmental information phase space vectors;

[0031] The forward information transfer amount of the environmental information variable in the set of environmental information phase space vectors to the motion state variable in the set of motion state phase space vectors and the reverse information transfer amount of the motion state variable in the set of motion state phase space vectors to the environmental information variable in the set of environmental information phase space vectors are calculated by the transfer entropy algorithm, and a bidirectional transfer entropy matrix is constructed based on the forward information transfer amount and the reverse information transfer amount;

[0032] The significant test algorithm is used to statistically test the transfer entropy values in the bidirectional transfer entropy matrix and determine effective transfer strength, and the starting variable and target variable of information transfer are marked according to the directionality of the effective transfer strength, so as to construct a causal influence matrix with motion state variables and environmental information variables as nodes and effective transfer strength and transfer direction as directed edges.

[0033] In an optional implementation,

[0034] The causal influence matrix is subjected to spectral decomposition to obtain a dominant eigenvalue and judge the dominant action degree of environmental change on motion state change, and the candidate abnormal time period is judged based on the dominant action degree to obtain a confirmed abnormal event set, including:

[0035] The causal influence matrix is subjected to spectral decomposition to extract all eigenvalues and corresponding eigenvectors, the eigenvalues are arranged in descending order according to the modulus value, and the first eigenvalue is taken as the dominant eigenvalue, and the eigenvector corresponding to the dominant eigenvalue is extracted as the dominant eigenvector;

[0036] The absolute value of each element in the dominant eigenvector is calculated, and the element with the maximum absolute value is determined, the index position of the element with the maximum absolute value in the dominant eigenvector is obtained and mapped to the corresponding variable node in the causal influence matrix, the variable type corresponding to the variable node is determined, and the dominant direction of causal influence is determined according to the variable type;

[0037] Based on the dominant eigenvalue, the dominant action strength is calculated by a spectral radius measurement method, and the dominant action strength is compared with a preset dominant action strength threshold, and the dominant action degree of environmental change on motion state change is determined according to the comparison result;

[0038] The type of the candidate abnormal time period is determined based on the dominant action degree and the dominant direction, if the dominant action degree exceeds a preset action degree threshold and the dominant direction is that of environmental information variables, the candidate abnormal time period is classified as an environment-driven abnormal event, otherwise, the candidate abnormal time period is classified as an autonomous behavior type abnormal event, and the confirmed abnormal event set is obtained by repeating the determination.

[0039] In an optional implementation,

[0040] The confirmed abnormal event set is traversed in time sequence and sustained abnormality and burst abnormality are identified, hierarchical warning information is generated according to the sustained time length of the sustained abnormality and the deviation degree of the burst abnormality, and output includes:

[0041] arranging all the abnormal events in the confirmed abnormal event set in ascending order according to timestamps and traversing, determining a time interval between a start time and an end time corresponding to each abnormal event as a duration, and extracting a numerical change feature of the motion state sequence in the abnormal event period;

[0042] calculating a time sequence change rate of the motion state sequence in the abnormal event period by a sliding window method, marking an abnormal event as a sustained abnormality if the time sequence change rate is lower than a preset time sequence change rate threshold and the duration exceeds a preset duration threshold, otherwise, marking the abnormal event as a burst abnormality;

[0043] calculating an evolution rate corresponding to the sustained abnormality by a trend fitting method and determining a duration level in combination with the duration, extracting a peak deviation and a change gradient corresponding to the burst abnormality and determining a deviation level, generating a warning level based on the duration level and the deviation level, and combining the warning level with a pre-determined abnormality type and a pre-acquired dominant direction to obtain graded warning information and output.

[0044] In a second aspect of the embodiment of the present application, a lake area ship trajectory abnormality recognition and warning system based on data fusion is provided, comprising:

[0045] a data fusion module configured to acquire multi-source trajectory data and perform spatio-temporal alignment processing to obtain fusion trajectory data in combination with pre-set timestamp association and spatial coordinate mapping rules, divide the fusion trajectory data into a plurality of time windows and construct a spatio-temporal data block corresponding to each time window;

[0046] an abnormality recognition module configured to map a position sequence in the spatio-temporal data block to a high-dimensional feature space by embedding spatial mapping, compare a distribution position of a current spatio-temporal data block and a historical normal navigation corresponding spatio-temporal data block in a topological space in the high-dimensional feature space, and recognize a candidate abnormal period;

[0047] a causal analysis module configured to extract a motion state sequence and an environment information sequence corresponding to the candidate abnormal period, and perform phase space reconstruction in a time dimension, calculate mutual information transfer strength and transfer direction between the motion state sequence and the environment information sequence in the reconstructed phase space, and construct a causal influence matrix, perform spectral decomposition on the causal influence matrix, obtain a dominant eigenvalue, and judge a dominant action degree of environment change on motion state change, judge the candidate abnormal period based on the dominant action degree to obtain a confirmed abnormal event set;

[0048] an abnormality warning module configured to traverse the confirmed abnormal event set in chronological order and recognize sustained abnormality and burst abnormality, generate graded warning information according to a duration of the sustained abnormality and a deviation degree of the burst abnormality, and output.

[0049] In a third aspect, the present application provides an electronic device, comprising:

[0050] a processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke the instructions stored in the memory to perform the method described above.

[0051] In a fourth aspect, the present application provides a computer-readable storage medium having stored thereon computer program instructions, which, when executed by a processor, implement the method described above.

[0052] In the present application, through multi-source trajectory data fusion and spatio-temporal alignment processing, comprehensive collection and integration of lake ship movement information are realized, the integrity and accuracy of the trajectory data are improved, and a high-quality data foundation is provided for subsequent anomaly identification. The spatial mapping embedding technology is used to map the position sequence to a high-dimensional feature space, which can more effectively capture the distribution characteristics of the ship trajectory in the topological space, improve the sensitivity and accuracy of candidate abnormal period identification, introduce phase space reconstruction and causal influence matrix analysis method, realize the quantitative characterization of the complex correlation between the movement state sequence and the environment information sequence, effectively distinguish the normal deviation caused by environmental factors and the real abnormal behavior through spectral decomposition and dominant eigenvalue judgment, significantly reduce the false positive rate, realize the fine classification of abnormal events through the analysis of the persistence and suddenness of the confirmed abnormal events, and generate hierarchical warning information based on the persistence and deviation degree, improve the pertinence and timeliness of the warning, and provide scientific and effective decision support for lake ship safety management. BRIEF DESCRIPTION OF DRAWINGS

[0053] Figure 1 FIG. 1 is a flowchart of a lake ship trajectory anomaly identification and warning method based on data fusion according to an embodiment of the present application;

[0054] Figure 2 FIG. 4 is an abnormal event analysis flowchart of the lake ship trajectory anomaly identification and warning method based on data fusion according to an embodiment of the present application. DETAILED DESCRIPTION

[0055] In order to make the purpose, technical scheme and advantages of the embodiments of the present application clearer, the technical scheme of the embodiments of the present application will be described clearly and completely below with reference to the drawings of the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0056] The technical solutions of the present application are described in detail below with specific examples. The following specific examples can be combined with each other, and the same or similar concepts or processes may not be described in detail in some examples.

[0057] Figure 1 The flowchart of the lake area ship trajectory anomaly identification and early warning method based on data fusion of the embodiments of the present application is shown in FIG. 1, which comprises the following steps: Figure 1

[0058] Obtaining multi-source trajectory data and combining pre-set timestamp association and spatial coordinate mapping rules to perform spatio-temporal alignment processing to obtain fusion trajectory data, dividing the fusion trajectory data into multiple time windows and constructing a spatio-temporal data block corresponding to each time window;

[0059] Mapping the position sequence in the spatio-temporal data block to a high-dimensional feature space through embedded space mapping, comparing the distribution positions of the current spatio-temporal data block and the spatio-temporal data block corresponding to the historical normal navigation in the topological space in the high-dimensional feature space, and identifying a candidate abnormal period;

[0060] Extracting the motion state sequence and the environment information sequence corresponding to the candidate abnormal period and performing phase space reconstruction in the time dimension, calculating the mutual information transfer strength and transfer direction between the motion state sequence and the environment information sequence in the reconstructed phase space, and constructing a causal influence matrix, performing spectral decomposition on the causal influence matrix, obtaining a dominant eigenvalue, and judging the dominant action degree of the environment change on the motion state change, judging the candidate abnormal period based on the dominant action degree to obtain a confirmed abnormal event set;

[0061] Iterating the confirmed abnormal event set in chronological order and identifying sustained abnormality and sudden abnormality, generating hierarchical warning information according to the duration of the sustained abnormality and the deviation degree of the sudden abnormality, and outputting.

[0062] In an optional implementation,

[0063] Obtaining multi-source trajectory data and combining pre-set timestamp association and spatial coordinate mapping rules to perform spatio-temporal alignment processing to obtain fusion trajectory data, dividing the fusion trajectory data into multiple time windows and constructing a spatio-temporal data block corresponding to each time window comprises:

[0064] Obtaining the position sequence and the motion state sequence corresponding to the ship and the environment information sequence corresponding to the lake area, and combining to obtain multi-source trajectory data, extracting the timestamp corresponding to the multi-source trajectory data as an association key value to establish a synchronization index table, and grouping and associating data records with timestamp difference within a tolerance range through the synchronization index table to obtain an associated data group;

[0065] ​performing coordinate system conversion on the position sequence in the associated data set to a standard coordinate system, generating an aligned position sequence by combining interpolation correction, performing resampling on the motion state sequence in the associated data set to a standard sampling frequency, generating an aligned motion sequence, performing resampling on the environment information sequence in the associated data set to a standard sampling frequency, generating an aligned environment sequence, and merging the aligned position sequence, the aligned motion sequence, and the aligned environment sequence according to timestamps to generate fused trajectory data;

[0066] dividing the fused trajectory data into a preliminary window sequence according to a fixed step length, extracting a fused trajectory data segment corresponding to each preliminary window in the preliminary window sequence, calculating a complexity index of the aligned position sequence in the fused trajectory data segment, and adjusting a time span of the preliminary window according to the complexity index to generate an adaptive window sequence;

[0067] extracting a fused trajectory data segment corresponding to each adaptive window in the adaptive window sequence and performing difference calculation to obtain a plurality of spatiotemporal data blocks.

[0068] A position sequence and a motion state sequence corresponding to a ship and an environment information sequence corresponding to a lake area are obtained. The position sequence includes longitude, latitude, and a collection timestamp. The motion state sequence includes ship speed, heading angle, sailing speed, and a collection timestamp. The environment information sequence includes wind speed, wind direction, lake water level, visibility, and a collection timestamp. The obtained sequences are combined to obtain multi-source trajectory data, and a timestamp corresponding to the multi-source trajectory data is extracted as an associated key value to establish a synchronization index table. For example, in a lake area ship monitoring system, the position sequence data of a certain ship is recorded as "north latitude 30.526 degrees, east longitude 114.357 degrees, 2025-01-05 08:15:30", the motion state is "heading angle 75 degrees, sailing speed 8.5 knots, 2025-01-05 08:15:32", and the environment information is "wind speed 5.2 meters / second, northeast wind, water level 125 meters, visibility 800 meters, 2025-01-05 08:15:27".

[0069] The data records with a timestamp difference within a tolerance range are grouped and associated by the synchronization index table to obtain an associated data set. The tolerance range is set to 5 seconds, which means that data with a timestamp difference of not more than 5 seconds will be grouped into the same group. In the foregoing example, the timestamps of the position, motion state, and environment information are 08:15:30, 08:15:32, and 08:15:27, respectively, and the difference is not more than 5 seconds, so they are grouped into the same associated data set.

[0070] The position sequence in the associated data set is converted into a coordinate system, and the latitude and longitude coordinates are uniformly converted into UTM projection coordinates. If there are position points with incomplete timestamp alignment, linear interpolation method is used for correction to generate aligned position sequence. For example, the latitude and longitude coordinates "30.526 degrees north latitude, 114.357 degrees east longitude" are converted into UTM coordinates "east coordinate 548256.75 meters, north coordinate 3375890.22 meters". The motion state sequence is resampled to a standard sampling frequency of 1 Hz to generate an aligned motion sequence. The environment information sequence is also similarly resampled to a standard sampling frequency of 1 Hz to generate an aligned environment sequence. The aligned position sequence, the aligned motion sequence and the aligned environment sequence are merged according to the timestamp to generate the fusion trajectory data.

[0071] The fusion trajectory data is divided into a preliminary window sequence according to a fixed step length of 60 seconds. The fusion trajectory data segment corresponding to each preliminary window in the preliminary window sequence is extracted. The complexity indicators of the aligned position sequence in the fusion trajectory data segment are calculated, including the trajectory curvature and the heading rate of change. For curvature calculation, three consecutive position points P1, P2 and P3 are selected, and the reciprocal of the radius of the circle formed by the three position points is calculated as the curvature value; for the heading rate of change, the change amount of the heading angle of the adjacent two time points is divided by the time interval. For example, in a certain 60-second window, the average curvature is 0.0025, the maximum curvature is 0.0075, the average heading rate of change is 1.8 degrees per second, and the maximum heading rate of change is 5.2 degrees per second.

[0072] The time span of the preliminary window is adjusted according to the complexity indicators to generate an adaptive window sequence. The specific adjustment rules are as follows: when the average curvature is greater than 0.005 or the average heading rate of change is greater than 3 degrees per second, the window time span is reduced to 30 seconds; when the average curvature is less than 0.001 and the average heading rate of change is less than 0.5 degrees per second, the window time span is expanded to 120 seconds; otherwise, the original 60-second window remains unchanged. In the above example, since the average curvature 0.0025 is lower than 0.005 and the average heading rate of change 1.8 degrees per second is lower than 3 degrees per second, and at the same time the average curvature is higher than 0.001 or the average heading rate of change is higher than 0.5 degrees per second, the original 60-second window remains unchanged.

[0073] The fusion track data segment corresponding to each adaptive window in the adaptive window sequence is extracted, and a plurality of spatio-temporal data blocks are obtained by difference calculation. The difference calculation includes position difference, speed difference, heading angle difference and environmental parameter difference. For the position difference, the difference between the position coordinates of two adjacent time points is calculated to obtain a displacement vector; for the speed difference, the difference between the speeds of two adjacent time points is calculated; for the heading angle difference, the difference between the heading angles of two adjacent time points is calculated, and the case where the angle crosses 360 degrees is handled; for the environmental parameter difference, the difference between the environmental parameters of two adjacent time points is calculated. For example, in a certain adaptive window, the UTM coordinates of two adjacent time points are "548256.75 meters east, 3375890.22 meters north" and "548258.95 meters east, 3375893.42 meters north", and the displacement vector is calculated as "2.2 meters east, 3.2 meters north"; the speed changes from 8.5 knots to 8.7 knots, and the difference is 0.2 knots; the heading angle changes from 75 degrees to 78 degrees, and the difference is 3 degrees; the wind speed changes from 5.2 m / s to 5.5 m / s, and the difference is 0.3 m / s.

[0074] Each spatio-temporal data block contains the position coordinates, displacement vector, speed, heading angle, environmental parameter and their difference values of all time points in the window, forming a complete feature description.

[0075] In this embodiment, by time stamp association and synchronous indexing construction of the ship position sequence, the motion state sequence and the lake area environmental information sequence, accurate alignment of multi-source heterogeneous track data in the time dimension is realized, the information mismatch problem caused by inconsistent sampling frequency or time drift is avoided, and the integrity and consistency of track data fusion are significantly improved.

[0076] In an optional implementation,

[0077] The position sequence in the spatio-temporal data block is mapped to a high-dimensional feature space by embedding space mapping, the distribution positions of the current spatio-temporal data block and the spatio-temporal data block corresponding to the historical normal navigation in the topological space are compared in the high-dimensional feature space, and a candidate abnormal period is identified, including:

[0078] For the position sequence in the spatio-temporal data block, a neighborhood matrix is constructed by a local linear embedding algorithm and a reconstruction weight is calculated, the position sequence is projected to a high-dimensional feature space to obtain a current embedding vector based on the reconstruction weight, and a metric tensor is constructed by a Riemannian geometry algorithm based on the distribution of the current embedding vector in the high-dimensional feature space.

[0079] The spatiotemporal data blocks corresponding to historical normal navigation are obtained, and the position sequences in the historical spatiotemporal data blocks are projected onto a high-dimensional feature space to obtain a set of historical embedding vectors. Based on the metric tensor, the topological distance between each historical embedding vector in the set of historical embedding vectors and its nearest neighbor is calculated by an adaptive density clustering algorithm, and the local density radius is determined by combining the preset quantile position. With the historical embedding vector as the center, the set of hypersphere surface points is constructed by combining the local density radius to obtain a set of normal distribution boundary points.

[0080] The topological distance between the current embedded vector and each boundary point in the set of normal distribution boundary points is calculated based on the metric tensor, and the minimum value is extracted as the deviation metric. Anomaly score is obtained based on the statistical characteristics of the topological distance between the deviation metric and the set of historical embedded vectors, and the corresponding time window is determined as a candidate anomaly period based on the anomaly score.

[0081] For a location sequence in a spatiotemporal data block, a neighborhood matrix is ​​constructed and reconstruction weights are calculated using a local linear embedding algorithm. The location sequence is then extracted from the spatiotemporal data block, with each location point containing eastward and northward coordinates. For example, a ship's location sequence within a 60-second window consists of 60 coordinate points, where the first point is "eastward coordinate 548256.75 meters, northward coordinate 3375890.22 meters". For each location point, its 15 nearest neighbors are selected to form a neighborhood set. The Euclidean distance from the current point to each point in the neighborhood is calculated, forming a distance matrix. A neighborhood covariance matrix is ​​constructed based on the distance matrix, and the reconstruction weights are obtained by solving the eigenvalue decomposition of the covariance matrix. The reconstruction weights represent how a linear combination of neighborhood points can approximate the current point. For example, for the location point "eastward coordinate 548262.35 meters, northward coordinate 3375895.68 meters", the reconstruction weights corresponding to the 15 neighboring points are 0.12, 0.08, and 0.05, respectively, summing to 1.

[0082] Based on the calculated reconstruction weights, the position sequence is projected onto a high-dimensional feature space to obtain the current embedding vector. The projection process is achieved by constructing a global covariance matrix and solving for its eigenvectors. The eigenvectors with the smallest non-zero eigenvalues ​​are selected to form the embedding vector. For example, projecting 60 position points onto an 8-dimensional feature space yields an 8-dimensional current embedding vector with values ​​[0.235, -0.178, 0.053, -0.092, 0.126, -0.045, 0.067, 0.021]. Based on the distribution of the current embedding vector in the high-dimensional feature space, a metric tensor is constructed using the Riemannian geometric algorithm. The metric tensor describes the local scaling properties in different directions in the high-dimensional feature space and is obtained by calculating the distribution covariance of points around the embedding vector and performing singular value decomposition. For the aforementioned 8-dimensional embedding vector, the constructed metric tensor is an 8×8 matrix, with diagonal elements of 1.25, 0.98, and 1.12, representing the weights of each dimension, and off-diagonal elements representing the correlations between different dimensions.

[0083] Obtain spatiotemporal data blocks corresponding to historical normal navigation, for example, selecting 1000 spatiotemporal data blocks of normal navigation from the past month's lake area ship navigation data. Project the position sequences in the historical spatiotemporal data blocks onto a high-dimensional feature space to obtain a set of historical embedding vectors. Each historical spatiotemporal data block corresponds to an 8-dimensional historical embedding vector, such as [0.221, -0.165, 0.048, -0.085, 0.119, -0.042, 0.063, 0.020]. Based on the previously constructed metric tensor, calculate the topological distance between each historical embedding vector in the set and its nearest neighbor using an adaptive density clustering algorithm. The topological distance calculation uses Mahalanobis distance, which considers the influence of the metric tensor and reflects the true distribution characteristics in high-dimensional space better than Euclidean distance. For example, the topological distance between two historical embedding vectors [0.221, -0.165, 0.048, -0.085, 0.119, -0.042, 0.063, 0.020] and [0.228, -0.171, 0.050, -0.088, 0.122, -0.044, 0.065, 0.020] is calculated to be 0.026.

[0084] The local density radius is determined by combining the preset quantile locations. For example, the 95th percentile of the topological distance between historical embedding vectors is taken as the local density radius, which is 0.085 in this embodiment. A set of hypersphere surface points is constructed with each historical embedding vector as the center, combined with the local density radius, to obtain a set of normally distributed boundary points. Specifically, for each historical embedding vector, 50 points located on the surface of the hypersphere are uniformly sampled in 8-dimensional space, and the topological distance from these points to the center is 0.085. All surface points generated from historical embedding vectors are merged to form a set of normally distributed boundary points containing approximately 50,000 points.

[0085] The topological distance between the current embedded vector [0.235, -0.178, 0.053, -0.092, 0.126, -0.045, 0.067, 0.021] and each boundary point in the set of boundary points in the normal distribution is calculated using a metric tensor, and the minimum value is extracted as the deviation metric. For example, the calculated topological distance from the current embedded vector to the nearest boundary point is 0.023, which is the deviation metric. Anomaly scores are obtained based on the statistical characteristics of the topological distances between the deviation metric and the set of historical embedded vectors. The mean and standard deviation of the topological distance distribution between the deviation metric and historical embedded vectors are standardized. Assuming the mean of the topological distances between historical embedded vectors is 0.055 and the standard deviation is 0.018, the anomaly score of the current embedded vector is (0.023-0.055) / 0.018 = -1.778, indicating that the current embedded vector is within the normal distribution and has a small degree of deviation.

[0086] The anomaly score determines the corresponding time window as a candidate anomaly period. An anomaly score threshold of 2.0 is set; that is, when the anomaly score is greater than 2.0, the corresponding time window is determined to be a candidate anomaly period. In this embodiment, the anomaly score is -1.778, which is less than the threshold of 2.0, therefore this time window is determined to be a normal period. When encountering abnormal behaviors such as sharp turns or sudden deceleration, the corresponding embedding vector will deviate from the normal distribution, resulting in a larger deviation metric and anomaly score. For example, if a ship suddenly deviates from its course, the calculated anomaly score is 2.87, which is greater than the threshold of 2.0, and this time window is determined to be a candidate anomaly period.

[0087] In this embodiment, by performing nonlinear manifold modeling on the position sequence in the spatiotemporal data block, the original low-dimensional trajectory information is mapped to a high-dimensional feature space that can preserve the local geometric structure. This effectively overcomes the problem of difficulty in characterizing the intrinsic structure of complex navigation trajectories, significantly improving the ability to express and distinguish subtle trajectory changes. By introducing a metric tensor construction method based on Riemannian geometry into the high-dimensional feature space, distance calculation can adaptively reflect the directionality and curvature characteristics of the trajectory distribution, avoiding the defect of fixed metrics in inaccurate characterization of true similarity under complex trajectory morphology, and improving the accuracy and stability of anomaly detection. By adaptively constructing normal distribution boundaries using historical normal navigation data, the range of normal behavior can be dynamically adjusted with changes in trajectory distribution density, effectively reducing the risk of misjudgment caused by threshold migration under different lake areas and different navigation conditions, and enhancing environmental adaptability and generalization ability.

[0088] In one alternative implementation,

[0089] Based on the metric tensor, the topological distance between each historical embedding vector in the historical embedding vector set and its nearest neighbor historical embedding vector is calculated using an adaptive density clustering algorithm, and the local density radius is determined by combining this with the preset quantile positions, including:

[0090] Based on the metric tensor, the topological distance is obtained by calculating the geodesic path integral between different historical embedding vectors in the historical embedding vector set. For each historical embedding vector, the historical embedding vector corresponding to the minimum topological distance is extracted as the nearest neighbor historical embedding vector. The topological distance from each historical embedding vector to the nearest neighbor historical embedding vector is recorded to obtain the nearest neighbor distance set.

[0091] The topological distances in the nearest neighbor distance set are arranged in ascending order to obtain an ordered distance sequence. Based on the preset quantile position, the corresponding index position is located in the ordered distance sequence and the topological distance at the index position is extracted as the global quantile distance.

[0092] The neighborhood density value is obtained by calculating the number of historical embedding vectors contained in the hypersphere with the global quantile distance as the radius for each historical embedding vector in the historical embedding vector set. The local density radius is obtained by weighted fusion of the nearest neighbor topological distance corresponding to each historical embedding vector and the reciprocal of the neighborhood density value.

[0093] Based on the metric tensor, the topological distance is obtained by calculating the geodesic path integral between different historical embedding vectors in the historical embedding vector set. The geodesic path integral takes into account the curvature characteristics in the high-dimensional feature space, and can more accurately reflect the true distance of the ship trajectory in the feature space. For any two historical embedding vectors, such as vectors [0.221, -0.165, 0.048, -0.085, 0.119, -0.042, 0.063, 0.020] and vectors [0.228, -0.171, 0.050, -0.088, 0.122, -0.044, 0.065, 0.020], a parameterized path connecting the two points is constructed. The path is discretized into 100 small segments. The metric tensor is applied to each small segment to calculate the local distance. The geodesic distance is obtained by accumulating the distances of all small segments. For example, the geodesic topological distance between the two vectors mentioned above is 0.031, which is slightly larger than the Euclidean distance of 0.026, because the geodesic takes into account the curvature changes in the feature space.

[0094] For each historical embedding vector, after calculation, the historical embedding vector corresponding to the minimum topological distance is extracted as the nearest neighbor historical embedding vector. For example, for the historical embedding vector [0.221, -0.165, 0.048, -0.085, 0.119, -0.042, 0.063, 0.020], the minimum topological distance calculated among 1000 historical embedding vectors is 0.031, and the corresponding nearest neighbor historical embedding vectors are [0.228, -0.171, 0.050, -0.088, 0.122, -0.044, 0.065, 0.020]. The topological distance from each historical embedding vector to the nearest neighbor historical embedding vector is recorded to obtain the nearest neighbor distance set. For 1000 historical embedding vectors, 1000 nearest neighbor topological distance values ​​are obtained, forming the nearest neighbor distance set.

[0095] The topological distances in the nearest neighbor distance set are arranged in ascending order to obtain an ordered distance sequence. For example, after sorting 1000 nearest neighbor topological distance values, the ordered distance sequence starts with [0.015, 0.018, 0.022, 0.025, 0.027, 0.031, ...] and ends with [..., 0.103, 0.112, 0.125]. Based on a preset quantile position, the corresponding index position in the ordered distance sequence is located, and the topological distance at that index position is extracted as the global quantile distance. Assuming the preset quantile position is 95%, the corresponding index position is 950, and the global quantile distance is the value of the 950th element in the ordered distance sequence, which is 0.085. The global quantile distance represents that the nearest neighbor distance between 95% of the vector pairs in the historical embedding vector set is less than this value.

[0096] The neighborhood density value is obtained by calculating the number of historical embedding vectors contained within a hypersphere with a global quantile distance as its radius for each historical embedding vector in the historical embedding vector set. Taking the historical embedding vector [0.221, -0.165, 0.048, -0.085, 0.119, -0.042, 0.063, 0.020] as an example, a hypersphere is constructed with a global quantile distance of 0.085 as its radius, and the number of historical embedding vectors contained within the sphere is calculated. Assuming the calculated result is 45, meaning there are 45 historical embedding vectors with a topological distance less than 0.085 from the current vector, the neighborhood density value is 45. Historical embedding vectors located in high-density regions have larger neighborhood density values; historical embedding vectors located in sparse regions have smaller neighborhood density values.

[0097] The local density radius is obtained by weighting and fusing the nearest neighbor topological distance and the reciprocal of the neighborhood density value for each historical embedding vector. The weighted fusion uses a linear combination method, with the weight coefficients adjusted according to the actual application scenario. Assuming the weight coefficients are 0.7 and 0.3, for the aforementioned historical embedding vector, the nearest neighbor topological distance is 0.031, the neighborhood density value is 45, and the reciprocal of the neighborhood density value is 1 / 45 = 0.022. The weighted fusion calculation is: 0.7 × 0.031 + 0.3 × 0.022 = 0.0217 + 0.0066 = 0.0283, resulting in a local density radius of 0.0283.

[0098] The local density radius of each historical embedding vector was calculated using the aforementioned method, resulting in 1000 local density radius values. These local density radius values ​​are adapted to the density distribution characteristics of different regions in the feature space. In the analysis of ship trajectories in the lake area, this means that a smaller local density radius is used for high-density trajectories in regular channel areas, while a larger local density radius is used for low-density trajectories in the occasionally traversed lake edge areas. For example, the average local density radius of historical embedding vectors in the main channel area is 0.032, while the average local density radius of historical embedding vectors in the lake edge area is 0.078.

[0099] For specific ship types and application scenarios, such as passenger ships, experiments determined that weighting coefficients of 0.65 and 0.35 yielded the best results; while for cargo ships, weighting coefficients of 0.75 and 0.25 yielded the best results. The adaptive nature of the local density radius allows the anomaly detection system to dynamically adjust its judgment criteria based on the trajectory density characteristics of different areas, reducing false alarm rates. For example, in a dense waterway area, the historical embedding vector density value corresponding to normal navigation of a passenger ship is as high as 68, and the local density radius is calculated as 0.65×0.028+0.35×(1 / 68)=0.0182+0.0018=0.02. The smaller local density radius makes the system more sensitive to trajectory anomalies in this area. In contrast, in remote waters of a lake area, the historical embedding vector density value corresponding to normal navigation of a cargo ship is only 12, and the local density radius is calculated as 0.75×0.075+0.35×(1 / 12)=0.0563+0.0292=0.0855. The larger local density radius allows for a greater range of trajectory variation.

[0100] In this embodiment, by introducing a geodesic-based topological distance calculation method in the high-dimensional embedding space, the similarity measurement between historical trajectory samples can be expanded along the real structure of the data manifold, avoiding the problem of distortion in the characterization of proximity relationships in nonlinear distribution scenarios, and improving the accuracy and stability of trajectory distribution modeling. By taking the overall distribution of nearest neighbor topological distance as the basis and combining it with the extraction of global quantile distance from quantile points, an adaptive characterization of the spatial scale of historical normal navigation samples is achieved, which can reduce the risk of boundary shift caused by uneven sample distribution. By calculating the neighborhood density of each historical embedding vector under a unified global quantile scale and fusing the nearest neighbor topological distance with the neighborhood density to determine the local density radius, the normal distribution boundary has differentiated scaling capabilities in dense and sparse regions, significantly enhancing the ability to characterize complex distribution patterns.

[0101] In one alternative implementation,

[0102] Extracting the motion state sequence and environmental information sequence corresponding to the candidate abnormal time period and reconstructing the phase space in the time dimension, calculating the mutual information transmission strength and direction between the motion state sequence and the environmental information sequence in the reconstructed phase space, and constructing a causal influence matrix includes:

[0103] Based on the start and end times of the candidate abnormal time period, the motion state sequence and environmental information sequence within the range from the start time to the end time are extracted from the spatiotemporal data block.

[0104] For each subsequence in the motion state sequence and the environment information sequence, the optimal delay time is calculated using a time delay embedding algorithm and a pre-set mutual information function, the optimal embedding dimension is calculated based on the false nearest neighbor algorithm, and the phase space is reconstructed based on the optimal delay time and the optimal embedding dimension to obtain the motion state phase space vector set and the environment information phase space vector set.

[0105] The forward information transfer amount from the environmental information variables in the environmental information phase space vector set to the motion state variables in the motion state phase space vector set and the reverse information transfer amount from the motion state variables in the motion state phase space vector set to the environmental information variables in the environmental information phase space vector set are calculated using the transfer entropy algorithm. A bidirectional transfer entropy matrix is ​​then constructed based on the forward and reverse information transfer amounts.

[0106] The transfer entropy values ​​in the bidirectional transfer entropy matrix are statistically tested based on the significance test algorithm to determine the effective transfer intensity. Based on the directional marker information of the effective transfer intensity, the starting variable and target variable of the transfer are used to construct a causal influence matrix with motion state variables and environmental information variables as nodes and effective transfer intensity and transfer direction as directed edges.

[0107] Based on the start and end times of candidate abnormal time periods, motion state sequences and environmental information sequences within the range from the start to the end time are extracted from the spatiotemporal data block. For example, for the identified candidate abnormal time period from 09:45:20 to 09:47:50 on January 5, 2025, motion state sequences and environmental information sequences within this time range are extracted from the corresponding spatiotemporal data block. The motion state sequences include ship speed, heading angle, and speed data, while the environmental information sequences include wind speed, wind direction, lake level, and visibility data. The sampling interval is 1 second, totaling 150 data points.

[0108] For each subsequence in the motion state sequence and environmental information sequence, the optimal delay time is calculated using a time delay embedding algorithm and a pre-set mutual information function. Taking the speed subsequence as an example, the mutual information value between the sequence and itself is calculated under different time delays: for a delay time of 1 second, the mutual information value is 2.85; for a delay time of 2 seconds, the mutual information value is 2.63; for a delay time of 3 seconds, the mutual information value is 2.42; and so on. The delay time corresponding to the first local minimum of the mutual information function is selected as the optimal delay time, which is 5 seconds in this embodiment. For the heading angle subsequence, the optimal delay time is calculated to be 4 seconds; for the wind speed subsequence, the optimal delay time is 8 seconds; for the wind direction subsequence, the optimal delay time is 7 seconds; for the water level subsequence, the optimal delay time is 12 seconds; and for the visibility subsequence, the optimal delay time is 10 seconds.

[0109] The optimal embedding dimension is calculated based on the false nearest neighbor (FNN) algorithm. The FNN algorithm determines the optimal embedding dimension by calculating the proportion of false nearest neighbors (FNNs) for different embedding dimensions. Taking the airspeed subsequence as an example, for embedding dimension 2, the proportion of FNNs is 0.35; for embedding dimension 3, it is 0.18; for embedding dimension 4, it is 0.06; and for embedding dimension 5, it is 0.05. When the proportion of FNNs is lower than a preset threshold of 0.1, the corresponding embedding dimension is selected as the optimal embedding dimension, which is 4 in this embodiment. For the heading angle subsequence, the optimal embedding dimension is 3; for the wind speed subsequence, it is 4; for the wind direction subsequence, it is 3; for the water level subsequence, it is 2; and for the visibility subsequence, it is 3.

[0110] Phase space reconstruction is performed based on the optimal delay time and optimal embedding dimension to obtain the set of motion state phase space vectors and the set of environmental information phase space vectors. Taking the speed subsequence as an example, the optimal delay time is 5 seconds and the optimal embedding dimension is 4. For the speed value v at time point t, the constructed phase space vector is [v(t), v(t+5), v(t+10), v(t+15)]. For example, the speed value is 8.5 knots at t=10 seconds, 8.7 knots at t=15 seconds, 8.3 knots at t=20 seconds, and 7.9 knots at t=25 seconds, then the phase space vector is [8.5, 8.7, 8.3, 7.9]. Similarly, phase space reconstruction is performed on all subsequences to obtain the set of motion state phase space vectors and the set of environmental information phase space vectors.

[0111] The transfer entropy algorithm is used to calculate the positive information transfer from environmental information variables in the environmental information phase space vector set to motion state variables in the motion state phase space vector set, and the negative information transfer from motion state variables in the motion state phase space vector set to environmental information variables in the environmental information phase space vector set. Based on conditional probability distribution, the transfer entropy algorithm calculates the contribution of source variables to the prediction of the future state of the target variable. Taking the influence of wind speed on flight speed as an example, the transfer entropy value from the wind speed phase space vector to the flight speed phase space vector is calculated, and the difference between the joint probability distribution and conditional probability distribution of the wind speed and flight speed phase space vectors is calculated. In this embodiment, the positive information transfer from wind speed to flight speed is 0.42, indicating that wind speed changes provide 42% of the information about future flight speed changes; the negative information transfer from flight speed to wind speed is 0.08, indicating that flight speed changes only provide 8% of the information about future wind speed changes.

[0112] The transfer entropy between all environmental information variables and motion state variables is calculated to construct a two-way transfer entropy matrix. Rows in the matrix represent source variables, and columns represent target variables. For example, a value of 0.42 in the wind speed row and ship speed column indicates the positive information transfer from wind speed to ship speed; a value of 0.08 in the ship speed row and wind speed column indicates the negative information transfer from ship speed to wind speed. The complete two-way transfer entropy matrix contains all transfer entropy values ​​between the six variables: wind speed, wind direction, water level, visibility, ship speed, and heading angle.

[0113] The significance test algorithm is used to statistically test the transfer entropy values ​​in the two-way transfer entropy matrix and determine the effective transfer strength. The significance test employs a permutation test method, generating 5000 randomized sequences by randomly shuffling the time series order, and calculating the distribution of transfer entropy values ​​in the randomized sequences. If the original transfer entropy value exceeds the 95th percentile of the randomized transfer entropy value, the transfer entropy value is considered statistically significant, and the corresponding information transfer relationship is valid. Taking the transfer entropy value from wind speed to ship speed of 0.42 as an example, the p-value obtained through the permutation test is 0.002, which is less than the significance level of 0.05, indicating that the information transfer relationship from wind speed to ship speed is valid, and the effective transfer strength is 0.42. However, the p-value corresponding to the transfer entropy value from ship speed to wind speed of 0.08 is 0.35, which is greater than the significance level of 0.05, indicating that the information transfer relationship from ship speed to wind speed is not significant, and the effective transfer strength is 0.

[0114] Based on the directional markers of the effective transmission intensity, indicating the starting and target variables of the information transmission, a causal influence matrix is ​​constructed, with motion state variables and environmental information variables as nodes and effective transmission intensity and transmission direction as directed edges. In this embodiment, the effective information transmission relationships include: wind speed to ship speed, intensity 0.42; wind direction to ship heading angle, intensity 0.38; wind direction to ship speed, intensity 0.25; water level to ship speed, intensity 0.15; and ship speed to ship heading angle, intensity 0.33. The constructed causal influence matrix reflects the causal relationship network between the variables and is visualized as a directed graph, where nodes represent variables, directed edges represent information transmission directions, and the thickness of the edges represents the magnitude of the transmission intensity.

[0115] In this embodiment, by selectively extracting the motion state sequence and environmental information sequence within the candidate anomaly period, subsequent analysis focuses on the key time range of anomaly occurrence. This effectively reduces interference from irrelevant data and improves the relevance and timeliness of anomaly cause analysis. By adaptively determining the time delay and embedding dimension and reconstructing the phase space, the dynamic evolution process of motion state and environmental information is fully unfolded in a high-dimensional state space. This avoids the shortcomings of fixed parameter settings, which are difficult to adapt to the dynamic characteristics of different variables. It can more realistically reflect the intrinsic dynamic behavior of the system and provide a reliable state representation basis for causal relationship analysis. By introducing bidirectional transfer entropy to quantify the direction and intensity of information transfer between environmental information and motion state, the influence paths between various factors within the anomaly period can be clearly characterized, thereby significantly improving the explanatory power of anomaly mechanism analysis.

[0116] In one alternative implementation,

[0117] The causal influence matrix is ​​subjected to spectral decomposition to obtain dominant eigenvalues ​​and determine the degree of dominance of environmental changes on changes in motion state. Based on the degree of dominance, the candidate abnormal time periods are judged to obtain a confirmed abnormal event set, including:

[0118] The causal influence matrix is ​​subjected to spectral decomposition to extract all eigenvalues ​​and corresponding eigenvectors. The eigenvalues ​​are sorted in descending order according to their magnitudes, and the first eigenvalue is taken as the dominant eigenvalue. The eigenvectors corresponding to the dominant eigenvalues ​​are extracted as the dominant eigenvectors.

[0119] Calculate the absolute value of each element in the dominant feature vector and determine the element with the largest absolute value. Obtain the index position of the element with the largest absolute value in the dominant feature vector and map it to the corresponding variable node in the causal influence matrix. Determine the variable type corresponding to the variable node and determine the dominant direction of causal influence based on the variable type.

[0120] Based on the dominant feature value, the intensity of the dominant effect is calculated by the spectral radius measurement method and compared with the intensity of the dominant effect with a preset threshold. The degree of the dominant effect of environmental change on motion state change is determined according to the comparison result.

[0121] Based on the degree of dominance and the direction of dominance, the candidate abnormal time period is classified into different types. If the degree of dominance exceeds a preset threshold and the direction of dominance is dominated by environmental information variables, the candidate abnormal time period is classified into an environment-driven abnormal event. Otherwise, the candidate abnormal time period is classified into an autonomous behavior abnormal event. The determination is repeated to obtain the confirmed abnormal event set.

[0122] The causal influence matrix is ​​subjected to spectral decomposition to extract all eigenvalues ​​and corresponding eigenvectors. This matrix reflects the information transmission relationship between six variables: wind speed, wind direction, water level, visibility, navigation speed, and heading angle, and is represented as a 6×6 asymmetric matrix. Spectral decomposition uses eigenvalue decomposition to calculate the matrix's eigenvalues ​​and corresponding eigenvectors. For the causal influence matrix in the previous example, the matrix eigenvalue decomposition algorithm yields six eigenvalues: 1.42, 0.85, 0.53, 0.28, 0.12, and 0.05. Each eigenvalue corresponds to an eigenvector; for example, the eigenvector corresponding to eigenvalue 1.42 is [0.55, 0.63, 0.18, 0.08, -0.42, -0.31], and the eigenvector corresponding to eigenvalue 0.85 is [0.32, 0.27, 0.15, 0.12, 0.65, 0.58].

[0123] The eigenvalues ​​are sorted in descending order of their moduli, and the first eigenvalue is taken as the dominant eigenvalue. The eigenvector corresponding to the dominant eigenvalue is then extracted as the dominant eigenvector. In this embodiment, the eigenvalues ​​are sorted by moduli as 1.42, 0.85, 0.53, 0.28, 0.12, and 0.05. Therefore, the dominant eigenvalue is 1.42, and the dominant eigenvector is [0.55, 0.63, 0.18, 0.08, -0.42, -0.31]. The dominant eigenvalue represents the most important pattern in the matrix, and the dominant eigenvector reflects the contribution or influence of each variable on this main pattern.

[0124] Calculate the absolute value of each element in the dominant eigenvector and determine the element with the largest absolute value. The absolute values ​​of the elements in the dominant eigenvector [0.55, 0.63, 0.18, 0.08, -0.42, -0.31] are 0.55, 0.63, 0.18, 0.08, 0.42, and 0.31, respectively. The element with the largest absolute value is 0.63, located in the second position of the eigenvector. Obtain the index position of the element with the largest absolute value in the dominant eigenvector and map it to the corresponding variable node in the causal influence matrix. The index position is 1, which maps to the second variable node in the causal influence matrix, i.e., the wind direction variable.

[0125] The variable type corresponding to each variable node is determined, and the dominant direction of causal influence is identified based on the variable type. In the causal influence matrix, variable nodes are divided into environmental information variables and motion state variables. Environmental information variables include wind speed, wind direction, water level, and visibility, while motion state variables include ship speed and heading angle. The second variable node corresponds to wind direction, which is an environmental information variable. Determining that the dominant direction of causal influence is dominated by environmental information variables indicates that environmental factors, especially changes in wind direction, have a significant impact on ship motion state during the current candidate anomaly period.

[0126] Based on the dominant eigenvalue, the intensity of the dominant effect is calculated using the spectral radius measurement method and compared with a preset dominant effect intensity threshold. The spectral radius measurement method calculates the ratio of the dominant eigenvalue to the sum of the magnitudes of all eigenvalues ​​as the dominant effect intensity. In this embodiment, the sum of the magnitudes of all eigenvalues ​​is 1.42 + 0.85 + 0.53 + 0.28 + 0.12 + 0.05 = 3.25, and the dominant eigenvalue is 1.42. Therefore, the dominant effect intensity is 1.42 / 3.25 = 0.437, indicating that the dominant mode accounts for 43.7% of the overall influence. Comparing the dominant effect intensity of 0.437 with the preset dominant effect intensity threshold of 0.4, 0.437 is greater than 0.4, indicating that the dominant mode has a strong influence.

[0127] The degree of dominance of environmental changes on changes in motion state was determined based on the comparison results. Since the dominance intensity of 0.437 exceeded the preset dominance intensity threshold of 0.4, and the dominant direction was dominated by environmental information variables, it was determined that environmental changes had a strong dominant effect on changes in motion state. During the candidate abnormal period, the changes in the ship's motion state were mainly affected by environmental factors, especially wind direction changes, rather than by the ship's autonomous behavior.

[0128] The candidate abnormal time period was classified based on the degree and direction of dominance. The degree of dominance was 0.437, exceeding the preset threshold of 0.4, and the direction of dominance was dominated by environmental information variables. Therefore, the candidate abnormal time period was classified as an environment-driven abnormal event. Environment-driven abnormal events indicate that the anomaly in the ship's trajectory is mainly caused by changes in environmental factors, such as strong winds or surges. The ship's response is within the normal operating range and does not require a high-level warning.

[0129] For another candidate anomaly period, from 10:35:10 on January 5, 2025 to 10:37:45 on January 5, 2025, the aforementioned analysis process was repeated. A causal influence matrix was constructed and spectral decomposition was performed, yielding eigenvalues ​​of 1.05, 0.88, 0.67, 0.45, 0.23, and 0.11, with a dominant eigenvalue of 1.05 and a dominant eigenvector of [0.25, 0.18, 0.12, 0.05, 0.68, 0.65]. The absolute values ​​of each element in the dominant eigenvector were calculated to be 0.25, 0.18, 0.12, 0.05, 0.68, and 0.65, respectively. The element with the largest absolute value was 0.68, corresponding to index position 4, which maps to the fifth variable in the causal influence matrix, flight speed, belonging to the motion state variable. The ratio of the dominant eigenvalue to the sum of the magnitudes of all eigenvalues ​​is 1.05 / 3.39 = 0.31, which is less than the preset dominant influence strength threshold of 0.4. Furthermore, the dominant direction is determined by the motion state variable. Therefore, this candidate abnormal period is classified as an autonomous behavior-type abnormal event. Autonomous behavior-type abnormal events indicate that the anomaly in the ship's trajectory is mainly caused by the ship's autonomous operation, potentially involving violations or emergency avoidance, requiring a high-level early warning.

[0130] Repeat the above judgment process to determine the type of all candidate abnormal time periods, thus obtaining a confirmed set of abnormal events. In the lake area vessel monitoring system, the confirmed set of abnormal events will be used for subsequent early warning processing. Environment-driven abnormal events usually only need to be recorded and given a low-level alert, while autonomous behavior-related abnormal events require an immediate high-level alert to regulatory personnel and may trigger an emergency response mechanism.

[0131] In this embodiment, by performing spectral decomposition on the causal influence matrix and extracting the dominant eigenvalues ​​and corresponding eigenvectors, complex multivariate causal relationships can be compressed into a small amount of representative dominant structural information. This allows for a more intuitive and stable depiction of the dominant pattern of the overall causal effect of the system during abnormal periods, improving the robustness and interpretability of the causal analysis results. By using the variable with the largest weight in the dominant eigenvector to identify key influencing nodes and further determining their variable type and influence direction, the problem of difficulty in identifying dominant factors due to the mutual coupling of multiple variables is avoided. This enables the factors that play a dominant role in the anomaly formation process to be clearly distinguished, improving the accuracy of anomaly tracing and classification. By introducing a dominant effect intensity measure based on spectral radius, the causal influence is transformed from a qualitative judgment into a quantifiable intensity index, which is compared with a preset threshold to achieve a graded judgment of the degree of influence of environmental factors on the motion state, enhancing the accuracy of anomaly type classification.

[0132] Figure 2 This is a flowchart illustrating the abnormal event analysis process of the data fusion-based method for identifying and warning of abnormal ship trajectories in lake areas, as described in this embodiment of the invention.

[0133] In one alternative implementation,

[0134] The confirmed abnormal event set is traversed in chronological order to identify persistent and sudden abnormalities. Based on the duration of persistent abnormalities and the degree of deviation of sudden abnormalities, tiered early warning information is generated and output, including:

[0135] All abnormal events in the confirmed abnormal event set are sorted in ascending order by timestamp and traversed. The time interval between the start time and end time of each abnormal event is determined as the duration. The numerical change features of the motion state sequence within the abnormal event period are extracted.

[0136] The temporal change rate of the motion state sequence within the period of the abnormal event is calculated by the sliding window method. Abnormal events with a temporal change rate lower than a preset temporal change rate threshold and a duration exceeding a preset duration threshold are marked as continuous abnormal events; otherwise, the abnormal events are marked as sudden abnormal events.

[0137] The evolution rate corresponding to the persistent anomaly is calculated by a trend fitting method, and the duration level is determined by combining the duration. The peak deviation and change gradient corresponding to the sudden anomaly are extracted and the deviation level is determined. Based on the duration level and the deviation level, an early warning level is generated and combined with a pre-determined anomaly type and a pre-acquired dominant direction to obtain graded early warning information and output it.

[0138] All anomalous events in the confirmed anomalous event set are sorted in ascending order by timestamp and traversed. The confirmed anomalous event set contains multiple anomalous events that have been confirmed through causal impact analysis. The anomalous events are then sorted by their occurrence time. For example, the top three anomalous events after sorting are: Anomalous event 1 occurred from 09:45:20 to 09:47:50 on January 5, 2025; Anomalous event 2 occurred from 10:35:10 to 10:37:45 on January 5, 2025; and Anomalous event 3 occurred from 13:25:30 to 13:42:15 on January 5, 2025.

[0139] The time interval between the start and end times of each anomalous event is defined as the duration. Numerical variation characteristics of the motion state sequence within the anomalous event period are extracted. For anomalous event 1, the start time is 2025-01-05 09:45:20, the end time is 2025-01-05 09:47:50, and the duration is 2 minutes and 30 seconds. The motion state sequence within the anomalous event 1 period is extracted, including airspeed and heading angle data. Numerical variation characteristics include mean, standard deviation, maximum value, minimum value, and variation range. In anomalous event 1, the mean airspeed is 6.8 knots, the standard deviation is 1.2 knots, the maximum value is 8.5 knots, the minimum value is 5.3 knots, and the variation range is 3.2 knots; the mean heading angle is 125 degrees, the standard deviation is 8.5 degrees, the maximum value is 140 degrees, the minimum value is 110 degrees, and the variation range is 30 degrees. The duration of anomalous event 2 is 2 minutes and 35 seconds, and the duration of anomalous event 3 is 16 minutes and 45 seconds.

[0140] The sliding window method was used to calculate the temporal rate of change of motion state sequences within the period of an abnormal event. The sliding window method uses a fixed-size time window, which slides gradually within the period of the abnormal event, calculating the rate of change of the motion state sequence within each window. The window size was set to 20 seconds, and the sliding step size was 5 seconds. For each window, the average rate of change of the data within the window was calculated, which is the ratio of the average difference between adjacent sampling points to the sampling interval. Taking abnormal event 1 as an example, the average rate of change of air speed in the first window from 09:45:20 to 09:45:40 was -0.11 knots / second, and the average rate of change of heading angle was 0.85 degrees / second; the average rate of change of air speed in the second window from 09:45:25 to 09:45:45 was -0.09 knots / second, and the average rate of change of heading angle was 0.78 degrees / second. The rate of change for all sliding windows was calculated. The average of the rates of change of all sliding windows during the entire period of the abnormal event is calculated as the time-series rate of change. The time-series rate of change for abnormal event 1 is 0.08 knots / second for air speed and 0.65 degrees / second for heading angle.

[0141] Anomalies with a time-series change rate lower than a preset time-series change rate threshold and a duration exceeding a preset duration threshold are marked as persistent anomalies; otherwise, they are marked as sudden anomalies. The preset time-series change rate thresholds are set to airspeed 0.1 knots / second and heading angle 1.0 degrees / second; the preset duration threshold is set to 5 minutes. For anomaly 1, the time-series change rate is 0.08 knots / second for airspeed and 0.65 degrees / second for heading angle, both lower than the preset thresholds, but the duration of 2 minutes and 30 seconds is less than the preset duration threshold of 5 minutes; therefore, it is marked as a sudden anomaly. For anomaly 2, the time-series change rate is 0.12 knots / second for airspeed and 1.2 degrees / second for heading angle, both higher than the preset thresholds, and the duration of 2 minutes and 35 seconds is less than the preset duration threshold; therefore, it is also marked as a sudden anomaly. For abnormal event 3, the time-series change rate is 0.05 knots / second for air speed and 0.35 degrees / second for heading angle, both of which are lower than the preset threshold, and the duration of 16 minutes and 45 seconds is greater than the preset duration threshold, so it is marked as a continuous abnormality.

[0142] The evolution rate corresponding to the sustained anomaly was calculated using a trend fitting method, and the duration was combined with the duration to determine the duration level. The trend fitting method uses linear regression to fit the motion state sequence within the sustained anomaly period, obtaining the slope of the trend line as the evolution rate. Taking anomaly event 3 as an example, linear regression fitting of the airspeed data yielded a slope of -0.03 knots / second, indicating an average decrease in airspeed of 1.8 knots per minute; linear regression fitting of the heading angle data yielded a slope of 0.2 degrees / second, indicating an average change in heading angle of 12 degrees per minute. The duration level was determined based on the duration, divided into three levels: low, medium, and high. Duration less than 10 minutes was classified as low, 10-20 minutes as medium, and greater than 20 minutes as high. Anomaly event 3 lasted 16 minutes and 45 seconds, belonging to the medium duration level.

[0143] Extract the peak deviation and gradient of sudden anomalies and determine the deviation level. The peak deviation is calculated as the difference between the maximum or minimum value of the motion state sequence during the anomaly event period and the normal behavior pattern. The gradient is calculated as the ratio of the peak deviation to the time required to reach the peak. Taking anomaly 1 as an example, the airspeed suddenly drops from the normal value of 8.5 knots to the minimum value of 5.3 knots, with a peak deviation of 3.2 knots; the time to reach the peak is 55 seconds, and the gradient is 0.058 knots / second. The heading angle suddenly changes from the normal value of 110 degrees to the maximum value of 140 degrees, with a peak deviation of 30 degrees; the time to reach the peak is 45 seconds, and the gradient is 0.67 degrees / second. The deviation level is determined based on the peak deviation and gradient, and is also divided into three levels: low, medium, and high. For speed, a peak deviation of less than 2 knots is classified as low, 2-5 knots as medium, and greater than 5 knots as high. For heading angle, a peak deviation of less than 20 degrees is classified as low, 20-45 degrees as medium, and greater than 45 degrees as high. Anomaly 1 has a speed deviation of 3.2 knots, which falls under the medium deviation category; its heading angle deviation is 30 degrees, also falling under the medium deviation category. Considering both speed and heading angle deviations, the overall deviation level for anomaly 1 is determined to be medium.

[0144] Warning levels are generated based on duration and deviation levels and combined with pre-determined anomaly types and pre-acquired dominant directions to obtain graded warning information, which is then output. Warning levels are determined according to either duration or deviation levels and are also divided into three levels: low, medium, and high. For persistent anomalies, the warning level is consistent with the duration level; for sudden anomalies, the warning level is consistent with the deviation level. Anomaly 3 is a persistent anomaly with a medium duration level, therefore its warning level is medium; anomaly 1 is a sudden anomaly with a medium deviation level, therefore its warning level is medium. Combining the pre-determined anomaly type and dominant direction, complete graded warning information is generated. Anomaly 1 is an environment-driven anomaly, with wind direction as the dominant direction, and its warning level is medium. Therefore, the generated warning information is: "Medium Warning: An environment-driven sudden anomaly has been detected. Changes in wind direction have caused abnormal changes in the ship's course and speed. Enhanced monitoring is recommended." Anomaly 3 is classified as an autonomous behavior anomaly, with speed as the dominant factor and a medium warning level. Therefore, the generated warning message is: "Medium warning: Autonomous behavior anomaly detected. The ship's speed is abnormally and continuously decreasing, which may indicate a malfunction or violation of regulations. It is recommended to contact the ship immediately and dispatch a patrol boat to investigate."

[0145] In this embodiment, by uniformly organizing confirmed abnormal events in chronological order and quantifying their duration, abnormalities are no longer merely discrete alarm points but are characterized as complete event processes. This more accurately reflects the actual impact range and duration of the abnormalities on system operation. By introducing a joint judgment mechanism of time-series change rate and duration, the automatic distinction between continuous and sudden abnormalities is achieved. This overcomes the problem of treating abnormalities with different evolutionary characteristics indiscriminately and making it difficult to distinguish between slowly accumulating risks and instantaneous violent disturbances. It improves the rationality and refinement of abnormal event classification. By using differentiated indicators such as evolution rate, duration level, peak deviation, and change gradient for quantitative evaluation of continuous and sudden abnormalities, the characterization of abnormal intensity is more consistent with its actual evolutionary characteristics. This avoids the inadequacy of a single indicator in taking into account different abnormal forms and significantly enhances the accuracy and stability of abnormality severity assessment.

[0146] A second aspect of this invention provides a data fusion-based system for identifying and warning of anomalies in ship trajectories in lake areas, comprising:

[0147] The data fusion module is used to acquire multi-source trajectory data and perform spatiotemporal alignment processing by combining it with pre-set timestamp association and spatial coordinate mapping rules to obtain fused trajectory data. The fused trajectory data is divided into multiple time windows and a spatiotemporal data block corresponding to each time window is constructed.

[0148] An anomaly identification module is used to map the position sequence in the spatiotemporal data block to a high-dimensional feature space by embedding spatial mapping, and compare the distribution position of the current spatiotemporal data block with the spatiotemporal data block corresponding to historical normal navigation in the topological space in the high-dimensional feature space to identify candidate abnormal time periods.

[0149] The causal analysis module is used to extract the motion state sequence and environmental information sequence corresponding to the candidate abnormal time period and reconstruct the phase space in the time dimension. In the reconstructed phase space, the mutual information transmission strength and transmission direction between the motion state sequence and the environmental information sequence are calculated and a causal influence matrix is ​​constructed. The causal influence matrix is ​​subjected to spectral decomposition to obtain the dominant feature value and to determine the degree of dominance of environmental changes on motion state changes. Based on the degree of dominance, the candidate abnormal time period is judged to obtain a set of confirmed abnormal events.

[0150] The anomaly warning module is used to traverse the confirmed abnormal event set in chronological order and identify persistent anomalies and sudden anomalies. It generates and outputs graded warning information based on the duration of the persistent anomaly and the degree of deviation of the sudden anomaly.

[0151] A third aspect of the present invention provides an electronic device, comprising:

[0152] A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke instructions stored in the memory to perform the aforementioned method.

[0153] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0154] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

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

Claims

1. A method for identifying and issuing early warnings of anomalies in vessel trajectories in lake areas based on data fusion, characterized in that: include: Multi-source trajectory data is acquired and spatiotemporally aligned by combining it with pre-set timestamp association and spatial coordinate mapping rules to obtain fused trajectory data. The fused trajectory data is divided into multiple time windows and spatiotemporal data blocks corresponding to each time window are constructed. By embedding spatial mapping, the position sequence in the spatiotemporal data block is mapped to a high-dimensional feature space. In the high-dimensional feature space, the distribution positions of the current spatiotemporal data block and the spatiotemporal data blocks corresponding to historical normal navigation are compared in the topological space, and candidate abnormal time periods are identified. Extract the motion state sequence and environmental information sequence corresponding to the candidate abnormal time period and reconstruct the phase space in the time dimension. Calculate the mutual information transmission strength and direction between the motion state sequence and the environmental information sequence in the reconstructed phase space and construct a causal influence matrix. Perform spectral decomposition on the causal influence matrix to obtain the dominant feature value and determine the degree of dominance of environmental changes on motion state changes. Based on the degree of dominance, determine the candidate abnormal time period to obtain a confirmed abnormal event set. The confirmed abnormal event set is traversed in chronological order to identify persistent and sudden abnormalities. Based on the duration of the persistent abnormalities and the degree of deviation of the sudden abnormalities, hierarchical early warning information is generated and output.

2. The method according to claim 1, characterized in that, Multi-source trajectory data is acquired and spatiotemporally aligned using pre-set timestamp association and spatial coordinate mapping rules to obtain fused trajectory data. The fused trajectory data is then divided into multiple time windows, and a spatiotemporal data block corresponding to each time window is constructed, including: The position sequence and motion state sequence of the ship and the environmental information sequence of the lake area are obtained and combined to obtain multi-source trajectory data. The timestamps corresponding to the multi-source trajectory data are extracted as association keys to establish a synchronization index table. Data records with timestamp differences within the tolerance range are grouped and associated through the synchronization index table to obtain associated data groups. The position sequence in the associated data group is transformed to a standard coordinate system and then interpolated to generate an aligned position sequence. The motion state sequence in the associated data group is resampled and unified to a standard sampling frequency to generate an aligned motion sequence. The environmental information sequence in the associated data group is resampled and unified to a standard sampling frequency to generate an aligned environment sequence. The aligned position sequence, the aligned motion sequence, and the aligned environment sequence are merged one by one according to the timestamp to generate fused trajectory data. The fused trajectory data is divided into preliminary window sequences according to a fixed step size. The fused trajectory data segments corresponding to each preliminary window in the preliminary window sequence are extracted. The complexity index of the alignment position sequence in the fused trajectory data segments is calculated. The time span of the preliminary window is adjusted according to the complexity index to generate an adaptive window sequence. Extract the fused trajectory data fragments corresponding to each adaptive window in the adaptive window sequence and perform differential calculations to obtain multiple spatiotemporal data blocks.

3. The method according to claim 1, characterized in that, The location sequence in the spatiotemporal data block is mapped to a high-dimensional feature space by embedding spatial mapping. In the high-dimensional feature space, the distribution positions of the current spatiotemporal data block and the spatiotemporal data blocks corresponding to historical normal navigation are compared in the topological space, and candidate abnormal time periods are identified, including: For the position sequence in the spatiotemporal data block, a neighborhood matrix is ​​constructed and a reconstruction weight is calculated using the local linear embedding algorithm. Based on the reconstruction weight, the position sequence is projected onto a high-dimensional feature space to obtain the current embedding vector. Based on the distribution of the current embedding vector in the high-dimensional feature space, a metric tensor is constructed using the Riemannian geometry algorithm. The spatiotemporal data blocks corresponding to historical normal navigation are obtained, and the position sequences in the historical spatiotemporal data blocks are projected onto a high-dimensional feature space to obtain a set of historical embedding vectors. Based on the metric tensor, the topological distance between each historical embedding vector in the set of historical embedding vectors and its nearest neighbor is calculated by an adaptive density clustering algorithm, and the local density radius is determined by combining the preset quantile position. With the historical embedding vector as the center, the set of hypersphere surface points is constructed by combining the local density radius to obtain a set of normal distribution boundary points. The topological distance between the current embedded vector and each boundary point in the set of normal distribution boundary points is calculated based on the metric tensor, and the minimum value is extracted as the deviation metric. Anomaly score is obtained based on the statistical characteristics of the topological distance between the deviation metric and the set of historical embedded vectors, and the corresponding time window is determined as a candidate anomaly period based on the anomaly score.

4. The method according to claim 3, characterized in that, Based on the metric tensor, the topological distance between each historical embedding vector in the historical embedding vector set and its nearest neighbor historical embedding vector is calculated using an adaptive density clustering algorithm, and the local density radius is determined by combining this with the preset quantile positions, including: Based on the metric tensor, the topological distance is obtained by calculating the geodesic path integral between different historical embedding vectors in the historical embedding vector set. For each historical embedding vector, the historical embedding vector corresponding to the minimum topological distance is extracted as the nearest neighbor historical embedding vector. The topological distance from each historical embedding vector to the nearest neighbor historical embedding vector is recorded to obtain the nearest neighbor distance set. The topological distances in the nearest neighbor distance set are arranged in ascending order to obtain an ordered distance sequence. Based on the preset quantile position, the corresponding index position is located in the ordered distance sequence and the topological distance at the index position is extracted as the global quantile distance. The neighborhood density value is obtained by calculating the number of historical embedding vectors contained in the hypersphere with the global quantile distance as the radius for each historical embedding vector in the historical embedding vector set. The local density radius is obtained by weighted fusion of the nearest neighbor topological distance corresponding to each historical embedding vector and the reciprocal of the neighborhood density value.

5. The method according to claim 1, characterized in that, Extracting the motion state sequence and environmental information sequence corresponding to the candidate abnormal time period and reconstructing the phase space in the time dimension, calculating the mutual information transmission strength and transmission direction between the motion state sequence and the environmental information sequence in the reconstructed phase space, and constructing a causal influence matrix includes: Based on the start and end times of the candidate abnormal time period, the motion state sequence and environmental information sequence within the range from the start time to the end time are extracted from the spatiotemporal data block. For each subsequence in the motion state sequence and the environment information sequence, the optimal delay time is calculated using a time delay embedding algorithm and a pre-set mutual information function, the optimal embedding dimension is calculated based on the false nearest neighbor algorithm, and the phase space is reconstructed based on the optimal delay time and the optimal embedding dimension to obtain the motion state phase space vector set and the environment information phase space vector set. The forward information transfer amount from the environmental information variables in the environmental information phase space vector set to the motion state variables in the motion state phase space vector set, and the reverse information transfer amount from the motion state variables in the motion state phase space vector set to the environmental information variables in the environmental information phase space vector set are calculated using the transfer entropy algorithm. A bidirectional transfer entropy matrix is ​​then constructed based on the forward and reverse information transfer amounts. The transfer entropy values ​​in the bidirectional transfer entropy matrix are statistically tested based on the significance test algorithm to determine the effective transfer intensity. Based on the directional marker information of the effective transfer intensity, the starting variable and target variable of the transfer are used to construct a causal influence matrix with motion state variables and environmental information variables as nodes and effective transfer intensity and transfer direction as directed edges.

6. The method according to claim 1, characterized in that, The causal influence matrix is ​​subjected to spectral decomposition to obtain dominant eigenvalues ​​and determine the degree of dominance of environmental changes on changes in motion state. Based on the degree of dominance, the candidate abnormal time periods are judged to obtain a confirmed abnormal event set, including: The causal influence matrix is ​​subjected to spectral decomposition to extract all eigenvalues ​​and corresponding eigenvectors. The eigenvalues ​​are sorted in descending order according to their magnitudes, and the first eigenvalue is taken as the dominant eigenvalue. The eigenvectors corresponding to the dominant eigenvalues ​​are extracted as the dominant eigenvectors. Calculate the absolute value of each element in the dominant feature vector and determine the element with the largest absolute value. Obtain the index position of the element with the largest absolute value in the dominant feature vector and map it to the corresponding variable node in the causal influence matrix. Determine the variable type corresponding to the variable node and determine the dominant direction of causal influence based on the variable type. Based on the dominant feature value, the intensity of the dominant effect is calculated by the spectral radius measurement method and compared with the intensity of the dominant effect with a preset threshold. The degree of the dominant effect of environmental change on motion state change is determined according to the comparison result. Based on the degree of dominance and the direction of dominance, the candidate abnormal time period is classified into different types. If the degree of dominance exceeds a preset threshold and the direction of dominance is dominated by environmental information variables, the candidate abnormal time period is classified into an environment-driven abnormal event. Otherwise, the candidate abnormal time period is classified into an autonomous behavior abnormal event. The determination is repeated to obtain the confirmed abnormal event set.

7. The method according to claim 1, characterized in that, The confirmed abnormal event set is traversed in chronological order to identify persistent and sudden abnormalities. Based on the duration of persistent abnormalities and the degree of deviation of sudden abnormalities, tiered early warning information is generated and output, including: All abnormal events in the confirmed abnormal event set are sorted in ascending order by timestamp and traversed. The time interval between the start time and end time of each abnormal event is determined as the duration. The numerical change features of the motion state sequence within the abnormal event period are extracted. The temporal change rate of the motion state sequence within the period of the abnormal event is calculated by the sliding window method. Abnormal events with a temporal change rate lower than a preset temporal change rate threshold and a duration exceeding a preset duration threshold are marked as continuous abnormal events; otherwise, the abnormal events are marked as sudden abnormal events. The evolution rate corresponding to the persistent anomaly is calculated by a trend fitting method, and the duration level is determined by combining the duration. The peak deviation and change gradient corresponding to the sudden anomaly are extracted and the deviation level is determined. Based on the duration level and the deviation level, an early warning level is generated and combined with a pre-determined anomaly type and a pre-acquired dominant direction to obtain graded early warning information and output it.

8. A data fusion-based system for identifying and warning of anomalies in vessel trajectories in lake areas, used to implement the method described in any one of claims 1-7, characterized in that, include: The data fusion module is used to acquire multi-source trajectory data and perform spatiotemporal alignment processing by combining it with pre-set timestamp association and spatial coordinate mapping rules to obtain fused trajectory data. The fused trajectory data is divided into multiple time windows and a spatiotemporal data block corresponding to each time window is constructed. An anomaly identification module is used to map the position sequence in the spatiotemporal data block to a high-dimensional feature space by embedding spatial mapping, and compare the distribution position of the current spatiotemporal data block with the spatiotemporal data block corresponding to historical normal navigation in the topological space in the high-dimensional feature space to identify candidate abnormal time periods. The causal analysis module is used to extract the motion state sequence and environmental information sequence corresponding to the candidate abnormal time period and reconstruct the phase space in the time dimension. In the reconstructed phase space, the mutual information transmission strength and transmission direction between the motion state sequence and the environmental information sequence are calculated and a causal influence matrix is ​​constructed. The causal influence matrix is ​​subjected to spectral decomposition to obtain the dominant feature value and to determine the degree of dominance of environmental changes on motion state changes. Based on the degree of dominance, the candidate abnormal time period is judged to obtain a set of confirmed abnormal events. The anomaly warning module is used to traverse the confirmed abnormal event set in chronological order and identify persistent anomalies and sudden anomalies. It generates and outputs graded warning information based on the duration of the persistent anomaly and the degree of deviation of the sudden anomaly.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images