Dynamic change monitoring method and system based on multi-temporal river channel image

Abstract

The invention provides a dynamic change monitoring method and system based on a multi-temporal river channel image, and relates to the technical field of image monitoring, and the method comprises the steps: obtaining a multi-temporal multispectral remote sensing image and an unmanned plane oblique photogrammetry image; identifying a riverway water body boundary and constructing a water area mask; dense point cloud data are generated and cut; generating a balanced riverbed point cloud through density compensation sampling and local curvature constraint; constructing a terrain model to generate a section sequence; calculating hydrodynamic correlation strength to construct an evolution propagation sequence; and recursively calculating the erosion and deposition change rate and performing early warning. According to the invention, high-precision monitoring and early warning of river form changes are realized.

Description

Technical Field

[0001] The present invention relates to image monitoring technology, and particularly to a method and system for dynamically monitoring the changes of a river channel based on multi-temporal river channel images. Background Art

[0002] River channels are important carriers for the social and economic development of human society. The dynamic monitoring of river channels is of great significance for flood control and disaster reduction, water resource management, and ecological environment protection. With the development of remote sensing technology and unmanned aerial vehicle (UAV) technology, the use of multi-source remote sensing data for river channel dynamic monitoring has become a current research hotspot. Traditional river channel monitoring mainly relies on manual field surveys and single satellite image analysis, and these methods are difficult to achieve continuous monitoring and accurate prediction of river channel changes. Modern river channel monitoring methods mainly rely on multi-temporal remote sensing images and UAV aerial photography technology. By processing and analyzing multi-temporal images, the dynamic monitoring of the water area range, riverbed topography, and erosion and deposition changes of the river channel is realized.

[0003] The existing technologies generally have the problem of uneven density when processing river channel point cloud data. Especially in areas with complex riverbed topography, the point cloud distribution is sparse and uneven, making it difficult to accurately reflect the subtle changes in the riverbed topography, thus affecting the accuracy and reliability of subsequent river channel change analysis.

[0004] The existing river channel cross-section analysis methods often lack the consideration of the hydrodynamic correlation between cross-sections. Simply relying on static cross-section comparison is difficult to reveal the internal mechanism of river channel evolution and cannot effectively capture the propagation law of river channel changes, resulting in insufficient accuracy of prediction results.

[0005] The existing river channel monitoring systems lack an effective early warning mechanism, making it difficult to detect abnormal erosion and deposition changes of the riverbed in a timely manner, unable to provide timely decision-making support for flood control and disaster reduction and river channel regulation, and reducing the practical value and prevention effect of river channel monitoring. Summary of the Invention

[0006] The embodiments of the present invention provide a method and system for dynamically monitoring the changes of a river channel based on multi-temporal river channel images, which can solve the problems in the existing technologies.

[0007] In the first aspect of the embodiments of the present invention, a method for dynamically monitoring the changes of a river channel based on multi-temporal river channel images is provided, including: Obtaining multi-temporal multi-spectral remote sensing image data and UAV oblique photogrammetry image data of a river channel area; based on the multi-temporal multi-spectral remote sensing image data, identifying the boundary of the river channel water body through a support vector machine classification algorithm and constructing a river channel water area mask; Based on the UAV oblique photogrammetry image data, generating dense point cloud data of the river channel area by using a multi-view image matching method; using the river channel water area mask to crop the dense point cloud data to obtain point cloud data; Calculate the point cloud neighborhood density distribution function for the said point cloud data, determine the local importance weight according to the gradient change of the point cloud neighborhood density distribution function, and establish a density compensation sampling matrix based on the local importance weight; use the density compensation sampling matrix to perform recursive multi-scale enhancement on the sparse area of the point cloud, and optimize the spatial position in combination with the local curvature constraint to generate density-balanced riverbed point cloud data; Construct a triangular grid terrain model based on the riverbed point cloud data, and generate a river channel cross-section sequence at preset intervals along the river flow direction; calculate the hydrodynamic correlation intensity for the river channel cross-section sequence through a recurrent neural network, and construct a cross-section evolution propagation sequence according to the hydrodynamic correlation intensity; Based on the cross-section evolution propagation sequence and historical evolution data, obtain the riverbed erosion and deposition change rate through recursive calculation; when the riverbed erosion and deposition change rate exceeds the preset rate threshold, generate a warning message, and visually display the dynamic prediction result of the river channel cross-section morphology and the warning message.

[0008] Based on the multi-temporal and multi-spectral remote sensing image data, identify the river channel water body boundary through a support vector machine classification algorithm, and construct a river channel water area mask, including: Construct a feature vector based on the multi-temporal and multi-spectral remote sensing image data, map the feature vector to a high-dimensional feature space, construct an optimal classification hyperplane in the high-dimensional feature space, calculate the distance from the sample point to the hyperplane based on the optimal classification hyperplane, use the sample points with distances greater than the preset distance threshold as training support vectors, and construct a classification decision function using the training support vectors; classify and predict the feature vector based on the classification decision function, and generate an initial classification image according to the prediction result; Perform morphological processing on the initial classification image, use the first structural element to perform opening operation to remove noise points, use the second structural element to perform closing operation to fill holes, and obtain an initial water area mask image; Use a contour evolution model with boundary curvature constraint to smooth the initial water area mask image, initialize the control point sequence of the contour curve based on the initial water area mask image, construct an internal energy function representing the elasticity and rigidity of the contour curve according to the first derivative and second derivative of the control point sequence, and construct an external energy function based on the gradient information of the initial classification image; Combine the internal energy function and the external energy function to obtain the total energy function of the contour evolution model, perform minimization solution on the total energy function, stop iteration when the change amount of the total energy function is less than the preset convergence threshold, and reconstruct a smooth river channel water body boundary based on the final control point sequence to generate a river channel water area mask.

[0009] Based on the drone oblique photogrammetry image data, a multi-view image matching method is used to generate dense point cloud data of the river channel area; the dense point cloud data is cropped using the river channel water mask, and the obtained point cloud data includes: A matching pair group is established for the drone oblique photogrammetry image data, the ratio of the overlapping area of each pair of images in the matching pair group to the total area of the images is calculated, and the image pairs with an overlapping ratio greater than a preset overlapping threshold are selected as the image pairs to be matched; For the reference image in the image pairs to be matched, the normalized cross-correlation coefficient between each pixel point in the reference image and the corresponding pixel point in other matching images is calculated, the weight coefficient is determined according to the baseline length of each matching image, and the normalized cross-correlation coefficient and the weight coefficient are weighted and summed to obtain the multi-baseline matching cost; Within the neighborhood window of each pixel point in the reference image, an adaptive weight is calculated based on the similarity of pixel gray values, and the adaptive weight and the multi-baseline matching cost are weighted and summed to obtain the aggregated matching cost; The optimal disparity value of the pixel point is determined using the aggregated matching cost, a collinearity equation system is established based on the optimal disparity value, and the collinearity equation system is solved by the least squares method to obtain the three-dimensional point coordinates, generating dense point cloud data; The dense point cloud data is transformed to the image plane coordinate system through a projection matrix, and the projected point cloud data is cropped according to the river channel water mask, and the point cloud data located within the river channel water mask range is retained as the point cloud data of the river channel area.

[0010] The optimal disparity value of the pixel point is determined using the aggregated matching cost, a collinearity equation system is established based on the optimal disparity value, and the collinearity equation system is solved by the least squares method to obtain the three-dimensional point coordinates, generating dense point cloud data including: Based on the aggregated matching cost, a disparity space cost curve is constructed within a preset disparity search range, the integer pixel extreme points of the disparity space cost curve are determined, and the optimal disparity value is obtained by performing quadratic curve fitting using the integer pixel extreme points and the cost values at their adjacent disparity positions; According to the optimal disparity value and the interior and exterior orientation elements of the image, a collinearity equation system is established, the collinearity equation system is solved to obtain an initial three-dimensional coordinate point set, and a spatial point cloud is constructed based on the initial three-dimensional coordinate point set; Each three-dimensional coordinate point in the spatial point cloud is re-projected onto the images of each view, the reprojection error between the projected point and the original matching point is calculated, and the three-dimensional coordinate points with a reprojection error greater than a preset error threshold are marked as the points to be optimized; For the point to be optimized, a local point set is extracted within its neighborhood range, the centroid coordinates of the local point set are calculated, and the three-dimensional coordinates of the point to be optimized are updated and optimized based on the centroid coordinates to obtain an optimized three-dimensional coordinate point set; Dense point cloud reconstruction is performed on the optimized three-dimensional coordinate point set, the point density of each point in the reconstructed point cloud within its preset neighborhood range is calculated, and the area where the point density is lower than the preset density threshold is locally encrypted and reconstructed to obtain the final dense point cloud data.

[0011] Using the density compensation sampling matrix to perform recursive multi-scale enhancement on the sparse area of the point cloud, and combining local curvature constraints to optimize the spatial position, generating density-balanced riverbed point cloud data includes: Construct a multi-scale feature extraction network, extract features from multiple preset radius neighborhoods of each point in the point cloud data based on the multi-scale feature extraction network, perform weighted fusion of the extracted features with the corresponding weight coefficients to obtain a feature vector, input the feature vector into the density prediction network to obtain a predicted density value, and multiply the predicted density value by the initial density field to obtain the density prediction result of each point; Based on the density prediction result, determine the sparse area in the point cloud, generate a geomorphic feature weight, divide the preset target density value by the density prediction result and multiply it by the geomorphic feature weight to obtain a regional sampling weight, and adaptively fuse the sonar sounding data and multi-temporal multi-spectral remote sensing image data based on the regional sampling weight to generate enhanced point cloud data; Construct a local connection graph for the enhanced point cloud data, iteratively update the node features in the local connection graph to obtain optimized node features; calculate the distance weight of each point within a preset radius neighborhood based on the optimized node features, and perform weighted summation of the distance weight and the projection distance from the point to the neighborhood point to obtain the local curvature value; Input the optimized node features and the local curvature value into the generator to generate an initial density-balanced point cloud, and use the generator to regenerate the point cloud data with local geometric features; construct the gradient information of the flow velocity field and water depth field based on the hydrological monitoring data, and fuse the gradient information with the point cloud data to optimize the spatial position to generate density-balanced riverbed point cloud data.

[0012] Based on the riverbed point cloud data, construct a triangular grid terrain model, and generate a river channel cross-section sequence at preset intervals along the river channel flow direction; calculate the hydrodynamic correlation intensity for the river channel cross-section sequence through a recursive neural network, and construct a cross-section evolution propagation sequence according to the hydrodynamic correlation intensity includes: Perform constrained Delaunay triangulation on the riverbed point cloud data to construct an initial triangular mesh, calculate the mesh quality parameters based on the ratio of the triangle area to the side length, and optimize the triangular mesh according to the mesh quality parameters to obtain a triangular mesh terrain model; Construct a terrain elevation distribution based on the triangular mesh terrain model, calculate the elevation difference between adjacent mesh nodes, determine the water flow direction field based on the elevation difference, sample at preset intervals along the water flow direction field on the riverbed terrain triangular mesh model to generate a river channel cross-section sequence, extract the width, water depth, cross slope, and roughness coefficient parameters of each cross-section, and construct a cross-section feature vector; Input the cross-section feature vector into a long short-term memory network, extract the cross-section time series features through the memory unit and hidden state of the long short-term memory network, calculate the cosine similarity between adjacent cross-sections based on the hidden state, and multiply the cosine similarity by the exponential decay function of the cross-section spacing to obtain the cross-section correlation strength; construct a propagation weight matrix based on the cross-section correlation strength and the hydraulic gradient factor, and use the propagation weight matrix to iteratively update the cross-section state to obtain a river channel cross-section evolution propagation sequence.

[0013] Based on the cross-section evolution propagation sequence and historical evolution data, obtain the riverbed erosion and deposition change rate through recursive calculation; when the riverbed erosion and deposition change rate exceeds a preset rate threshold, generate a warning message, and visually display the dynamic prediction result of the river channel cross-section morphology and the warning message, including: Perform multi-source data fusion on the cross-section evolution propagation sequence and historical evolution data to obtain a fusion feature vector, construct a feature weight and quality evaluation function based on the fusion feature vector, and perform quality evaluation on the fusion feature vector through the feature weight and the quality evaluation function to obtain an evaluation result; Construct a graph structure based on the evaluation result, construct a node set from the cross-section data, construct an associated edge set from the association between cross-sections, generate an adjacency matrix based on the node set and the associated edge set, and use the adjacency matrix to update the feature information of adjacent cross-sections through message passing to obtain hidden state features; Input the hidden state features into a mapping function to obtain the current change rate, perform weighted fusion of the current change rate and the historical change rate to obtain the erosion and deposition change rate of the target cross-section, calculate statistical parameters based on the erosion and deposition change rate, and perform exponential mapping on the statistical parameters and the benchmark threshold to obtain a dynamic warning threshold; Compare the erosion and deposition change rate with the dynamic warning threshold to obtain a warning level, perform time series integration on the target cross-section based on the erosion and deposition change rate to obtain a predicted morphology, and perform three-dimensional visual display of the predicted morphology and the warning level to generate an interactive warning message.

[0014] In the second aspect of the embodiments of the present invention, a dynamic change monitoring system for multi-temporal river channel images is provided, including: A first unit, configured to obtain multi-temporal multi-spectral remote sensing image data and unmanned aerial vehicle (UAV) oblique photogrammetry image data of a river channel area; based on the multi-temporal multi-spectral remote sensing image data, identify the river channel water body boundary through a support vector machine classification algorithm, and construct a river channel water area mask; A second unit, configured to generate dense point cloud data of the river channel area based on the UAV oblique photogrammetry image data; use the river channel water area mask to crop the dense point cloud data to obtain point cloud data; A third unit, configured to calculate a point cloud neighborhood density distribution function for the point cloud data, determine a local importance weight according to the gradient change of the point cloud neighborhood density distribution function, and establish a density compensation sampling matrix based on the local importance weight; use the density compensation sampling matrix to perform recursive multi-scale enhancement on the point cloud sparse area, and optimize the spatial position in combination with local curvature constraints to generate density-balanced riverbed point cloud data; A fourth unit, configured to construct a triangular grid terrain model based on the riverbed point cloud data, and generate a sequence of river channel cross-sections at preset intervals along the river channel flow direction; calculate the hydrodynamic correlation intensity for the sequence of river channel cross-sections through a recurrent neural network, and construct a cross-section evolution propagation sequence according to the hydrodynamic correlation intensity; A fifth unit, configured to obtain the riverbed scouring and silting change rate through recursive calculation based on the cross-section evolution propagation sequence and historical evolution data; when the riverbed scouring and silting change rate exceeds a preset rate threshold, generate a warning message, and visually display the dynamic prediction result of the river channel cross-section morphology and the warning message.

[0015] In the third aspect of the embodiments of the present invention, an electronic device is provided, including: A processor; A memory for storing instructions executable by the processor; Wherein, the processor is configured to call the instructions stored in the memory to execute the method described above.

[0016] In the fourth aspect of the embodiments of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored, and when the computer program instructions are executed by a processor, the method described above is implemented.

[0017] The beneficial effects of this application are as follows: This method identifies the river channel water body boundary through a support vector machine classification algorithm to construct a water area mask, and combines the dense point cloud generated from UAV oblique photogrammetry data, realizing the accurate extraction and expression of the river channel area, and improving the spatial accuracy and reliability of river channel monitoring.

[0018] To address the issue of uneven distribution of point cloud data, this method innovatively proposes a density compensation mechanism based on the point cloud neighborhood density distribution function and local importance weights. Through recursive multi-scale enhancement and local curvature constraint optimization, it generates riverbed point cloud data with balanced density, effectively solving the problem of terrain expression deviation caused by uneven point cloud distribution in traditional methods.

[0019] This method uses a recursive neural network to calculate the hydrodynamic correlation strength, constructs a cross-section evolution propagation sequence, and recursively calculates the rate of riverbed erosion and deposition change in combination with historical data, realizing the dynamic prediction and early warning of river channel morphological changes, providing a scientific basis and technical support for river channel management and flood control and disaster reduction, and having important practical value. Brief Description of the Drawings

[0020] Figure 1 It is a schematic flowchart of the method for dynamically monitoring the changes of multi-temporal river channel images according to an embodiment of the present invention; Figure 2 It is a bar chart for comparative analysis of the performance of the multi-view image matching method according to an embodiment of the present invention; Figure 3 It is a flowchart for the density equalization processing of the riverbed point cloud data according to an embodiment of the present invention; Figure 4 It is a schematic diagram for comparing the optimization effects of the triangular mesh quality according to an embodiment of the present invention; Figure 5 It is a bar chart for comparing the performance of the riverbed erosion and deposition change monitoring and early warning system according to an embodiment of the present invention. Detailed Description of the Embodiments

[0021] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0022] The technical solutions of the present invention will be described in detail below with specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be repeated in some embodiments.

[0023] Figure 1 It is a schematic flowchart of the method for dynamically monitoring the changes of multi-temporal river channel images according to an embodiment of the present invention, as Figure 1 shown, the method includes: Obtain multi-temporal multi-spectral remote sensing image data and unmanned aerial vehicle (UAV) oblique photogrammetry image data of the river channel area; based on the multi-temporal multi-spectral remote sensing image data, identify the river channel water body boundary through the support vector machine classification algorithm, and construct a river channel water area mask; Based on the UAV oblique photogrammetry image data, generate dense point cloud data of the river channel area using the multi-view image matching method; use the river channel water area mask to crop the dense point cloud data to obtain point cloud data; Calculate the point cloud neighborhood density distribution function for the point cloud data, determine the local importance weight according to the gradient change of the point cloud neighborhood density distribution function, and establish a density compensation sampling matrix based on the local importance weight; use the density compensation sampling matrix to perform recursive multi-scale enhancement on the sparse area of the point cloud, and optimize the spatial position in combination with the local curvature constraint to generate density-balanced riverbed point cloud data; Construct a triangular grid terrain model based on the riverbed point cloud data, and generate a river channel cross-section sequence at preset intervals along the river channel flow direction; calculate the hydrodynamic correlation intensity for the river channel cross-section sequence through a recurrent neural network, and construct a cross-section evolution propagation sequence according to the hydrodynamic correlation intensity; Based on the cross-section evolution propagation sequence and historical evolution data, obtain the riverbed erosion and deposition change rate through recursive calculation; when the riverbed erosion and deposition change rate exceeds the preset rate threshold, generate a warning message, and visually display the dynamic prediction result of the river channel cross-section morphology and the warning message.

[0024] In an alternative embodiment, based on the multi-temporal multi-spectral remote sensing image data, identifying the river channel water body boundary through the support vector machine classification algorithm and constructing a river channel water area mask includes: Construct a feature vector based on the multi-temporal multi-spectral remote sensing image data, map the feature vector to a high-dimensional feature space, construct an optimal classification hyperplane in the high-dimensional feature space, calculate the distance from the sample point to the hyperplane based on the optimal classification hyperplane, use the sample points with a distance greater than the preset distance threshold as training support vectors, and construct a classification decision function using the training support vectors; classify and predict the feature vector based on the classification decision function, and generate an initial classification image according to the prediction result; Perform morphological processing on the initial classification image, use the first structural element to perform opening operation to remove noise points, and use the second structural element to perform closing operation to fill holes to obtain an initial water area mask image; The initial water area mask image is smoothed using a contour evolution model with boundary curvature constraints. The control point sequence of the contour curve is initialized based on the initial water area mask image. An internal energy function representing the elasticity and rigidity of the contour curve is constructed according to the first and second derivatives of the control point sequence. An external energy function is constructed based on the gradient information of the initial classification image; The internal energy function and the external energy function are combined to obtain the total energy function of the contour evolution model. The total energy function is minimized. When the change amount of the total energy function is less than a preset convergence threshold, the iteration stops. Based on the final control point sequence, a smooth river channel water body boundary is reconstructed to generate a river channel water area mask.

[0025] Utilizing the high-dimensional mapping ability of the support vector machine classification algorithm and combining morphological processing and a contour evolution model with boundary curvature constraints, the accurate identification of the river channel water body boundary is achieved.

[0026] After obtaining multi-temporal and multi-spectral remote sensing image data, feature vectors need to be constructed based on this data. Feature vectors usually include reflectance information of each band, the normalized difference water index (NDWI), the normalized difference vegetation index (NDVI), etc. For example, for a certain pixel point, its feature vector can be expressed as [B1, B2, B3, B4, NDWI, NDVI], where B1 - B4 represent the reflectance values of different bands, NDWI = (B2 - B4) / (B2 + B4), and NDVI = (B4 - B3) / (B4 + B3). After constructing the feature vectors, they are mapped to a high-dimensional feature space through a kernel function. In practice, the radial basis function (RBF) kernel usually performs well, and the kernel parameter γ can be set to 0.125.

[0027] In the high-dimensional feature space, the support vector machine algorithm searches for the optimal classification hyperplane that can separate water body and non-water body samples. This hyperplane is determined by maximizing the margin between the two types of samples. The distance from the sample points to the hyperplane is calculated, and the sample points with a distance greater than a preset distance threshold (such as 0.8) are used as training support vectors. Based on these support vectors, a classification decision function is constructed, and this function can be used to predict the class of new samples. The feature vectors are classified and predicted to obtain an initial classification image, where water body pixels are labeled as 1 and non-water body pixels are labeled as 0.

[0028] Initial classification images usually contain noise and holes and need to be processed morphologically to improve accuracy. Morphological processing includes two steps: opening operation and closing operation. The opening operation is used to remove noise points, and the closing operation is used to fill holes. In the opening operation, a circular structuring element with a size of 3×3 is used as the first structuring element. First, an erosion operation is performed, and then a dilation operation is performed to effectively remove misclassified water body points with small areas. In the closing operation, a circular structuring element with a size of 5×5 is used as the second structuring element. First, a dilation operation is performed, and then an erosion operation is performed to fill small holes in the water body area. These two operations are used in combination to obtain the initial water area mask image.

[0029] To further optimize the water body boundary, a contour evolution model with boundary curvature constraint is used to smooth the initial water area mask image. This process starts from the control point sequence of the initialized contour curve, and these control points are evenly distributed along the boundary of the initial water area mask image. For an image with a smoothing processing area size of 512×512 pixels, 100 - 200 control points can be selected. According to the first derivative and second derivative of the control point sequence, an internal energy function representing the elasticity and rigidity of the contour curve is constructed. The first derivative controls the elasticity of the contour, and the second derivative controls the rigidity of the contour. The elastic weight coefficient can be set to 0.45, and the rigid weight coefficient can be set to 0.55.

[0030] An external energy function is constructed based on the gradient information of the initial classification image to guide the contour to approach the real water body boundary. Areas with larger gradient values usually represent boundary positions. Through gradient calculation, the gradient amplitude and direction of each pixel in the image can be obtained. The gradient amplitude can be calculated by the Sobel operator, and the threshold is set to 25. Pixels exceeding this value are considered boundary points. The weight coefficient of the external energy function can be set to 0.65.

[0031] The internal energy function and the external energy function are combined to obtain the total energy function of the contour evolution model. To minimize the total energy function, a greedy algorithm can be used. In each iteration, the control points are moved in the direction of decreasing energy. The iteration step size is set to 0.5 pixels to prevent the contour from deforming too quickly. When the change amount of the total energy function is less than the preset convergence threshold (such as 0.001), the iteration stops. Usually, convergence can be achieved when the number of iterations is between 50 - 100 times. Based on the final control point sequence, a smooth river channel water body boundary is reconstructed by methods such as B-spline interpolation or polynomial fitting to generate a river channel water area mask.

[0032] Experiments show that when this method is applied to river channel remote sensing images in different seasons, the average accuracy can reach 93.5%, which is 7.8 percentage points higher than the traditional threshold segmentation method, and the boundary smoothness is increased by 65%. It is especially suitable for river channel extraction tasks under complex backgrounds.

[0033] In an alternative embodiment, based on the drone oblique photogrammetry image data, a multi-view image matching method is used to generate dense point cloud data of the river channel area; the dense point cloud data is cropped using the river channel water mask, and the obtained point cloud data includes: A matching pair group is established for the drone oblique photogrammetry image data, the ratio of the overlapping area of each pair of images in the matching pair group to the total area of the images is calculated, and the image pairs with an overlapping ratio greater than a preset overlapping threshold are selected as the image pairs to be matched; For the reference image in the image pairs to be matched, the normalized cross-correlation coefficient between each pixel point in the reference image and the corresponding pixel point in other matching images is calculated, the weight coefficient is determined according to the baseline length of each matching image, and the normalized cross-correlation coefficient and the weight coefficient are weighted and summed to obtain the multi-baseline matching cost; Within the neighborhood window of each pixel point in the reference image, an adaptive weight is calculated based on the similarity of pixel gray values, and the adaptive weight and the multi-baseline matching cost are weighted and summed to obtain the aggregated matching cost; The optimal disparity value of the pixel point is determined using the aggregated matching cost, a collinearity equation system is established based on the optimal disparity value, and the collinearity equation system is solved by the least squares method to obtain the three-dimensional point coordinates, generating dense point cloud data; The dense point cloud data is transformed to the image plane coordinate system through a projection matrix, and the projected point cloud data is cropped according to the river channel water mask, and the point cloud data located within the river channel water mask range is retained as the point cloud data of the river channel area.

[0034] Drone oblique photogrammetry image data of the river channel area is collected, including images in the vertical and oblique directions. At the same time, a river channel water mask is generated through image segmentation technology for subsequent cropping of point cloud data.

[0035] After obtaining the drone oblique photogrammetry image data, a matching pair group is established for these data. Specifically, the ratio of the overlapping area of each pair of images in the matching pair group to the total area of the images is calculated. For example, for two images A and B, the overlapping area is calculated to be 500 square meters, the total area of image A is 1000 square meters, the total area of image B is 1200 square meters, and the average value of the two is taken as 1100 square meters, then the overlapping ratio is 500 / 1100≈0.45. The preset overlapping threshold is set to 0.3. Since 0.45 is greater than 0.3, this pair of images is selected as the image pair to be matched.

[0036] For the selected image pairs to be matched, dense point cloud generation is performed. In the image pairs to be matched, one is selected as the reference image, and the normalized cross-correlation coefficient between each pixel point in the reference image and the corresponding pixel point in other matching images is calculated. The normalized cross-correlation coefficient is calculated by comparing the similarity of the gray-scale distributions within the neighborhood window of the pixel points, and the window size is set to 11×11 pixels.

[0037] The weight coefficient is determined according to the baseline length of each matching image. The baseline length refers to the distance between the shooting positions of two cameras. For example, if the baseline length between the reference image and matching image A is 20 meters, between the reference image and matching image B is 10 meters, and between the reference image and matching image C is 30 meters, then the weight coefficients of these three matching images can be set to 0.3, 0.2, and 0.5. The normalized cross-correlation coefficient and the weight coefficient are weighted and summed to obtain the multi-baseline matching cost.

[0038] To improve the matching accuracy, an adaptive weight is calculated based on the similarity of pixel gray-scale values within the neighborhood window of each pixel point in the reference image. The neighborhood window size is set to 5×5 pixels. The higher the similarity of the pixel gray-scale values, the greater the corresponding adaptive weight. The adaptive weight and the multi-baseline matching cost are weighted and summed to obtain the aggregated matching cost.

[0039] The optimal disparity value of the pixel point is determined using the aggregated matching cost. For each pixel point in the reference image, within a certain disparity range (such as 0 - 64 pixels), the disparity value with the minimum aggregated matching cost is searched for as the optimal disparity value. In practical applications, dynamic programming or semi-global matching algorithms can be used to optimize the disparity search process and improve the robustness of the matching.

[0040] Based on the optimal disparity value, a collinearity equation set is established, that is, the relationship equation between the object point coordinates and the image point coordinates is established according to the internal and external parameters of the camera and the image point coordinates. The collinearity equation set is solved by the least squares method to obtain the three-dimensional point coordinates. For example, for a pixel point in the river channel area, its coordinates in the reference image are (1024, 768), the calculated optimal disparity value is 32 pixels, and the three-dimensional coordinates obtained by solving the collinearity equation are (521368.45, 3876215.78, 125.62). By performing the above processing on all pixel points in the reference image, dense point cloud data is generated.

[0041] After generating the dense point cloud data, it is transformed to the image plane coordinate system through the projection matrix. The projection matrix consists of the internal parameters of the camera (such as focal length, principal point coordinates, distortion coefficients) and external parameters (such as camera position, attitude). For example, assuming the camera focal length is 35 mm, the principal point coordinates are (1024, 768), there is no distortion, the camera position is (521000, 3876000, 200), and the attitude angles are (0°, 0°, 0°), then the corresponding projection matrix can be constructed to project the three-dimensional point cloud onto the two-dimensional image plane.

[0042] The projected point cloud data is cropped according to the river channel water area mask. The river channel water area mask is a binary image, where the pixels with a value of 1 represent the river channel water area, and the pixels with a value of 0 represent non-river channel areas. The projected point cloud data is compared with the river channel water area mask, and the point cloud data located within the range of the river channel water area mask (i.e., the area with a mask value of 1) is retained as the point cloud data of the river channel area.

[0043] For example, for a three-dimensional point, its projected coordinates on the image plane are (800, 600). If the value of the river channel water area mask at this position is 1, then this point is retained; if the mask value is 0, then this point is discarded. In this way, the point cloud data of the river channel area can be effectively extracted, providing data support for subsequent river channel analysis and monitoring.

[0044] After the above processing, the finally obtained point cloud data of the river channel area contains rich three-dimensional information and can be used for applications such as river channel cross-section analysis, water surface elevation measurement, and river channel deformation monitoring. In practical applications, according to different river channel characteristics and application requirements, the settings of various parameters, such as overlap threshold, window size, disparity search range, etc., can be adjusted to obtain the best point cloud extraction effect.

[0045] Figure 2 This is the bar chart for the performance comparison and analysis of the multi-view image matching method in the embodiments of the present invention: The image shows the performance comparison of three different matching methods in various scenarios: the traditional NCC matching method, the multi-baseline weighted matching method, and the adaptive weight aggregation matching method. From different scenarios, in the low-texture area, the matching success rates of the three methods are 34.0%, 42.0%, and 50.0% respectively; in the specular reflection area on the water surface, the matching success rates are 28.0%, 38.0%, and 44.0% respectively; in the vegetation-covered area, the matching success rates are 32.0%, 40.0%, and 46.0% respectively. From the average level of the whole scenario, the overall matching success rates of the three methods are 34.0%, 42.0%, and 48.0% respectively. The data shows that the adaptive weight aggregation matching method has achieved the best results in all scenarios, with an average increase of 14 percentage points compared to the traditional NCC matching method and an average increase of 6 percentage points compared to the multi-baseline weighted matching method. Especially in the low-texture area, the advantage of the adaptive weight aggregation matching method is the most obvious, reaching a matching success rate of 50.0%.

[0046] In an optional implementation manner, the optimal disparity value of a pixel point is determined by using the aggregated matching cost, a collinearity equation set is established based on the optimal disparity value, and the three-dimensional point coordinates are obtained by solving the collinearity equation set through the least squares method. Generating the dense point cloud data includes: Construct a disparity space cost curve within a preset disparity search range based on the aggregated matching cost, determine the integer-pixel extreme points of the disparity space cost curve, and perform quadratic curve fitting using the cost values at the integer-pixel extreme points and their adjacent disparity positions to obtain the optimal disparity value; Construct a collinearity equation set according to the optimal disparity value and the internal and external orientation elements of the image, solve the collinearity equation set to obtain an initial three-dimensional coordinate point set, and construct a spatial point cloud based on the initial three-dimensional coordinate point set; Reproject each three-dimensional coordinate point in the spatial point cloud to the images of each view, calculate the reprojection error between the projection point and the original matching point, and mark the three-dimensional coordinate points with a reprojection error greater than the preset error threshold as points to be optimized; For the points to be optimized, extract a local point set within their neighborhood range, calculate the centroid coordinates of the local point set, and update and optimize the three-dimensional coordinates of the points to be optimized based on the centroid coordinates to obtain an optimized three-dimensional coordinate point set; Perform dense point cloud reconstruction on the optimized three-dimensional coordinate point set, calculate the point density of each point in the reconstructed point cloud within its preset neighborhood range, and perform local encryption reconstruction on the area where the point density is lower than the preset density threshold to obtain the final dense point cloud data.

[0047] When constructing a disparity space cost curve based on the aggregated matching cost within a preset disparity search range, the search can be conducted within a disparity range of 0 to 64 pixels. The matching cost values are arranged in the order of disparity values to form a cost curve. By analyzing this curve, the integer-pixel extreme points can be determined. For example, when the aggregated matching cost of a certain pixel point is the smallest at a disparity value of 28 pixels and the value is 0.15, then this point is an integer-pixel extreme point. Subsequently, a quadratic curve equation is constructed using the cost values at this extreme point and its adjacent disparity positions 27 and 29 (0.18 and 0.17 respectively). By solving the extreme point of this quadratic equation, a more accurate disparity value, such as 28.12 pixels, can be obtained as the optimal disparity value for this pixel point.

[0048] When constructing a collinearity equation system based on the obtained optimal disparity value and the interior and exterior orientation elements of the image, known camera internal parameters (such as focal length f = 35mm, principal point coordinates x p = y p = 0mm) and external parameters (such as projection center coordinates X s = 108.25m, Y s = 215.36m, Z s = 1520.85m, rotation matrix elements r 11 = 0.9986, r 12 = 0.0012, etc.) are used. For each image point, a collinearity relationship from the image point to the object point is established to form an equation system. For example, for the image point coordinates (x = 12.5mm, y = 8.2mm) and its corresponding disparity value of 28.12 pixels, the initial three-dimensional coordinates (X = 156.32m, Y = 278.45m, Z = 32.16m) are calculated through the collinearity equation system. This process is repeated for all image points to form a set of initial three-dimensional coordinate points and construct a spatial point cloud.

[0049] When re-projecting each three-dimensional coordinate point in the spatial point cloud to the images of each perspective, using the same interior and exterior orientation elements, the three-dimensional coordinate point (X = 156.32m, Y = 278.45m, Z = 32.16m) is projected back to the original image, and the projected point coordinates (x' = 12.48mm, y' = 8.23mm) are obtained. Calculate the reprojection error between the projected point and the original matching point, that is . Set the preset error threshold to 0.1mm. If the reprojection error is greater than this threshold, then this point is marked as a point to be optimized.

[0050] For the points to be optimized, local point sets are extracted within their neighborhood ranges. All points within a spherical region centered at this point with a radius of 2 meters can be taken as the local point sets. For example, for the three-dimensional coordinate point (X = 186.45m, Y = 321.78m, Z = 28.92m), its reprojection error is 0.15mm, exceeding the threshold of 0.1mm, and it is marked as a point to be optimized. A local point set containing 52 points is extracted within its 2-meter neighborhood range, and the centroid coordinates are calculated as (X c = 186.38m, Y c = 321.82m, Z c = 28.95m). The points to be optimized are weighted and updated based on the centroid coordinates. The weight depends on the distance from the point to the centroid. The updated three-dimensional coordinates become (X' = 186.40m, Y' = 321.80m, Z' = 28.94m), and the reprojection error drops to 0.08mm, which is lower than the threshold, completing the optimization.

[0051] When performing dense point cloud reconstruction on the optimized three-dimensional coordinate point set, calculate the point density of each point in the reconstructed point cloud within a preset neighborhood range. Set the preset neighborhood range to 1 cubic meter, and the preset density threshold to 10 points per cubic meter. For example, for the three-dimensional point (X = 245.67m, Y = 418.93m, Z = 35.24m), calculate the number of points within its 1-cubic-meter neighborhood to be 8, which is lower than the density threshold of 10 points / cubic meter. Mark this area as an area that needs to be encrypted and reconstructed. In this area, by increasing the matching point density on the original image, such as from taking one matching point per 5 pixels originally to taking one matching point per 2 pixels, re-perform the above steps to generate more three-dimensional points. After encrypted reconstruction, the point density in this area increases to 15 points / cubic meter, meeting the density requirements.

[0052] Through the above steps, dense point cloud data generation for the entire image area is completed. The final point cloud data contains approximately 5 million three-dimensional points, with an average point density reaching 25 points per square meter, the mean reprojection error being 0.06mm, and the standard deviation being 0.03mm, meeting the requirements for high-precision three-dimensional reconstruction. The generated dense point cloud data can be used for subsequent applications such as digital surface model generation and three-dimensional model reconstruction.

[0053] In an alternative embodiment, recursively multi-scale enhancing the sparse point cloud region using the density compensation sampling matrix, and optimizing the spatial position in combination with local curvature constraints to generate density-balanced riverbed point cloud data includes: Construct a multi-scale feature extraction network. Based on the multi-scale feature extraction network, extract features from multiple preset radius neighborhoods of each point in the point cloud data, perform weighted fusion on the extracted features and corresponding weight coefficients to obtain feature vectors, input the feature vectors into a density prediction network to obtain predicted density values, and multiply the predicted density values by the initial density field to obtain the density prediction result for each point; Determine the sparse regions in the point cloud based on the density prediction results, generate landform feature weights, divide the preset target density value by the density prediction results and multiply by the landform feature weights to obtain regional sampling weights, and adaptively fuse sonar sounding data and multi-temporal multi-spectral remote sensing image data based on the regional sampling weights to generate enhanced point cloud data; Construct a local connection graph for the enhanced point cloud data, iteratively update the node features in the local connection graph to obtain optimized node features; calculate the distance weights of each point within a preset radius neighborhood based on the optimized node features, and perform weighted summation of the distance weights and the projection distances from the point to its neighborhood points to obtain local curvature values; Input the optimized node features and the local curvature values into a generator to generate an initial density-balanced point cloud, and use the generator to regenerate point cloud data with local geometric features; construct gradient information of the flow velocity field and water depth field based on hydrological monitoring data, and fuse the gradient information with the point cloud data to optimize the spatial position to generate density-balanced riverbed point cloud data.

[0054] As Figure 3 shown, the method further includes: Construct a multi-scale feature extraction network, which includes three convolutional layers with convolutional kernel sizes of 32, 64, and 128 respectively, for capturing local and global features of the point cloud. For each point in the point cloud, select points within the neighborhood ranges with preset radii of 0.5 meters, 1.0 meters, and 2.0 meters respectively for feature extraction. During the feature extraction process, use max pooling operations to aggregate neighborhood information, and set the dimension of the feature vector extracted at each scale to 128 dimensions. Perform weighted fusion on the features extracted at three different scales with corresponding weight coefficients of 0.3, 0.3, and 0.4 to obtain a 384-dimensional feature vector.

[0055] This feature vector is then input into a density prediction network composed of three fully connected layers. The number of neurons in the fully connected layers is 256, 128, and 1 respectively. The activation function uses the ReLU function, and the last layer uses the Sigmoid function to map the output to between 0 and 1 to obtain the predicted density value. Multiply the predicted density value by the initial density field (the initial density field is obtained by calculating the number of points within a 1.5-meter radius neighborhood of each point) to obtain the density prediction result for each point.

[0056] Determine the sparse regions in the point cloud based on the density prediction results. The specific method is to set the density threshold to 20 points per cubic meter, and the regions below this threshold are marked as sparse regions. Meanwhile, generate the geomorphic feature weights, which are calculated based on the curvature and slope information of the point cloud. The curvature is obtained by fitting a quadratic surface to the local neighborhood, and the slope is obtained by calculating the angle between the normal vector of the local plane and the vertical direction. After normalizing the curvature value and the slope value, they are weighted with a ratio of 0.6 and 0.4 to obtain the geomorphic feature weights. Divide the preset target density value (set to 50 points per cubic meter) by the density prediction result, and then multiply it by the geomorphic feature weights to obtain the regional sampling weights.

[0057] Perform adaptive fusion on the sonar sounding data and multi-temporal multi-spectral remote sensing image data based on this sampling weight. The sonar data provides accurate depth information with a weight coefficient of 0.7; the remote sensing image provides planar position information through edge detection and texture analysis with a weight coefficient of 0.3. The fusion process uses a spatial index structure based on an octree, and fuses the sonar point cloud and the point cloud derived from the remote sensing image in the overlapping area through weights, and directly merges them in the non-overlapping area to generate enhanced point cloud data.

[0058] Construct a local connection graph for the enhanced point cloud data. Each node in the connection graph connects to its neighboring points within a range of 1.2 meters, and on average each node connects to about 12 neighbors. Use a graph convolutional network to iteratively update the node features in the connection graph for three rounds. In each round of iteration, aggregate the information of the neighboring nodes and update the features of the central node. The update formula takes into account the distance weights of the edges, and the closer neighbors have a greater impact. After the iteration is completed, the optimized node features are obtained, and the feature dimension is 256.

[0059] Based on the optimized node features, calculate the distance weights of each point within a neighborhood with a radius of 1.0 meter. The distance weights are calculated using a Gaussian kernel function with a standard deviation set to 0.3 meters. Perform a weighted sum of the distance weights and the projection distance from the point to the neighboring points. The projection distance refers to the distance from the point to the best-fitting plane formed by the neighboring points to obtain the local curvature value.

[0060] Concatenate the optimized node features and the local curvature values and input them into the generator. The generator adopts a three-layer encoder-decoder structure. The encoder maps the input features to the latent space, and the decoder reconstructs the point cloud coordinates from the latent space. The generator outputs an initial density-balanced point cloud with a point cloud density reaching 45 - 55 points per cubic meter. Use the same generator to regenerate the point cloud data with local geometric features, and introduce random noise (with an amplitude of 0.05 meters) during the regeneration process to increase the natural variation of the point cloud.

[0061] Construct the gradient information of the flow velocity field and water depth field based on hydrological monitoring data. The flow velocity field data comes from the measurement of an acoustic Doppler current profiler (ADCP) with a sampling interval of 5 meters; the water depth field data comes from an echo sounder with a sampling interval of 2 meters. The measured data is used to construct continuous flow velocity field and water depth field through trilinear interpolation, and the gradient information of the field is calculated to represent the changing trends of water flow and terrain.

[0062] Fuse the gradient information with the point cloud data, and slightly adjust the positions of the points along the gradient direction. The adjustment amplitude is proportional to the gradient magnitude, and the maximum adjustment does not exceed 0.3 meters to maintain the continuity and authenticity of the terrain. After adjustment, the final riverbed point cloud data with balanced density is generated, and the standard deviation of the point cloud density distribution is reduced to less than 30% of the original data, which can accurately express the micro-topography features of the riverbed.

[0063] In an alternative embodiment, a triangular grid terrain model is constructed based on the riverbed point cloud data, and a sequence of river channel cross-sections is generated at preset intervals along the river channel flow direction; for the sequence of river channel cross-sections, the hydrodynamic correlation intensity is calculated through a recurrent neural network, and the construction of the cross-section evolution propagation sequence based on the hydrodynamic correlation intensity includes: Perform constrained Delaunay triangulation on the riverbed point cloud data to construct an initial triangular grid. Calculate the grid quality parameter based on the ratio of the triangle area to the side length, and optimize the triangular grid according to the grid quality parameter to obtain a triangular grid terrain model; Construct a terrain elevation distribution based on the triangular grid terrain model, calculate the elevation difference between adjacent grid nodes, determine the water flow direction field based on the elevation difference, sample at preset intervals along the water flow direction field on the riverbed terrain triangular grid model to generate a sequence of river channel cross-sections, and extract the width, water depth, cross slope, and roughness parameters of each cross-section to construct a cross-section feature vector; Input the cross-section feature vector into a long short-term memory network. Extract the cross-section time series features through the memory unit and hidden state of the long short-term memory network. Calculate the cosine similarity between adjacent cross-sections based on the hidden state, and multiply the cosine similarity by the exponential decay function of the cross-section spacing to obtain the cross-section correlation intensity; construct a propagation weight matrix based on the cross-section correlation intensity and the hydraulic gradient factor, and use the propagation weight matrix to iteratively update the cross-section state to obtain the river channel cross-section evolution propagation sequence.

[0064] The constrained Delaunay triangulation algorithm is used to construct the initial triangular mesh. By considering the spatial position relationships of each point in the point cloud, this algorithm ensures that the generated triangles satisfy the Delaunay condition, that is, no other points are contained within the circumcircle of any triangle. In practical applications, for riverbed point cloud data containing 500,000 points, the system uses a quadtree spatial index structure to divide the space into grid cells of size 5m×5m, which improves the point query efficiency and reduces the time complexity of triangulation from O(n 2 ) to O(nlogn).

[0065] After constructing the initial triangular mesh, the system calculates mesh quality parameters to evaluate the mesh quality. The mesh quality parameter Q is defined as the ratio of the triangle area to the product of the squares of the three side lengths. For an ideal equilateral triangle, the Q value is close to 0.5; while for a slender triangle, the Q value is close to 0. In this embodiment, the system sets the quality threshold to 0.2, and optimizes the triangles with Q values lower than the threshold. The optimization process adopts two strategies: edge flipping and point insertion. When two adjacent triangles form a convex quadrilateral, the system determines whether flipping the shared edge can improve the overall quality; when the interior angle of a triangle is less than 20 degrees, a new point is inserted inside it and re-triangulation is performed. After 3 rounds of iterative optimization, the average value of the mesh quality parameter is increased from 0.35 to 0.42, forming a high-quality triangular mesh terrain model.

[0066] Based on the optimized triangular mesh terrain model, the system constructs the terrain elevation distribution and calculates the elevation differences between adjacent grid nodes. The elevation difference calculation method is the difference between the elevation values of two adjacent nodes divided by the horizontal distance between the two points. In a practical case, the riverbed terrain triangular mesh contains 78,562 nodes and 156,934 triangles. The system calculates the elevation differences between each node and its adjacent nodes, and determines the water flow direction field based on this. The water flow direction is defined as the direction of the steepest elevation drop, represented by the vector pointing from each node to its lowest adjacent node.

[0067] When generating a sequence of river cross-sections along the water flow direction field, the system sets the cross-section spacing to 50 meters, and a total of 200 cross-sections are generated on a 10-kilometer-long river. For each cross-section position, the system intersects with the triangular mesh model in a plane perpendicular to the water flow direction to obtain the cross-section contour. For each cross-section, the system extracts the cross-section width, average water depth, left and right bank slopes, and cross-section roughness parameters.

[0068] The cross-section width is defined as the horizontal distance from the leftmost point to the rightmost point of the cross-section; the average water depth is the difference between the average of the elevations of all points on the cross-section and the elevation of the lowest point; the left and right bank slopes are the slopes of the lines connecting the bank edge points and the lowest point of the riverbed; the roughness parameter is calculated from the undulation degree of the cross-section morphology and the distribution of the bed material. In an example, the extracted cross-section feature vector includes a cross-section width of 102 m, an average water depth of 2.8 m, a left bank slope of 0.32, a right bank slope of 0.28, and a roughness coefficient of 0.035.

[0069] The extracted cross-section feature vector is input into a long short-term memory network for processing. The network consists of an input layer, a forget gate, an input gate, an output gate, and a hidden layer. The network input is a sequence of cross-section feature vectors arranged in time series, and the cross-section time series features are extracted through a three-layer LSTM structure, with each layer containing 128 neurons. The system uses 10 cross-sections as a time window for sliding training, sets the learning rate to 0.001, the batch size to 32, and trains for 500 iterations. The memory unit of the network stores long-term dependence information, and the hidden state represents the comprehensive features of the current cross-section.

[0070] Based on the hidden state of the LSTM network, the cosine similarity between adjacent cross-sections is calculated. For each pair of adjacent cross-sections i and j, the dot product of their hidden state vectors is divided by the product of the norms of the two vectors to obtain the cosine similarity value. Then, the cosine similarity is multiplied by an exponential decay function of the cross-section spacing, with the decay coefficient set to 0.05, to obtain the cross-section correlation strength matrix. In the example, the average correlation strength between adjacent cross-sections is 0.78, and the correlation strength decreases exponentially as the cross-section spacing increases.

[0071] A propagation weight matrix is constructed based on the cross-section correlation strength and the hydraulic gradient factor. The hydraulic gradient factor takes into account the riverbed slope and the cross-section morphology change, and is calculated by the ratio of the average elevation difference of the cross-section to the cross-section spacing. The propagation weight matrix defines the influence weights of each adjacent cross-section in the cross-section state update process. The system uses an iterative method to update the cross-section state, sets the number of iterations to 30, and the convergence threshold to 0.001. Through the iterative update of the propagation weight matrix, a cross-section evolution propagation sequence reflecting the dynamic change process of the river channel is obtained, and this sequence can be used for river channel evolution analysis and hydrodynamic simulation.

[0072] Figure 4 Schematic diagram for comparing the optimization effect of triangular mesh quality in the embodiment of the present invention: This figure shows the comparison of network quality indicators for three different technical solutions at different grid density levels. The horizontal axis in the figure represents the grid density levels from extremely low to extremely high, and the vertical axis represents the network quality indicators. This technical solution (circular markers) shows optimal performance at each density level: starting from 0.16 at extremely low density, passing through low density (0.30), medium-low density (0.40), medium density (0.50), medium-high density (0.58), high density (0.64) to extremely high density (0.68), showing a stable upward trend. The triangulation scheme based on edge constraints (diamond markers) comes next, gradually increasing from 0.12 at extremely low density to 0.52 at extremely high density. The traditional Delaunay triangulation (square markers) performs the worst, only increasing from 0.08 at extremely low density to 0.34 at extremely high density. The data shows that the network quality indicator of this technical solution reaches 0.68 under extremely high density conditions, which is 0.34 units higher than the traditional method and 0.16 units higher than the method based on edge constraints, demonstrating significant technical advantages.

[0073] In an alternative embodiment, based on the cross-section evolution propagation sequence and historical evolution data, the riverbed scouring and silting change rate is obtained through recursive calculation; when the riverbed scouring and silting change rate exceeds a preset rate threshold, a warning message is generated, and the dynamic prediction result of the river channel cross-section morphology and the warning message are visually displayed, including: Perform multi-source data fusion on the cross-section evolution propagation sequence and historical evolution data to obtain a fusion feature vector, construct a feature weight and quality evaluation function based on the fusion feature vector, and perform quality evaluation on the fusion feature vector through the feature weight and the quality evaluation function to obtain an evaluation result; Construct a graph structure based on the evaluation result, construct the cross-section data as a node set, construct the association between cross-sections as an association edge set, generate an adjacency matrix based on the node set and the association edge set, and use the adjacency matrix to update the feature information of adjacent cross-sections through message passing to obtain a hidden state feature; Input the hidden state feature into a mapping function to obtain the current change rate, perform weighted fusion of the current change rate and the historical change rate to obtain the scouring and silting change rate of the target cross-section, calculate statistical parameters based on the scouring and silting change rate, and perform exponential mapping of the statistical parameters and a reference threshold to obtain a dynamic warning threshold; Compare the scouring and silting change rate with the dynamic warning threshold to obtain a warning level, perform time-series integration on the target cross-section based on the scouring and silting change rate to obtain a predicted morphology, and perform three-dimensional visual display of the predicted morphology and the warning level to generate an interactive warning message.

[0074] Based on the cross-section evolution propagation sequence and historical evolution data, the riverbed erosion and deposition change rate is obtained through recursive calculation. When the riverbed erosion and deposition change rate exceeds the preset rate threshold, a warning message will be generated, and the dynamic prediction results of the river channel cross-section morphology and the warning message will be visually displayed.

[0075] Multi-source data fusion is performed on the cross-section evolution propagation sequence and historical evolution data to obtain a fused feature vector. The multi-source data includes hydrological station monitoring data, cross-section measurement data, remote sensing image data, etc. Taking hydrological stations as an example, assume there are 6 hydrological stations, and each station contains 10 features such as flow rate, sediment concentration, water level, etc., then a 6×10 feature matrix is formed. At the same time, temporal features are constructed through historical evolution data. For example, the erosion and deposition change amounts of a certain cross-section in the past 5 years are +0.2 m, -0.15 m, +0.1 m, -0.05 m, +0.3 m respectively. These heterogeneous data are converted into vectors of a unified dimension through feature extraction, and then a weighted average method is used for fusion to obtain a fused feature vector with a dimension of 128.

[0076] Based on the fused feature vector, a feature weight and a quality evaluation function are constructed. The feature weight is determined according to the contribution degree of each feature to the erosion and deposition change. For example, the weight of the flow rate is 0.4, the weight of the sediment concentration is 0.3, the weight of the water level is 0.2, and the weights of other features are 0.1. The quality evaluation function comprehensively considers data integrity, timeliness, and consistency, and conducts a comprehensive score by setting different evaluation indicators. For example, the data completeness rate of 95% gets a score of 0.95, the data update time within 3 days gets a score of 0.9, the data consistency test passing rate of 98% gets a score of 0.98, and the overall score is 0.95×0.4 + 0.9×0.3 + 0.98×0.3 = 0.947.

[0077] The quality of the fused feature vector is evaluated through the feature weight and the quality evaluation function to obtain an evaluation result. The evaluation result includes a reliability index, an uncertainty range, etc. For example, the reliability index of a certain cross-section data is 0.92, indicating that the data has a credibility of 92%; the uncertainty range is ±0.05 m, indicating that there is an error of ±0.05 m in the prediction result.

[0078] Based on the evaluation result, a graph structure is constructed. The cross-section data is constructed as a node set, and the association between cross-sections is constructed as an association edge set. Assume there are 10 cross-section points, then 10 nodes are formed. Association edges are established between adjacent cross-sections. For example, there are association edges between cross-section 1 and cross-section 2, and between cross-section 2 and cross-section 3, while cross-section 1 and cross-section 3 are not directly connected. Based on the node set and the association edge set, an adjacency matrix is generated. The adjacency matrix is a 10×10 matrix, where the adjacency matrix A[i][j]=1 indicates that cross-section i and cross-section j are adjacent, and A[i][j]=0 indicates that they are not adjacent.

[0079] The hidden state features are obtained by message passing and updating the characteristic information of adjacent cross-sections using the adjacency matrix. During the message passing process, each node aggregates its own features with the features of adjacent nodes through weighting. For example, the updated features of cross-section 2 are composed of 0.6 weight of its own features, 0.2 weight of the features of cross-section 1, and 0.2 weight of the features of cross-section 3. Through 3 rounds of iterative updates, the hidden state features containing global information are finally obtained.

[0080] The hidden state features are input into a mapping function to obtain the current change rate. The mapping function uses a piecewise linear function to convert the hidden state features into the scouring and silting rate. For example, when the eigenvalue is in the range of [-0.5, 0.5], it is mapped to a change rate of [-0.1 m / year, 0.1 m / year]; when the eigenvalue is greater than 0.5, it is mapped to a change rate of 0.1 + 0.2×(eigenvalue - 0.5) m / year.

[0081] The current change rate and the historical change rate are weighted and fused to obtain the scouring and silting change rate of the target cross-section. The weight of the current change rate is 0.7, and the weight of the historical change rate is 0.3. For example, if the currently calculated change rate of a certain cross-section is 0.15 m / year and the historical average change rate is 0.1 m / year, then the fused change rate is 0.15×0.7 + 0.1×0.3 = 0.135 m / year.

[0082] Statistical parameters are calculated based on the scouring and silting change rate, including the average value, standard deviation, etc. For example, the average scouring and silting change rate of 10 cross-sections in a certain river section is 0.08 m / year, and the standard deviation is 0.03 m / year. The statistical parameters are exponentially mapped with a reference threshold to obtain a dynamic warning threshold. The reference threshold is set to 0.2 m / year and is mapped to a dynamic threshold through an exponential mapping function. For example, when the river is in the flood season, the dynamic threshold is adjusted to 0.2×1.5 = 0.3 m / year; when the river is in the dry season, the dynamic threshold is adjusted to 0.2×0.8 = 0.16 m / year.

[0083] The scouring and silting change rate is compared with the dynamic warning threshold to obtain the warning level. The warning level is divided into four levels: normal, attention, warning, and danger. When the change rate is less than 50% of the threshold, it is at the normal level; when the change rate is between 50% and 80% of the threshold, it is at the attention level; when the change rate is between 80% and 100% of the threshold, it is at the warning level; when the change rate exceeds the threshold, it is at the danger level. For example, if the scouring and silting change rate of a certain cross-section is 0.25 m / year and the dynamic threshold is 0.3 m / year, and the change rate accounts for 83.3% of the threshold, it is determined to be at the warning level.

[0084] The predicted morphology is obtained by performing time-series integration on the target cross-section based on the erosion and deposition change rate. The time-series integration adopts a step-by-step accumulation method to predict the future cross-section morphology. For example, if the current elevation of a certain cross-section is 100 meters and the erosion and deposition change rate is 0.135 meters per year, the predicted elevation after one year is 100 + 0.135 = 100.135 meters, and the predicted elevation after two years is 100.135 + 0.135 = 100.27 meters.

[0085] The predicted morphology and the warning level are visually displayed in three dimensions to generate interactive warning information. The three-dimensional visualization includes a plan view, a cross-sectional view, and a time-series change graph. The plan view shows the river channel layout and the warning area; the cross-sectional view shows the current and predicted morphologies of the cross-section; the time-series change graph shows the changing trend of the cross-section elevation over time. The warning information is displayed through color coding, green for normal, yellow for attention, orange for warning, and red for danger. Users can view the predicted results at different time points through the interactive interface and receive warning prompt information.

[0086] Figure 5 This is a bar chart comparing the performance of the riverbed erosion and deposition change monitoring and warning system according to the embodiments of the present invention: This figure shows a comparative analysis of the warning performance of three different technical solutions, including three performance indicators: warning accuracy, response time, and change prediction accuracy. In terms of warning accuracy, the technical solution of the present invention reaches 48.0%, significantly better than 36.0% of the traditional statistical adjustment method and 40.0% of the simple neural network method; in terms of the response time index, the technical solution of the present invention performs best with an advantage of 52.0%, while the traditional statistical adjustment method is 32.0% and the simple neural network method is 36.0%; in terms of the change prediction accuracy, the technical solution of the present invention also maintains a leading advantage, reaching 44.0%. In contrast, the traditional statistical adjustment method is 28.0% and the simple neural network method is 32.0%. From the overall performance, the technical solution of the present invention has significant advantages in all performance indicators, leading the other two methods by about 12 - 16 percentage points on average, especially in the key indicator of response time, where the advantage is more obvious, reflecting the outstanding value of this solution in practical applications.

[0087] In the second aspect of the embodiments of the present invention, a dynamic change monitoring system for multi-temporal river channel images is provided, including: A first unit for acquiring multi-temporal multi-spectral remote sensing image data and unmanned aerial vehicle (UAV) oblique photogrammetry image data of the river channel area; based on the multi-temporal multi-spectral remote sensing image data, identifying the river channel water body boundary through a support vector machine classification algorithm and constructing a river channel water area mask; A second unit for generating dense point cloud data of the river channel area by using a multi-view image matching method based on the UAV oblique photogrammetry image data; cropping the dense point cloud data by using the river channel water area mask to obtain point cloud data; A third unit is configured to calculate a point cloud neighborhood density distribution function for the point cloud data, determine local importance weights according to the gradient change of the point cloud neighborhood density distribution function, and establish a density compensation sampling matrix based on the local importance weights; use the density compensation sampling matrix to perform recursive multi-scale enhancement on the sparse area of the point cloud, and optimize the spatial position in combination with local curvature constraints to generate density-balanced riverbed point cloud data; A fourth unit is configured to construct a triangular grid terrain model based on the riverbed point cloud data, and generate a sequence of river channel cross-sections at preset intervals along the river flow direction; calculate the hydrodynamic correlation intensity for the sequence of river channel cross-sections through a recurrent neural network, and construct a cross-section evolution propagation sequence according to the hydrodynamic correlation intensity; A fifth unit is configured to obtain the rate of riverbed erosion and deposition change through recursive calculation based on the cross-section evolution propagation sequence and historical evolution data; when the rate of riverbed erosion and deposition change exceeds a preset rate threshold, generate a warning message, and visually display the dynamic prediction result of the river channel cross-section morphology and the warning message.

[0088] In a third aspect of the embodiments of the present invention, an electronic device is provided, including: A processor; A memory for storing instructions executable by the processor; Wherein, the processor is configured to call the instructions stored in the memory to execute the method described above.

[0089] In a fourth aspect of the embodiments of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored, and when the computer program instructions are executed by a processor, the method described above is implemented.

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

[0091] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that: they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for dynamically monitoring the changes of a river channel based on multi-temporal river channel images, characterized in that, Including: Obtain multi-temporal multi-spectral remote sensing image data and unmanned aerial vehicle (UAV) oblique photogrammetry image data of the river channel area; Based on the multi-temporal multi-spectral remote sensing image data, identify the river channel water body boundary through the support vector machine classification algorithm, and construct a river channel water area mask; Based on the UAV oblique photogrammetry image data, use the multi-view image matching method to generate dense point cloud data of the river channel area; use the river channel water area mask to crop the dense point cloud data to obtain point cloud data; Calculate the point cloud neighborhood density distribution function for the point cloud data, determine the local importance weight according to the gradient change of the point cloud neighborhood density distribution function, and establish a density compensation sampling matrix based on the local importance weight; use the density compensation sampling matrix to perform recursive multi-scale enhancement on the point cloud sparse area, and optimize the spatial position in combination with the local curvature constraint to generate density-balanced riverbed point cloud data; Construct a triangular grid terrain model based on the riverbed point cloud data, and generate a river channel cross-section sequence at preset intervals along the river channel flow direction; calculate the hydrodynamic correlation intensity for the river channel cross-section sequence through a recurrent neural network, and construct a cross-section evolution propagation sequence according to the hydrodynamic correlation intensity; Based on the cross-section evolution propagation sequence and historical evolution data, obtain the riverbed erosion and deposition change rate through recursive calculation; when the riverbed erosion and deposition change rate exceeds the preset rate threshold, generate a warning message, and visually display the dynamic prediction result of the river channel cross-section morphology and the warning message.

2. The method according to claim 1, wherein Based on the multi-temporal multi-spectral remote sensing image data, identifying the river channel water body boundary through the support vector machine classification algorithm and constructing a river channel water area mask includes: Construct a feature vector based on the multi-temporal multi-spectral remote sensing image data, map the feature vector to a high-dimensional feature space, construct an optimal classification hyperplane in the high-dimensional feature space, calculate the distance from the sample point to the hyperplane based on the optimal classification hyperplane, use the sample points with a distance greater than the preset distance threshold as training support vectors, and construct a classification decision function using the training support vectors; classify and predict the feature vector based on the classification decision function, and generate an initial classification image according to the prediction result; Perform morphological processing on the initial classification image, use the first structural element to perform opening operation to remove noise points, use the second structural element to perform closing operation to fill holes, and obtain an initial water area mask image; Use a contour evolution model with boundary curvature constraint to smooth the initial water area mask image, initialize the control point sequence of the contour curve based on the initial water area mask image, construct an internal energy function representing the elasticity and rigidity of the contour curve according to the first-order derivative and second-order derivative of the control point sequence, and construct an external energy function based on the gradient information of the initial classification image; Combine the internal energy function and the external energy function to obtain the total energy function of the contour evolution model, perform minimization solution on the total energy function, stop the iteration when the change amount of the total energy function is less than the preset convergence threshold, and reconstruct a smooth river channel water body boundary based on the final control point sequence to generate a river channel water area mask.

3. The method according to claim 1, wherein Based on the drone oblique photogrammetry image data, a multi-view image matching method is used to generate dense point cloud data of the river channel area; the dense point cloud data is cropped by using the river channel water mask, and the obtained point cloud data includes: A matching pair group is established for the drone oblique photogrammetry image data, the ratio of the overlapping area of each pair of images in the matching pair group to the total area of the images is calculated, and the image pairs with an overlapping ratio greater than a preset overlapping threshold are selected as the image pairs to be matched; For the reference image in the image pairs to be matched, the normalized cross-correlation coefficient between each pixel point in the reference image and the corresponding pixel points in other matching images is calculated, the weight coefficient is determined according to the baseline length of each matching image, and the normalized cross-correlation coefficient and the weight coefficient are weighted and summed to obtain the multi-baseline matching cost; Within the neighborhood window of each pixel point in the reference image, an adaptive weight is calculated based on the similarity of pixel gray values, and the adaptive weight and the multi-baseline matching cost are weighted and summed to obtain the aggregated matching cost; The optimal disparity value of the pixel point is determined by using the aggregated matching cost, a collinearity equation set is established based on the optimal disparity value, and the collinearity equation set is solved by the least squares method to obtain the three-dimensional point coordinates, and dense point cloud data is generated; The dense point cloud data is transformed to the image plane coordinate system through a projection matrix, and the projected point cloud data is cropped according to the river channel water mask, and the point cloud data located within the river channel water mask range is retained as the point cloud data of the river channel area.

4. The method according to claim 3, wherein The optimal disparity value of the pixel point is determined by using the aggregated matching cost, a collinearity equation set is established based on the optimal disparity value, and the collinearity equation set is solved by the least squares method to obtain the three-dimensional point coordinates, and the generated dense point cloud data includes: Based on the aggregated matching cost, a disparity space cost curve is constructed within a preset disparity search range, the integer pixel extreme points of the disparity space cost curve are determined, and the optimal disparity value is obtained by performing quadratic curve fitting on the integer pixel extreme points and the cost values at their adjacent disparity positions; According to the optimal disparity value and the interior and exterior orientation elements of the image, a collinearity equation set is constructed, the collinearity equation set is solved to obtain an initial three-dimensional coordinate point set, and a spatial point cloud is constructed based on the initial three-dimensional coordinate point set; Each three-dimensional coordinate point in the spatial point cloud is re-projected onto the images of each view, the reprojection error between the projected point and the original matching point is calculated, and the three-dimensional coordinate points with a reprojection error greater than a preset error threshold are marked as points to be optimized; For the points to be optimized, a local point set is extracted within their neighborhood range, the centroid coordinates of the local point set are calculated, and the three-dimensional coordinates of the points to be optimized are updated and optimized based on the centroid coordinates to obtain an optimized three-dimensional coordinate point set; Dense point cloud reconstruction is performed on the optimized three-dimensional coordinate point set, the point density of each point in the reconstructed point cloud within its preset neighborhood range is calculated, and the areas with a point density lower than the preset density threshold are locally encrypted and reconstructed to obtain the final dense point cloud data.

5. The method according to claim 1, characterized in that, Recursively multi-scale enhancing the sparse area of the point cloud using the density compensation sampling matrix, and optimizing the spatial position in combination with local curvature constraints to generate density-balanced riverbed point cloud data includes: Construct a multi-scale feature extraction network, extract features from multiple preset radius neighborhoods of each point in the point cloud data based on the multi-scale feature extraction network, perform weighted fusion of the extracted features with corresponding weight coefficients to obtain feature vectors, input the feature vectors into a density prediction network to obtain predicted density values, and multiply the predicted density values by the initial density field to obtain the density prediction results for each point; Determine the sparse area in the point cloud based on the density prediction results, generate geomorphic feature weights, divide the preset target density value by the density prediction results and multiply by the geomorphic feature weights to obtain regional sampling weights, and adaptively fuse sonar sounding data and multi-temporal multi-spectral remote sensing image data based on the regional sampling weights to generate enhanced point cloud data; Construct a local connection graph for the enhanced point cloud data, iteratively update the node features in the local connection graph to obtain optimized node features; calculate the distance weights of each point within a preset radius neighborhood based on the optimized node features, and perform weighted summation of the distance weights and the projection distance from the point to the neighborhood points to obtain local curvature values; Input the optimized node features and the local curvature values into a generator to generate initial density-balanced point clouds, and use the generator to regenerate point cloud data with local geometric features; construct gradient information of the flow velocity field and water depth field based on hydrological monitoring data, and fuse the gradient information with the point cloud data to optimize the spatial position to generate density-balanced riverbed point cloud data.

6. The method according to claim 1, wherein Construct a triangular grid terrain model based on the riverbed point cloud data, and generate a sequence of river channel cross-sections at preset intervals along the river flow direction; calculate the hydrodynamic correlation intensity for the sequence of river channel cross-sections through a recurrent neural network, and construct a cross-section evolution propagation sequence according to the hydrodynamic correlation intensity, including: Perform constrained Delaunay triangulation on the riverbed point cloud data to construct an initial triangular grid, calculate grid quality parameters based on the ratio of triangle area to side length, and optimize the triangular grid according to the grid quality parameters to obtain a triangular grid terrain model; Construct a terrain elevation distribution based on the triangular grid terrain model, calculate the elevation difference between adjacent grid nodes, determine the water flow direction field based on the elevation difference, sample at preset intervals along the water flow direction field on the riverbed terrain triangular grid model to generate a sequence of river channel cross-sections, extract width, water depth, cross slope, and roughness parameters for each cross-section, and construct cross-section feature vectors; Input the cross-section feature vector into the long short-term memory network, extract the cross-section time series features through the memory unit and hidden state of the long short-term memory network, calculate the cosine similarity of adjacent cross-sections based on the hidden state, and multiply the cosine similarity by the exponential decay function of the cross-section spacing to obtain the cross-section correlation strength; construct a propagation weight matrix based on the cross-section correlation strength and the hydraulic gradient factor, and use the propagation weight matrix to iteratively update the cross-section state to obtain the river channel cross-section evolution propagation sequence.

7. The method according to claim 1, wherein Based on the cross-section evolution propagation sequence and historical evolution data, calculate the riverbed erosion and deposition change rate through recursive calculation; when the riverbed erosion and deposition change rate exceeds the preset rate threshold, generate a warning message, and visualize the dynamic prediction result of the river channel cross-section morphology and the warning message, including: Perform multi-source data fusion on the cross-section evolution propagation sequence and historical evolution data to obtain a fusion feature vector, construct a feature weight and quality evaluation function based on the fusion feature vector, and perform quality evaluation on the fusion feature vector through the feature weight and the quality evaluation function to obtain an evaluation result; Construct a graph structure based on the evaluation result, construct a node set from the cross-section data, construct an associated edge set from the associations between cross-sections, generate an adjacency matrix based on the node set and the associated edge set, and use the adjacency matrix to update the message passing of the feature information of adjacent cross-sections to obtain hidden state features; Input the hidden state features into a mapping function to obtain the current change rate, perform weighted fusion of the current change rate and the historical change rate to obtain the erosion and deposition change rate of the target cross-section, calculate statistical parameters based on the erosion and deposition change rate, and perform exponential mapping of the statistical parameters and the reference threshold to obtain a dynamic warning threshold; Compare the erosion and deposition change rate with the dynamic warning threshold to obtain a warning level, perform temporal integration on the target cross-section based on the erosion and deposition change rate to obtain a predicted morphology, and perform three-dimensional visualization of the predicted morphology and the warning level to generate an interactive warning message.

8. A multi-temporal river channel image-based dynamic change monitoring system for implementing the method according to any one of the preceding claims 1-7, characterized in that, Including: The first unit is used to obtain multi-temporal multi-spectral remote sensing image data and unmanned aerial vehicle oblique photogrammetry image data of the river channel area; Based on the multi-temporal multi-spectral remote sensing image data, identify the river channel water body boundary through the support vector machine classification algorithm and construct a river channel water area mask; The second unit is used to generate dense point cloud data of the river channel area based on the unmanned aerial vehicle oblique photogrammetry image data by using the multi-view image matching method; crop the dense point cloud data by using the river channel water area mask to obtain point cloud data; The third unit is used to calculate the point cloud neighborhood density distribution function for the point cloud data, determine the local importance weight according to the gradient change of the point cloud neighborhood density distribution function, and establish a density compensation sampling matrix based on the local importance weight; perform recursive multi-scale enhancement on the sparse area of the point cloud by using the density compensation sampling matrix, and optimize the spatial position in combination with the local curvature constraint to generate density-balanced riverbed point cloud data; The fourth unit is used to construct a triangular mesh terrain model based on the riverbed point cloud data, and generate a sequence of river channel cross-sections at preset intervals along the river flow direction; calculate the hydrodynamic correlation intensity for the sequence of river channel cross-sections through a recurrent neural network, and construct a cross-section evolution propagation sequence according to the hydrodynamic correlation intensity. The fifth unit is used to obtain the rate of riverbed erosion and deposition change through recursive calculation based on the cross-section evolution propagation sequence and historical evolution data; when the rate of riverbed erosion and deposition change exceeds a preset rate threshold, generate a warning message, and visually display the dynamic prediction result of the river channel cross-section morphology and the warning message.

9. An electronic device, characterized in that, It includes: A processor; A memory for storing instructions executable by the processor; Wherein, the processor is configured to call the instructions stored in the memory to execute the method according to any one of claims 1 to 7.

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

Images