Traffic network data modeling method and system based on graph neural network

Abstract

The application discloses a traffic network data modeling method and system based on a graph neural network, and belongs to the field of intelligent transportation systems.The method comprises the following steps: 1, splitting into fast and slow subgraphs according to road types and modal activity, assigning a space-time scale to each subgraph, and setting an asymmetric fissure operator for the boundary; 2, converting meteorological events into a negative volume permeation field irrelevant to topology, calculating the phase difference between the negative volume permeation field and the fast and slow subgraphs according to the space-time scale, and taking the ratio as an opening prediction coefficient; and 3, projecting trajectory data to an equipotential surface in a reverse direction according to the opening prediction coefficient to generate slip anchor points.The application can split fast and slow subgraphs according to road types and modal activity, set an asymmetric fissure operator, and convert meteorological events into a negative volume permeation field, thereby solving the problem that existing models ignore the heterogeneity of a road network and the space-time mismatch of multi-modal data, accurately capturing the dynamic evolution law of traffic flow, and reducing the traffic state prediction error in a complex scene.

Description

Technical Field

[0001] This invention relates to the field of intelligent transportation systems, and more specifically, to a method and system for modeling traffic network data based on graph neural networks. Background Technology

[0002] With the rapid development of intelligent transportation, traffic network data modeling has become the management foundation for traffic condition prediction and congestion management. Graph neural networks, due to their powerful spatial correlation modeling capabilities, have been widely used in the field of traffic network data modeling. Most existing modeling models use a unified graph structure to model all the correlations in the traffic network, failing to fully consider the complex heterogeneity within traffic data and the interaction characteristics between multimodal traffic information. This leads to deviations in the model's understanding of complex traffic conditions and makes it difficult to accurately capture the dynamic evolution of traffic flow.

[0003] The transportation network itself is significantly heterogeneous, with fundamental differences in the physical dynamic characteristics of different types of roads: highways have high vehicle speeds, fast traffic congestion propagation speeds, and weak spatial correlation decay; ordinary roads have low vehicle speeds, slow congestion propagation speeds, and are more significantly affected by local events such as intersections and pedestrians. At the same time, the multimodal data involved in the modeling of the transportation network, such as trajectory, meteorology, and events, have a natural mismatch in spatiotemporal scales: trajectory data provides discrete and sparse vehicle position sampling information, meteorological data presents continuous and slowly changing field distribution characteristics, and event data are instantaneous and sudden markers that are unrelated to the road network topology.

[0004] In the transitional zones between urban and suburban areas where highways and ordinary roads frequently intersect, and in complex scenarios involving both sparse trajectory data and sudden localized rainstorm events, a one-sided pursuit of graph structure uniformity and computational efficiency can lead to the simplification of the dynamic coupling relationship between heterogeneous roads into a uniform propagation model. At the same time, the slow infiltration effect of weather events and the instantaneous response of trajectory data can cause a phase mismatch, resulting in a sharp amplification of traffic speed prediction errors near weather change points. For example, within 5 minutes after the start of a rainstorm, the predicted speed on highways is more than 20% higher, while that on ordinary roads is more than 15% lower, seriously affecting the scientific nature and effectiveness of traffic control decisions. Summary of the Invention

[0005] To address the problems existing in the prior art, the present invention aims to provide a traffic network data modeling method and system based on graph neural networks. This method can solve the problems of neglecting the heterogeneity of the road network and the spatiotemporal mismatch of multimodal data in existing models by splitting fast and slow subgraphs according to road type and modal activity, setting asymmetric crack operators, and combining negative volume infiltration field to transform meteorological events. It can accurately capture the dynamic evolution law of traffic flow and reduce the prediction error of traffic state in complex scenarios.

[0006] To solve the above problems, the present invention adopts the following technical solution:

[0007] The first aspect is a traffic network data modeling method based on graph neural networks, including:

[0008] Step 1: Divide the graph into fast and slow subgraphs according to road type and modal activity. Assign a spatiotemporal scale to each subgraph and set an asymmetric fracture operator at the boundary.

[0009] Step 2: The meteorological event is converted into a negative volume infiltration field that is independent of topology. The phase difference between the field and the fast and slow subgraphs is calculated using a spatiotemporal scale, and the ratio is used as the opening prediction coefficient.

[0010] Step 3: Project the trajectory data back onto the equipotential surface using the opening prediction coefficient to generate sliding anchor points. The offset is determined by the local gradient of the phase difference.

[0011] Step 4: Construct a hinge matrix at the boundary nodes using a set of sliding anchor points, with the element table showing the maximum sliding amount; define the interlocking coefficient by the permeation field gradient, and when the gradient exceeds a preset gradient threshold, the interlocking coefficient approaches zero, and the coupling in the corresponding region becomes intermittent jump transmission;

[0012] Step 5: Using the sliding anchor point as an elastic constraint, construct a pre-distortion space with the hinge matrix. Each road segment is a time-varying second-order tensor, and the deformation rate is controlled by the spatiotemporal scale and the biting coefficient. Apply a no-reference correction to the tensor field and force the sum of the sub-traces to be zero.

[0013] Step 6: Using the tensor field as the stress source and the hinge matrix as the boundary condition, solve the Poisson equation for the road network to obtain the stress and strain, and mark the stress concentration area as an uncertainty mask.

[0014] Further, step 1 includes:

[0015] Step 11: For each road segment, the physical type features and modal activity are mapped to the boundary ambiguity coefficient through boundary entropy. When the value is greater than the threshold, it is marked as a transition zone. The rest are divided into fast and slow sub-graph core areas. The mapping rule is: the ratio of features to activity is compressed into a base value through nonlinear compression, and then modulated by the local gradient of adjacent boundaries.

[0016] Step 12: Only the road sections within the core area are assigned a spatiotemporal scale, with the scale scale being proportional to the product of the reciprocal of the design speed and the modal lag time. The road sections in the transition zone are not assigned a scale and are used as the sole attachment objects of the asymmetric crack operator. The operator is a directional gating matrix, with rows corresponding to the side boundary nodes of the fast subgraph and columns corresponding to the side boundary nodes of the slow subgraph. The elements are -1, 0, and +1, which correspond to slow to fast jumps, no coupling, and fast to slow jumps, respectively. The three-state values ​​are determined by the average ambiguity coefficient within the transition zone. The ambiguity coefficient is monotonically positively correlated with the proportion of 0, and the sign of non-zero values ​​is determined by the positive or negative sign of the ambiguity gradient direction.

[0017] Further, step 2 includes:

[0018] Step 21: Discretize the meteorological event into a topology-free grid, assign an initial negative scalar field value to each grid point, and define the negative volume as the integral of the negative scalar field values ​​of all grid points; propagate the field according to the anisotropic permeability tensor, the principal direction of the permeability tensor is jointly modulated by the terrain slope and road density, and the secondary direction is an orthogonal rotation of the principal direction, generating a negative volume permeability field that is independent of the road network topology.

[0019] Step 22: Relying on the negative volume permeation field and the spatiotemporal scales of the fast and slow subgraphs, the permeation field values ​​are inversely transformed along the time axis using the Laplace transform to obtain the Green's function response spectrum under the spatiotemporal scale of each subgraph; the extreme point time of the rate of change of the permeation field is extracted from the response spectrum, and time offset is performed with the extreme point time of the rate of change of the modal activity in Step 1 to obtain the phase difference field of the spatial distribution of each subgraph.

[0020] Furthermore, step 2 also includes:

[0021] Step 23: Relying on the phase difference field of the fast and slow subgraphs, smooth the two anisotropically along the direction indicated by the asymmetric fracture operator, then perform logarithmic compression on the smoothed ratio and retain the original sign to obtain the initial prediction coefficient field; multiply the prediction coefficient field of the previous time step by the attenuation factor and then superimpose them to form a fracture opening prediction coefficient field with hysteresis.

[0022] Step 24: Relying on the fracture opening prediction coefficient field, use it as a modulation signal to act in reverse on the negative volume permeation field: when the prediction coefficient is greater than 1, the local diffusion coefficient is multiplied by the reciprocal of the coefficient, and when it is less than 1, it is multiplied by the square of the coefficient; re-input the modulated permeation field into step 22, and repeat steps 22 to 24 until the rate of change of the prediction coefficient field is lower than the threshold, and output the stable fracture opening prediction coefficient field.

[0023] Further, step 3 includes:

[0024] Step 31: Relying on the fracture opening prediction coefficient field and the spatiotemporal scale of the fast and slow sub-maps, determine the equipotential surface level of each point in the trajectory data in the negative volume permeation field, calculate the local gradient vector of the phase difference field of the sub-map, multiply the gradient vector by the projection direction deflection angle of the prediction coefficient field at that point, and project the trajectory point onto the equipotential surface along the deflection angle to generate a temporary offset vector.

[0025] Step 32: Relying on the temporary offset vector and the direction of the asymmetric crack operator, if the included angle is less than 90 degrees, the offset vector is retained and multiplied by the scale value as the anchor point displacement; otherwise, it is applied after rotating 180 degrees around the equipotential surface normal. The gated anchor point displacement is applied to the original trajectory point to generate a sliding anchor point, and the relative displacement between adjacent anchor points is forced not to exceed the preset neighborhood consistency tolerance.

[0026] Further, step 4 includes:

[0027] Step 41: Based on the set of sliding anchor points and the set of boundary nodes, calculate the relative displacement vector magnitude between the two sliding anchor points for each pair of boundary nodes, subtract the lower limit of the neighborhood consistency tolerance value as the initial allowable sliding margin, and arrange them in matrix form to obtain the original hinge matrix.

[0028] Step 42: Relying on the original hinge matrix and the crack opening prediction coefficient field, the prediction coefficient field value is used as the modulation factor: when it is greater than 1, it is divided by the square of the modulation factor; when it is less than 1, it is multiplied by the reciprocal of the modulation factor; when it is equal to 1, it remains unchanged, and the dynamic hinge matrix is ​​output.

[0029] Step 43: Relying on the negative volume permeation field, extract the local gradient magnitude at the boundary node and output the engagement coefficient through the nonlinear response function: the engagement coefficient is 1 when the local gradient magnitude is lower than the first threshold, decays exponentially to less than 1% of the initial value when it is between the first and second thresholds, and is 0 when it is higher than the second threshold; multiply it element-wise with the dynamic hinge matrix to obtain the hinge matrix with engagement control.

[0030] Furthermore, step 4 also includes:

[0031] Step 44: Relying on the hinge matrix with engagement control, set an accumulated slip counter and a jump trigger flag for each pair of boundary nodes. Accumulate the actual relative displacement increment in each virtual time step. When the counter value reaches or exceeds the hinge matrix element value, set the trigger flag to 1 and retain the excess amount. If the element value is 0, force it to be set to 0 and clear the counter.

[0032] Step 45: Relying on the jump trigger flag and the direction of the asymmetric crack operator, when the flag is 1, perform a one-way exchange of the difference of the elastic tensor trace according to the operator direction. After the exchange, the counter is deducted from the corresponding hinge matrix element value and the flag is reset. The region with a biting coefficient of 0 is forcibly prohibited from triggering until the coefficient returns to a positive number.

[0033] Further, step 5 includes:

[0034] Step 51: Relying on the set of sliding anchor points and the hinge matrix with biting control, the anchor point displacement vector is decomposed into tangential and normal components as boundary conditions; the initial zero second-order tensor of each road segment is defined as follows: tangential strain is positively correlated with the tangential anchor point displacement gradient, and normal strain is positively correlated with the normal anchor point displacement gradient; the deformation rate is controlled by using the reciprocal of the spatiotemporal scale as the relaxation rate and the biting coefficient as the boundary stiffness: the tangential deformation rate is the ratio of tangential strain to relaxation rate multiplied by the boundary stiffness, and the normal deformation rate is the ratio of normal strain to relaxation rate divided by the boundary stiffness, thus obtaining the time-varying second-order tensor field;

[0035] Step 52: Relying on the time-varying second-order tensor field, sum all tensor traces in each subgraph. If the sum is not zero, calculate the correction amount for each tensor. The correction amount is proportional to the deviation between the trace value and the average trace value of the subgraph multiplied by a weighting factor that is inversely proportional to the current value of the cumulative slip counter. Subtract the correction amount from each tensor to make the sum of traces in the subgraph zero. If a tensor trace value crosses the zero point, scale its tangential and normal components proportionally to maintain the original sign. Output the corrected elastic tensor field.

[0036] Further, step 6 includes:

[0037] Step 61: Relying on the corrected elastic tensor field and dynamic hinge matrix, perform an asymmetric rotation transformation on the tensor field of each road segment, extract the first derivative along the road direction and the second derivative in the vertical direction to form an equivalent body force distribution, and obtain the stress source term.

[0038] Step 62: Relying on the equivalent body force distribution and dynamic hinge matrix, the non-zero elements in the dynamic hinge matrix are used as crack constraints to construct Lagrange multipliers. The Poisson equation is transformed into a linear complementary problem with inequality boundary conditions. The problem is solved by successive over-relaxation iterations with a fixed number of iterations equal to twice the road network diameter. The road network diameter is defined as the maximum value of the shortest path between all node pairs. The stress and strain field is then output.

[0039] Step 63: Based on the stress-strain field and the region with zero interlocking coefficient in step 45, calculate the strain energy density of each road segment, obtain the stress concentration index by second-order difference between the strain energy density of the road segment and the adjacent road segment, mark the road segment with the index exceeding three times the global median as a high uncertainty mask, and automatically include the region with zero interlocking coefficient in the mask, and output the mask label.

[0040] Secondly, the present invention also provides a traffic network data modeling system based on graph neural networks, comprising:

[0041] The split map scaling module is used to split the map into fast and slow sub-maps according to road type and modal activity. Each sub-map is assigned a spatiotemporal scale and the boundary is set with an asymmetric split operator.

[0042] The infiltration field phase judgment module is used to convert meteorological events into negative volume infiltration fields that are independent of topology. It calculates the phase difference between the field and the fast and slow subgraphs using a spatiotemporal scale, and the ratio is used as the opening prediction coefficient.

[0043] The reverse projection anchor point module is used to project the trajectory data onto the equipotential surface with the opening prediction coefficient to generate a sliding anchor point. The offset is determined by the local gradient of the phase difference.

[0044] The hinge jump module is used to construct a hinge matrix at the boundary nodes using a set of sliding anchor points. The element table shows the maximum sliding amount. The engagement coefficient is defined by the permeation field gradient. When the gradient exceeds a preset gradient threshold, the engagement coefficient approaches zero, and the coupling in the corresponding region becomes intermittent jump transmission.

[0045] The tensor zeroing module is used to construct a pre-distortion space with the sliding anchor point as an elastic constraint and the hinge matrix. Each road segment is a time-varying second-order tensor, and the deformation rate is controlled by the spatiotemporal scale and the biting coefficient. The tensor field is subjected to a referenceless correction, which forces the sum of the sub-traces to be zero.

[0046] The stress mask module is used to solve the Poisson equation of the road network with the tensor field as the stress source and the hinge matrix as the boundary condition to obtain the stress and strain, and to mark the stress concentration area as an uncertainty mask.

[0047] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0048] (1) This scheme solves the problem of existing models ignoring road network heterogeneity and spatiotemporal mismatch of multimodal data by splitting fast and slow subgraphs according to road type and modal activity, setting asymmetric crack operators, and combining negative volume infiltration field to transform meteorological events. It accurately captures the dynamic evolution law of traffic flow and reduces the prediction error of traffic state in complex scenarios.

[0049] (2) This scheme generates sliding anchor points by back-projecting the trajectory data through the opening prediction coefficient, and controls the offset by the phase difference local gradient, which effectively solves the defect of discrete and sparse trajectory data, makes the trajectory data adaptable to meteorological events and sub-graph dynamic characteristics, provides a reliable position benchmark for boundary coupling control, and improves the utilization and adaptability of modeling data.

[0050] (3) This scheme constructs a hinge matrix by sliding anchor points and defines the biting coefficient by combining the permeation field gradient, so as to realize the dynamic control of the coupling state of the subgraph boundary, adapt to the intermittent jump transmission requirements when the weather events change drastically, avoid the deviation caused by the simplification of the boundary coupling relationship, and improve the model's adaptability to complex traffic scenarios.

[0051] (4) This scheme uses time-varying second-order tensor field construction and reference-free correction, combined with the solution of road network Poisson equation and stress concentration area marking, to accurately identify high uncertainty areas, provide clear risk guidance for traffic control decision-making, improve the scientificity and effectiveness of traffic control decision-making, and help congestion management and accurate prediction of traffic conditions. Attached Figure Description

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

[0053] Figure 1 This is a flowchart illustrating the overall method of the present invention;

[0054] Figure 2 This is a schematic diagram of the fast and slow subgraph splitting and asymmetric crack operator of the present invention;

[0055] Figure 3 This is a schematic diagram illustrating the generation of the negative volume permeation field and sliding anchor point in this invention.

[0056] Figure 4 This is a schematic diagram of the data flow between the various modules of the present invention. Detailed Implementation

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

[0058] Example 1

[0059] Please see Figures 1 to 3 A traffic network data modeling method based on graph neural networks is proposed. This method includes: Step 1, splitting the data into fast and slow subgraphs according to road type and modal activity, assigning a spatiotemporal scale to each subgraph, and setting an asymmetric fracture operator at the boundary. The specific operations are as follows:

[0060] First, by combining the differences in road types and modal activity characteristics of each road segment in the road network, the overall traffic network map is decomposed into two independent subgraph structures: a fast subgraph and a slow subgraph, which are related by boundaries. The fast subgraph mainly corresponds to road types with high traffic efficiency and weak spatial correlation decay characteristics, while the slow subgraph mainly corresponds to road types with low traffic efficiency and strong local event influence characteristics. After the decomposition, a spatiotemporal scale is assigned to each subgraph. This scale is used to quantify the differences in spatiotemporal characteristics between different subgraphs, achieving accurate matching of spatiotemporal scales in subsequent modeling. At the same time, an asymmetric gap operator is set in the boundary region between the fast and slow subgraphs. This operator is used to regulate the coupling relationship between the two subgraphs, adapt to the dynamic evolution characteristics of different subgraphs, avoid the simplification of coupling relationships caused by a unified graph structure, and ensure that the interaction between subgraphs can truly reflect the actual dynamics of the road network.

[0061] Step 1 further includes the following steps:

[0062] Step 11: For each road segment, map the physical type features and modal activity using boundary entropy to obtain the boundary ambiguity coefficient. If the coefficient exceeds a threshold, it is marked as a transition zone; the rest are divided into fast and slow sub-map core regions. The mapping rule is: the feature-to-activity ratio is non-linearly compressed to a base value, and then modulated using local gradients at adjacent boundaries. The specific operations are as follows:

[0063] For each road segment in the road network, two key parameters need to be extracted: physical type features and modal activity. Physical type features reflect the inherent traffic attributes of the road segment, while modal activity reflects the response intensity of multimodal traffic data on the road segment. Using a boundary entropy mapping method, these two parameters are transformed into boundary ambiguity coefficients that characterize the boundary features of the road segment. This mapping process follows specific rules: first, the ratio of physical type features to modal activity is calculated; this ratio is then subjected to nonlinear compression to obtain a mapping base value; finally, the base value is modulated using the local gradients of adjacent road segment boundaries. To eliminate biases caused by a single parameter, the boundary ambiguity coefficient of each road segment is obtained. Then, a preset threshold is set, which is determined based on the actual road network scenario and modeling accuracy requirements. Road segments with boundary ambiguity coefficients greater than this threshold are marked as transition zones. Due to the ambiguity of their boundary characteristics, these road segments cannot be clearly classified into fast or slow subgraphs. Road segments with boundary ambiguity coefficients not greater than this threshold are split into fast subgraph core areas and slow subgraph core areas. During the splitting process, the physical type characteristics of the road segments and the specific values ​​of modal activity are combined to ensure that the road segments in the core areas have similar dynamic characteristics.

[0064] Step 12: Only the road sections within the core area are assigned a spatiotemporal scale. The scale scale is proportional to the product of the reciprocal of the design speed and the modal lag time. The road sections in the transition zone are not assigned a scale and are the only attachment objects of the asymmetric crack operator. The operator is a directional gating matrix, with rows corresponding to the side boundary nodes of the fast subgraph and columns corresponding to the side boundary nodes of the slow subgraph. The elements are -1, 0, and +1, which correspond to slow to fast jumps, no coupling, and fast to slow jumps, respectively. The three-state values ​​are determined by the average ambiguity coefficient within the transition zone. The ambiguity coefficient is monotonically positively correlated with 0. The sign of non-zero values ​​is determined by the positive or negative sign of the ambiguity gradient direction. The specific operation is as follows:

[0065] The assignment of the spatiotemporal scale has a clear range limitation. The road segments in the core areas of the fast subgraph and the core areas of the slow subgraph determined in step 11 are assigned values. The road segments in the transition zone are not assigned a spatiotemporal scale, and the transition zone is the only area that can be attached to the asymmetric crack operator. This design ensures the coupling and control function of the transition zone as the boundary of the subgraph. The scale setting of the spatiotemporal scale follows a strict quantification relationship. Its scale size is proportional to the product of the reciprocal of the road segment's design speed and the modal lag time. The design speed is an inherent parameter of the road segment, which determines the basic traffic efficiency of the road segment. The modal lag time reflects the response delay of multimodal traffic data in the road segment. The product of the two can accurately quantify the spatiotemporal characteristics of different core area road segments, so that the time and space calculations in the subsequent modeling process can be adapted to the actual dynamics of each subgraph.

[0066] The essence of the asymmetric fracture operator is a directional gating matrix. The rows of this matrix correspond to the boundary nodes on the fast subgraph side, and the columns correspond to the boundary nodes on the slow subgraph side. The matrix elements take only three values: -1, 0, and +1, each corresponding to a different subgraph coupling state: -1 represents a sudden jump coupling from the slow subgraph to the fast subgraph, meaning the dynamic characteristics of the slow subgraph are transferred to the fast subgraph; 0 represents no coupling between the fast and slow subgraphs, and they evolve independently; +1 represents a sudden jump coupling from the fast subgraph to the slow subgraph, meaning the dynamic characteristics of the fast subgraph are transferred to the slow subgraph. The graph is transmitted; the three-state values ​​of the matrix elements are determined by the average ambiguity coefficient within the transition zone, where the ambiguity coefficient is monotonically positively correlated with the proportion of 0, that is, the larger the average ambiguity coefficient, the higher the proportion of non-zero elements in the matrix, and the more active the coupling relationship between subgraphs; the sign of the non-zero elements is determined by the direction of the ambiguity gradient within the transition zone. When the ambiguity gradient is positive, the non-zero elements take a positive sign, and when the ambiguity gradient is negative, the non-zero elements take a negative sign. In this way, the dynamic control of the coupling relationship between the fast and slow subgraph boundaries is achieved, ensuring that the coupling state matches the actual boundary characteristics of the transition zone.

[0067] In a preferred embodiment of the present invention, step 2 is included: the meteorological event is converted into a topology-independent negative volume permeation field, and its phase difference with the fast and slow subgraphs is calculated using a spatiotemporal scale. The ratio is used as the opening prediction coefficient. The specific operation is as follows:

[0068] First, the meteorological event is transformed into a negative volume permeation field independent of the road network topology. This eliminates the interference of the road network topology on the propagation characteristics of the meteorological event, ensuring that the propagation law of the meteorological event can truly reflect its own physical evolution characteristics. Then, using the spatiotemporal scales assigned to the fast and slow subgraphs in step 1, the spatiotemporal phase difference between the negative volume permeation field and the two subgraphs is quantified. This phase difference can accurately reflect the time offset characteristics of the dynamic evolution of the meteorological event and the subgraphs. Finally, the phase difference between the two subgraphs is calculated as a ratio, which is the opening prediction coefficient used to predict the degree of coupling opening of the fast and slow subgraph boundaries.

[0069] Step 2 further includes the following steps:

[0070] Step 21: Discretize the meteorological event into a topology-free grid, assign an initial negative scalar field value to each grid point, and define the negative volume as the integral of the negative scalar field values ​​of all grid points; propagate this field according to the anisotropic permeability tensor, where the principal direction of the permeability tensor is jointly modulated by the terrain slope and road density, and the secondary direction is an orthogonal rotation of the principal direction, generating a negative volume permeability field independent of the road network topology. The specific operations are as follows:

[0071] The area covered by the meteorological event is discretized to construct a topologically unrelated grid. This grid is independent of the node and road segment distribution of the transportation network and is used solely to accurately characterize the spatial distribution features of the meteorological event. Each grid point is assigned an initial negative scalar field value, the magnitude of which is determined by the intensity of the meteorological event; the greater the intensity of the meteorological event, the larger the absolute value of the initial negative scalar field value at the corresponding grid point. Next, a negative volume parameter is defined, which is the integral of the initial negative scalar field values ​​of all grid points. Through integration, the discrete grid point field values ​​are transformed into a representation of the overall intensity of the meteorological event. The quantitative indicators are used to achieve a comprehensive characterization of the intensity of meteorological events. Then, the initial negative scalar field is propagated according to the anisotropic permeability tensor. The anisotropic permeability tensor is used to simulate the differences in the propagation of meteorological events in different directions. Its principal direction is jointly modulated by the terrain slope and road density. The terrain slope determines the natural propagation tendency of meteorological events, while the road density affects the propagation speed and range of meteorological events in the road network area. The two work together to determine the principal direction of the permeability tensor. The secondary direction is an orthogonal rotation of the principal direction to ensure that the permeability tensor can fully cover the propagation characteristics in all directions of space.

[0072] The local diffusion coefficient used here for infiltration field propagation is a parameter of infiltration field evolution. It is obtained by fitting historical observation data in combination with meteorological event type, terrain features, and road network density. The specific value needs to match the actual scenario: if the meteorological event is a severe convective weather such as heavy rain or heavy snow, the local diffusion coefficient ranges from 0.8 to 1.2; if it is a weak weather event such as light rain or light wind, the value ranges from 0.3 to 0.7; for areas with plains and high road network density, the diffusion coefficient takes the upper limit of the corresponding interval; for areas with mountainous terrain and low road network density, the lower limit of the corresponding interval is taken. The magnitude of the diffusion coefficient directly determines the propagation speed of the infiltration field. The larger the coefficient, the faster the propagation speed. Through the above anisotropic propagation process, the initial negative scalar field gradually evolves and finally generates a negative volume infiltration field that is independent of the road network topology. This infiltration field can truly reflect the spatial propagation law and intensity distribution of meteorological events.

[0073] Step 22: Relying on the negative volume permeation field and the spatiotemporal scales of the fast and slow subgraphs, the permeation field values ​​are inversely Laplaced along the time axis to obtain the Green's function response spectrum under the spatiotemporal scale of each subgraph; the extreme point of the rate of change of the permeation field is extracted from the response spectrum, and time-shifted to the extreme point of the rate of change of the modal activity in Step 1 to obtain the spatial distribution phase difference field of each subgraph. The specific operation is as follows:

[0074] Using the negative volume permeation field generated in step 21 and the spatiotemporal scales assigned to the fast and slow subplots in step 1 as core dependent parameters, the field value of the negative volume permeation field is first subjected to an inverse Laplace transform along the time axis. The purpose of the inverse Laplace transform is to convert the frequency domain characteristics of the permeation field into time domain response characteristics, thereby obtaining the Green's function response spectrum under the spatiotemporal scale of each subplot. Since the spatiotemporal scales of the fast and slow subplots are different, their corresponding Green's function response spectra are also different. This response spectrum can accurately reflect the dynamic response law of the negative volume permeation field under different spatiotemporal scales of the subplots. When extracting the extreme points of the rate of change of the permeation field from the Green's function response spectrum... The moment corresponds to the time point when the infiltration field intensity changes most drastically, which can characterize the key time node of the impact of meteorological events on the road network. At the same time, the extreme point of the modal activity change rate in step 1 is extracted. This moment corresponds to the key node of the subgraph dynamic evolution and reflects the peak moment of the subgraph's response to multimodal data. By calculating the time offset between these two extreme point moments, the phase difference value corresponding to each subgraph is obtained. The phase difference values ​​of all subgraphs are integrated according to the spatial distribution law to finally form the spatial distribution phase difference field of each subgraph. This phase difference field can quantify the degree of spatiotemporal misalignment between the negative volume infiltration field and the dynamic evolution of each subgraph.

[0075] Step 23: Relying on the phase difference fields of the fast and slow subgraphs, smooth them anisotropically along the direction indicated by the asymmetric fracture operator. Then, logarithmically compress the smoothed ratio while retaining the original sign to obtain the initial prediction coefficient field. Multiply the prediction coefficient field of the previous time step by the attenuation factor and then superimpose them to form a fracture opening prediction coefficient field with hysteresis. The specific operation is as follows:

[0076] Based on the fast and slow subgraph phase difference fields obtained in step 22, anisotropic smoothing is performed along the direction indicated by the asymmetric fracture operator set in step 1. The direction of anisotropic smoothing is consistent with the directionality of the asymmetric fracture operator. The purpose is to eliminate local noise in the phase difference field, making the spatial distribution of the phase difference field smoother, while preserving the phase difference characteristics of the subgraph boundary region. After smoothing, the ratio of the fast and slow subgraph phase difference fields is calculated, and this ratio is logarithmically compressed. The purpose of logarithmic compression is to adjust the numerical range of the ratio to a reasonable range, avoiding the interference of extreme values ​​on subsequent calculations, while preserving the original sign of the ratio to ensure that the positive and negative characteristics of the phase difference can be accurately transmitted. Through the above processing, an initial prediction coefficient field is obtained, which can... The initial analysis reflects the widening tendency of the fast and slow subgraph boundaries. Subsequently, the prediction coefficient field from the previous time step is introduced and multiplied by an attenuation factor. The attenuation factor is obtained by combining the dynamic change frequency of the road network and the modal data update cycle, and is calibrated through historical data. The value range is fixed at 0.1 to 0.5, and the specific value is adjusted according to the dynamic characteristics of the road network: when the road network traffic flow fluctuates greatly and the modal data is updated frequently, such as once every 1 minute, the attenuation factor is 0.1 to 0.2 to reduce the influence of historical prediction coefficients and ensure that the prediction coefficients can quickly respond to real-time changes; when the road network traffic flow is stable and the modal data update cycle is long, such as once every 5 minutes, the attenuation factor is 0.3 to 0.5 to enhance the influence of historical prediction coefficients and improve the stability of the prediction results.

[0077] The attenuation factor always takes a value between 0 and 1, and will not be greater than 1 or less than 0. This avoids excessive interference or no impact from historical data on the current prediction coefficient. The attenuated prediction coefficient field of the previous time step is superimposed with the current initial prediction coefficient field to finally form a hysteresis-induced crack opening prediction coefficient field. The hysteresis characteristic can prevent the prediction coefficient from abruptly changing due to instantaneous phase changes, ensuring the stability of the prediction results and conforming to the continuous characteristics of the dynamic evolution of the traffic network.

[0078] Step 24: Relying on the fracture opening prediction coefficient field, use it as a modulation signal to act inversely on the negative volume permeation field: when the prediction coefficient is greater than 1, multiply the local diffusion coefficient by the reciprocal of the coefficient; when it is less than 1, multiply it by the square of the coefficient. Re-input the modulated permeation field into step 22, and repeat steps 22 to 24 until the rate of change of the prediction coefficient field is lower than the threshold. Output the stable fracture opening prediction coefficient field. The specific operation is as follows:

[0079] The crack opening prediction coefficient field obtained in step 23 is used as the modulation signal and acts in the opposite direction on the negative volume permeation field generated in step 21. The modulation rule is determined according to the magnitude of the prediction coefficient: when the prediction coefficient is greater than 1, the local diffusion coefficient of the negative volume permeation field is multiplied by the reciprocal of the prediction coefficient to slow down the diffusion speed of the permeation field and avoid excessive boundary opening due to an excessively large prediction coefficient; when the prediction coefficient is less than 1, the local diffusion coefficient is multiplied by the square of the prediction coefficient to accelerate the local diffusion of the permeation field and compensate for insufficient boundary opening due to an excessively small prediction coefficient; when the prediction coefficient is equal to 1, the local diffusion coefficient remains unchanged and the permeation field propagates normally. The local diffusion coefficient used here is the same parameter as the local diffusion coefficient used for the initial propagation of the permeation field in step 21, and the value standard is completely consistent. It does not need to be obtained again. Dynamic modulation can be performed only according to the magnitude of the prediction coefficient.

[0080] After modulation, the new negative volume permeation field is re-inputted into step 22, and the entire process from step 22 to step 24 is repeated for iterative iteration. During each iteration, the phase difference field and the prediction coefficient field are recalculated based on the new permeation field, and the deviation of the prediction coefficient is gradually corrected. The iteration terminates when the rate of change of the prediction coefficient field is lower than a preset threshold. This threshold is determined according to the modeling accuracy requirements and is usually between 0.001 and 0.01. The higher the modeling accuracy requirements, the smaller the threshold value. When the rate of change is lower than the threshold, it indicates that the prediction coefficient field has reached a stable state and no longer changes significantly with the iteration process. At this time, the stable fracture opening prediction coefficient field is output.

[0081] In a preferred embodiment of the present invention, step 3 is included, in which the trajectory data is back-projected onto an equipotential surface using the opening prediction coefficient to generate a sliding anchor point, the offset of which is determined by the local gradient of the phase difference. The specific operation is as follows:

[0082] Using the stable fracture opening prediction coefficient output in step 2 as the core control parameter, the discrete trajectory data is projected inversely onto the equipotential surface of the negative volume infiltration field. This eliminates the bias caused by the sparsity and discreteness of the trajectory data, enabling the trajectory data to accurately reflect the dynamic characteristics of the road network under the influence of meteorological events. During the projection process, the offset of the sliding anchor point is not a fixed value, but is determined by the local gradient of the phase difference field obtained in step 2. The local gradient can quantify the spatial rate of change of the phase difference, thereby accurately controlling the magnitude and direction of the offset, ensuring that the generated sliding anchor point can truly match the dynamic evolution law of the sub-map and the infiltration influence of meteorological events.

[0083] Step 3 further includes the following steps:

[0084] Step 31: Relying on the fracture opening prediction coefficient field and the spatiotemporal scales of the fast and slow sub-maps, determine the equipotential surface level of each point in the trajectory data in the negative volume permeation field, calculate the local gradient vector of the phase difference field of the sub-map, multiply the gradient vector by the projection direction deflection angle of the prediction coefficient field at that point, and project the trajectory point onto the equipotential surface along this deflection angle to generate a temporary offset vector. The specific operations are as follows:

[0085] For each trajectory point in the trajectory data, its equipotential surface level in the negative volume permeation field needs to be determined first. An equipotential surface is a spatial surface formed by points with equal field values ​​in the negative volume permeation field. The equipotential surface level is determined based on the magnitude of the permeation field value; trajectory points with similar field values ​​are grouped into the same level. This division process must consider the overall field value distribution of the negative volume permeation field to ensure that the level division accurately reflects the spatial intensity differences of the permeation field. After determining the equipotential surface level, the local gradient vector of the phase difference field in the subgraph where the trajectory point is located is calculated. The calculation of the local gradient vector is based on the phase difference field data around the trajectory point, obtained by performing a difference operation on the phase difference values ​​of adjacent grid points, which accurately reflects the spatial direction and rate of phase difference change at that location. The calculated local gradient vector is then multiplied by the value of the fracture opening prediction coefficient field at that trajectory point. The result of the multiplication is used to determine the projection direction deflection angle. The core logic of this operation is to quantify the boundary opening tendency using the prediction coefficient, and by combining it with the gradient vector, calibrate the projection direction so that the projection can adapt to the coupling state of the subgraph boundary. The multiplication formula is: ;

[0086] This formula is derived from the geometric definition of the vector dot product. The result of the vector dot product is equal to the product of the magnitudes of the two vectors and the cosine of the angle between the two vectors. By transformation, the cosine of the projection direction deflection angle θ can be obtained, and thus the magnitude of the deflection angle can be determined. In the formula, It is the local gradient vector of the phase difference field. The cosine of the spatial angle between the scalar k and the crack opening prediction coefficient is directly derived from the geometric definition of vector dot product and is used to quantify the degree of directional agreement between the two vectors in space. The local gradient vector representing the phase difference field characterizes the spatial direction and rate of phase difference change at this trajectory point; k represents the value of the fracture opening prediction coefficient field at this trajectory point, characterizing the boundary opening tendency at this location. θ is the magnitude of the local gradient vector, |k| is the absolute value of the prediction coefficient, and θ is the projection direction deflection angle. After determining the deflection angle, the trajectory points are projected in the opposite direction onto the corresponding equipotential surface along the deflection angle. During the projection process, it is necessary to ensure the spatial correlation between the projected position of the trajectory points on the equipotential surface and the original trajectory points. After the projection is completed, a temporary offset vector is generated. The magnitude of this vector is equal to the spatial distance between the original trajectory points and the projected points, and the direction is from the original trajectory points to the projected points. This vector is used for further processing of the subsequent anchor point displacement.

[0087] Step 32: Relying on the temporary offset vector and the direction of the asymmetric fracture operator, if the included angle is less than 90 degrees, retain the offset vector multiplied by the scale value as the anchor point displacement; otherwise, apply it after rotating 180 degrees around the equipotential surface normal. Apply the gated anchor point displacement to the original trajectory point to generate a sliding anchor point, and force the relative displacement between adjacent anchor points not to exceed the preset neighborhood consistency tolerance. The specific operation is as follows:

[0088] The angle between the temporary offset vector and the direction of the asymmetric fracture operator is calculated. This angle calculation is based on the geometric principle of vector angles, determined by the ratio of the dot product of the two vectors to the product of their magnitudes. The temporary offset vector is then gated based on this angle: if the angle is less than 90 degrees, it indicates that the direction of the temporary offset vector matches the direction indicated by the asymmetric fracture operator, reflecting the coupling tendency of the subgraph boundary. In this case, the temporary offset vector is retained, and it is multiplied by the spatiotemporal scale value of the corresponding subgraph core area to obtain the anchor point displacement. The scale value is the same as in step 12 and does not need to be obtained again. Its function is to calibrate the temporary offset vector to the spatiotemporal scale of the subgraph to ensure that the anchor point displacement matches the dynamic characteristics of the subgraph. If the included angle is not less than 90 degrees, it means that the direction of the temporary offset vector is contrary to the indication direction of the asymmetric fracture operator. In this case, the temporary offset vector needs to be rotated 180 degrees around the normal of the equipotential surface and its direction needs to be adjusted to match the direction of the operator. Then, it is multiplied by the spatiotemporal scale value to obtain the anchor point displacement. During the rotation process, the magnitude of the temporary offset vector needs to remain unchanged, and only the direction needs to be adjusted.

[0089] After gating, the obtained anchor point displacement is applied to the original trajectory points. The final spatial position of the sliding anchor point is obtained through spatial coordinate superposition calculation. Each original trajectory point corresponds to a sliding anchor point, achieving accurate calibration of the trajectory data. At the same time, a preset neighborhood consistency tolerance needs to be set. This tolerance is obtained by combining the trajectory data sampling frequency, road network segment length, and modeling accuracy requirements, and is calibrated through historical trajectory data. The value ranges from 0.5 to 2 meters, and the specific value is adjusted according to the actual scenario: when the trajectory sampling frequency is high, such as once every 10 seconds, the road segment length is short, and the modeling accuracy requirement is high, the tolerance is 0.5 to 1 meter; when the sampling frequency is low, the road segment length is long, and the modeling accuracy requirement is moderate, the tolerance is 1 to 2 meters. The neighborhood consistency tolerance is used to constrain the relative displacement of adjacent sliding anchor points, forcing the spatial relative displacement between adjacent anchor points to not exceed the tolerance. If it exceeds the tolerance, the position of adjacent anchor points is fine-tuned to ensure that the spatial distribution of sliding anchor points is continuous and uniform, conforming to the actual topological characteristics and dynamic evolution law of the road network.

[0090] In a preferred embodiment of the present invention, step 4 is included: constructing a hinge matrix at the boundary nodes using a set of sliding anchor points, with the element table representing the maximum slip; defining the engagement coefficient by the permeation field gradient, where the engagement coefficient approaches zero when the gradient exceeds a preset gradient threshold, and the corresponding region coupling becomes intermittent jump transmission. The specific operation is as follows:

[0091] At the boundary nodes of the fast and slow subgraphs, a hinge matrix is ​​constructed using the set of sliding anchor points generated in step 3. Each element in the matrix is ​​used to characterize the maximum amount of sliding between the corresponding pair of boundary nodes, providing a quantitative constraint for boundary sliding. The biting coefficient is defined based on the gradient change of the negative volume permeation field. The biting coefficient is used to characterize the tightness of coupling at the boundary nodes. When the gradient of the permeation field exceeds the preset gradient threshold, the biting coefficient approaches zero. At this time, the subgraph coupling state of the corresponding region changes to intermittent jump transmission, which adapts to the boundary dynamic characteristics when meteorological events change drastically, ensuring that the boundary coupling can accurately respond to the dynamic interaction between meteorological events and the subgraph.

[0092] Step 4 further includes the following steps:

[0093] Step 41: Relying on the set of sliding anchor points and the set of boundary nodes, calculate the relative displacement vector magnitude between the two sliding anchor points for each pair of boundary nodes, subtract the lower limit of the neighborhood consistency tolerance value as the initial allowable sliding margin, and arrange them into matrix form to obtain the original hinge matrix. The specific operations are as follows:

[0094] First, the composition of the boundary node set is defined. Boundary nodes are all nodes at the boundary between the fast and slow subgraphs. Each pair of boundary nodes corresponds to one boundary node on the fast subgraph side and one boundary node on the slow subgraph side. The matching rule for node pairs is consistent with the row-column correspondence of the asymmetric fracture operator in step 1, ensuring that the matrix structure is adapted to the boundary coupling direction. For each pair of boundary nodes, the corresponding sliding anchor points on both sides need to be extracted, and the relative displacement vector magnitude between these two sliding anchor points needs to be calculated. The calculation of the relative displacement vector magnitude is based on the spatial coordinate difference and is obtained through the vector magnitude formula. It is used to quantify the spatial offset of the two anchor points. Then, subtract... The preset lower limit of neighborhood consistency tolerance is used to obtain the initial allowable slip margin. The lower limit of neighborhood consistency tolerance corresponds to the range of neighborhood consistency tolerance in step 32 and is fixed at 0.5 meters. Its function is to reserve reasonable slip space and avoid excessive slip constraint caused by small fluctuations in anchor points. The initial allowable slip margins of all boundary node pairs are arranged in matrix form. The rows of the matrix correspond to the boundary nodes on the fast subgraph side and the columns correspond to the boundary nodes on the slow subgraph side. The value of each element in the matrix is ​​the initial allowable slip margin of the corresponding boundary node pair. Thus, the original hinge matrix is ​​obtained. This matrix initially quantifies the maximum slip capability of each boundary node pair.

[0095] Step 42: Relying on the original hinge matrix and the crack opening prediction coefficient field, the prediction coefficient field value is used as the modulation factor: if it is greater than 1, it is divided by the square of the modulation factor; if it is less than 1, it is multiplied by the reciprocal of the modulation factor; if it is equal to 1, it remains unchanged. The dynamic hinge matrix is ​​then output. The specific operation is as follows:

[0096] The modulation process uses the value of the crack opening prediction coefficient field as the modulation factor. Different modulation rules are adopted according to different values ​​of the prediction coefficient to ensure that the hinge matrix can accurately adapt to the boundary opening tendency. When the crack opening prediction coefficient is greater than 1, it indicates that the boundary tends to open excessively. In this case, the corresponding element in the original hinge matrix needs to be attenuated by dividing the element by the square of the modulation factor to reduce the maximum slip and avoid excessive boundary slip. When the crack opening prediction coefficient is less than 1, it indicates that the boundary opening is insufficient. In this case, the corresponding element in the original hinge matrix needs to be enhanced. The element is multiplied by the reciprocal of the modulation factor to increase the maximum slip and improve the flexibility of boundary slip. When the crack opening prediction coefficient is equal to 1, it means that the boundary opening state is within a reasonable range. At this time, the corresponding element in the original hinge matrix remains unchanged and no modulation is required. The modulation process requires matching the prediction coefficient value of each element in the original hinge matrix one by one and performing the corresponding modulation operation. After all elements are modulated, a dynamic hinge matrix is ​​output. This matrix can dynamically adjust the maximum slip according to the boundary opening prediction state, making the boundary coupling control more targeted.

[0097] Step 43: Relying on the negative volume permeation field, extract the local gradient magnitude at the boundary nodes, and output the engagement coefficient through the nonlinear response function: the engagement coefficient is 1 when the local gradient magnitude is below the first threshold, decays exponentially to less than 1% of the initial value when it is between the first and second thresholds, and is 0 when it is above the second threshold; multiply it element-wise with the dynamic hinge matrix to obtain the hinge matrix with engagement control, the specific operation is as follows:

[0098] The local gradient magnitude of the negative volume infiltration field at the boundary node is extracted. The calculation method of the local gradient magnitude is the same as that of the local gradient vector of the phase difference field in step 31. It is obtained by performing a difference operation on the infiltration field values ​​of adjacent grid points around the boundary node, and is used to quantify the spatial rate of change of meteorological event infiltration intensity at the boundary node. The local gradient magnitude is input into the nonlinear response function, and the interlocking coefficient is output through this function. The design of the nonlinear response function is based on the logic of the influence of meteorological event infiltration intensity on boundary coupling. The more drastic the change in infiltration intensity, the smaller the interlocking coefficient and the looser the boundary coupling. The piecewise calculation formula of the interlocking coefficient is as follows: ;

[0099] This formula is derived from the infiltration characteristics of meteorological events and the boundary coupling law. When the infiltration field gradient is small, the boundary coupling is tight, and the interlocking coefficient is 1; as the gradient increases, the coupling gradually loosens, and the interlocking coefficient decreases exponentially; when the gradient is too large, the coupling is completely released, and the interlocking coefficient is 0; in the formula, The interlocking coefficient represents the degree of coupling at the boundary nodes; The local gradient magnitude represents the negative volume permeation field at the boundary node; The first threshold, The second threshold is used, and both are obtained by calibrating historical permeation field data with boundary coupling effects. The value ranges from 0.2 to 0.5. The value ranges from 1.0 to 1.5; the higher the required modeling accuracy, the smaller the threshold value. The attenuation coefficient is fixed at a value of 0.8 to 1.2 and is used to control the rate of exponential decay. After obtaining the engagement coefficient, it is multiplied element by element with the dynamic hinge matrix. Each matrix element is multiplied by the engagement coefficient of the corresponding boundary node to obtain a hinge matrix with engagement control. This matrix combines the maximum slip amount with the tightness of coupling to achieve precise control of boundary coupling.

[0100] Step 44: Relying on the hinge matrix with engagement control, set an accumulated slip counter and a jump trigger flag for each pair of boundary nodes. Accumulate the actual relative displacement increment every virtual time step. When the counter value reaches or exceeds the hinge matrix element value, set the trigger flag to 1 and retain the excess amount. If the element value is 0, force it to be set to 0 and clear the counter. The specific operation is as follows:

[0101] For each pair of boundary nodes, a cumulative slip counter and a jump trigger flag are set. The cumulative slip counter records the sum of the actual relative displacement increments of the boundary node pair within each virtual time step. The jump trigger flag indicates whether the jump trigger condition has been met. Both flags and counters are initially set to 0. The virtual time step is determined by combining the road network dynamic update cycle and modeling accuracy, ranging from 1 to 5 seconds. When road network traffic flow fluctuates greatly and high modeling accuracy is required, 1 to 2 seconds are used; when traffic flow is stable and modeling accuracy requirements are moderate, 3 to 5 seconds are used. Within each virtual time step, the actual relative displacement increment of the boundary node pair is calculated. The value is accumulated in the cumulative slip counter to achieve real-time accumulation of slip amount. When the value of the cumulative slip counter reaches or exceeds the value of the corresponding element in the hinge matrix with engagement control, it indicates that the boundary slip amount has reached the maximum allowable value. At this time, the jump trigger flag is set to 1, and the part of the counter value that exceeds the matrix element value is retained as the initial accumulation amount for the next virtual time step. If the value of the corresponding element in the hinge matrix is ​​0, it means that the boundary node pair is not allowed to slip. At this time, the cumulative slip counter is forcibly set to 0, and the jump trigger flag is reset to 0 to ensure that the boundary slip does not exceed the constraint range, thereby achieving accurate monitoring and trigger control of slip amount.

[0102] Step 45: Relying on the sudden jump trigger flag and the direction of the asymmetric crack operator, when the flag is 1, perform a unidirectional exchange of the difference of the elastic tensor traces according to the operator direction. After the exchange, the counter is decremented by the corresponding hinge matrix element value and the flag is reset. The region with a bite coefficient of 0 is forcibly prohibited from triggering until the coefficient returns to a positive value. The specific operation is as follows:

[0103] Determine the status of the jump trigger flag for each pair of boundary nodes. When the flag is 1, it means that the pair of boundary nodes has reached the maximum slip and a one-way exchange operation needs to be performed in the direction of the asymmetric crack operator. The exchanged content is the difference between the elastic tensor traces at the corresponding boundary nodes of the fast and slow subgraphs. The exchange process strictly follows the directionality of the operator. If the operator element is +1, a one-way exchange from the fast subgraph to the slow subgraph is performed; if the operator element is -1, a one-way exchange from the slow subgraph to the fast subgraph is performed; if the operator element is 0, no exchange operation is performed, ensuring that the exchange direction is consistent with the boundary coupling tendency. After the exchange is completed, the value of the accumulated slip counter is subtracted from the value of the corresponding element in the hinge matrix with biting control, deducting the slip amount that has been triggered to jump, and the jump trigger flag is reset to 0, preparing for the slip accumulation and trigger judgment of the next virtual time step. At the same time, for the region with a biting coefficient of 0, the boundary coupling of the region has been completely released, and the jump trigger operation must be forcibly prohibited until the biting coefficient returns to a positive number. At this time, the boundary coupling is reactivated, and the jump trigger judgment can be performed again to avoid invalid jumps in the state of decoupling, ensuring the rationality and stability of the boundary coupling control and adapting to the boundary dynamic characteristics when meteorological events change drastically.

[0104] In a preferred embodiment of the present invention, step 5 is included: using the sliding anchor point as an elastic constraint, a pre-distortion space is constructed with the hinge matrix, each road segment is a time-varying second-order tensor, and the deformation rate is controlled by the spatiotemporal scale and the engagement coefficient; a referenceless correction is applied to the tensor field to force the sum of the sub-traces to be zero. The specific operation is as follows:

[0105] Using the sliding anchor points generated in step 3 as elastic constraints, and combined with the hinge matrix with biting control obtained in step 4, a pre-distortion space is constructed. This space can simulate the subtle distortion state of the road network under the coupling effect of meteorological events and subgraphs, which fits the dynamic evolution characteristics of the actual road network. Each road segment in the road network is defined as a time-varying second-order tensor. The deformation rate of the tensor is jointly controlled by the spatiotemporal scale given in step 1 and the biting coefficient defined in step 4, ensuring that the deformation rate can adapt to the spatiotemporal characteristics and boundary coupling tightness of different road segments. A no-reference correction operation is applied to the generated tensor field. By forcing the sum of the traces of all tensors in each subgraph to be zero, the systematic bias in the tensor field is eliminated, ensuring that the tensor field can accurately reflect the real deformation state of the road segment.

[0106] Step 5 further includes the following steps:

[0107] Step 51: Relying on the set of sliding anchor points and the hinge matrix with engagement control, the anchor point displacement vector is decomposed into tangential and normal components as boundary conditions. The initial zero-second-order tensor for each road segment is defined as follows: tangential strain is positively correlated with the tangential anchor point displacement gradient, and normal strain is positively correlated with the normal anchor point displacement gradient. Using the reciprocal of the spatiotemporal scale as the relaxation rate and the engagement coefficient as the boundary stiffness, the deformation rate is controlled: the tangential deformation rate is the ratio of tangential strain to relaxation rate multiplied by the boundary stiffness, and the normal deformation rate is the ratio of normal strain to relaxation rate divided by the boundary stiffness, resulting in a time-varying second-order tensor field. The specific operations are as follows:

[0108] The displacement vector of each sliding anchor point is decomposed into tangential and normal components. The tangential component is parallel to the extension direction of the corresponding road segment, and the normal component is perpendicular to the extension direction of the road segment. The decomposition process uses a vector projection method, which projects the displacement vector onto the tangential and normal unit vectors of the road segment. The two decomposed components together serve as the boundary conditions for constructing the pre-distortion space, ensuring that the boundary constraints can accurately adapt to the spatial orientation of the road segment. Next, an initial zero-second-order tensor is defined for each road segment in the road network. In the initial state, all components of the tensor are zero. Subsequently, the tensor will be adjusted according to the road segment... The deformation characteristics evolve gradually. The strain characteristics of the initial zero-second tensor are defined. Tangential strain is positively correlated with the tangential anchor point displacement gradient. The larger the tangential anchor point displacement gradient, the larger the tangential strain, which represents the more obvious deformation of the road segment along the extension direction. Normal strain is positively correlated with the normal anchor point displacement gradient. The larger the normal anchor point displacement gradient, the larger the normal strain, which represents the more obvious deformation of the road segment perpendicular to the extension direction. The positive correlation coefficient here is obtained by calibrating historical deformation data and modeling accuracy, and the value ranges from 0.8 to 1.2 to ensure the matching degree between strain and displacement gradient.

[0109] The relaxation rate is represented by the reciprocal of the spatiotemporal scale. The larger the spatiotemporal scale, the smaller the relaxation rate, indicating a slower response speed to road segment deformation. The bite coefficient is represented by the boundary stiffness. The larger the bite coefficient, the greater the boundary stiffness, indicating a tighter boundary coupling and stronger constraint on road segment deformation. The deformation rate of the tensor is controlled by the above parameters. The tangential deformation rate is the ratio of tangential strain to relaxation rate multiplied by the boundary stiffness. This calculation logic allows the tangential deformation to adapt to the spatiotemporal characteristics of the road segment while being constrained by the boundary coupling state. The normal deformation rate is the ratio of normal strain to relaxation rate divided by the boundary stiffness, taking into account both the flexibility of normal deformation and the influence of boundary constraints, avoiding excessive normal deformation. After the deformation rate of all road segments is calculated, a time-varying second-order tensor field is obtained by combining the time evolution process. This tensor field can reflect the deformation state of each road segment at different times in real time.

[0110] Step 52: Relying on the time-varying second-order tensor field, sum all tensor traces within each subgraph. If the sum is not zero, calculate the correction for each tensor. The correction is proportional to the deviation between the trace value and the average trace value of the subgraph multiplied by a weighting factor inversely proportional to the current value of the cumulative slip counter. Subtract the correction for each tensor until the sum of traces within the subgraph is zero. If a tensor trace value crosses zero, scale its tangential and normal components proportionally to maintain its original sign. Output the corrected elastic tensor field. The specific operations are as follows:

[0111] For each subgraph, calculate and sum the traces of all time-varying second-order tensors within that subgraph. The trace of a tensor is the sum of its main diagonal elements, which quantifies the overall deformation intensity of the tensor. The summation operation must cover the tensors of all road segments within the subgraph to ensure that the calculation result reflects the overall deviation of the tensor field of the subgraph. If the summation result is not zero, it indicates that there is a systematic deviation in the tensor field of the subgraph, and the correction amount for each tensor needs to be calculated. The calculation of the correction amount follows a specific logic, being directly proportional to the deviation of the trace value of that tensor and the average deviation of the trace values ​​of all tensors within the subgraph, and inversely proportional to the current value of the cumulative slip counter. The formula for calculating the correction amount is as follows: ;

[0112] The formula is derived from the requirement of no-reference correction, aiming to achieve precise correction by quantifying the deviation and combining it with the slip accumulation state, ensuring that the sum of the sub-traces is zero while adapting to the boundary slip characteristics; in the formula, The correction amount represents a single tensor; This is a correction factor, obtained through calibration using historical correction data, with a value ranging from 0.1 to 0.3. It is used to control the magnitude of the correction and avoid over- or under-correction. The trace value represents a single tensor; The value represents the average of all tensor traces within the subgraph; C represents the current value of the cumulative slip counter. The larger the cumulative slip, the smaller the correction weight, to avoid over-correction affecting the slip coupling effect. After calculating the correction amount for each tensor, each tensor is substituted into the correction one by one, that is, each tensor is subtracted from the corresponding correction amount, until the sum of the traces of all tensors within the subgraph is zero. If the trace value of a tensor crosses the zero point during the correction process, that is, the trace value changes from positive to negative or from negative to positive, the tangential and normal components of the tensor need to be scaled proportionally. During the scaling process, the ratio of the two components remains unchanged to ensure that the deformation direction of the tensor is consistent with the original sign, and to avoid distortion of the deformation direction due to correction. After all correction operations are completed, the corrected elastic tensor field is output. This tensor field eliminates system bias and can accurately reflect the true elastic deformation state of the road segment.

[0113] In a preferred embodiment of the present invention, step 6 is included: using the tensor field as the stress source and the hinge matrix as the boundary condition, solving the Poisson equation of the road network to obtain the stress and strain, and marking the stress concentration region as an uncertainty mask. The specific operation is as follows:

[0114] The corrected elastic tensor field output in step 5 is used as the stress source. This tensor field has eliminated system bias and can accurately reflect the true elastic deformation state of the road segment. The dynamic hinge matrix output in step 4 is used as the boundary condition to ensure that the solution process can adapt to the coupling control characteristics of the subgraph boundary. By solving the Poisson equation of the road network, the stress-strain field of the entire road network is obtained. This field can quantify the stress and strain distribution of each road segment and reflect the stress and deformation state of the road network under the coupling effect of meteorological events and subgraphs. Subsequently, stress concentration areas are identified based on the stress-strain field and these areas are marked as uncertainty masks. At the same time, the areas with a zero interlocking coefficient in step 4 are automatically included in the mask range. The coupling state of these areas is unstable and the prediction bias risk is high. By marking them with masks, clear risk guidance can be provided for subsequent traffic prediction and control, improving the scientificity and effectiveness of decision-making.

[0115] Step 6 further includes the following steps:

[0116] Step 61: Relying on the corrected elastic tensor field and dynamic hinge matrix, perform an asymmetric rotation transformation on the tensor field of each road segment, extract the first derivative along the road direction and the second derivative in the vertical direction to form an equivalent body force distribution, and obtain the stress source term. The specific operation is as follows:

[0117] An asymmetric rotation transformation is performed on the corrected elastic tensor field of each road segment. The direction of the rotation transformation is consistent with the extension direction of the road segment and the directionality of the asymmetric fracture operator. The purpose is to adjust the components of the tensor field to match the actual stress direction of the road network and eliminate stress calculation deviations caused by the mismatch between the tensor field direction and the road network orientation. After the rotation transformation, the first derivative of the tensor field along the road direction and the second derivative perpendicular to the road direction are extracted. The first derivative along the road direction can quantify the rate of change of the tensor field along the extension direction of the road segment, reflecting the longitudinal stress change of the road segment. The second derivative perpendicular to the road direction can quantify the degree of curvature of the tensor field in the lateral direction, reflecting the lateral stress difference of the road segment. These two derivatives are combined to form an equivalent body force distribution. The equivalent body force distribution can transform the deformation characteristics of the tensor field into an equivalent load of the road network stress, accurately characterizing the stress intensity at each location of the road network. Through the above processing, the stress source term is finally obtained. This stress source term is the core input for solving the Poisson equation of the road network, which can ensure that the solution of the equation closely matches the actual stress and deformation state of the road network.

[0118] Step 62: Relying on the equivalent body force distribution and dynamic hinge matrix, the non-zero elements in the dynamic hinge matrix are used as crack constraints to construct Lagrange multipliers. The Poisson equation is transformed into a linear complementarity problem with inequality boundary conditions. This problem is solved through successive over-relaxation iterations with a fixed number of iterations equal to twice the road network diameter, where the road network diameter is defined as the maximum value of the shortest path between all node pairs. The stress-strain field is then output. The specific operations are as follows:

[0119] Using a dynamic hinge matrix as boundary constraints, the non-zero elements in the dynamic hinge matrix represent the allowed slip coupling between corresponding boundary node pairs. These non-zero elements are used as gap constraints to construct Lagrange multipliers. The role of the Lagrange multipliers is to integrate the boundary constraints into the Poisson equation, ensuring that the equation solution process meets the boundary coupling requirements. By introducing Lagrange multipliers, the original road network Poisson equation is transformed into a linear complementarity problem with inequality boundary conditions. This transformation logic is based on the principle of constraint optimization, which can combine boundary constraints with equation solution and avoid the solution result from exceeding the boundary constraint range. Subsequently, a successive over-relaxation iterative method is used to solve the linear complementarity problem. The number of iterations is fixed at twice the network diameter, where the network diameter is defined as the maximum value of the shortest path between all pairs of nodes. The calculation method uses the shortest path algorithm in graph theory, traversing all pairs of nodes, calculating the shortest path length between each pair, and taking the maximum value as the network diameter. The value of the network diameter is determined according to the network size; for small-scale networks, it is typically 5 to 10, and for large-scale networks, it is typically 10 to 20. The iteration count of twice the network diameter ensures convergence of the iteration results. During the iteration process, the over-relaxation factor is calibrated using historical solution data, with a value range of 0.8 to 1.2, used to control the iteration convergence speed. A larger factor results in faster convergence, but excessively large factors should be avoided to prevent iteration oscillations. After the iteration is completed, the stress-strain field of the entire network is output, which accurately reflects the stress magnitude, direction, and strain state of each road segment.

[0120] Step 63: Based on the stress-strain field and the region with zero interlocking coefficient in Step 45, calculate the strain energy density of each road segment. Obtain the stress concentration index by second-order difference between the strain energy density of each segment and that of adjacent road segments. Mark road segments with an index exceeding three times the global median as high uncertainty masks. Automatically include regions with zero interlocking coefficient in the mask and output the mask labels. The specific operations are as follows:

[0121] Calculate the strain energy density for each road segment. Strain energy density is a physical quantity characterizing the energy stored in the deformation of a road segment, and can indirectly reflect the degree of stress concentration in the road segment. Its calculation formula is as follows: ;

[0122] This formula is derived from the definition of strain energy in elasticity. Strain energy is equal to half the product of stress and strain. The strain energy per unit volume is the strain energy density, which can accurately quantify the local stress and energy storage state of a road section. In the formula, w represents the strain energy density. The stress vector representing the road segment is taken from the stress-strain field; The strain vector representing the road segment is also taken from the stress-strain field. The dot product of the two is used to quantify the coupling effect of stress and strain. After calculating the strain energy density of all road segments, the second difference between the strain energy density of each road segment and its adjacent road segments is calculated. The second difference can amplify the local abrupt change in strain energy density. The more obvious the abrupt change, the more drastic the stress change in the area and the higher the stress concentration. Thus, the stress concentration index of each road segment is obtained. The global median of the stress concentration index of all road segments is calculated. The global median is calculated by arranging all stress concentration indices in ascending order and taking the value of the middle position. If the number of indices is even, the average of the two middle values ​​is taken. Road segments with stress concentration indices exceeding three times the global median are marked as high-uncertainty masks. These road segments exhibit significant stress concentration and have a high risk of traffic condition prediction deviation. Simultaneously, areas with a zero interlocking coefficient from step 45 are automatically included in the mask range. These areas have completely decoupled boundaries, exhibit unstable dynamic characteristics, and also have a high risk of prediction deviation. After all markings are completed, the final mask markings are output, providing clear risk area guidance for subsequent traffic condition predictions and ensuring that control decisions can specifically avoid prediction deviations in high-risk areas.

[0123] Example 2

[0124] Please see Figure 4 Based on Example 1, this example provides a traffic network data modeling system based on graph neural networks, including:

[0125] The split map scaling module is used to split the map into fast and slow sub-maps according to road type and modal activity. Each sub-map is assigned a spatiotemporal scale and the boundary is set with an asymmetric split operator.

[0126] The infiltration field phase judgment module is used to convert meteorological events into topology-independent negative volume infiltration fields, calculate their phase difference with the fast and slow subgraphs using a spatiotemporal scale, and use the ratio as the opening prediction coefficient.

[0127] The reverse projection anchor point module is used to project the trajectory data onto the equipotential surface with the opening prediction coefficient to generate a sliding anchor point. The offset is determined by the local gradient of the phase difference.

[0128] The hinge jump module is used to construct a hinge matrix at the boundary nodes using a set of sliding anchor points. The element table shows the maximum sliding amount. The engagement coefficient is defined by the permeation field gradient. When the gradient exceeds a preset gradient threshold, the engagement coefficient approaches zero, and the coupling in the corresponding region becomes intermittent jump transmission.

[0129] The tensor zeroing module is used to construct a pre-distortion space with the sliding anchor point as an elastic constraint and the hinge matrix. Each road segment is a time-varying second-order tensor, and the deformation rate is controlled by the spatiotemporal scale and the biting coefficient. The tensor field is subjected to a referenceless correction, which forces the sum of the sub-traces to be zero.

[0130] The stress mask module is used to solve the Poisson equation of the road network with the tensor field as the stress source and the hinge matrix as the boundary condition to obtain the stress and strain, and to mark the stress concentration area as an uncertainty mask.

[0131] The above description is merely a preferred embodiment of the present invention; however, the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and its improved concepts, should be covered within the scope of protection of the present invention.

Claims

1. A traffic network data modeling method based on graph neural networks, characterized in that, include: Step 1: Divide the graph into fast and slow subgraphs according to road type and modal activity. Assign a spatiotemporal scale to each subgraph and set an asymmetric fracture operator at the boundary. Step 2: The meteorological event is converted into a negative volume infiltration field that is independent of topology. The phase difference between the field and the fast and slow subgraphs is calculated using a spatiotemporal scale, and the ratio is used as the opening prediction coefficient. Step 3: Project the trajectory data back onto the equipotential surface using the opening prediction coefficient to generate sliding anchor points. The offset is determined by the local gradient of the phase difference. Step 4: Construct a hinge matrix at the boundary nodes using a set of sliding anchor points, with the element table showing the maximum sliding amount; define the interlocking coefficient by the permeation field gradient, and when the gradient exceeds a preset gradient threshold, the interlocking coefficient approaches zero, and the coupling in the corresponding region becomes intermittent jump transmission; Step 5: Using the sliding anchor point as an elastic constraint, construct a pre-distortion space with the hinge matrix. Each road segment is a time-varying second-order tensor, and the deformation rate is controlled by the spatiotemporal scale and the biting coefficient. Apply a no-reference correction to the tensor field and force the sum of the sub-traces to be zero. Step 6: Using the tensor field as the stress source and the hinge matrix as the boundary condition, solve the Poisson equation for the road network to obtain the stress and strain, and mark the stress concentration area as an uncertainty mask.

2. The traffic network data modeling method based on graph neural networks according to claim 1, characterized in that, Step 1 includes: Step 11: For each road segment, the physical type features and modal activity are mapped to the boundary ambiguity coefficient through boundary entropy. When the value is greater than the threshold, it is marked as a transition zone. The rest are divided into fast and slow sub-graph core areas. The mapping rule is: the ratio of features to activity is compressed into a base value through nonlinear compression, and then modulated by the local gradient of adjacent boundaries. Step 12: Only the road sections within the core area are assigned a spatiotemporal scale, with the scale scale being proportional to the product of the reciprocal of the design speed and the modal lag time. The road sections in the transition zone are not assigned a scale and are used as the sole attachment objects of the asymmetric crack operator. The operator is a directional gating matrix, with rows corresponding to the side boundary nodes of the fast subgraph and columns corresponding to the side boundary nodes of the slow subgraph. The elements are -1, 0, and +1, which correspond to slow to fast jumps, no coupling, and fast to slow jumps, respectively. The three-state values ​​are determined by the average ambiguity coefficient within the transition zone. The ambiguity coefficient is monotonically positively correlated with the proportion of 0, and the sign of non-zero values ​​is determined by the positive or negative sign of the ambiguity gradient direction.

3. The traffic network data modeling method based on graph neural networks according to claim 2, characterized in that, Step 2 includes: Step 21: Discretize the meteorological event into a topology-free grid, assign an initial negative scalar field value to each grid point, and define the negative volume as the integral of the negative scalar field values ​​of all grid points; propagate the field according to the anisotropic permeability tensor, the principal direction of the permeability tensor is jointly modulated by the terrain slope and road density, and the secondary direction is an orthogonal rotation of the principal direction, generating a negative volume permeability field that is independent of the road network topology. Step 22: Relying on the negative volume permeation field and the spatiotemporal scales of the fast and slow subgraphs, the permeation field values ​​are inversely transformed along the time axis using the Laplace transform to obtain the Green's function response spectrum under the spatiotemporal scale of each subgraph; the extreme point time of the rate of change of the permeation field is extracted from the response spectrum, and time offset is performed with the extreme point time of the rate of change of the modal activity in Step 1 to obtain the phase difference field of the spatial distribution of each subgraph.

4. The traffic network data modeling method based on graph neural networks according to claim 3, characterized in that, Step 2 also includes: Step 23: Relying on the phase difference field of the fast and slow subgraphs, smooth the two anisotropically along the direction indicated by the asymmetric fracture operator, then perform logarithmic compression on the smoothed ratio and retain the original sign to obtain the initial prediction coefficient field; multiply the prediction coefficient field of the previous time step by the attenuation factor and then superimpose them to form a fracture opening prediction coefficient field with hysteresis. Step 24: Relying on the fracture opening prediction coefficient field, use it as a modulation signal to act in reverse on the negative volume permeation field: when the prediction coefficient is greater than 1, the local diffusion coefficient is multiplied by the reciprocal of the coefficient, and when it is less than 1, it is multiplied by the square of the coefficient; re-input the modulated permeation field into step 22, and repeat steps 22 to 24 until the rate of change of the prediction coefficient field is lower than the threshold, and output the stable fracture opening prediction coefficient field.

5. The traffic network data modeling method based on graph neural networks according to claim 4, characterized in that, Step 3 includes: Step 31: Relying on the fracture opening prediction coefficient field and the spatiotemporal scale of the fast and slow sub-maps, determine the equipotential surface level of each point in the trajectory data in the negative volume permeation field, calculate the local gradient vector of the phase difference field of the sub-map, multiply the gradient vector by the projection direction deflection angle of the prediction coefficient field at that point, and project the trajectory point onto the equipotential surface along the deflection angle to generate a temporary offset vector. Step 32: Relying on the temporary offset vector and the direction of the asymmetric crack operator, if the included angle is less than 90 degrees, the offset vector is retained and multiplied by the scale value as the anchor point displacement; otherwise, it is applied after rotating 180 degrees around the equipotential surface normal. The gated anchor point displacement is applied to the original trajectory point to generate a sliding anchor point, and the relative displacement between adjacent anchor points is forced not to exceed the preset neighborhood consistency tolerance.

6. The traffic network data modeling method based on graph neural networks according to claim 5, characterized in that, Step 4 includes: Step 41: Based on the set of sliding anchor points and the set of boundary nodes, calculate the relative displacement vector magnitude between the two sliding anchor points for each pair of boundary nodes, subtract the lower limit of the neighborhood consistency tolerance value as the initial allowable sliding margin, and arrange them in matrix form to obtain the original hinge matrix. Step 42: Relying on the original hinge matrix and the crack opening prediction coefficient field, the prediction coefficient field value is used as the modulation factor: when it is greater than 1, it is divided by the square of the modulation factor; when it is less than 1, it is multiplied by the reciprocal of the modulation factor; when it is equal to 1, it remains unchanged, and the dynamic hinge matrix is ​​output. Step 43: Relying on the negative volume permeation field, extract the local gradient magnitude at the boundary node and output the engagement coefficient through the nonlinear response function: the engagement coefficient is 1 when the local gradient magnitude is lower than the first threshold, decays exponentially to less than 1% of the initial value when it is between the first and second thresholds, and is 0 when it is higher than the second threshold; multiply it element-wise with the dynamic hinge matrix to obtain the hinge matrix with engagement control.

7. The traffic network data modeling method based on graph neural networks according to claim 6, characterized in that, Step 4 also includes: Step 44: Relying on the hinge matrix with engagement control, set an accumulated slip counter and a jump trigger flag for each pair of boundary nodes. Accumulate the actual relative displacement increment in each virtual time step. When the counter value reaches or exceeds the hinge matrix element value, set the trigger flag to 1 and retain the excess amount. If the element value is 0, force it to be set to 0 and clear the counter. Step 45: Relying on the jump trigger flag and the direction of the asymmetric crack operator, when the flag is 1, perform a one-way exchange of the difference of the elastic tensor trace according to the operator direction. After the exchange, the counter is deducted from the corresponding hinge matrix element value and the flag is reset. The region with a biting coefficient of 0 is forcibly prohibited from triggering until the coefficient returns to a positive number.

8. The traffic network data modeling method based on graph neural networks according to claim 7, characterized in that, Step 5 includes: Step 51: Relying on the set of sliding anchor points and the hinge matrix with biting control, the anchor point displacement vector is decomposed into tangential and normal components as boundary conditions; the initial zero second-order tensor of each road segment is defined as follows: tangential strain is positively correlated with the tangential anchor point displacement gradient, and normal strain is positively correlated with the normal anchor point displacement gradient; the deformation rate is controlled by using the reciprocal of the spatiotemporal scale as the relaxation rate and the biting coefficient as the boundary stiffness: the tangential deformation rate is the ratio of tangential strain to relaxation rate multiplied by the boundary stiffness, and the normal deformation rate is the ratio of normal strain to relaxation rate divided by the boundary stiffness, thus obtaining the time-varying second-order tensor field; Step 52: Relying on the time-varying second-order tensor field, sum all tensor traces in each subgraph. If the sum is not zero, calculate the correction amount for each tensor. The correction amount is proportional to the deviation between the trace value and the average trace value of the subgraph multiplied by a weighting factor that is inversely proportional to the current value of the cumulative slip counter. Subtract the correction amount from each tensor to make the sum of traces in the subgraph zero. If a tensor trace value crosses the zero point, scale its tangential and normal components proportionally to maintain the original sign. Output the corrected elastic tensor field.

9. The traffic network data modeling method based on graph neural networks according to claim 8, characterized in that, Step 6 includes: Step 61: Relying on the corrected elastic tensor field and dynamic hinge matrix, perform an asymmetric rotation transformation on the tensor field of each road segment, extract the first derivative along the road direction and the second derivative in the vertical direction to form an equivalent body force distribution, and obtain the stress source term. Step 62: Relying on the equivalent body force distribution and dynamic hinge matrix, the non-zero elements in the dynamic hinge matrix are used as crack constraints to construct Lagrange multipliers. The Poisson equation is transformed into a linear complementary problem with inequality boundary conditions. The problem is solved by successive over-relaxation iterations with a fixed number of iterations equal to twice the road network diameter. The road network diameter is defined as the maximum value of the shortest path between all node pairs. The stress and strain field is then output. Step 63: Based on the stress-strain field and the region with zero interlocking coefficient in step 45, calculate the strain energy density of each road segment, obtain the stress concentration index by second-order difference between the strain energy density of the road segment and the adjacent road segment, mark the road segment with the index exceeding three times the global median as a high uncertainty mask, and automatically include the region with zero interlocking coefficient in the mask, and output the mask label.

10. A traffic network data modeling system based on graph neural networks, applied to the traffic network data modeling method based on graph neural networks as described in any one of claims 1-9, characterized in that, include: The split map scaling module is used to split the map into fast and slow sub-maps according to road type and modal activity. Each sub-map is assigned a spatiotemporal scale and the boundary is set with an asymmetric split operator. The infiltration field phase judgment module is used to convert meteorological events into negative volume infiltration fields that are independent of topology. It calculates the phase difference between the field and the fast and slow subgraphs using a spatiotemporal scale, and the ratio is used as the opening prediction coefficient. The reverse projection anchor point module is used to project the trajectory data onto the equipotential surface with the opening prediction coefficient to generate a sliding anchor point. The offset is determined by the local gradient of the phase difference. The hinge jump module is used to construct a hinge matrix at boundary nodes using a set of sliding anchor points, with the element table showing the maximum sliding amount. The interlocking coefficient is defined by the permeation field gradient. When the gradient exceeds a preset gradient threshold, the interlocking coefficient approaches zero, and the coupling in the corresponding region becomes intermittent jump transmission. The tensor zeroing module is used to construct a pre-distortion space with the sliding anchor point as an elastic constraint and the hinge matrix. Each road segment is a time-varying second-order tensor, and the deformation rate is controlled by the spatiotemporal scale and the biting coefficient. The tensor field is subjected to a referenceless correction, which forces the sum of the sub-traces to be zero. The stress mask module is used to solve the Poisson equation of the road network with the tensor field as the stress source and the hinge matrix as the boundary condition to obtain the stress and strain, and to mark the stress concentration area as an uncertainty mask.

Images