Method and system for monitoring dynamic changes of river channels based on multi-temporal river images

Abstract

This invention provides a method and system for monitoring dynamic changes in river channels based on multi-temporal imagery. These methods involve acquiring multi-temporal multispectral remote sensing images and unmanned aerial vehicle (UAV) oblique photogrammetry images; identifying river water boundaries and constructing water masks; generating dense point cloud data and cropping it; generating a balanced riverbed point cloud through density-compensated sampling and local curvature constraints; constructing a terrain model to generate a cross-section sequence; calculating hydrodynamic correlation strength to construct an evolutionary propagation sequence; and recursively calculating the scouring and deposition change rate and issuing warnings. This invention enables high-precision monitoring and early warning of river channel morphological changes.

Description

Technical Field

[0001] The present invention relates to image monitoring technology, and in particular to a method and system for monitoring dynamic changes of multi-temporal river channel images. Background Art

[0002] Rivers are crucial vehicles for socioeconomic development, and monitoring their dynamic changes is crucial for flood prevention and disaster reduction, water resource management, and ecological and environmental protection. With the advancement of remote sensing and drone technologies, the use of multi-source remote sensing data for dynamic river monitoring has become a research hotspot. Traditional river monitoring relies primarily on manual field surveys and analysis of single satellite images, which struggle to achieve continuous monitoring and accurate prediction of river changes. Modern river monitoring methods, primarily based on multi-temporal remote sensing imagery and drone aerial photography, enable dynamic monitoring of river water extent, riverbed topography, and erosion and deposition changes through the processing and analysis of multi-temporal images.

[0003] Existing technologies generally have density imbalance problems 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 subtle changes in riverbed topography, thereby affecting the accuracy and reliability of subsequent river channel change analysis.

[0004] Existing river section analysis methods often lack consideration of the hydrodynamic correlation between sections. Relying solely on static section comparisons makes it difficult to reveal the inherent mechanism of river evolution and cannot effectively capture the propagation laws of river changes, resulting in inaccurate prediction results.

[0005] The existing river monitoring system lacks an effective early warning mechanism, making it difficult to detect abnormal scouring and silting changes in the riverbed in a timely manner, and unable to provide timely decision-making support for flood control and disaster reduction and river management, which reduces the practical value and preventive effect of river monitoring. Summary of the Invention

[0006] The embodiments of the present invention provide a method and system for monitoring dynamic changes of river channels based on multi-temporal images, which can solve the problems in the prior art.

[0007] A first aspect of an embodiment of the present invention provides a method for monitoring dynamic changes of river channels based on multi-temporal river channel images, comprising:

[0008] Acquire multi-temporal multispectral remote sensing image data and unmanned aerial vehicle (UAV) oblique photogrammetry image data of the river area; identify the river water boundary using a support vector machine classification algorithm based on the multi-temporal multispectral remote sensing image data, and construct a river water mask;

[0009] Based on the UAV oblique photogrammetry image data, a multi-view image matching method is used to generate dense point cloud data of the river area; the dense point cloud data is cropped using the river water mask to obtain point cloud data;

[0010] Calculating a point cloud neighborhood density distribution function for the point cloud data, determining a local importance weight based on a gradient change of the point cloud neighborhood density distribution function, and establishing a density compensation sampling matrix based on the local importance weight; recursively performing multi-scale enhancement on sparse regions of the point cloud using the density compensation sampling matrix, and optimizing spatial positions in combination with local curvature constraints to generate density-balanced riverbed point cloud data;

[0011] A triangular mesh terrain model is constructed based on the riverbed point cloud data, and a river section sequence is generated at preset intervals along the river flow direction; a hydrodynamic correlation strength is calculated for the river section sequence using a recursive neural network, and a section evolution propagation sequence is constructed based on the hydrodynamic correlation strength;

[0012] 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, an early warning message is generated, and the dynamic prediction results of the river section morphology and the early warning message are visualized.

[0013] Based on the multi-temporal and multispectral remote sensing image data, the river water body boundary is identified by using a support vector machine classification algorithm, and a river water mask is constructed, including:

[0014] Constructing a feature vector based on the multi-temporal and multispectral remote sensing image data, mapping the feature vector to a high-dimensional feature space, constructing an optimal classification hyperplane in the high-dimensional feature space, calculating the distance between a sample point and the hyperplane based on the optimal classification hyperplane, using the sample points whose distance is greater than a preset distance threshold as training support vectors, and constructing a classification decision function using the training support vectors; performing classification prediction on the feature vector based on the classification decision function, and generating an initial classification image according to the prediction result;

[0015] Performing morphological processing on the initial classified image, using the first structuring element to perform an opening operation to remove noise points, and using the second structuring element to perform a closing operation to fill holes, to obtain an initial water area mask image;

[0016] The initial water area mask image is smoothed using a contour evolution model constrained by boundary curvature, a control point sequence of a contour curve is initialized based on the initial water area mask image, an internal energy function characterizing the elasticity and rigidity of the contour curve is constructed based on the first-order derivatives and second-order derivatives of the control point sequence, and an external energy function is constructed based on the gradient information of the initial classification image;

[0017] The internal energy function and the external energy function are combined to obtain the total energy function of the contour evolution model, and the total energy function is minimized. When the change of the total energy function is less than the preset convergence threshold, the iteration is stopped, and a smooth river water body boundary is reconstructed based on the final control point sequence to generate a river water area mask.

[0018] Based on the UAV oblique photogrammetry image data, a multi-view image matching method is used to generate dense point cloud data of the river area; the dense point cloud data is cropped using the river water mask to obtain point cloud data including:

[0019] Establishing a matching pair group for the UAV oblique photogrammetry image data, calculating the ratio of the overlapping area of ​​each pair of images in the matching pair group to the total area of ​​the images, and selecting the image pairs whose overlapping ratio is greater than a preset overlapping threshold as the image pairs to be matched;

[0020] For the reference image in the pair of images to be matched, a normalized cross-correlation coefficient is calculated between each pixel in the reference image and the corresponding pixel in the other matching images, a weight coefficient is determined according to the baseline length of each matching image, and a weighted sum of the normalized cross-correlation coefficient and the weight coefficient is performed to obtain a multi-baseline matching cost;

[0021] In a neighborhood window of each pixel point in the reference image, an adaptive weight is calculated based on the similarity of the pixel grayscale values, and an aggregated matching cost is obtained by weighted summing the adaptive weight and the multi-baseline matching cost;

[0022] Determining an optimal disparity value of a pixel point using the aggregate matching cost, establishing a collinearity equation system based on the optimal disparity value, solving the collinearity equation system using a least squares method to obtain three-dimensional point coordinates, and generating dense point cloud data;

[0023] The dense point cloud data is transformed into an image plane coordinate system through a projection matrix, the projected point cloud data is cropped according to the river water mask, and the point cloud data within the river water mask is retained as the point cloud data of the river area.

[0024] Determining the optimal disparity value of the pixel point using the aggregate matching cost, establishing a collinearity equation system based on the optimal disparity value, solving the collinearity equation system by the least squares method to obtain the three-dimensional point coordinates, and generating dense point cloud data includes:

[0025] constructing a disparity space cost curve within a preset disparity search range based on the aggregate matching cost, determining an integer pixel extreme point of the disparity space cost curve, and performing quadratic curve fitting using the cost value of the integer pixel extreme point and its adjacent disparity positions to obtain an optimal disparity value;

[0026] Constructing a collinearity equation group according to the optimal disparity value and the internal and external orientation elements of the image, solving the collinearity equation group to obtain an initial three-dimensional coordinate point set, and constructing a spatial point cloud based on the initial three-dimensional coordinate point set;

[0027] Reprojecting each 3D coordinate point in the spatial point cloud to each perspective image, calculating the reprojection error between the projected point and the original matching point, and marking the 3D coordinate point with a reprojection error greater than a preset error threshold as a point to be optimized;

[0028] For the point to be optimized, a local point set is extracted within its neighborhood, 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;

[0029] The optimized three-dimensional coordinate point set is reconstructed into a dense point cloud, the point density of each point in the reconstructed point cloud within its preset neighborhood 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.

[0030] The density compensation sampling matrix is ​​used to perform recursive multi-scale enhancement on the sparse area of ​​the point cloud, and the spatial position is optimized in combination with the local curvature constraint to generate density-balanced riverbed point cloud data, including:

[0031] Constructing a multi-scale feature extraction network, performing feature extraction on multiple preset radius neighborhoods of each point in the point cloud data based on the multi-scale feature extraction network, performing weighted fusion of the extracted features with corresponding weight coefficients to obtain a feature vector, inputting the feature vector into a density prediction network to obtain a predicted density value, and multiplying the predicted density value with the initial density field to obtain a density prediction result for each point;

[0032] Determine sparse areas in the point cloud based on the density prediction result, generate landform feature weights, divide a preset target density value by the density prediction result, and multiply the result by the landform feature weight to obtain a regional sampling weight, and adaptively fuse the sonar bathymetric data and the multi-temporal and multispectral remote sensing image data based on the regional sampling weight to generate enhanced point cloud data;

[0033] Constructing a local connection graph for the enhanced point cloud data, iteratively updating node features in the local connection graph to obtain optimized node features; calculating the distance weight of each point within a preset radius neighborhood based on the optimized node features, and performing a weighted summation of the distance weight and the projection distance from the point to the neighboring points to obtain a local curvature value;

[0034] The optimized node features and the local curvature values ​​are input into a generator to generate an initial density-balanced point cloud, and the generator is used to regenerate point cloud data with local geometric features; based on the hydrological monitoring data, the gradient information of the velocity field and the water depth field is constructed, and the gradient information is fused with the point cloud data to optimize the spatial position and generate density-balanced riverbed point cloud data.

[0035] Constructing a triangular mesh terrain model based on the riverbed point cloud data, generating a river section sequence at preset intervals along the river flow direction; calculating the hydrodynamic correlation strength of the river section sequence using a recursive neural network, and constructing a section evolution propagation sequence based on the hydrodynamic correlation strength includes:

[0036] Performing constrained Delaunay triangulation on the riverbed point cloud data to construct an initial triangular mesh, calculating a mesh quality parameter based on a ratio of triangle area to side length, optimizing the triangular mesh according to the mesh quality parameter, and obtaining a triangular mesh terrain model;

[0037] constructing a terrain elevation distribution based on the triangular mesh terrain model, calculating the elevation differences between adjacent mesh nodes, determining a water flow direction field based on the elevation differences, sampling the riverbed terrain triangular mesh model at preset intervals along the water flow direction field to generate a river channel section sequence, extracting the width, water depth, cross slope, and roughness parameters of each section, and constructing a section feature vector;

[0038] The section feature vector is input into a long short-term memory network, and the section time series features are extracted through the memory units and hidden states of the long short-term memory network. The cosine similarity of adjacent sections is calculated based on the hidden state, and the section association strength is obtained by multiplying the cosine similarity with the exponential decay function of the section distance. A propagation weight matrix is ​​constructed based on the section association strength and the hydraulic gradient factor, and the section state is iteratively updated using the propagation weight matrix to obtain the river section evolution propagation sequence.

[0039] 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, an early warning message is generated, and the dynamic prediction results of the river channel cross-section morphology and the early warning message are visually displayed, including:

[0040] Performing multi-source data fusion on the section evolution propagation sequence and the historical evolution data to obtain a fused feature vector, constructing a feature weight and a quality assessment function based on the fused feature vector, and performing a quality assessment on the fused feature vector using the feature weight and the quality assessment function to obtain an assessment result;

[0041] Constructing a graph structure based on the evaluation results, constructing the section data into a node set, constructing the associations between sections into an associated edge set, generating an adjacency matrix based on the node set and the associated edge set, and using the adjacency matrix to perform message passing and update feature information of adjacent sections to obtain hidden state features;

[0042] Inputting the latent state feature into a mapping function to obtain a current change rate, performing weighted fusion on the current change rate and the historical change rate to obtain a scouring and silting change rate of the target section, calculating a statistical parameter based on the scouring and silting change rate, and performing exponential mapping on the statistical parameter and a benchmark threshold to obtain a dynamic warning threshold;

[0043] The scouring and silting change rate is compared with the dynamic warning threshold to obtain a warning level, the target section is time-series integrated based on the scouring and silting change rate to obtain a predicted shape, the predicted shape and the warning level are three-dimensionally visualized to generate interactive warning information.

[0044] A second aspect of an embodiment of the present invention provides a system for monitoring dynamic changes in river channels based on multi-temporal river channel images, comprising:

[0045] The first unit is used to obtain multi-temporal multispectral remote sensing image data and unmanned aerial vehicle oblique photogrammetry image data of the river area; based on the multi-temporal multispectral remote sensing image data, the support vector machine classification algorithm is used to identify the river water body boundary and construct the river water area mask;

[0046] The second unit is configured to generate dense point cloud data of the river area based on the UAV oblique photogrammetry image data using a multi-view image matching method; and to crop the dense point cloud data using the river water mask to obtain point cloud data;

[0047] A third unit is configured to calculate a point cloud neighborhood density distribution function for the point cloud data, determine a local importance weight based on a 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 sparse areas of the point cloud using the density compensation sampling matrix, and optimize spatial positions in combination with local curvature constraints to generate density-balanced riverbed point cloud data;

[0048] A fourth unit is configured to construct a triangular mesh terrain model based on the riverbed point cloud data, generate a river section sequence at preset intervals along the river flow direction, calculate the hydrodynamic correlation strength of the river section sequence using a recursive neural network, and construct a section evolution propagation sequence based on the hydrodynamic correlation strength;

[0049] The fifth unit is used 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, an early warning message is generated, and the dynamic prediction results of the river section morphology and the early warning message are visualized.

[0050] According to a third aspect of an embodiment of the present invention, an electronic device is provided, including:

[0051] processor;

[0052] a memory for storing processor-executable instructions;

[0053] The processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0054] According to a fourth aspect of an embodiment of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the method described above is implemented.

[0055] The beneficial effects of this application are as follows:

[0056] This method uses the support vector machine classification algorithm to identify the river water body boundary and construct a water area mask. Combined with the dense point cloud generated by UAV oblique photogrammetry data, it achieves accurate extraction and expression of the river area, improving the spatial accuracy and reliability of river monitoring.

[0057] To address the uneven distribution problem 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 constrained optimization, it generates density-balanced riverbed point cloud data, effectively solving the terrain expression deviation problem caused by uneven point cloud distribution in traditional methods.

[0058] This method uses a recursive neural network to calculate the hydrodynamic correlation strength, constructs a cross-section evolution propagation sequence, and combines historical data to recursively calculate the riverbed scouring and deposition change rate, achieving dynamic prediction and early warning of river morphological changes. It provides a scientific basis and technical support for river management and flood prevention and disaster reduction, and has important practical value. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 Schematic diagram of the process of a method for monitoring dynamic changes of multi-temporal river channel images according to an embodiment of the present invention;

[0060] Figure 2 A histogram showing the performance comparison and analysis of the multi-view image matching method according to an embodiment of the present invention;

[0061] Figure 3This is a flow chart of density equalization processing of riverbed point cloud data according to an embodiment of the present invention;

[0062] Figure 4 Schematic diagram comparing the triangular mesh quality optimization effects according to an embodiment of the present invention;

[0063] Figure 5 This is a bar chart comparing the performance of the riverbed scour and sedimentation change monitoring and early warning system according to an embodiment of the present invention. DETAILED DESCRIPTION

[0064] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0065] The technical solution of the present invention is described in detail below with reference to specific embodiments. The following specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described in detail in some embodiments.

[0066] Figure 1 FIG. 1 is a flow chart of a method for monitoring dynamic changes of multi-temporal river channel images according to an embodiment of the present invention. Figure 1 As shown, the method includes:

[0067] Acquire multi-temporal multispectral remote sensing image data and unmanned aerial vehicle (UAV) oblique photogrammetry image data of the river area; identify the river water boundary using a support vector machine classification algorithm based on the multi-temporal multispectral remote sensing image data, and construct a river water mask;

[0068] Based on the UAV oblique photogrammetry image data, a multi-view image matching method is used to generate dense point cloud data of the river area; the dense point cloud data is cropped using the river water mask to obtain point cloud data;

[0069] Calculating a point cloud neighborhood density distribution function for the point cloud data, determining a local importance weight based on a gradient change of the point cloud neighborhood density distribution function, and establishing a density compensation sampling matrix based on the local importance weight; recursively performing multi-scale enhancement on sparse regions of the point cloud using the density compensation sampling matrix, and optimizing spatial positions in combination with local curvature constraints to generate density-balanced riverbed point cloud data;

[0070] A triangular mesh terrain model is constructed based on the riverbed point cloud data, and a river section sequence is generated at preset intervals along the river flow direction; a hydrodynamic correlation strength is calculated for the river section sequence using a recursive neural network, and a section evolution propagation sequence is constructed based on the hydrodynamic correlation strength;

[0071] 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, an early warning message is generated, and the dynamic prediction results of the river section morphology and the early warning message are visualized.

[0072] In an optional embodiment, based on the multi-temporal and multispectral remote sensing image data, identifying the river water boundary by a support vector machine classification algorithm and constructing the river water mask includes:

[0073] Constructing a feature vector based on the multi-temporal and multispectral remote sensing image data, mapping the feature vector to a high-dimensional feature space, constructing an optimal classification hyperplane in the high-dimensional feature space, calculating the distance between a sample point and the hyperplane based on the optimal classification hyperplane, using the sample points whose distance is greater than a preset distance threshold as training support vectors, and constructing a classification decision function using the training support vectors; performing classification prediction on the feature vector based on the classification decision function, and generating an initial classification image according to the prediction result;

[0074] Performing morphological processing on the initial classified image, using the first structuring element to perform an opening operation to remove noise points, and using the second structuring element to perform a closing operation to fill holes, to obtain an initial water area mask image;

[0075] The initial water area mask image is smoothed using a contour evolution model constrained by boundary curvature, a control point sequence of a contour curve is initialized based on the initial water area mask image, an internal energy function characterizing the elasticity and rigidity of the contour curve is constructed based on the first-order derivatives and second-order derivatives of the control point sequence, and an external energy function is constructed based on the gradient information of the initial classification image;

[0076] The internal energy function and the external energy function are combined to obtain the total energy function of the contour evolution model, and the total energy function is minimized. When the change of the total energy function is less than the preset convergence threshold, the iteration is stopped, and a smooth river water body boundary is reconstructed based on the final control point sequence to generate a river water area mask.

[0077] By utilizing the high-dimensional mapping capability of the support vector machine classification algorithm, combined with morphological processing and a contour evolution model constrained by boundary curvature, accurate identification of river water body boundaries is achieved.

[0078] After acquiring multi-temporal and multispectral remote sensing imagery data, a feature vector must be constructed based on this data. This feature vector typically contains information about the reflectance of each band, the Normalized Difference Water Width Index (NDWI), the Normalized Difference Vegetation Index (NDVI), and other parameters. For example, for a given pixel, 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 vector, it is mapped to a high-dimensional feature space using a kernel function. In practice, the radial basis function (RBF) kernel generally performs well, and the kernel parameter γ can be set to 0.125.

[0079] In a high-dimensional feature space, the support vector machine algorithm searches for the optimal classification hyperplane that separates water samples from non-water samples. This hyperplane is determined by maximizing the separation between the two classes of samples. The distance from the sample points to the hyperplane is calculated, and sample points with a distance greater than a preset distance threshold (e.g., 0.8) are used as training support vectors. Based on these support vectors, a classification decision function is constructed, which can be used to predict the category of new samples. Classification predictions are performed on the feature vectors to obtain an initial classification image, in which water pixels are labeled 1 and non-water pixels are labeled 0.

[0080] Initial classified images often contain noise and holes, requiring morphological processing to improve accuracy. Morphological processing consists of two steps: opening and closing. The opening operation removes noise points, while the closing operation fills holes. In the opening operation, a 3×3 circular structuring element is used as the first structuring element. An erosion operation is performed first, followed by a dilation operation, effectively removing small misclassified water points. In the closing operation, a 5×5 circular structuring element is used as the second structuring element. A dilation operation is performed first, followed by an erosion operation, filling small holes within the water area. These two steps are combined to produce an initial water area mask image.

[0081] To further optimize the water boundary, a contour evolution model with boundary curvature constraints is used to smooth the initial water mask image. This process begins with initializing a sequence of control points for the contour curve, evenly distributed along the boundary of the initial water mask image. For an image with a smoothing region of 512 × 512 pixels, 100–200 control points can be selected. An internal energy function representing the elasticity and rigidity of the contour curve is constructed based on the first and second derivatives of the control point sequence. The first-order derivative controls the elasticity of the contour, while the second-order derivative controls the rigidity. The elasticity weight coefficient can be set to 0.45, and the rigidity weight coefficient can be set to 0.55.

[0082] An external energy function is constructed based on the gradient information of the initial classified image to guide the contour toward the actual water boundary. Regions with large gradient values ​​typically represent boundaries. Gradient calculations provide the gradient magnitude and direction for each pixel in the image. The gradient magnitude is calculated using the Sobel operator, with a threshold of 25. Pixels exceeding this threshold are considered boundary points. The weight coefficient of the external energy function can be set to 0.65.

[0083] 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 is set to 0.5 pixels to prevent the contour from deforming too quickly. The iteration is stopped when the change in the total energy function is less than the preset convergence threshold (such as 0.001). Usually, convergence can be achieved after 50-100 iterations. Based on the final control point sequence, the smooth river water boundary is reconstructed through methods such as B-spline interpolation or polynomial fitting to generate a river water mask.

[0084] Experiments show that this method can achieve an average accuracy of 93.5% when applied to river remote sensing images of different seasons, which is 7.8 percentage points higher than the traditional threshold segmentation method, and the boundary smoothness is improved by 65%. It is particularly suitable for river extraction tasks under complex backgrounds.

[0085] In an optional embodiment, based on the UAV oblique photogrammetry image data, a multi-view image matching method is used to generate dense point cloud data of the river area; the dense point cloud data is cropped using the river water mask to obtain point cloud data including:

[0086] Establishing a matching pair group for the UAV oblique photogrammetry image data, calculating the ratio of the overlapping area of ​​each pair of images in the matching pair group to the total area of ​​the images, and selecting the image pairs whose overlapping ratio is greater than a preset overlapping threshold as the image pairs to be matched;

[0087] For the reference image in the pair of images to be matched, a normalized cross-correlation coefficient is calculated between each pixel in the reference image and the corresponding pixel in the other matching images, a weight coefficient is determined according to the baseline length of each matching image, and a weighted sum of the normalized cross-correlation coefficient and the weight coefficient is performed to obtain a multi-baseline matching cost;

[0088] In a neighborhood window of each pixel point in the reference image, an adaptive weight is calculated based on the similarity of the pixel grayscale values, and an aggregated matching cost is obtained by weighted summing the adaptive weight and the multi-baseline matching cost;

[0089] Determining an optimal disparity value of a pixel point using the aggregate matching cost, establishing a collinearity equation system based on the optimal disparity value, solving the collinearity equation system using a least squares method to obtain three-dimensional point coordinates, and generating dense point cloud data;

[0090] The dense point cloud data is transformed into an image plane coordinate system through a projection matrix, the projected point cloud data is cropped according to the river water mask, and the point cloud data within the river water mask is retained as the point cloud data of the river area.

[0091] Collect oblique photogrammetric image data from the river channel using drones, including both vertical and oblique images. Simultaneously, use image segmentation techniques to generate a river channel mask for subsequent point cloud data cropping.

[0092] After acquiring drone oblique photogrammetry image data, matching pairs are created for this data. Specifically, the ratio of the overlapping area of ​​each pair of images in the matching pair to the total area of ​​the images is calculated. For example, for two images A and B, their overlapping area is calculated to be 500 square meters, the total area of ​​image A is 1000 square meters, and the total area of ​​image B is 1200 square meters. Taking the average of the two, 1100 square meters, the overlap ratio is 500 / 1100 ≈ 0.45. The preset overlap 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.

[0093] For the selected image pairs to be matched, a dense point cloud is generated. One image is selected as the reference image in the pair, and the normalized cross-correlation coefficient (NCC) is calculated between each pixel in the reference image and the corresponding pixel in the other matching image. The NCC is calculated by comparing the similarity of the grayscale distribution within a pixel neighborhood window, with the window size set to 11×11 pixels.

[0094] A weighting factor is determined based on the baseline length of each matching image. The baseline length refers to the distance between the two camera positions. For example, if the baseline length between the reference image and matching image A is 20 meters, the baseline length between the reference image and matching image B is 10 meters, and the baseline length between the reference image and matching image C is 30 meters, then the weighting factors for these three matching images can be set to 0.3, 0.2, and 0.5, respectively. The normalized cross-correlation coefficient and the weighting factors are weighted and summed to obtain the multi-baseline matching cost.

[0095] To improve matching accuracy, an adaptive weight is calculated based on the similarity of the pixel grayscale values ​​within a neighborhood window of each pixel in the reference image. The neighborhood window size is set to 5×5 pixels. The higher the pixel grayscale similarity, the larger the corresponding adaptive weight. The adaptive weight is weighted and summed with the multi-baseline matching cost to obtain the aggregate matching cost.

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

[0097] Based on the optimal disparity value, a set of collinearity equations is established. Specifically, the relationship between the object point coordinates and the image point coordinates is established based on the camera's intrinsic and extrinsic parameters and the image point coordinates. The collinearity equations are solved using the least squares method to obtain the three-dimensional coordinates of the point. For example, for a pixel in the river channel, 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 equations are (521368.45, 3876215.78, 125.62). This processing is performed on all pixels in the reference image to generate dense point cloud data.

[0098] After generating the dense point cloud data, it is transformed into the image plane coordinate system using a projection matrix. The projection matrix consists of camera intrinsic parameters (such as focal length, principal point coordinates, and distortion coefficients) and extrinsic parameters (such as camera position and attitude). For example, assuming a camera focal length of 35 mm, principal point coordinates of (1024, 768), no distortion, camera position of (521000, 3876000, 200), and attitude angle of (0°, 0°, 0°), the corresponding projection matrix can be constructed to project the 3D point cloud onto the 2D image plane.

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

[0100] For example, for a 3D point with coordinates (800, 600) after projection on the image plane, if the value of the river water mask at that location is 1, the point is retained; if the mask value is 0, the point is discarded. This method can effectively extract point cloud data of the river area, providing data support for subsequent river analysis and monitoring.

[0101] After the above processing, the resulting point cloud data of the river area contains rich 3D information, which can be used for applications such as river cross-section analysis, water surface elevation measurement, and river deformation monitoring. In practice, according to different river characteristics and application requirements, various parameters such as overlap threshold, window size, and disparity search range can be adjusted to achieve the best point cloud extraction results.

[0102] Figure 2 This is a bar chart comparing the performance of the multi-view image matching method according to an embodiment of the present invention:

[0103] The image shows a performance comparison of three different matching methods across various scenarios: traditional NCC matching, multi-baseline weighted matching, and adaptive weighted aggregation matching. Across different scenarios, the matching success rates of the three methods were 34.0%, 42.0%, and 50.0%, respectively, in low-texture areas; 28.0%, 38.0%, and 44.0%, respectively, in reflective areas; and 32.0%, 40.0%, and 46.0%, respectively, in areas covered by vegetation. Averagely across all scenarios, the overall matching success rates of the three methods were 34.0%, 42.0%, and 48.0%, respectively. The data shows that the adaptive weighted aggregation matching method achieved the best results in all scenarios, improving by an average of 14 percentage points compared to the traditional NCC matching method and by an average of 6 percentage points compared to the multi-baseline weighted matching method. The adaptive weighted aggregation matching method was particularly advantageous in low-texture areas, achieving a matching success rate of 50.0%.

[0104] In an optional embodiment, the aggregate matching cost is used to determine the optimal disparity value of the pixel point, 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. The generation of dense point cloud data includes:

[0105] constructing a disparity space cost curve within a preset disparity search range based on the aggregate matching cost, determining an integer pixel extreme point of the disparity space cost curve, and performing quadratic curve fitting using the cost value of the integer pixel extreme point and its adjacent disparity positions to obtain an optimal disparity value;

[0106] Constructing a collinearity equation group according to the optimal disparity value and the internal and external orientation elements of the image, solving the collinearity equation group to obtain an initial three-dimensional coordinate point set, and constructing a spatial point cloud based on the initial three-dimensional coordinate point set;

[0107] Reprojecting each 3D coordinate point in the spatial point cloud to each perspective image, calculating the reprojection error between the projected point and the original matching point, and marking the 3D coordinate point with a reprojection error greater than a preset error threshold as a point to be optimized;

[0108] For the point to be optimized, a local point set is extracted within its neighborhood, 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;

[0109] The optimized three-dimensional coordinate point set is reconstructed into a dense point cloud, the point density of each point in the reconstructed point cloud within its preset neighborhood 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.

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

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

[0112] When reprojecting each 3D coordinate point in the spatial point cloud to each view image, use the same internal and external orientation elements to project the 3D coordinate point (X=156.32m, Y=278.45m, Z=32.16m) back to the original image to obtain the projected point coordinates (x'=12.48mm, y'=8.23mm). Calculate the reprojection error between the projected point and the original matching point, that is, The preset error threshold is set to 0.1 mm. If the reprojection error is greater than this threshold, the point is marked as a point to be optimized.

[0113] For the point to be optimized, a local point set is extracted within its neighborhood. All points within a spherical area with a radius of 2 meters and centered at the point can be taken as the local point set. For example, for a 3D coordinate point (X=186.45m, Y=321.78m, Z=28.92m), its reprojection error is 0.15mm, which exceeds the threshold of 0.1mm and is marked as a point to be optimized. A local point set containing 52 points is extracted within its 2-meter neighborhood, and the centroid coordinates are calculated as (X c =186.38m, Y c =321.82m, Z c =28.95m). A weighted update is performed on the points to be optimized based on the centroid coordinates, with the weight depending on the distance from the point to the centroid. The updated 3D coordinates become (X'=186.40m, Y'=321.80m, Z'=28.94m). The reprojection error is reduced to 0.08mm, which is below the threshold, completing the optimization.

[0114] When reconstructing a dense point cloud for the optimized three-dimensional coordinate point set, the point density of each point in the reconstructed point cloud within the preset neighborhood is calculated. The preset neighborhood range is set to 1 cubic meter, and the preset density threshold is 10 points per cubic meter. For example, for a three-dimensional point (X=245.67m, Y=418.93m, Z=35.24m), the number of points within its 1 cubic meter neighborhood is calculated to be 8, which is lower than the density threshold of 10 / cubic meter. This area is marked as an area that requires encrypted reconstruction. In this area, by increasing the density of matching points on the original image, such as from taking a matching point every 5 pixels to taking a matching point every 2 pixels, the above steps are repeated 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.

[0115] Through these steps, dense point cloud data was generated for the entire image area. The final point cloud data contained approximately 5 million 3D points, with an average point density of 25 points per square meter. The mean reprojection error was 0.06mm, and the standard deviation was 0.03mm, meeting the requirements for high-precision 3D reconstruction. The resulting dense point cloud data can be used for subsequent applications such as digital surface model generation and 3D model reconstruction.

[0116] In an optional embodiment, recursive multi-scale enhancement is performed on sparse regions of the point cloud using the density compensation sampling matrix, and spatial position is optimized in combination with local curvature constraints to generate density-balanced riverbed point cloud data, including:

[0117] Constructing a multi-scale feature extraction network, performing feature extraction on multiple preset radius neighborhoods of each point in the point cloud data based on the multi-scale feature extraction network, performing weighted fusion of the extracted features with corresponding weight coefficients to obtain a feature vector, inputting the feature vector into a density prediction network to obtain a predicted density value, and multiplying the predicted density value with the initial density field to obtain a density prediction result for each point;

[0118] Determine sparse areas in the point cloud based on the density prediction result, generate landform feature weights, divide a preset target density value by the density prediction result, and multiply the result by the landform feature weight to obtain a regional sampling weight, and adaptively fuse the sonar bathymetric data and the multi-temporal and multispectral remote sensing image data based on the regional sampling weight to generate enhanced point cloud data;

[0119] Constructing a local connection graph for the enhanced point cloud data, iteratively updating node features in the local connection graph to obtain optimized node features; calculating the distance weight of each point within a preset radius neighborhood based on the optimized node features, and performing a weighted summation of the distance weight and the projection distance from the point to the neighboring points to obtain a local curvature value;

[0120] The optimized node features and the local curvature values ​​are input into a generator to generate an initial density-balanced point cloud, and the generator is used to regenerate point cloud data with local geometric features; based on the hydrological monitoring data, the gradient information of the velocity field and the water depth field is constructed, and the gradient information is fused with the point cloud data to optimize the spatial position and generate density-balanced riverbed point cloud data.

[0121] like Figure 3 As shown, the method further includes:

[0122] A multi-scale feature extraction network was constructed, consisting of three convolutional layers with kernel sizes of 32, 64, and 128, respectively, to capture local and global features of the point cloud. For each point in the point cloud, features were extracted from neighborhoods with preset radii of 0.5, 1.0, and 2.0 meters. During feature extraction, a maximum pooling operation was used to aggregate neighborhood information, setting the dimension of the feature vector extracted at each scale to 128. The features extracted at the three different scales were weighted and fused with corresponding weight coefficients of 0.3, 0.3, and 0.4, resulting in a 384-dimensional feature vector.

[0123] This feature vector is then fed into a density prediction network consisting of three fully connected layers, with 256, 128, and 1 neurons, respectively. The activation function uses the Reluctant Unified Unit (ReLU) function, and the final layer uses a Sigmoid function to map the output to a range between 0 and 1, yielding the predicted density value. The predicted density value is then multiplied by the initial density field (obtained by calculating the number of points within a 1.5-meter radius neighborhood for each point) to produce the density prediction for each point.

[0124] The sparse areas in the point cloud are determined based on the density prediction results. The specific method is to set the density threshold to 20 points per cubic meter. The areas below the threshold are marked as sparse areas. At the same time, the landform feature weights are generated, which are calculated based on the curvature and slope information of the point cloud. The curvature is obtained by fitting the quadratic surface of 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 and slope values, the landform feature weights are obtained by weighting them at a ratio of 0.6 and 0.4. The preset target density value (set to 50 points per cubic meter) is divided by the density prediction result, and then multiplied by the landform feature weight to obtain the regional sampling weight.

[0125] Based on this sampling weight, sonar bathymetric data and multi-temporal and multispectral remote sensing imagery are adaptively fused. Sonar data provides precise depth information, with a weighting factor of 0.7; remote sensing imagery provides planar position information through edge detection and texture analysis, with a weighting factor of 0.3. The fusion process employs an octree-based spatial index structure, merging the sonar point cloud and the point cloud derived from the remote sensing imagery using weighted fusion in overlapping areas and directly merging non-overlapping areas to generate enhanced point cloud data.

[0126] A local connection graph is constructed for the augmented point cloud data. Each node in the connection graph connects to its neighboring points within a 1.2-meter radius, with an average of approximately 12 connections per node. A graph convolutional network is used to iterate and update the node features in the connection graph over three rounds. In each round, information from neighboring nodes is aggregated and the features of the central node are updated. The update formula takes into account edge distance weights, with closer neighbors having a greater influence. After the iterations are complete, optimized node features are obtained, with a feature dimension of 256.

[0127] Based on the optimized node features, the distance weight of each point within a 1.0-meter radius neighborhood is calculated using a Gaussian kernel function with a standard deviation of 0.3 meters. The distance weight is then weighted and summed with the projected distance from the point to its neighbors. The projected distance is the distance from the point to the best-fit plane formed by the neighboring points to obtain the local curvature value.

[0128] The optimized node features and local curvature values ​​are concatenated and fed into a generator. The generator employs a three-layer encoder-decoder architecture. The encoder maps the input features into a latent space, and the decoder reconstructs the point cloud coordinates from this latent space. The generator outputs an initial density-balanced point cloud with a density of 45-55 points per cubic meter. The same generator is then used to regenerate point cloud data with local geometric features. Random noise (with an amplitude of 0.05 meters) is introduced during the regeneration process to increase the natural variation of the point cloud.

[0129] Gradient information for velocity and depth fields was constructed based on hydrological monitoring data. Velocity data was obtained from an acoustic Doppler current profiler (ADCP) with a sampling interval of 5 meters, while depth data was obtained from an echo sounder with a sampling interval of 2 meters. Trilinear interpolation was used to construct continuous velocity and depth fields, and the gradient information of the calculated fields was used to represent the changing trends of water flow and topography.

[0130] Gradient information is integrated with point cloud data, and point positions are fine-tuned along the gradient direction. The adjustment is proportional to the gradient magnitude, with a maximum adjustment of no more than 0.3 meters to maintain the continuity and authenticity of the terrain. This adjustment generates a final, density-balanced riverbed point cloud, with the standard deviation of the point cloud density distribution reduced to less than 30% of the original data, accurately representing the riverbed's micro-topography.

[0131] In an optional embodiment, a triangular mesh terrain model is constructed based on the riverbed point cloud data, and a river section sequence is generated at preset intervals along the river flow direction; a hydrodynamic correlation strength is calculated for the river section sequence using a recursive neural network, and a section evolution propagation sequence is constructed based on the hydrodynamic correlation strength, including:

[0132] Performing constrained Delaunay triangulation on the riverbed point cloud data to construct an initial triangular mesh, calculating a mesh quality parameter based on a ratio of triangle area to side length, optimizing the triangular mesh according to the mesh quality parameter, and obtaining a triangular mesh terrain model;

[0133] constructing a terrain elevation distribution based on the triangular mesh terrain model, calculating the elevation differences between adjacent mesh nodes, determining a water flow direction field based on the elevation differences, sampling the riverbed terrain triangular mesh model at preset intervals along the water flow direction field to generate a river channel section sequence, extracting the width, water depth, cross slope, and roughness parameters of each section, and constructing a section feature vector;

[0134] The section feature vector is input into a long short-term memory network, and the section time series characteristics are extracted through the memory units and hidden states of the long short-term memory network. The cosine similarity of adjacent sections is calculated based on the hidden state, and the section association strength is obtained by multiplying the cosine similarity with the exponential decay function of the section distance. A propagation weight matrix is ​​constructed based on the section association strength and the hydraulic gradient factor, and the section state is iteratively updated using the propagation weight matrix to obtain the river section evolution propagation sequence.

[0135] The constrained Delaunay triangulation algorithm is used to construct the initial triangular mesh. This algorithm ensures that the generated triangles meet the Delaunay condition by considering the spatial position relationship of each point in the point cloud, that is, the circumcircle of any triangle does not contain other points. 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 efficiency of point query and reduces the time complexity of triangulation from O(n) to O(n). 2 ) is reduced to O(nlogn).

[0136] After constructing the initial triangular mesh, the system calculates the 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 lengths of the three sides. For an ideal equilateral triangle, the Q value is close to 0.5; 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 triangles with Q values ​​below 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 the overall quality can be improved by flipping the shared edge; when the internal angle of the triangle is less than 20 degrees, new points are inserted inside it and re-divided. After three rounds of iterative optimization, the average value of the mesh quality parameter increased from 0.35 to 0.42, forming a high-quality triangular mesh terrain model.

[0137] Based on the optimized triangular mesh terrain model, the system constructs the terrain elevation distribution and calculates the elevation difference between adjacent mesh nodes. The elevation difference is calculated by dividing the difference in the elevation values ​​of two adjacent nodes by the horizontal distance between the two points. In an actual case, the riverbed terrain triangular mesh contains 78,562 nodes and 156,934 triangles. The system calculates the elevation difference between each node and its adjacent nodes and uses this to determine the water flow direction field. The water flow direction is defined as the direction of the steepest drop in elevation and is represented by the vector from each node to its lowest adjacent node.

[0138] When generating a sequence of river channel cross sections along the flow direction, the system sets the cross-section spacing to 50 meters, generating a total of 200 cross sections along a 10-kilometer river. For each cross-section, the system intersects the triangular mesh model on a plane perpendicular to the flow direction to obtain the cross-section profile. For each cross section, the system extracts parameters such as cross-section width, average water depth, left and right bank slopes, and cross-section roughness.

[0139] The cross-section width is defined as the horizontal distance from the leftmost to the rightmost point of the cross-section; the average water depth is the difference between the average elevation 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 and the lowest point of the riverbed; and the roughness parameter is calculated from the cross-section morphology and the distribution of the bottom material. In one example, the extracted cross-section feature vector includes a cross-section width of 102 meters, an average water depth of 2.8 meters, a left bank slope of 0.32, a right bank slope of 0.28, and a roughness coefficient of 0.035.

[0140] The extracted cross-section feature vectors are input into a long short-term memory (LSTM) network for processing. This network consists of an input layer, a forget gate, an input gate, an output gate, and a hidden layer. The network input is a time-sequential sequence of cross-section feature vectors. A three-layer LSTM architecture extracts cross-section temporal features, with each layer containing 128 neurons. The system uses 10 cross-sections as a time window for sliding training, with a learning rate of 0.001, a batch size of 32, and 500 training iterations. The network's memory cells store long-term dependency information, and the hidden state represents the comprehensive characteristics of the current cross-section.

[0141] The cosine similarity of adjacent sections is calculated based on the hidden state of the LSTM network. For each pair of adjacent sections i and j, the dot product of their hidden state vectors is divided by the product of the two vectors' moduli to obtain the cosine similarity value. The cosine similarity is then multiplied by the exponential decay function of the section spacing, with the decay coefficient set to 0.05, to obtain the section correlation strength matrix. In this example, the average correlation strength of adjacent sections is 0.78, and the correlation strength decreases exponentially with increasing section spacing.

[0142] A propagation weight matrix is ​​constructed based on the cross-sectional correlation strength and hydraulic gradient factors. The hydraulic gradient factor accounts for variations in riverbed slope and cross-sectional morphology and is calculated as the ratio of the average elevation difference between cross-sectional areas to the distance between cross-sectional areas. The propagation weight matrix defines the influence weight of each adjacent cross-sectional area during the cross-sectional state update process. The system updates cross-sectional states using an iterative method, with a set number of iterations of 30 and a convergence threshold of 0.001. Through iterative updates of the propagation weight matrix, a cross-sectional evolution propagation sequence is generated that reflects the dynamic changes in the river channel. This sequence can be used for river channel evolution analysis and hydrodynamic simulation.

[0143] Figure 4 This is a schematic diagram comparing the triangular mesh quality optimization effects of an embodiment of the present invention:

[0144] This figure compares the network quality indicators of three different technical solutions at different mesh density levels. The horizontal axis represents the mesh density level from very low to very high, and the vertical axis represents the network quality indicator. This technical solution (circular markers) shows the best performance at all density levels: starting from 0.16 at very 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), and finally to very high density (0.68), showing a steady upward trend. The edge-constrained triangulation solution (diamond markers) is second, gradually improving from 0.12 at very low density to 0.52 at very high density. Traditional Delaunay triangulation (square markers) performs the worst, improving only from 0.08 at very low density to 0.34 at very high density. Data shows that the network quality index 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 edge constraint-based method, reflecting significant technical advantages.

[0145] In an optional embodiment, based on the cross-section evolution propagation sequence and historical evolution data, a 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 results of the river channel cross-section morphology and the warning message are visually displayed, including:

[0146] Performing multi-source data fusion on the section evolution propagation sequence and the historical evolution data to obtain a fused feature vector, constructing a feature weight and a quality assessment function based on the fused feature vector, and performing a quality assessment on the fused feature vector using the feature weight and the quality assessment function to obtain an assessment result;

[0147] Constructing a graph structure based on the evaluation results, constructing the section data into a node set, constructing the associations between sections into an associated edge set, generating an adjacency matrix based on the node set and the associated edge set, and using the adjacency matrix to perform message passing and update feature information of adjacent sections to obtain hidden state features;

[0148] Inputting the latent state feature into a mapping function to obtain a current change rate, performing weighted fusion on the current change rate and the historical change rate to obtain a scouring and silting change rate of the target section, calculating a statistical parameter based on the scouring and silting change rate, and performing exponential mapping on the statistical parameter and a benchmark threshold to obtain a dynamic warning threshold;

[0149] The scouring and silting change rate is compared with the dynamic warning threshold to obtain a warning level, the target section is time-series integrated based on the scouring and silting change rate to obtain a predicted shape, the predicted shape and the warning level are three-dimensionally visualized to generate interactive warning information.

[0150] Based on the cross-section evolution propagation sequence and historical evolution data, the rate of change of riverbed erosion and deposition is recursively calculated. When the rate of change exceeds a preset rate threshold, an early warning message is generated, and the dynamic prediction results and early warning information of the river cross-section morphology are visualized.

[0151] Multi-source data fusion is performed on the cross-section evolution propagation sequence and historical evolution data to generate a fused feature vector. This multi-source data includes hydrological station monitoring data, cross-section measurement data, and remote sensing imagery. For example, assuming there are six hydrological stations, each containing 10 features such as flow, sediment content, and water level, this results in a 6×10 feature matrix. Furthermore, temporal features are constructed using historical evolution data. For example, the changes in erosion and sedimentation at a certain cross-section over the past five years were +0.2 meters, -0.15 meters, +0.1 meters, -0.05 meters, and +0.3 meters, respectively. These heterogeneous data are converted into a vector of uniform dimension through feature extraction and then fused using a weighted average method, resulting in a fused feature vector of 128 dimensions.

[0152] Based on the fused feature vectors, feature weights and a quality assessment function are constructed. Feature weights are determined based on the contribution of each feature to changes in erosion and deposition. For example, flow is weighted 0.4, sediment content is weighted 0.3, water level is weighted 0.2, and other features are weighted 0.1. The quality assessment function comprehensively considers data completeness, timeliness, and consistency, and assigns a comprehensive score based on various evaluation indicators. For example, a data completeness rate of 95% scores 0.95, a data update time of within 3 days scores 0.9, and a data consistency check pass rate of 98% scores 0.98, resulting in an overall score of 0.95 × 0.4 + 0.9 × 0.3 + 0.98 × 0.3 = 0.947.

[0153] The fused feature vector is evaluated using feature weights and a quality assessment function to generate a quality assessment result. This assessment includes reliability indexes and uncertainty ranges. For example, a reliability index of 0.92 for a particular cross-section indicates a 92% confidence level; an uncertainty range of ±0.05 meters indicates a ±0.05-meter error in the predicted result.

[0154] Based on the evaluation results, a graph structure is constructed. The section data is structured as a node set, and the associations between sections are structured as an edge set. For example, if there are 10 section points, 10 nodes are formed. Edges are established between adjacent sections. For example, there are edges between sections 1 and 2, and between sections 2 and 3, while sections 1 and 3 are not directly connected. An adjacency matrix is ​​generated based on the node set and edge set. The adjacency matrix is ​​a 10×10 matrix, where A[i][j]=1 indicates that sections i and j are adjacent, and A[i][j]=0 indicates that they are not adjacent.

[0155] The adjacency matrix is ​​used to update the feature information of adjacent sections through message passing to obtain hidden state features. During message passing, each node performs a weighted aggregation of its own features with those of adjacent nodes. For example, the updated features of Section 2 are composed of a 0.6 weight for its own features, a 0.2 weight for the features of Section 1, and a 0.2 weight for the features of Section 3. After three rounds of iterative updates, the hidden state features containing global information are finally obtained.

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

[0157] The current rate of change and the historical rate of change are weighted together to obtain the scouring and deposition rate of the target section. The current rate of change is weighted as 0.7, and the historical rate of change is weighted as 0.3. For example, if the current calculated rate of change for a section is 0.15 m / year and the historical average rate of change is 0.1 m / year, the fused rate of change is 0.15 × 0.7 + 0.1 × 0.3 = 0.135 m / year.

[0158] Statistical parameters, including mean and standard deviation, are calculated based on the rate of change of erosion and deposition. For example, the mean rate of change of erosion and deposition across 10 sections of a river is 0.08 meters per year, with a standard deviation of 0.03 meters per year. Dynamic warning thresholds are obtained by exponentially mapping statistical parameters to baseline thresholds. The baseline threshold is set at 0.2 meters per year, and then mapped to dynamic thresholds using an exponential mapping function. For example, during the flood season, the dynamic threshold is adjusted to 0.2 × 1.5 = 0.3 meters per year; during the dry season, the dynamic threshold is adjusted to 0.2 × 0.8 = 0.16 meters per year.

[0159] The rate of change in erosion and deposition is compared with the dynamic warning threshold to determine the warning level. Warning levels are categorized as normal, caution, warning, and danger. A rate of change less than 50% of the threshold is considered normal; between 50% and 80% of the threshold is considered caution; between 80% and 100% of the threshold is considered warning; and exceeding the threshold is considered danger. For example, if the rate of change in erosion and deposition at a section is 0.25 meters per year and the dynamic threshold is 0.3 meters per year, the rate of change accounts for 83.3% of the threshold, resulting in a warning level.

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

[0161] The forecast and warning levels are visualized in 3D, generating interactive warning information. The 3D visualization includes plan views, cross-sectional views, and time-series graphs. The plan view shows the river layout and warning areas; the cross-sectional view displays the current and predicted cross-sectional shapes; and the time-series graph displays the elevation trends over time. Warning information is color-coded: normal (green), caution (yellow), warning (orange), and danger (red). Users can view forecast results at different time points and receive warning alerts through an interactive interface.

[0162] Figure 5 This is a bar chart comparing the performance of the riverbed scouring and silting change monitoring and early warning system according to an embodiment of the present invention:

[0163] The figure shows a comparative analysis of the early warning performance of three different technical solutions, including three performance indicators: early warning accuracy, response time, and change prediction accuracy. In terms of early warning accuracy, this technical solution achieved 48.0%, significantly outperforming the traditional statistical adjustment method (36.0%) and the simple neural network method (40.0%). In terms of response time, this technical solution performed best, with a 52.0% advantage, compared to 32.0% for the traditional statistical adjustment method and 36.0% for the simple neural network method. In terms of change prediction accuracy, this technical solution also maintained a leading position, reaching 44.0%, compared to 28.0% for the traditional statistical adjustment method and 32.0% for the simple neural network method. Overall, this technical solution achieved significant advantages across all performance indicators, with an average lead of approximately 12-16 percentage points over the other two methods. The advantage was particularly pronounced in the key indicator of response time, demonstrating the outstanding value of this solution in practical applications.

[0164] A second aspect of an embodiment of the present invention provides a system for monitoring dynamic changes in river channels based on multi-temporal river channel images, comprising:

[0165] The first unit is used to obtain multi-temporal multispectral remote sensing image data and unmanned aerial vehicle oblique photogrammetry image data of the river area; based on the multi-temporal multispectral remote sensing image data, the support vector machine classification algorithm is used to identify the river water body boundary and construct the river water area mask;

[0166] The second unit is configured to generate dense point cloud data of the river area based on the UAV oblique photogrammetry image data using a multi-view image matching method; and to crop the dense point cloud data using the river water mask to obtain point cloud data;

[0167] A third unit is configured to calculate a point cloud neighborhood density distribution function for the point cloud data, determine a local importance weight based on a 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 sparse areas of the point cloud using the density compensation sampling matrix, and optimize spatial positions in combination with local curvature constraints to generate density-balanced riverbed point cloud data;

[0168] A fourth unit is configured to construct a triangular mesh terrain model based on the riverbed point cloud data, generate a river section sequence at preset intervals along the river flow direction, calculate the hydrodynamic correlation strength of the river section sequence using a recursive neural network, and construct a section evolution propagation sequence based on the hydrodynamic correlation strength;

[0169] The fifth unit is used 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, an early warning message is generated, and the dynamic prediction results of the river section morphology and the early warning message are visualized.

[0170] According to a third aspect of an embodiment of the present invention, an electronic device is provided, including:

[0171] processor;

[0172] a memory for storing processor-executable instructions;

[0173] The processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0174] According to a fourth aspect of an embodiment of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the method described above is implemented.

[0175] 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 carrying computer-readable program instructions for executing various aspects of the present invention.

[0176] 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 it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the above embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for monitoring dynamic changes in river channels based on multi-temporal river channel images, characterized in that: include: Acquire multi-temporal and multi-spectral remote sensing image data and UAV oblique photogrammetry image data of river areas; Based on the multi-temporal and multispectral remote sensing image data, identifying the river water body boundary by a support vector machine classification algorithm and constructing a river water mask; Based on the UAV oblique photogrammetry image data, a multi-view image matching method is used to generate dense point cloud data of the river area; the dense point cloud data is cropped using the river water mask to obtain point cloud data; Calculating a point cloud neighborhood density distribution function for the point cloud data, determining a local importance weight based on a gradient change of the point cloud neighborhood density distribution function, and establishing a density compensation sampling matrix based on the local importance weight; recursively performing multi-scale enhancement on sparse regions of the point cloud using the density compensation sampling matrix, and optimizing spatial positions in combination with local curvature constraints to generate density-balanced riverbed point cloud data; A triangular mesh terrain model is constructed based on the riverbed point cloud data, and a river section sequence is generated at preset intervals along the river flow direction; a hydrodynamic correlation strength is calculated for the river section sequence using a recursive neural network, and a section evolution propagation sequence is constructed based on the hydrodynamic correlation strength; 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, an early warning message is generated, and the dynamic prediction results of the river section morphology and the early warning message are visualized.

2. The method according to claim 1, characterized in that Based on the multi-temporal and multispectral remote sensing image data, the river water body boundary is identified by using a support vector machine classification algorithm, and a river water mask is constructed, including: Constructing a feature vector based on the multi-temporal and multispectral remote sensing image data, mapping the feature vector to a high-dimensional feature space, constructing an optimal classification hyperplane in the high-dimensional feature space, calculating the distance between a sample point and the hyperplane based on the optimal classification hyperplane, using the sample points whose distance is greater than a preset distance threshold as training support vectors, and constructing a classification decision function using the training support vectors; performing classification prediction on the feature vector based on the classification decision function, and generating an initial classification image according to the prediction result; Performing morphological processing on the initial classified image, using the first structuring element to perform an opening operation to remove noise points, and using the second structuring element to perform a 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 constrained by boundary curvature, a control point sequence of a contour curve is initialized based on the initial water area mask image, an internal energy function characterizing the elasticity and rigidity of the contour curve is constructed based on the first-order derivatives and second-order derivatives of the control point sequence, and 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, and the total energy function is minimized. When the change of the total energy function is less than the preset convergence threshold, the iteration is stopped, and a smooth river water body boundary is reconstructed based on the final control point sequence to generate a river water area mask.

3. The method according to claim 1, characterized in that Based on the UAV oblique photogrammetry image data, a multi-view image matching method is used to generate dense point cloud data of the river area; the dense point cloud data is cropped using the river water mask to obtain point cloud data including: Establishing a matching pair group for the UAV oblique photogrammetry image data, calculating the ratio of the overlapping area of ​​each pair of images in the matching pair group to the total area of ​​the images, and selecting the image pairs whose overlapping ratio is greater than a preset overlapping threshold as the image pairs to be matched; For the reference image in the pair of images to be matched, a normalized cross-correlation coefficient is calculated between each pixel in the reference image and the corresponding pixel in the other matching images, a weight coefficient is determined according to the baseline length of each matching image, and a weighted sum of the normalized cross-correlation coefficient and the weight coefficient is performed to obtain a multi-baseline matching cost; In a neighborhood window of each pixel point in the reference image, an adaptive weight is calculated based on the similarity of the pixel grayscale values, and an aggregated matching cost is obtained by weighted summing the adaptive weight and the multi-baseline matching cost; Determining an optimal disparity value of a pixel point using the aggregate matching cost, establishing a collinearity equation system based on the optimal disparity value, solving the collinearity equation system using a least squares method to obtain three-dimensional point coordinates, and generating dense point cloud data; The dense point cloud data is transformed into an image plane coordinate system through a projection matrix, the projected point cloud data is cropped according to the river water mask, and the point cloud data within the river water mask is retained as the point cloud data of the river area.

4. The method according to claim 3, characterized in that Determining the optimal disparity value of the pixel point using the aggregate matching cost, establishing a collinearity equation system based on the optimal disparity value, solving the collinearity equation system by the least squares method to obtain the three-dimensional point coordinates, and generating dense point cloud data includes: constructing a disparity space cost curve within a preset disparity search range based on the aggregate matching cost, determining an integer pixel extreme point of the disparity space cost curve, and performing quadratic curve fitting using the cost value of the integer pixel extreme point and its adjacent disparity positions to obtain an optimal disparity value; Constructing a collinearity equation group according to the optimal disparity value and the internal and external orientation elements of the image, solving the collinearity equation group to obtain an initial three-dimensional coordinate point set, and constructing a spatial point cloud based on the initial three-dimensional coordinate point set; Reprojecting each 3D coordinate point in the spatial point cloud to each perspective image, calculating the reprojection error between the projected point and the original matching point, and marking the 3D coordinate point with a reprojection error greater than a preset error threshold as a point to be optimized; For the point to be optimized, a local point set is extracted within its neighborhood, 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; The optimized three-dimensional coordinate point set is reconstructed into a dense point cloud, the point density of each point in the reconstructed point cloud within its preset neighborhood 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.

5. The method according to claim 1, wherein The density compensation sampling matrix is ​​used to perform recursive multi-scale enhancement on the sparse area of ​​the point cloud, and the spatial position is optimized in combination with the local curvature constraint to generate density-balanced riverbed point cloud data, including: Constructing a multi-scale feature extraction network, performing feature extraction on multiple preset radius neighborhoods of each point in the point cloud data based on the multi-scale feature extraction network, performing weighted fusion of the extracted features with corresponding weight coefficients to obtain a feature vector, inputting the feature vector into a density prediction network to obtain a predicted density value, and multiplying the predicted density value with the initial density field to obtain a density prediction result for each point; Determine sparse areas in the point cloud based on the density prediction result, generate landform feature weights, divide a preset target density value by the density prediction result, and multiply the result by the landform feature weight to obtain a regional sampling weight, and adaptively fuse the sonar bathymetric data and the multi-temporal and multispectral remote sensing image data based on the regional sampling weight to generate enhanced point cloud data; Constructing a local connection graph for the enhanced point cloud data, iteratively updating node features in the local connection graph to obtain optimized node features; calculating the distance weight of each point within a preset radius neighborhood based on the optimized node features, and performing a weighted summation of the distance weight and the projection distance from the point to the neighboring points to obtain a local curvature value; The optimized node features and the local curvature values ​​are input into a generator to generate an initial density-balanced point cloud, and the generator is used to regenerate point cloud data with local geometric features; based on the hydrological monitoring data, the gradient information of the velocity field and the water depth field is constructed, and the gradient information is fused with the point cloud data to optimize the spatial position and generate density-balanced riverbed point cloud data.

6. The method according to claim 1, characterized in that Constructing a triangular mesh terrain model based on the riverbed point cloud data, generating a river section sequence at preset intervals along the river flow direction; calculating the hydrodynamic correlation strength of the river section sequence using a recursive neural network, and constructing a section evolution propagation sequence based on the hydrodynamic correlation strength includes: Performing constrained Delaunay triangulation on the riverbed point cloud data to construct an initial triangular mesh, calculating a mesh quality parameter based on a ratio of triangle area to side length, optimizing the triangular mesh according to the mesh quality parameter, and obtaining a triangular mesh terrain model; constructing a terrain elevation distribution based on the triangular mesh terrain model, calculating the elevation differences between adjacent mesh nodes, determining a water flow direction field based on the elevation differences, sampling a river channel section sequence on the triangular mesh terrain model at preset intervals along the water flow direction field, extracting the width, water depth, cross slope, and roughness parameters of each section, and constructing a section feature vector; The section feature vector is input into a long short-term memory network, and the section time series characteristics are extracted through the memory units and hidden states of the long short-term memory network. The cosine similarity of adjacent sections is calculated based on the hidden state, and the section association strength is obtained by multiplying the cosine similarity with the exponential decay function of the section distance. A propagation weight matrix is ​​constructed based on the section association strength and the hydraulic gradient factor, and the section state is iteratively updated using the propagation weight matrix to obtain the river section evolution propagation sequence.

7. The method according to claim 1, characterized in that 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, an early warning message is generated, and the dynamic prediction results of the river channel cross-section morphology and the early warning message are visually displayed, including: Performing multi-source data fusion on the section evolution propagation sequence and the historical evolution data to obtain a fused feature vector, constructing a feature weight and a quality assessment function based on the fused feature vector, and performing a quality assessment on the fused feature vector using the feature weight and the quality assessment function to obtain an assessment result; Constructing a graph structure based on the evaluation results, constructing the section data into a node set, constructing the associations between sections into an associated edge set, generating an adjacency matrix based on the node set and the associated edge set, and using the adjacency matrix to perform message passing and update feature information of adjacent sections to obtain hidden state features; Inputting the latent state feature into a mapping function to obtain a current change rate, performing weighted fusion on the current change rate and the historical change rate to obtain a scouring and silting change rate of the target section, calculating a statistical parameter based on the scouring and silting change rate, and performing exponential mapping on the statistical parameter and a benchmark threshold to obtain a dynamic warning threshold; The scouring and silting change rate is compared with the dynamic warning threshold to obtain a warning level, the target section is time-series integrated based on the scouring and silting change rate to obtain a predicted shape, the predicted shape and the warning level are three-dimensionally visualized to generate interactive warning information.

8. A multi-temporal river channel image dynamic change monitoring system for implementing the method according to any one of claims 1 to 7, characterized in that: include: The first unit is used to obtain multi-temporal and multi-spectral remote sensing image data and UAV oblique photogrammetry image data of the river area; Based on the multi-temporal and multispectral remote sensing image data, identifying the river water body boundary by a support vector machine classification algorithm and constructing a river water mask; The second unit is configured to generate dense point cloud data of the river area based on the UAV oblique photogrammetry image data using a multi-view image matching method; and to crop the dense point cloud data using the river water 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 a local importance weight based on a 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 sparse areas of the point cloud using the density compensation sampling matrix, and optimize spatial positions in combination with local curvature constraints to generate density-balanced riverbed point cloud data; A fourth unit is configured to construct a triangular mesh terrain model based on the riverbed point cloud data, generate a river section sequence at preset intervals along the river flow direction, calculate the hydrodynamic correlation strength of the river section sequence using a recursive neural network, and construct a section evolution propagation sequence based on the hydrodynamic correlation strength; The fifth unit is used 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, an early warning message is generated, and the dynamic prediction results of the river section morphology and the early warning message are visualized.

9. An electronic device, characterized in that: include: processor; a memory for storing processor-executable instructions; 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 a processor, the method according to any one of claims 1 to 7 is implemented.

Images