Navigation risk dynamic assessment method based on parallel acceleration

Abstract

The invention belongs to the technical field of navigation risk assessment and flow field analysis, and particularly provides a parallel acceleration-based navigation risk dynamic assessment method, which comprises the following steps of: acquiring data, drawing a two-dimensional unstructured grid, and performing terrain interpolation and numerical interpolation; building a two-dimensional hydrodynamic model based on CUDA parallel calculation, and setting boundary conditions; performing model calculation on the established two-dimensional hydrodynamic model, exporting calculation data, drawing a virtual structure grid, and performing spatial data interpolation; establishing a navigation risk assessment model, and performing parallel calculation of space risks to obtain a space risk rating; establishing a risk assessment display platform based on a Cesium framework, and dynamically updating and displaying navigation risk assessment results under different water flow conditions; and according to the real-time updating of the flow and water level boundary, the space risk rating after the flow and water level boundary conditions are updated. According to the method, the efficiency of flow field simulation and navigation risk assessment is improved through a parallel computing technology, a dynamic risk assessment result is displayed in combination with a three-dimensional visualization platform, and powerful support is provided for safe navigation in a complex water area.

Description

Technical Field

[0001] The invention belongs to the field of flow field simulation and navigation risk assessment, and in particular, relates to a navigation risk dynamic assessment method based on parallel acceleration. Background Art

[0002] In the field of navigation safety assessment, the influence of water flow field characteristics on ship navigation is crucial. With the development of the shipping industry, the density of ship traffic continues to increase, and the prediction and management of navigation risks have become particularly important. However, the existing navigation risk assessment methods still have certain limitations when dealing with flow field information in complex waters, resulting in the accuracy and timeliness of the assessment results being difficult to meet actual needs.

[0003] Most current assessment methods focus on the ship's own operating parameters, such as speed, heading and meteorological conditions, but do not adequately consider real-time flow field information. Although some improved methods have attempted to introduce flow field information, they are mostly based on simplified flow velocity or flow direction, and fail to fully reflect key elements in hydrodynamic characteristics, such as flow velocity distribution, backflow phenomenon, water surface slope, etc. This limitation makes it difficult for existing methods to accurately assess navigation risks in complex waters or variable water flow conditions, limiting their practical application effects.

[0004] In addition, the lack of computational efficiency is another major challenge faced by existing methods. The calculation of hydrodynamic models usually requires the processing of large amounts of data, especially in scenarios with complex terrain or unstructured grids, where the amount of computation increases significantly. Each step of flow field simulation, risk assessment, and visual rendering takes a lot of time, which not only prolongs the duration of the entire evaluation process, but also reduces the real-time response capability of the system. Traditional computing methods are inefficient in dealing with these large-scale computing tasks, resulting in a long time to generate risk assessment results, especially when flow and water level change frequently and dynamically. The assessment results lag behind the actual situation, thus affecting the system's support for real-time navigation decisions.

[0005] Therefore, the existing navigation risk assessment technology has room for improvement in terms of the breadth of flow field information and computational efficiency. In order to improve the accuracy and real-time performance of the assessment, a solution is urgently needed that can not only fully integrate flow field information but also achieve fast calculations through efficient parallel computing technology, so as to better support navigation safety management. Summary of the invention

[0006] The technical problem to be solved by the present invention is to provide a method for dynamic assessment of navigation risk based on parallel acceleration, improve the efficiency of flow field simulation and navigation risk assessment through parallel computing technology, and display the dynamic risk assessment results in combination with a three-dimensional visualization platform.

[0007] In order to solve the above technical problems, the technical solution adopted by the present invention is: a navigation risk dynamic assessment method based on parallel acceleration, comprising the following steps: Step 1: Data collection, drawing two-dimensional unstructured grids, terrain interpolation and numerical interpolation; Step 2: Build a two-dimensional hydrodynamic model based on CUDA parallel computing, set boundary conditions, and calculate the two-dimensional hydrodynamic model; Step 3: Export the calculated data, draw the virtual structure grid, and perform spatial data interpolation; Step 4: Establish a navigation risk assessment model, perform parallel calculation of spatial risks, and obtain a spatial risk rating; Step 5: Establish a risk assessment display platform based on the Cesium framework to dynamically update and display the navigation risk assessment results under different water flow conditions; Step 6: Based on the real-time update of flow and water level boundaries, repeat steps 3 to 5 to update the spatial risk rating after flow and water level boundary conditions.

[0008] In the preferred solution, in step 1, a two-dimensional unstructured grid is constructed under complex terrain and boundary conditions in the water area, and the grid density is controlled by the flow velocity gradient and the local characteristic length. The refinement formula is: (1); in, is the flow velocity gradient, is the local characteristic length, is the refinement threshold.

[0009] In a preferred solution, in step 1, the terrain interpolation operation method is: block processing is performed on the space, the entire interpolation space is divided into n×n sub-blocks, and the terrain interpolation is performed in each sub-block using the inverse distance weighted method, and the difference formula is: (2); in, z ( x ) is the elevation of the node to be interpolated, is the elevation of the known nodes in the sub-block, is the distance between the node to be interpolated and the known node, p is the weight index, k The maximum number of searches to set.

[0010] In a preferred solution, in step 1, the numerical interpolation adopts a linear interpolation method to fill in the data missing and error problems that may occur during the data acquisition process. The linear interpolation expression is: (3); Where y is the estimated value of the interpolation point, x is the independent variable of the interpolation point, x1 and x2 are the independent variables of the known data points, and y1 and y2 are the dependent variables of the known data points.

[0011] In a preferred solution, in step 2, the control equation for solving the two-dimensional hydrodynamic model based on CUDA parallel computing is: (4); ; ; ; ; ; ; ; ; ; Where h is the water depth, U is the conservation vector, G is the flux vector, E and F are x and y Directional flux; S is the source term vector; u and v Respectively represent x and y The velocity component in the direction, g represents the acceleration due to gravity, and Respectively x and y The bottom slope term in the direction of S r is the source and sink term, Z b represents the bottom elevation, and Respectively x and y The friction slope in the direction, n is the Manning coefficient.

[0012] In a preferred embodiment, in step 2, the two-dimensional hydrodynamic model is calculated, and within a single time step, the calculation process of a single grid includes the following steps: S201, parameter update: Replace the result of the previous time step calculation with the initial conditions of the current calculation moment, and update the boundary conditions input at the corresponding moment; S202, state division: According to the updated water depth of each grid unit, the state of the corresponding grid is updated through the processing of the minimum water depth method; S203, flux calculation: Combined with the grid space discretization, the grid area is integrated, and the Gaussian method is applied to convert the volume integral into the grid perimeter line integral, and after discretization, the expression is obtained: (5); In the formula, is the length of each side, is the numerical flux of each edge, is the unit external normal vector passing through the cell boundary j; According to the Riemann problem, the interfaces of adjacent control bodies can be regarded as discontinuities, and a local Riemann problem is generated at any adjacent unit interface. For the local Riemann problem here, the Roe format approximate Riemann solver is used to replace the exact solution. The expression is as follows: (6); In the formula, , is the Jacobian matrix of the Roe format average, and are the conserved variables on both sides of the grid unit; and are the x- and y-direction fluxes of the left and right grids, respectively; S204, bottom slope source term calculation: The discrete equation for solving the bottom slope source term is: (7); ; In the formula, (j=1,2,3) is The components of are the eigenvectors corresponding to the three eigenvalues ​​of the Jacobian matrix, sign is the sign function, is the length of the corresponding side i, is the average wave speed, is the elevation difference between the bottom slopes of the left and right units; S205, friction source term calculation: The discrete conservation equation of the friction source term is: (8); In the formula, is the influence weight coefficient of the next moment, when When , the discrete format is semi-implicit. When it is fully explicit; It is fully implicit, n is the nth moment, and n+1 is the next calculation moment after n.

[0013] In a preferred solution, in step 3, a virtual structure grid is drawn and a spatial data interpolation method is performed as follows: S301, grid generation: constructing a virtual structured grid according to the unstructured grid of the two-dimensional hydrodynamic model, and generating a rectangular structured grid by dividing the simulation area in a regular manner, the area of ​​which is required to cover all nodes and units of the unstructured grid; S302, node data inverse distance weighted interpolation: For each node P in the virtual structured grid i , search and determine the nearest n unstructured grid nodes 𝑃 around it j The position and conservation value of S303, using the inverse distance weighted interpolation method to calculate the virtual structured grid node P i The data value D i , the interpolation formula is: (9); Where d(Pi, 𝑃j) is the distance from the virtual structured grid node Pi to the unstructured grid node 𝑃j, e is the weight index, and n is the number of neighboring unstructured grid nodes used for interpolation.

[0014] In a preferred embodiment, the operation steps of step 4 are as follows: S401. Establish risk assessment indicators, including water surface slope, backflow degree, flow rate and water depth, and establish a risk assessment table; S402. Calculation of risk assessment index: Through the grey correlation analysis method, the correlation is determined according to the set similarity and the discrete sequence is sorted. The closer the spatial set curve is to the standard sequence, the higher the correlation is. The risk criterion sequence level with the highest correlation with the actual indicator sequence is used as the navigation risk level of the node.

[0015] In a preferred solution, the operation process of step S402 is as follows: 1) Determine the risk assessment standard matrix R by taking the average values ​​of the three indicators of water surface slope, backflow degree and flow rate according to the risk assessment table. 5×3 ; 2) Risk assessment standard matrix R 5×3 Normalize, the normalization method is: (10); Among them, N i (k) represents the normalized value of the i-th navigation risk criterion of the k-th indicator, Represents the value of the risk assessment standard of the i-th item in the k-th indicator of the original matrix; 3) Calculate the index sequence of all nodes according to the index calculation method, and all index sequences finally constitute Y n×3 matrix, n is the number of nodes that make up all the grids of the hydrodynamic model, Y n×3The matrix is ​​normalized, and the normalization method is as follows: (11); In the formula, y i (k) is the value of k-item index at the i-th grid node; Finally, the normalized actual index matrix NY is obtained n×3 , as shown below: (12); 4) Calculate the correlation For any grid node i, the three index number sequences of the node constitute a vector , the standard matrix Each rating index in the vector (i=1, 2, ..., 5), calculate the correlation between the node index and the standard risk level, and obtain the absolute difference between the actual index sequence and each rating index sequence , through the relevant discrete function Calculate the correlation as follows: (13); (14); In the formula, Reflects the absolute difference between the k index in node i and the k index in the standard matrix risk level j. When , it means that the indicator k is consistent with the risk criterion, and , indicating the maximum correlation, when When , it means that the node index k does not meet the standard of risk level j, and the correlation is the smallest. , indicating that there is a certain correlation between the index k of node i and the risk level, The larger the correlation, the greater the Calculate all required As the associated discrete function: ; Finally, the objective weight is calculated by super-weight method, and the weighted correlation is calculated according to the weight result and discrete function. The calculation method is as follows: (15); (16); In the formula, is the weighted value of the k index of node i, is the j-level correlation of node i, and the maximum value of the 5-level correlation of node i is taken as the temporary navigation risk level of the node. Finally, the final navigation risk level of the node is determined according to the water depth limit range. If the water depth exceeds the threshold, the navigation risk level remains unchanged; if it is less than the threshold, the navigation node risk level is directly defined as the highest risk.

[0016] In the preferred solution, in step 5, risk assessment visualization is established based on the Cesium framework: the three-dimensional earth browsing function provided by the Cesium framework is used to load terrain and river base map data, combined with real-time ship position and navigation path data, to achieve dynamic display and assessment of risk points; the specific steps include: S501, map the flow field simulation results to the three-dimensional scene, render the high-risk area in real time according to the risk assessment model, and intuitively present the navigation risk through color gradient and other forms. The risk assessment model mapping adopts the internal primitive module rendering, and the flow field rendering adopts the shader method. The two instances use their actual locations as reference points, calculate the transfer matrix through the reference points, and place the whole in the virtual earth. The overall coordinate transfer formula is: (17); Among them, C is represented by a four-dimensional vector, are local coordinates, is the global coordinate, is the fourth-order model conversion matrix; The coordinate transfer formula in the Shader shader is: (18); In the formula, , , , Respectively represent viewport coordinates, clipping coordinates, standardized coordinates and screen coordinates, (C) x 、(C) y 、(C) z 、(C) w Represents vector C=(x,y,z,w) respectively T The corresponding components in , where x, y, z, and w represent four-dimensional vectors, where x, y, and z represent positions, and w is used for mathematical operations and coordinate transformations. , They represent the view transformation matrix and projection transformation matrix respectively. Width and height represent the width and height of the screen resolution respectively.

[0017] The present invention provides a method for dynamic assessment of navigation risk based on parallel acceleration, which has the following beneficial effects: 1. By introducing CUDA parallel computing technology, the solution speed of the hydrodynamic model is significantly accelerated, especially when processing unstructured grids and large-scale data, which greatly reduces the calculation time.

[0018] 2. The present invention combines the Cesium framework for three-dimensional visualization, which realizes the dynamic and intuitive presentation of risk assessment results, and helps to monitor and evaluate risk changes in the navigation environment in real time.

[0019] 3. Compared with traditional risk assessment methods, the present invention can comprehensively consider flow field characteristics, including key factors such as flow velocity, backflow, and water depth, thus avoiding the problem of inaccurate risk assessment caused by ignoring flow field information in previous methods.

[0020] 4. The present invention overcomes the limitations of low computing efficiency and inability to update assessment results in a timely manner in the prior art through an efficient parallel computing architecture, thereby enabling faster acquisition of dynamic assessment results of navigation risks.

[0021] 5. The parallel computing framework adopted by the present invention also has a flexible workload distribution mechanism, which can dynamically adjust computing resources to adapt to changes in complex aquatic environments, further improving the stability and responsiveness of the system.

[0022] 6. This method can effectively improve the computational efficiency of complex water flow field simulation, and improve the real-time and accuracy of navigation risk assessment, significantly improve the accuracy and real-time of navigation safety assessment, and provide more reliable risk prediction and decision-making support for ship navigation. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] The present invention will be further described below in conjunction with the accompanying drawings and embodiments: Figure 1 It is the overall process framework diagram of the present invention; Figure 2 Virtual grid for risk assessment; Figure 3 Calculate flow charts for hydrodynamic and risk assessment in parallel; Figure 4 Schematic diagram of the cutting method; Figure 5 is an example river section grid; Figure 6 A schematic diagram showing the risk assessment results. DETAILED DESCRIPTION

[0024] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0025] Embodiment 1: like Figure 1 and Figure 3 As shown, a method for dynamic assessment of navigation risk based on parallel acceleration includes the following steps: Step 1: Data collection, drawing of two-dimensional unstructured grids, terrain interpolation and numerical interpolation. In actual operation, data verification can also be performed to ensure the accuracy and consistency of model input data.

[0026] Generate an unstructured triangular mesh through Delaunay triangulation, randomly generate point sets, construct initial triangles, and iterate and optimize. The formula is: ; For any triangle , if and only if, there are no other grid points in the circumscribed circle of P, the Delaunay condition is satisfied, and a two-dimensional unstructured grid is constructed under complex terrain and boundary conditions in the water area. A high-precision two-dimensional unstructured grid is constructed in complex waters through an adaptive grid generation algorithm. The algorithm is based on the Delaunay triangulation technology and automatically refines the grid in areas with large velocity changes or complex terrain.

[0027] A two-dimensional unstructured grid is constructed under complex terrain and boundary conditions in the water area. The grid density is controlled by the velocity gradient and the local characteristic length. The refinement formula is: (1); in, is the flow velocity gradient, is the local characteristic length, is the refinement threshold.

[0028] The operation method of terrain interpolation is: divide the space into blocks, divide the entire interpolation space into n×n sub-blocks, and use the inverse distance weighted method to perform terrain interpolation in each sub-block. The difference formula is: (2); in, z ( x ) is the elevation of the node to be interpolated, is the elevation of the known nodes in the sub-block, is the distance between the node to be interpolated and the known node, p is the weight index, which is set to 2 in this embodiment. k The maximum number of searches to set.

[0029] In order to improve the interpolation efficiency, the present invention limits the search range of interpolation points in each sub-block, and ensures that the interpolation calculation can be completed quickly in the case of a super-large grid by setting the maximum search number k and the search time threshold.

[0030] The numerical interpolation adopts a linear interpolation method to complete the data missing and error problems that may occur during the data collection process.

[0031] The linear interpolation expression is: (3); Where y is the estimated value of the interpolation point, x is the independent variable of the interpolation point, x1 and x2 are the independent variables of the known data points, and y1 and y2 are the dependent variables of the known data points.

[0032] The interpolation results are then verified through a data validation check mechanism, which includes deviation analysis and physical consistency check.

[0033] Step 2: Build a two-dimensional hydrodynamic model based on CUDA parallel computing, set boundary conditions, and calculate the two-dimensional hydrodynamic model.

[0034] Boundary types include open boundary, closed boundary (land boundary), inner boundary, and dynamic boundary wet and dry processing; there is no internal boundary problem in this patent. The closed boundary setting is determined according to the external connection of the grid division. The dynamic boundary wet and dry processing adopts the minimum water depth method. Here is the open boundary condition setting. The boundary conditions include common boundary water depth, flow, inflow and sink source items (rainfall, generalized tributary inflow and sink, etc.). The boundary conditions are given water level process, given flow process, and given water level flow relationship. At the inflow open boundary and the outflow open boundary, the corresponding compatibility relationship is determined according to the slow flow and rapid flow states, and the boundary condition setting formula for the corresponding state is further obtained.

[0035] The CUDA parallel computing framework is used to construct a two-dimensional hydrodynamic model, set reasonable boundary conditions, and perform efficient simulation of the flow field.

[0036] Based on the CUDA parallel acceleration framework, implicit global synchronization of multiple kernels is used in the association processing of the hydrodynamic model and the risk assessment model to ensure the consistency of the calculation. In the parallel calculation of the gradual update duration within the hydrodynamics, kernel implicit global synchronization is also used to coordinate the progress between each step. In the local task allocation, in order to avoid frequent thread bundle switching due to the large number of grids, a single thread bundle is responsible for multiple grid calculations, and grid calculation tasks are allocated at fixed intervals to achieve an approximately balanced workload. Furthermore, through the thread bundle monitoring method, combined with the calculation characteristics of the dry and wet boundaries of the model, static load adjustment is performed in advance to optimize the calculation efficiency.

[0037] The control equations for solving the two-dimensional hydrodynamic model based on CUDA parallel computing are: (4); ; ; ; ; ; ; ; ; ; Where h is the water depth, U is the conservation vector, G is the flux vector, E and F are x and y Directional flux; S is the source term vector; u and v Respectively represent x and y The velocity component in the direction, g represents the acceleration due to gravity, and Respectively x and y The bottom slope term in the direction of S r is the source and sink term, Z b represents the bottom elevation, and Respectively x and y The friction slope in the direction, n is the Manning coefficient.

[0038] The explicit finite volume method is used to discretize the integral equations. The discretization scheme uses unstructured triangular grids. The Roe format is used to approximate the flux of the Riemann discontinuity problem. The dry and wet boundary conditions are handled by the minimum water depth method.

[0039] Combined with the grid space discretization, the grid area is integrated, and the Gaussian method is applied to convert the volume integral into the grid perimeter line integral, and after discretization, the following is obtained: (5); In the formula, is the length of each side, is the numerical flux of each edge, is the unit external normal vector passing through the cell boundary j.

[0040] According to the Riemann problem, the interfaces of adjacent control bodies can be regarded as discontinuities, and a local Riemann problem is generated at any adjacent unit interface. For the local Riemann problem here, the Godunov format is used to form the initial value of the piecewise function U (x, t) in the discrete distribution {U} at time t, and the exact solution of the local Riemann problem of each discontinuity is solved, and then the unit interface change at time t+1 is solved. Since the exact solution of the Riemann problem requires too much computation, the present invention uses the Roe format approximate Riemann solver to replace the exact solution, efficiently and accurately capturing the shock wave interface, and its expression is as follows: (6); In the formula , is the Jacobian matrix of the Roe format average, and are the conserved variables on both sides of the grid cell, and are the x and y direction fluxes of the left and right grids respectively.

[0041] The two-dimensional hydrodynamic model is used for model calculation to provide basic data for subsequent risk assessment.

[0042] In the two-dimensional hydrodynamic model based on the explicit finite volume method, the calculation at the time level depends on all the results of the previous step, so only a step-by-step serial operation mode can be used. However, within a single time step, the calculation of the grid depends on the state of the surrounding grids at the previous moment, which makes it possible to parallelize the calculation at the spatial level. Within a single time step, the calculation process of a single grid includes five steps: parameter update, state partition, flux calculation, bottom slope source term calculation, and friction source term calculation. Among them, the calculation of flux and bottom slope source terms involves the interaction between the current grid and the adjacent grids. Although these tasks can be divided in parallel, they need to be quickly reduced after each step, and this reduction operation may take up a lot of computing time. The calculation of the friction source term is performed independently for each grid, and the calculation results of each grid in the next step are obtained after the flux of all source terms is summarized. According to the calculation logic, reduction and synchronization are required before and after the parameter update, so the kernel's implicit global synchronization method is used for processing. The state division, flux and source term calculations are handled by another kernel, and the result output is controlled by the CPU. When the result output time point is reached, the result is temporarily stored in the video memory to provide support for subsequent risk assessment.

[0043] The calculation process of a two-dimensional hydrodynamic model for a single grid in a single time step includes the following steps: S201, parameter update: Replace the result of the previous time step calculation with the initial conditions of the current calculation moment, and update the boundary conditions input at the corresponding moment; S202, state division: according to the updated water depth of each grid unit, the state of the corresponding grid is updated by processing the minimum water depth method; S203, flux calculation: Calculate according to the above formulas (5) and (6).

[0044] S204, bottom slope source term calculation: The discrete equation for solving the bottom slope source term is: (7); ; In the formula, (j=1,2,3) is The components are the eigenvectors corresponding to the three eigenvalues ​​of the Jacobian matrix, sign is the sign function, is the length of the corresponding side i, is the average wave speed, is the elevation difference between the bottom slopes of the left and right units; S205, friction source term calculation: The discrete conservation equation of the friction source term is: (8); In the formula, is the influence weight coefficient of the next moment, when When , the discrete format is semi-implicit. When it is fully explicit; It is fully implicit, n is the nth moment, and n+1 is the next calculation moment after n.

[0045] Step 3: Export the calculated data, draw the virtual structure grid, and perform spatial data interpolation.

[0046] Draw a virtual structure grid and perform spatial data interpolation method as follows: S301, grid generation: constructing a virtual structured grid according to the unstructured grid of the two-dimensional hydrodynamic model, and generating a rectangular structured grid by dividing the simulation area in a regular manner, the area of ​​which is required to cover all nodes and units of the unstructured grid; S302, node data inverse distance weighted interpolation: For each node P in the virtual structured grid i , search and determine the nearest n unstructured grid nodes 𝑃 around it j The position and conservation value of S303, using the inverse distance weighted interpolation method to calculate the virtual structured grid node P i The data value D i , the interpolation formula is: (9); Where d(Pi, 𝑃j) is the distance from the virtual structured grid node Pi to the unstructured grid node 𝑃j, e is the weight index, and n is the number of neighboring unstructured grid nodes used for interpolation.

[0047] Step 4: Establish a navigation risk assessment model, perform parallel calculation of spatial risks, obtain spatial risk ratings, and improve assessment efficiency.

[0048] Before parallel computing risk assessment, it is necessary to prepare the preconditions for virtual grid interpolation in advance and read the global variables stored in the video memory to avoid repeated calculations before each calculation. To facilitate the application of risk assessment models, such as the calculation of water surface slope, it is necessary to establish a virtual grid with a side length of 20m, corresponding to 24 virtual nodes around each real node, as shown in the attached figure. Figure 2 As shown in the figure, the virtual nodes and real nodes form a 5×5 square node matrix. Then, the virtual node coordinates are deduced, and the unstructured grid nodes around each virtual coordinate are retrieved in advance, and the sequence number and distance weight of the associated grid are stored in advance. Furthermore, after the hydrodynamic model calculation results after a certain number of time steps, the risk assessment parallel calculation kernel is connected to obtain the real risk assessment level.

[0049] The steps are as follows: S401. Establish risk assessment indicators and create a risk assessment table; Risk assessment indicators include water surface slope, backflow degree, flow rate and water depth, and each risk assessment indicator includes five risk levels.

[0050] The risk assessment indicators are based on the "Inland Waterway Navigation Standards" to establish four risk assessment indicators: water surface slope, backflow degree, speed and depth. The backflow degree represents the deviation angle between the flow direction of each node and the average flow direction at a certain spatial scale, where the angle is used as a unit of measurement to reflect the disorder of the flow field. In terms of the safety threshold of the indicators, only the water depth has a clear threshold range. According to the different waters in various places, find the local regulations for safe navigation in local waters, clarify the threshold range, determine five risk assessment levels, and establish a risk assessment table.

[0051] S402. Calculation of risk assessment index: Through the grey correlation analysis method, the correlation is determined according to the set similarity and the discrete sequence is sorted. The closer the spatial set curve is to the standard sequence, the higher the correlation is. The risk criterion sequence level with the highest correlation with the actual indicator sequence is used as the navigation risk level of the node.

[0052] The operation process is as follows: 1) Determine the risk assessment standard matrix R by taking the average values ​​of the three indicators of water surface slope, backflow degree and flow rate according to the risk assessment table. 5×3 .

[0053] Take the 25 nodes of the virtual grid as an example, The calculation formula for water surface slope and backflow degree is: ; ; Among them, h i , x i ,y i , and angle i They represent the water depth, abscissa, ordinate, and clockwise angle between the vector velocity direction and the north direction at node i, sl 13 With bf 13 is the water surface slope and backflow degree at the real grid node.

[0054] 2) Risk assessment standard matrix R 5×3 Normalize, the normalization method is: (10); Among them, N i (k) represents the normalized value of the i-th navigation risk criterion of the k-th indicator, Represents the value of the risk assessment standard for the i-th item in the k-th indicator of the original matrix.

[0055] 3) Calculate the index sequence of all nodes according to the index calculation method, and all index sequences finally constitute Y n×3 matrix, n is the number of nodes that make up all the grids of the hydrodynamic model, Y n×3 The matrix is ​​normalized, and the normalization method is as follows: (11); In the formula, y i (k) is the value of k-item index at the i-th grid node; Finally, the normalized actual index matrix NY is obtained n×3 , as shown below: (12); Where n is the number of nodes constituting all grids of the hydrodynamic model. Equation (12) shows that the normalized change method normalizes all indicators in each item (k) to obtain a normalized matrix of 0-1.

[0056] 4) Calculate the correlation For any grid node i, the three index number sequences of the node constitute a vector , the standard matrix Each rating index in the vector (i=1, 2, ..., 5), calculate the correlation between the node index and the standard risk level, and obtain the absolute difference between the actual index sequence and each rating index sequence , through the relevant discrete function Calculate the correlation as follows: (13); (14); In the formula, Reflects the absolute difference between the k index in node i and the k index in the standard matrix risk level j. When , it means that the indicator k is consistent with the risk criterion, and , indicating the maximum correlation, when When , it means that the node index k does not meet the standard of risk level j, and the correlation is the smallest. , indicating that there is a certain correlation between the index k of node i and the risk level, The larger the correlation, the greater the Calculate all required As the associated discrete function: ; Finally, the objective weight is calculated by super-weight method, and the weighted correlation is calculated according to the weight result and discrete function. The calculation method is as follows: (15); (16); In the formula, is the weighted value of the k index of node i, is the j-level correlation of node i, and the maximum value of the 5-level correlation of node i is taken as the temporary navigation risk level of the node. Finally, the final navigation risk level of the node is determined according to the water depth limit range. If the water depth exceeds the threshold, the navigation risk level remains unchanged; if it is less than the threshold, the navigation node risk level is directly defined as 5.

[0057] Step 5: Establish a risk assessment display platform based on the Cesium framework to dynamically update and display the navigation risk assessment results under different water flow conditions.

[0058] The Cesium framework is used to build a risk assessment display platform to achieve real-time dynamic update and visualization of navigation risks, and intuitively display the changes in risks under different water flow conditions.

[0059] Risk assessment visualization is established based on the Cesium framework: the three-dimensional earth browsing function provided by the Cesium framework is used to load terrain and river basemap data, combined with real-time ship position and navigation path data, to achieve dynamic display and assessment of risk points.

[0060] The risk assessment results are finally rendered in the scalar field mode, and the rendering of the scalar field is realized by cutting and coloring. Select the corresponding rendering method according to the target magnitude and grid type.

[0061] The implementation method of the cutting method is as follows: First, interpolate each grid, calculate the focus of the contour line at the grid boundary, and then connect these intersection points to extract the contour of color segmentation; this method is also known as the clipping method. The original method can be applied to both triangular grids and quadrilateral grids. Since the primitives in Cesium are triangular grids, if the two-dimensional hydrodynamic model uses quadrilateral grids for calculation, then before clipping and cutting, using the order formed by the nodes counterclockwise, take the front and back adjacent three points of the quadrilateral as a triangle respectively, and after separation, it is transformed into two triangular grids with a common side.

[0062] For the reconstituted triangle, there are water depth, flow velocity values (or flow velocity in any direction), risk ratings, etc. at the three nodes. Each cut shows one of these scalar values. Here, the risk assessment is taken as an example. The three nodes have corresponding node information, and the triangular grid has the serial numbers corresponding to the three nodes. Before rendering and displaying, preprocess and segment the triangle on the CPU side. The segmentation process and method are as follows: First, set the current triangle state to 0 to indicate the uncut state, and create a linked list of the new triangle topology structure, which is responsible for storing the point-line serial numbers and classification of the new triangles after cutting.

[0063] Set the scalar attribute value, traverse the scalar values of the isosurface. If the scalar attribute value exceeds the range between the minimum value and the maximum value, no operation is performed. Otherwise, taking the point with the minimum value among the three nodes of the triangle as the standard, search again from small to large, judge the step range where the scalar value is exceeded for the first time, make a containment judgment on the remaining two points. After the containment judgment, determine the position of the isoline according to the linear interpolation range, and establish new vertices according to the position, and cut the original triangle into two or three triangles. Here, the following triangle is used as an example. Taking triangle ABC as an example, where the scalar attribute fA < fB < fC; assume that the numerical range of the isoline layer is from fM to fQ. If the value of fM is greater than fC, no segmentation is performed. If fQ is less than fA, no segmentation is performed. When the numerical value fQ of the isoline layer is between fA and fB, according to the linear interpolation method, determine the position of point Q. Set the value of point Q to fQ, and the coordinates are (xQ, yQ). There is the following formula: Figure 4 ; ; Among them, the coordinates of point A and point B are known, and fQ is known. The coordinates of point Q can be obtained by linear interpolation. The same method is used to find point P on the AC side. In the contour line, fQ=fP, and the coordinates of point p are deduced in the same way. Finally, the two newly found vertices are pushed into the updated triangle topology structure vertex list to obtain new triangles AQP, APC, and BQC. Since the grid sequence number does not need to be considered in the display, the three newly obtained unordered triangles are pushed into the triangle topology structure triangle mesh respectively. Then, the cutting state is changed to 1 in the mesh attribute corresponding to AQP, indicating that the triangle cutting is completed. Since it has not been determined whether the PQC triangle and the QBC triangle are between the next contour layer, there is still the possibility of cutting. The cutting state attributes of these two triangles are assigned to 0. And loop through all triangles again, and cut the mesh with the triangle cutting state attribute of 0 again.

[0064] When fQ=fB, according to Figure 4 The cutting is performed in the way (b) in the figure. The two newly generated triangles ABP are assigned the cutting status 1, and BPC is assigned the status 0. When fQ is between fB and fC, Figure 4 As shown in (c), the triangles with cutting status 1 are ABQ and QPA, and PQC is set to 0 to further determine whether the master needs to cut.

[0065] The implementation process of the coloring method is to convert the values ​​in the grid nodes, such as pollution level or water depth, into corresponding colors, and then change the colors of the internal points through interpolation. The scalar values ​​are associated with colors through a step function for rendering.

[0066] The calculation results of the two-dimensional hydrodynamic model based on unstructured grids are displayed. The coloring method has the characteristics of efficient operation and simple programming. The coloring method can use the natural advantages of GPU to concurrently search for internal point colors, while the cutting method needs to find the intersection points of the grid boundaries one by one, and then further bring them into the GPU for color rendering, which is less efficient.

[0067] Specific implementation process: (1) The calculated mesh node data is passed to the CPU, and a native primitive is created for it. The vertex coordinates and sorting are passed to the vertex buffer and index buffer of the primitive. A buffer is created for the scalar data (water depth, risk), which is bound as vertexArray and passed to Cesium's drawCommand together with the vertex array.

[0068] (2) The GPU rendering pipeline is processed in parts, first through the vertex shader. After the scalar data is rasterized by the GPU, it is automatically interpolated from the vertex shader to the fragment shader. The vertex coordinate information is converted into clip coordinates through the Cesium shader built-in variable czm_modelViewProjection.

[0069] (3) Use multiple blending functions to calculate the step function and convert the rasterized scalar data into the scalar color of the corresponding pixel unit.

[0070] The evaluation results need to be placed in the Cesium map to perform spatial coordinate transformation on the rendered object. The overall coordinate transfer formula is: (17); Among them, C is represented by a four-dimensional vector, Local coordinates, is the global coordinate, is the fourth-order model transformation matrix.

[0071] The coordinate transfer formula in the Shader shader is: (18); In the formula, , , , Respectively represent viewport coordinates, clipping coordinates, standardized coordinates and screen coordinates, (C) x (C) y 、(C) z 、(C) w Respectively represent vector C=(x,y,z,w) T The corresponding components in , where x, y, z, and w represent four-dimensional vectors, where x, y, and z represent positions, and w is used for mathematical operations and coordinate transformations. , They represent the view transformation matrix and projection transformation matrix respectively. Width and height represent the width and height of the screen resolution respectively.

[0072] Step 6: Repeat steps 3 to 5 according to the real-time update of flow and water level boundaries, update the spatial risk rating after flow and water level boundary conditions, and dynamically update the visualization results of navigation risk to improve the real-time and accuracy of the assessment.

[0073] Embodiment 2: The parallel accelerated flow field simulation and navigation risk dynamic assessment are applied to the Beijiang waterway. The Beijiang River is the second largest water system in the Pearl River Basin, and its entrance into the Xijiang River is adjacent to the urban area of ​​Guangzhou. In 2019, the Beijiang waterway capacity expansion and upgrading project achieved a major breakthrough, including the official opening of five thousand-ton ship locks including the Feilaixia hub and the Qingyuan hub, and the Beijiang waterway was upgraded from a V-class waterway to a III-class waterway. Finally, the river section between Feilaixia and Qingyuan Water Conservancy Hub of the Beijiang River was the main research area, with a mileage of 44 kilometers. The water conservancy conditions there were calculated, and the three-dimensional model was built and the result risk analysis was performed.

[0074] According to the terrain characteristics of the Beijiang waterway, 253438 unstructured grids were divided. The grid division results are shown in the attached figure. Figure 5 As shown, the grid side length of the V-shaped canyon river channel is about 20m, the grid side length of the spur dike and the boundary section of the entrance and exit is 5m, and the flat area is about 50-100m. The model terrain data file is the water depth terrain from Feilaixia Hub to Qingyuan Hub Station measured in 2019, which is interpolated to the grid terrain according to the 27 mapped water levels in the measured section and the 19,888 water depth measurement points inside. The model is calibrated according to the boundary conditions provided by the hydrological information of Feilaixia and Qingyuan Hub from 2016 to 2020. After the roughness is obtained by calibration, it is brought into the other three working conditions for model verification. The model verification results are shown in Table 1.

[0075] Table 1. Model verification of the section from Feilaixia hub to Qingyuan hub

[0076] Based on the boundary conditions of the outflow from Feilaixia upstream and the water level in front of the Qingyuan Dam downstream, the water conditions in the entire route were simulated and compared with the actual measurement results of the Qingyuan (III) station passing through Qingyuan City. The measured information of this station was provided by the "Open Guangdong" government data unified open platform.

[0077] The verification results show that after simulating the real water conditions of the three types of Feilaixia discharge, small, medium and large, the difference between the model simulation results and the actual observation results is within 0.3m, which meets the calculation accuracy requirements of the model. The simulation results can provide a data basis for subsequent applications and visualization displays.

[0078] In terms of confirming the safety threshold of the evaluation indicators, the relevant research conclusions on the characteristics of the water area in the current case are summarized: for the slope of the water surface, referring to the maximum allowable speed and slope of a 10,000-ton fleet, when the water surface slope is <1 ‱, the impact on ship safety is small, and for navigation safety when it is >4 ‱. For the degree of backflow, according to the regulations that the navigation angle of the ship (fleet) is <10° and the turning angle of the entrance area of ​​the fairway is <20°, it is believed that when the backflow degree is >16°, the flow state is turbulent enough and seriously endangers navigation safety. On the contrary, when the backflow degree is <4°, there is almost no turbulence. In terms of flow velocity, based on traffic signs, when ships are sailing on the water, the flow velocity is required to be small, generally not exceeding 3 meters / second. According to the "Regulations on Safe Navigation of Ships in the Pearl River Estuary Waters" issued by the Guangdong Maritime Safety Administration in 2010, when sailing in the waters north of the Nansha Port Area, the surplus water depth retained by the ship should not be less than 10% of the draft. At the same time, according to the inland navigation standards, the minimum design draft of 1,000-ton ships is 2.0m. It can be determined that passage is prohibited when the water depth is less than 2.2m.

[0079] Combined with the general parameter thresholds of shipping safety, a navigation risk assessment standard for flow conditions is established, and 5 risk levels are divided to represent the level of dynamic changes in navigation risks. Level 1 determines that the flow state has no impact on safe navigation, and level 5 flow state is considered to pose a serious threat to navigation safety. The risk assessment table is established by focusing on safety thresholds such as water surface slope, backflow degree, flow velocity, and water depth, as shown in Table 2.

[0080] Table 2 Risk Assessment Table

[0081] The navigation risk standard includes four indicators, among which the water depth indicator limits the navigation limit of ships, so the water depth indicator can be used as the final boundary judgment condition for risk assessment. The remaining three indicators are determined according to the average value of the risk assessment table (Table 2) to determine the risk assessment standard matrix R 5×3 , as follows: ; According to the process in the above patent method, the risk assessment standard matrix is ​​normalized, and the normalized matrix N 5×3 as follows: ; Four operating conditions were set: the downstream flow rate of Feilaixia upstream was 180 m 3 / s, 2000 m 3 / s, 4000 m 3 / s、6400 m 3 / s, and the water level in Qingyuan was uniformly set to 9.5m. The water conditions in the simulated study area were used to calculate the risk assessment index and risk assessment level. The results were loaded onto the virtual earth using a scalar field visualization method, and the earth map was selected from Tiantu. The overall results are shown in the figure below. Figure 6 shown.

[0082] The present invention realizes the rapid risk assessment response of the Beijiang waterway under different water flow conditions by combining parallel accelerated flow field simulation and dynamic assessment technology. The comparative analysis results show that when the upstream flow is small, the impact of the upstream spur dike on the river section is more significant, and the navigation risk rating is mostly high risk (level 5), while the deep-line area is usually low risk (level 1). As the flow rate increases, the overall navigation risk rating tends to decrease, but after further increasing the flow rate, the risk rating of some areas in the middle and lower reaches increases again, which is manifested as a gradual increase from level 1 to level 2 to level 4. In particular, when the upstream flow reaches 6400m³ / s, the navigation risk rating of narrow rivers in mountainous areas generally rises to level 3 to 4.

[0083] The combination of visualization technology, risk assessment model and WebGIS technology effectively improves the understanding of complex hydrological conditions and provides efficient and accurate navigation risk information for the region. Overall, this method effectively solves the problem of insufficient consideration of water flow conditions in traditional spatial waterway risk assessment, and has a strong universal application value. It is suitable for water resources management, waterway risk assessment and other fields, and can also be applied to parallel computing and WebGIS three-dimensional display of other scientific data (such as meteorological and oceanographic data).

[0084] It will be easily understood by those skilled in the art that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for dynamic assessment of navigation risk based on parallel acceleration, characterized in that: The following steps are involved: Step 1: Data collection, drawing two-dimensional unstructured grids, terrain interpolation and numerical interpolation; Step 2: Build a two-dimensional hydrodynamic model based on CUDA parallel computing, set boundary conditions, and calculate the two-dimensional hydrodynamic model; Step 3: Export the calculated data, draw the virtual structure grid, and perform spatial data interpolation; Step 4: Establish a navigation risk assessment model, perform parallel calculation of spatial risks, and obtain a spatial risk rating; Step 5: Establish a risk assessment display platform based on the Cesium framework to dynamically update and display the navigation risk assessment results under different water flow conditions; Step 6: Based on the real-time update of flow and water level boundaries, repeat steps 3 to 5 to update the spatial risk rating after flow and water level boundary conditions.

2. The method for dynamic assessment of navigation risk based on parallel acceleration according to claim 1, characterized in that: In step 1, a two-dimensional unstructured grid is constructed under complex terrain and boundary conditions in the water area. The grid density is controlled by the velocity gradient and the local characteristic length. The refinement formula is: （1）； in, is the flow velocity gradient, is the local characteristic length, is the refinement threshold.

3. The method for dynamic assessment of navigation risk based on parallel acceleration according to claim 1, characterized in that: In step 1, the terrain interpolation operation method is: block processing is performed on the space, the entire interpolation space is divided into n×n sub-blocks, and the terrain interpolation is performed in each sub-block using the inverse distance weighted method. The difference formula is: （2）； in, z ( x ) is the elevation of the node to be interpolated, is the elevation of the known nodes in the sub-block, is the distance between the node to be interpolated and the known node, p is the weight index, k The maximum number of searches to set.

4. The method for dynamic assessment of navigation risk based on parallel acceleration according to claim 1, characterized in that: In step 1, the numerical interpolation adopts a linear interpolation method to fill in the data missing and error problems that may occur during the data collection process. The linear interpolation expression is: （3）； Where y is the estimated value of the interpolation point, x is the independent variable of the interpolation point, x1 and x2 are the independent variables of the known data points, and y1 and y2 are the dependent variables of the known data points.

5. The method for dynamic assessment of navigation risk based on parallel acceleration according to claim 1, characterized in that: In step 2, the control equation for solving the two-dimensional hydrodynamic model based on CUDA parallel computing is: （4）； ； ； ； ； ； ； ； ； ； Where h is the water depth, U is the conservation vector, G is the flux vector, E and F are x and y Directional flux; S is the source term vector; u and v Respectively represent x and y The velocity component in the direction, g represents the acceleration due to gravity, and Respectively x and y The bottom slope term in the direction of S r is the source and sink term, Z b represents the bottom elevation, and Respectively x and y The friction slope in the direction, n is the Manning coefficient.

6. The method for dynamic assessment of navigation risk based on parallel acceleration according to claim 1, characterized in that: In step 2, the two-dimensional hydrodynamic model is calculated. In a single time step, the calculation process of a single grid includes the following steps: S201, parameter update: Replace the result of the previous time step calculation with the initial conditions of the current calculation moment, and update the boundary conditions input at the corresponding moment; S202, state division: According to the updated water depth of each grid unit, the state of the corresponding grid is updated through the processing of the minimum water depth method; S203, flux calculation: Combined with the grid space discretization, the grid area is integrated, and the Gaussian method is applied to convert the volume integral into the grid perimeter line integral, and after discretization, the expression is obtained: （5）； In the formula, is the length of each side, is the numerical flux of each edge, is the unit external normal vector passing through the cell boundary j; According to the Riemann problem, the interfaces of adjacent control bodies can be regarded as discontinuities, and a local Riemann problem is generated at any adjacent unit interface. For the local Riemann problem here, the Roe format approximate Riemann solver is used to replace the exact solution. The expression is as follows: （6）； In the formula, , is the Jacobian matrix of the Roe format average, and are the conserved variables on both sides of the grid unit; and are the x- and y-direction fluxes of the left and right grids, respectively; S204, bottom slope source term calculation: The discrete equation for solving the bottom slope source term is: （7）； ； In the formula, (j=1,2,3) is The components of are the eigenvectors corresponding to the three eigenvalues ​​of the Jacobian matrix, sign is the sign function, is the length of the corresponding side i, is the average wave speed, is the elevation difference between the bottom slopes of the left and right units; S205, friction source term calculation: The discrete conservation equation of the friction source term is: （8）； In the formula, is the influence weight coefficient of the next moment, when When , the discrete format is semi-implicit. When it is fully explicit; It is fully implicit, n is the nth moment, and n+1 is the next calculation moment after n.

7. The method for dynamic assessment of navigation risk based on parallel acceleration according to claim 1, characterized in that: In step 3, a virtual structure grid is drawn and a spatial data interpolation method is performed as follows: S301, grid generation: constructing a virtual structured grid according to the unstructured grid of the two-dimensional hydrodynamic model, and generating a rectangular structured grid by dividing the simulation area in a regular manner, the area of ​​which is required to cover all nodes and units of the unstructured grid; S302, node data inverse distance weighted interpolation: For each node P in the virtual structured grid i , search and determine the nearest n unstructured grid nodes 𝑃 around it j The position and conservation value of S303, using the inverse distance weighted interpolation method to calculate the virtual structured grid node P i The data value D i , the interpolation formula is: （9）； Where d(Pi, 𝑃j) is the distance from the virtual structured grid node Pi to the unstructured grid node 𝑃j, e is the weight index, and n is the number of neighboring unstructured grid nodes used for interpolation.

8. The method for dynamic assessment of navigation risk based on parallel acceleration according to claim 1, characterized in that: The operation steps of step 4 are as follows: S401. Establish risk assessment indicators, including water surface slope, backflow degree, flow rate and water depth, and establish a risk assessment table; S402. Calculation of risk assessment index: Through the grey correlation analysis method, the correlation is determined according to the set similarity and the discrete sequence is sorted. The closer the spatial set curve is to the standard sequence, the higher the correlation is. The risk criterion sequence level with the highest correlation with the actual indicator sequence is used as the navigation risk level of the node.

9. The method for dynamic assessment of navigation risk based on parallel acceleration according to claim 8, characterized in that: The operation process of step S402 is as follows: 1) Determine the risk assessment standard matrix R by taking the average values ​​of the three indicators of water surface slope, backflow degree and flow rate according to the risk assessment table. 5×3 ; 2) Risk assessment standard matrix R 5×3 Normalize, the normalization method is: （10）； Among them, N i (k) represents the normalized value of the i-th navigation risk criterion of the k-th indicator, Represents the value of the risk assessment standard of the i-th item in the k-th indicator of the original matrix; 3) Calculate the index sequence of all nodes according to the index calculation method, and all index sequences finally constitute Y n×3 matrix, n is the number of nodes that make up all the grids of the hydrodynamic model, Y n×3 The matrix is ​​normalized, and the normalization method is as follows: （11）； In the formula, y i (k) is the value of k-item index at the i-th grid node; Finally, the normalized actual index matrix NY is obtained n×3 , as shown below: （12）； 4) Calculate the correlation For any grid node i, the three index number sequences of the node constitute a vector , the standard matrix Each rating index in the vector (i=1, 2, ..., 5), calculate the correlation between the node index and the standard risk level, and obtain the absolute difference between the actual index sequence and each rating index sequence , through the relevant discrete function Calculate the correlation as follows: (13)； (14)； In the formula, Reflects the absolute difference between the k index in node i and the k index in the standard matrix risk level j. When , it means that the indicator k is consistent with the risk criterion, and , indicating the maximum correlation, when When , it means that the node index k does not meet the standard of risk level j, and the correlation is the smallest. , indicating that there is a certain correlation between the index k of node i and the risk level, The larger the correlation, the greater the Calculate all required As the associated discrete function: ; Finally, the objective weight is calculated by super-weight method, and the weighted correlation is calculated according to the weight result and discrete function. The calculation method is as follows: (15)； (16)； In the formula, is the weighted value of the k index of node i, is the j-level correlation of node i, and the maximum value of the 5-level correlation of node i is taken as the temporary navigation risk level of the node. Finally, the final navigation risk level of the node is determined according to the water depth limit range. If the water depth exceeds the threshold, the navigation risk level remains unchanged; if it is less than the threshold, the navigation node risk level is directly defined as the highest risk.

10. The method for dynamic assessment of navigation risk based on parallel acceleration according to claim 8, characterized in that: In step 5, risk assessment visualization is established based on the Cesium framework: the three-dimensional earth browsing function provided by the Cesium framework is used to load terrain and river base map data, combined with real-time ship position and navigation path data, to achieve dynamic display and assessment of risk points; The specific steps include: S501, map the flow field simulation results to the three-dimensional scene, render the high-risk area in real time according to the risk assessment model, and intuitively present the navigation risk through color gradient and other forms. The risk assessment model mapping adopts the internal primitive module rendering, and the flow field rendering adopts the shader method. The two instances use their actual locations as reference points, calculate the transfer matrix through the reference points, and place the whole in the virtual earth. The overall coordinate transfer formula is: (17)； Among them, C is represented by a four-dimensional vector, are local coordinates, is the global coordinate, is the fourth-order model conversion matrix; The coordinate transfer formula in the Shader shader is: （18）； In the formula, , , , Respectively represent viewport coordinates, clipping coordinates, standardized coordinates and screen coordinates, (C) x 、(C) y 、(C) z (C) w Represents vector C=(x,y,z,w) respectively T The corresponding components in , where x, y, z, and w represent four-dimensional vectors, where x, y, and z represent positions, and w is used for mathematical operations and coordinate transformations. , They represent the view transformation matrix and projection transformation matrix respectively. Width and height represent the width and height of the screen resolution respectively.

Images