Lake form monitoring method based on satellite remote sensing technology

Abstract

The invention discloses a lake form monitoring method based on a satellite remote sensing technology, and particularly relates to the technical field of water environment monitoring. By acquiring optical, SAR and thermal infrared remote sensing images of multiple time nodes, radiation correction, atmospheric correction and geometric registration are completed; fusing a multi-source image to extract an initial boundary of the lake, and constructing a morphological characteristic spectrum; constructing a time sequence characteristic matrix based on indexes such as curvature, compactness, unstable factors and gravity center, and establishing a morphological variation path diagram model; a sudden change node and a sudden change area are identified through a dynamic time warping algorithm, and a sudden change weight is calculated in combination with a curvature amplitude and a spatial position; comparing with a historical form database to form a microcosmic variation risk probability graph; and finally, outputting a lake form stability index, and generating an early warning layer of a high-risk area. The method is high in precision, strong in timeliness and suitable for lake dynamic monitoring and risk early warning in a complex environment.

Description

Technical Field

[0001] This invention relates to the field of water environment monitoring technology, specifically to a method for monitoring lake morphology based on satellite remote sensing technology. Background Technology

[0002] With the intensification of global climate change and human activities, the spatiotemporal changes in the water morphology (including area, shoreline curvature, and center of gravity) of lakes, as important carriers of freshwater resources, have become key indicators for environmental monitoring and ecological protection. However, traditional lake monitoring methods mainly rely on ground measurements, drone aerial photography, or data from fixed monitoring points, which have many problems such as long measurement cycles, limited spatial coverage, and poor dynamic response capabilities, making it difficult to achieve high-frequency dynamic monitoring of lake morphology on short time scales.

[0003] While some current remote sensing technologies can acquire information about large-scale lake bodies, they still face significant technical bottlenecks, such as insufficient resolution, complex multi-source data fusion, severe cloud cover, and difficulty in extracting morphological details. In particular, when monitoring small lakes or rapidly evolving seasonal water body boundaries, it is even more difficult to capture local deformation and subtle fissions of lake shorelines and their corresponding evolutionary mechanisms, resulting in insufficient timeliness and accuracy of early warning of lake ecological risks.

[0004] Especially in arid and semi-arid regions, lake morphology is significantly affected by abrupt changes in precipitation, sediment transport and deposition, and human-induced water diversion and regulation. This can easily lead to nonlinear morphological abrupt changes, such as rare structural variations like "funnel-shaped" rupture and "swallowtail-shaped" confluence. These morphological abrupt changes not only directly reflect the evolutionary trends of lake functional zones but also easily induce environmental disasters such as water quality degradation, wetland shrinkage, and even lake disappearance. Summary of the Invention

[0005] The purpose of this invention is to provide a lake morphology monitoring method based on satellite remote sensing technology to address the shortcomings of the prior art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a lake morphology monitoring method based on satellite remote sensing technology, comprising:

[0007] Acquire multi-source satellite remote sensing image data of the target area at multiple time points, including optical images, SAR radar images, and thermal infrared images;

[0008] The acquired remote sensing images are subjected to radiometric correction, atmospheric correction and geometric registration to construct a multi-source image sequence under a unified time reference;

[0009] A model for extracting the initial lake boundary was constructed based on fused image sequences. The initial morphological contour of the lake boundary was extracted by multi-scale edge detection and water spectral feature matching.

[0010] The initial morphological contour is transformed into a lake morphological feature map, and the following are extracted: shoreline curvature distribution index, morphological compactness index, boundary instability factor and water body distribution centroid value.

[0011] Based on the morphological feature map, a time-series feature matrix was constructed, a lake morphological variation path map model was established, and the dynamic time warping algorithm was used to align and match the variation paths of adjacent time periods to identify abrupt change nodes and morphological change regions.

[0012] Based on the spatial distribution of mutation nodes and the amplitude of shoreline curvature change within the morphological transition area, the morphological mutation weight value of the morphological transition area is calculated and compared with the feature deviation of the historical morphological database to form a micro-variation risk probability map.

[0013] Based on the micro-variation risk probability map, the lake morphology stability index is output, and a warning marker layer for high-risk areas is generated.

[0014] Preferably, the lake initial boundary extraction model is constructed based on the fused image sequence:

[0015] The preprocessed optical, SAR and thermal infrared images are fused at multiple resolution scales to output a multi-scale fused image.

[0016] The improved Canny operator combined with the Laplacian operator is applied to the fused image at each scale to perform multi-scale edge detection, and the edge candidate set is cross-validated to generate a boundary confidence map.

[0017] Based on the boundary confidence map, edges with high confidence are used as seeds, and adaptive region growing and spline curve fitting are combined to form the final initial shape contour of the lake.

[0018] Preferably, the extracted parameters include the shoreline curvature distribution index, morphological compactness index, boundary instability factor, and water body distribution centroid value, including:

[0019] The extraction of the shoreline curvature distribution index includes: discretizing the initial boundary contour of the lake into a sequence of equally spaced sampling points, calculating the curvature value of each point using the three-point sliding method, and performing segmented statistics on the curvature values ​​of all points to form a curvature distribution histogram, where curvature is defined as the reciprocal of the arc formed by three consecutive points, and the curvature distribution index is the ratio of the maximum curvature segment to the average curvature of the entire range.

[0020] The extraction of the morphological compactness index includes: calculating the area A of the region enclosed by the lake boundary contour and its boundary length P, according to the formula: Compactness Index = Solve the problem.

[0021] The extraction of the boundary instability factor includes: performing point-to-point matching between the boundary contour of the current time node and the boundary contour of the previous time node, calculating the proportion of boundary segments with a boundary change rate higher than a set threshold, and defining them as boundary instability factors.

[0022] The extraction of the water body distribution centroid value includes: converting the lake area into a binary mask image, using the water body pixel coordinates as input, calculating the weighted average position of all water body pixels on a two-dimensional plane, and obtaining the water body distribution centroid value. , where n is the total number of water pixels, and xi and yi are the coordinates of the i-th pixel.

[0023] Preferably, the step of constructing a time-series feature matrix based on morphological feature maps and establishing a lake morphological variation path map model includes:

[0024] The morphological features extracted at each time point are arranged in chronological order to form the original time series feature matrix, and each feature column is normalized to zero mean and unit variance.

[0025] Principal component analysis was performed on the normalized time series matrix, and principal components with a cumulative variance contribution rate of not less than 90% were retained to obtain a low-dimensional time series embedding representation.

[0026] The low-dimensional representation of each time node is used as a node in the path graph, and the weighted distance of dynamic time warping similarity and cosine similarity is calculated for each pair of adjacent time nodes as the edge weight.

[0027] Threshold pruning is performed on the constructed weighted path graph, and shortest path search is used to identify the main path of morphological variation and mutation nodes.

[0028] Preferably, the step of aligning and matching mutation paths in adjacent time periods using a dynamic time warping algorithm to identify mutation nodes and morphological transition regions includes:

[0029] Based on the low-dimensional temporal features of each node in the lake morphological variation path graph model, feature sequence pairs between adjacent time nodes are extracted to construct a local variation sub-path set.

[0030] For each pair of sub-path sequences, a dynamic time warping algorithm is used for nonlinear alignment. The optimal pairing path and its alignment cost matrix are calculated as a quantitative expression of the impact of time scale changes on morphological evolution.

[0031] The DTW values ​​of all paired paths are normalized, and abnormal alignment regions are identified by setting a mutation threshold and marked as candidate mutation nodes.

[0032] Within the time window of the candidate node, the boundary change rate map and centroid offset vector field corresponding to the original image sequence are traced back. Combining the morphological jump intensity and spatial position overlap, the morphological jump region is finally determined and spatial annotation results are generated.

[0033] Preferably, the step of calculating the morphological abrupt change weight value of the morphological abrupt change region based on the spatial distribution location of the abrupt change node and the shoreline curvature variation amplitude within the morphological abrupt change region includes:

[0034] Curvature values ​​of all boundary points are extracted within the identified morphological transition region. Local curvature is calculated using the three-point sliding method, and curvature change curves are constructed.

[0035] The difference between the maximum and minimum curvature values ​​within the transition region is calculated and used as an index of the curvature amplitude of the region.

[0036] Based on the spatial distribution of mutation nodes throughout the lake region, the normalized distance from the mutation node to the lake centroid is calculated to reflect the central sensitivity of the mutation.

[0037] The curvature amplitude index and the normalized spatial distance are weighted and combined to construct the morphological change weight value.

[0038] Preferably, the step of comparing the feature deviation with the historical morphology database to form a micro-variation risk probability map includes:

[0039] Establish a historical morphological database covering the target lake, storing morphological feature maps and their corresponding spatial distribution information at multiple past time points;

[0040] The Euclidean distance between the morphological feature vector of the current mutation region and the average feature vector of historical regions at the same location is calculated to obtain the feature deviation value.

[0041] The morphological mutation weight value and the feature deviation value are jointly normalized and converted into a micro-variation risk probability value using a probability mapping function.

[0042] The risk probability values ​​of all abrupt change regions are mapped to a spatial coordinate system to form a micro-variation risk probability map that includes the degree of mutation and spatial location.

[0043] Preferably, based on the micro-variation risk probability map, the lake morphology stability index is output, including:

[0044] Statistical analysis was performed on the risk probability values ​​of all abrupt regions in the microvariation risk probability diagram, and their mean and standard deviation were calculated.

[0045] Based on the risk probability value and spatial area of ​​each transition region, an area-weighted risk mean is constructed and defined as the lake morphology stability index.

[0046] A high-risk threshold is set, and areas with a lake morphological stability index greater than 0.6 or a single area risk probability value exceeding 0.8 are judged as high-risk.

[0047] Areas that meet the high-risk criteria are highlighted as layers on the remote sensing image to generate a warning label layer, and the spatial coordinates and corresponding risk levels are output.

[0048] The technical effects and advantages provided by the present invention in the above technical solution are as follows:

[0049] 1. This invention provides a lake morphology monitoring method based on satellite remote sensing technology, which can achieve high-precision, automated, and time-series monitoring of lake morphology evolution. By fusing multi-source remote sensing data and combining multi-scale edge detection and morphological feature extraction, it not only improves the accuracy of identifying complex lake boundaries but also overcomes the problem of unstable morphological extraction under rainy, nighttime, or seasonal conditions, effectively enhancing the adaptability of the monitoring system in multi-climate and multi-regional environments.

[0050] 2. The morphological variation path map model and micro-variation risk probability map introduced in this invention can comprehensively characterize the abrupt change behavior of lakes from two dimensions: temporal evolution and spatial distribution. By combining mutation weights, historical deviation comparisons, and risk mapping, it can achieve dynamic identification of high-risk areas and generation of early warning layers. This method has the advantages of clear structure, well-defined indicators, and ease of integration and deployment, providing strong technical support for lake ecological security assessment, water resource regulation, and disaster risk management. Attached Figure Description

[0051] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0052] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0053] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0054] For examples, please refer to Figure 1As shown in this embodiment, a lake morphology monitoring method based on satellite remote sensing technology includes:

[0055] Acquire multi-source satellite remote sensing image data of the target area at multiple time points, including optical images, SAR radar images, and thermal infrared images;

[0056] The acquired remote sensing images are subjected to radiometric correction, atmospheric correction and geometric registration to construct a multi-source image sequence under a unified time reference;

[0057] A model for extracting the initial lake boundary was constructed based on fused image sequences. The initial morphological contour of the lake boundary was extracted by multi-scale edge detection and water spectral feature matching.

[0058] The initial morphological contour is transformed into a lake morphological feature map, and the following are extracted: shoreline curvature distribution index, morphological compactness index, boundary instability factor and water body distribution centroid value.

[0059] Based on the morphological feature map, a time-series feature matrix was constructed, a lake morphological variation path map model was established, and the dynamic time warping algorithm was used to align and match the variation paths of adjacent time periods to identify abrupt change nodes and morphological change regions.

[0060] Based on the spatial distribution of mutation nodes and the amplitude of shoreline curvature change within the morphological transition area, the morphological mutation weight value of the morphological transition area is calculated and compared with the feature deviation of the historical morphological database to form a micro-variation risk probability map.

[0061] Based on the micro-variation risk probability map, the lake morphology stability index is output, and a warning marker layer for high-risk areas is generated.

[0062] This invention begins by acquiring multi-source satellite remote sensing image data of the target area at multiple time points. This step aims to provide a stable, continuous, and multi-dimensional data foundation for subsequent lake morphology identification, trend analysis, and abrupt change detection.

[0063] Specifically, the multi-source satellite remote sensing image data includes, but is not limited to, optical remote sensing images, synthetic aperture radar (SAR) images, and thermal infrared remote sensing images. Wherein:

[0064] Optical remote sensing imagery possesses excellent spatial and spectral resolution, making it suitable for extracting water body boundaries, analyzing water spectral characteristics, and identifying the boundary between water and land. Under cloudless and fog-free observation conditions, optical imagery can clearly reflect the spatial distribution and shoreline structure of lakes.

[0065] SAR radar imagery has all-weather, day-and-night observation capabilities, enabling stable acquisition of surface information under complex weather conditions such as clouds, fog, or nighttime. SAR imagery is particularly suitable for identifying morphological changes in lakes during rainy seasons or extreme weather conditions, and can effectively assist in the completion and correction of invalid areas in optical imagery.

[0066] Thermal infrared remote sensing images are mainly used to monitor the difference in thermal radiation between water bodies and the surrounding land surface, indirectly reflecting physical properties such as water depth, flow rate and evaporation rate. They are highly sensitive for identifying seasonal shoals and dry-up lake edges.

[0067] During implementation, multi-orbit satellite data sources covering the target lake area are preferred. These may include publicly available remote sensing platforms such as the Landsat series, Sentinel series, Gaofen series, Terra / Aqua MODIS, RADARSAT, and ALOS PALSAR, to ensure the complementarity of imagery in terms of spatial coverage, temporal resolution, and band combination.

[0068] The multiple time points can be set according to the periodicity and suddenness of lake morphological changes. Preferably, a combination of seasonal scale (e.g., once per quarter) and monthly scale (e.g., twice per month) is used to ensure that long-term trend evolution can be captured and that morphological responses caused by short-term mutations or external disturbances (e.g., rainstorms or artificial water diversion) can be reflected in a timely manner.

[0069] After acquisition, the aforementioned multi-source remote sensing images will undergo unified preprocessing, including image radiometric calibration, atmospheric correction, geometric registration, and data cropping, to ensure consistency between different source images on the time and spatial axes.

[0070] After acquiring optical images, SAR radar images, and thermal infrared remote sensing images of the target area at multiple time points, in order to ensure the spatiotemporal consistency and information comparability between different data sources, it is necessary to perform unified preprocessing operations on the remote sensing image data, including radiometric correction, atmospheric correction, and geometric registration, so as to construct a multi-source image sequence under a unified time reference.

[0071] Radiometric correction is primarily used to eliminate grayscale errors in remote sensing images caused by factors such as sensor response characteristics, variations in solar altitude angle, and differences in surface reflectance. This step converts the digital values ​​(DN values) of the original image into physically meaningful radiance or surface reflectance to ensure the contrast in radiance intensity between different images. Preferably, for optical images, an absolute radiometric correction method based on sensor metadata can be used; for thermal infrared images, a thermal radiation transfer model can be introduced to correct for thermal reflectance discrepancies between water and land.

[0072] Atmospheric correction is used to eliminate distortions in remote sensing signals caused by atmospheric interference such as atmospheric molecular scattering and aerosol absorption. For optical remote sensing data, algorithms such as Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) or Dark Object Subtraction (DOS) can be used for correction. For thermal infrared data, atmospheric radiative transfer models such as MODTRAN can be introduced to retrieve path radiatives. SAR imagery is generally unaffected by visible light atmospheric interference, but if it contains differences between polarization channels, electromagnetic radiation consistency correction is also necessary.

[0073] Geometric registration is used to ensure the spatial consistency of various remote sensing images and avoid morphological extraction errors caused by image misalignment. Specifically, remote sensing images with high spatial resolution and high geographic accuracy (such as Landsat or Gaofen satellite data) are used as master images, and the remaining images are used as slave images. Registration is performed based on feature point matching algorithms (such as SIFT, ORB, etc.) combined with affine or polynomial geometric transformation models, with the error preferably controlled within 1 pixel.

[0074] After completing the radiometric correction, atmospheric correction, and geometric registration processes described above, all images are projected onto a unified coordinate system (such as WGS84 / UTM Zone) and organized into a multi-source image sequence in chronological order. This image sequence not only maintains the temporal logical relationship between images at different time points but also ensures spatial coverage consistency and spectral response comparability across different sources, providing a data foundation for subsequent lake boundary extraction, morphological change analysis, and abrupt change path modeling.

[0075] In this invention, a lake initial boundary extraction model is constructed based on a fused image sequence. The initial morphological contour of the lake boundary is extracted through multi-scale edge detection and water spectral feature matching. Specifically, the construction of the lake initial boundary extraction model includes the following steps:

[0076] Preprocessed optical remote sensing images, SAR radar images, and thermal infrared images were spatially aligned, and wavelet transform and principal component analysis (PCA) were used as the fusion framework. During the fusion process, each source image was first decomposed into three scale levels using discrete wavelet transform to extract high-frequency and low-frequency components. Then, based on PCA, principal component weights were calculated for the multi-source components at the same scale, and the fused image was reconstructed along the direction of maximum variance. The final multi-scale fused image not only preserved the spatial texture information of water body boundaries but also enhanced the spectral contrast of water bodies at different spectral bands.

[0077] The fused image is then processed at various scales using an improved Canny edge detection operator and a Laplacian operator for edge enhancement. The improved Canny algorithm incorporates a multi-directional non-maximum suppression strategy during the gradient magnitude calculation stage and dynamically adjusts the high and low threshold ranges during the dual-threshold decision stage. Preferably, the high threshold is set to 1.5 times the average gradient value of the image, and the low threshold is set to 0.5 times it, thus enhancing the response capability to weak edges. The Laplacian operator is used to supplement the second derivative change information at the boundaries, compensating for the Canny algorithm's shortcomings in edge detection in regions with complex textures. The detection results are then processed through morphological closing operations and edge concatenation to output an edge candidate set image.

[0078] Each edge pixel in the candidate edge set is cross-validated with the NDWI (Water Index) from the optical image and the backscattering coefficient from the SAR image. The NDWI is calculated by dividing the difference between the green band and the near-infrared band by its sum; the SAR backscattering value is processed by a Gamma-MAP filter, and low-value areas usually correspond to water bodies. If an edge pixel simultaneously satisfies an NDWI value greater than 0.2 and a SAR backscattering coefficient less than -15 dB, then that pixel is assigned a higher confidence weight. The confidence values ​​of all edge pixels are normalized within the range of 0 to 1 to form a boundary confidence map, which is used for determining the region growing seed for subsequent boundary extraction.

[0079] Edge pixels with a confidence value greater than 0.7 in the boundary confidence map were used as initial seed points, and an adaptive region growing algorithm based on spectral consistency and texture similarity was employed for region expansion. The growth conditions for adjacent pixels included: NDWI difference less than 0.1, grayscale standard deviation less than 5, and texture correlation coefficient greater than 0.85; pixels meeting these conditions were included in the boundary region. The grown boundary underwent edge tracking processing to obtain continuous closed boundary lines. Finally, the boundary contour was fitted using Catmull-Rom spline curves to form a smooth initial lake morphology contour curve with geometrically consistent representation. Simultaneously, a fitting error map was extracted, and the boundary deviation range was recorded as a reference for subsequent dynamic evolution detection.

[0080] The initial morphological contour is transformed into a lake morphological feature map, and features including shoreline curvature distribution index, morphological compactness index, boundary instability factor, and water body distribution centroid value are extracted.

[0081] Method for extracting shoreline curvature distribution index: The lake boundary outline is discretized into a sequence of ordered sampling points at equal intervals, with the sampling interval preferably set to 5 to 10 pixels. For each sampling point, the curvature value is calculated using the three-point sliding method, that is, based on the arc formed by the current point and the two adjacent points, the corresponding curvature is calculated, and the curvature is defined as the reciprocal of the radius of the arc. The resulting curvature value sequence is segmented and statistically analyzed according to the magnitude of curvature to construct a curvature distribution histogram. The curvature distribution index is defined as the ratio between the number of pixels in the segment with the maximum curvature and the average curvature of all sampling points. This value is used to reflect the degree of high curvature concentration of the lake shoreline in a local area, that is, the "meandering density" of the shoreline.

[0082] Method for extracting the morphological compactness index: The area A and total boundary length P of the region enclosed by the lake's boundary contour are calculated. The area A can be calculated using polygon face accumulation (e.g., the Shoelace algorithm), and the boundary length P can be obtained by accumulating the distances between points. Then, the compactness index is calculated using the compactness formula: Compactness Index = π is 3.1416. The compactness index is used to assess the similarity between the shape of a lake and a standard geometric circle. The closer the value is to 1, the more regular the lake boundary and the more compact the shape; the smaller the value, the more fragmented or extended the lake boundary.

[0083] Boundary instability factor extraction method: The lake boundary contour extracted at the current time node is registered one-to-one with the boundary contour at the previous time node. A point-to-point matching strategy based on minimum Euclidean distance is preferred; that is, for each current boundary point, its nearest corresponding point in the previous frame is found, forming a boundary difference sequence. The degree of change of each boundary segment is calculated based on the position offset rate (rate of change) between points. If the offset rate exceeds a set threshold (preferably 5%), it is considered an unstable boundary segment. Finally, the total length ratio of all unstable boundary segments is statistically analyzed as the boundary instability factor, used to measure the dynamic fluctuation of the lake shoreline over time.

[0084] Method for extracting the centroid value of water body distribution: The lake area is converted into a binary mask image, where the pixel value of water body is 1 and the value of other areas is 0. The coordinates of all water body pixels are iterated, and the row and column coordinates of each pixel are denoted as xi and yi, respectively. The total number of water body pixels is n. By averaging the coordinates of all water body pixels along the X and Y axes, the horizontal and vertical coordinates of the centroid are obtained, denoted as the centroid value G(x,y). Its calculation method is as follows: The water body's center of gravity value can be used to analyze the overall spatial shift trend of lake water bodies, and is especially suitable for tracking the contraction and expansion centers of lakes in arid or water diversion areas.

[0085] A temporal feature matrix was constructed based on morphological feature maps to establish a lake morphological variation path map model. A dynamic time warping algorithm was then used to align and match variation paths between adjacent time periods, identifying abrupt change nodes and morphological transition regions. Specifically, this includes:

[0086] The extracted lake morphological features, including shoreline curvature distribution index, morphological compactness index, boundary instability factor, and water body centroid value, are arranged in chronological order to construct an original temporal feature matrix. Each row of this matrix corresponds to a time node, and each column corresponds to a morphological feature index. Since direct comparison of the various feature indices will introduce bias due to differences in their dimensions and numerical scales, the matrix needs to be standardized.

[0087] Preferably, each column in the matrix is ​​processed using a zero-mean, unit-variance normalization method. Specifically, the mean of a column is subtracted from all eigenvalues ​​in that column, and then divided by the standard deviation of that column. This results in a normalized eigenvalue column with a statistical distribution of mean 0 and standard deviation of 1, ensuring that the relative importance of different eigenvalues ​​is consistent in the calculation.

[0088] To compress data dimensionality and extract the main patterns in lake morphological changes, Principal Component Analysis (PCA) was used to reduce the dimensionality of the normalized time-series feature matrix. Specific steps included calculating the covariance matrix, eigenvalue decomposition, and sorting and selecting the top principal components, ensuring their cumulative variance contribution rate was no less than 90%. This method yielded a low-dimensional embedding vector sequence, where each vector represented the projected coordinates of the morphological state at the corresponding time point in the principal component space, preserving the main trends while eliminating noise and redundant features.

[0089] Each vector in the aforementioned low-dimensional embedding vector sequence is treated as a node in the mutation path graph, with node numbers arranged chronologically. Any two adjacent time nodes are connected by an edge, representing the degree of morphological change between adjacent time periods. To more accurately characterize this change, the edge weight is defined as a weighted combination of the Dynamic Time Warping Distance (DTW distance) and the cosine similarity distance, calculated as follows: Edge Weight Distance = 0.6 × DTW Distance + 0.4 × Cosine Distance. Here, the DTW distance measures the matching difference of the time series at non-linear time scales, while the cosine distance measures the difference in feature space angles. The combined weighting coefficients can be fine-tuned according to the specific lake type, with 0.6 and 0.4 being preferred to balance the influence of temporal elasticity and spatial structure changes.

[0090] The constructed path graph is a weighted directed graph. To reduce the influence of interfering edges, a threshold pruning strategy is used to simplify the graph structure. A preferred edge weight threshold τ = 0.5 is set, deleting all edges with weights less than τ. This operation can be seen as filtering weakly correlated periods, reinforcing the main trend of change.

[0091] Subsequently, a shortest path search algorithm (such as Dijkstra's algorithm or A* search algorithm) is used to traverse all reachable paths starting from the initial time node, and the path with the minimum cumulative weight is identified and defined as the main path of lake morphological variation. At the same time, nodes with local maximum edge weights are identified on the main path and marked as candidate mutation nodes, which serve as input for subsequent mutation detection and analysis.

[0092] After constructing the mutation path map, this invention further uses the Dynamic Time Warping (DTW) algorithm to perform a fine comparison of the morphological evolution process within adjacent time periods to identify local mutation regions and abrupt change events. The processing includes the following steps:

[0093] Low-dimensional embedded feature vectors at adjacent time points are extracted from the lake morphological variation path map model and arranged into time series pairs using a sliding window approach. For example, by using three time points as a window, a set of local variation sub-paths containing the early, middle, and late periods is constructed, enabling the model to simultaneously analyze continuous change trends and local mutation behaviors.

[0094] For each pair of variant sub-path sequences, the DTW algorithm is applied for nonlinear alignment. DTW selects the path with the minimum total cost as the optimal match by finding all possible alignment paths between two time series sequences, allowing for flexible compression and expansion on the time axis, thus better adapting to actual non-uniform variation processes.

[0095] The cost matrix construction rules for DTW are as follows: the cost between each time point pair is defined as the Euclidean distance between that time point and the feature space. The total path cost is the sum of the costs of all matching point pairs. After calculation, this alignment cost is used as a quantification of the degree of change between the current sub-paths.

[0096] The DTW alignment costs of all subpaths are collected and normalized to obtain a standardized cost sequence. The mean μ and standard deviation σ are calculated, and a mutation identification threshold is set to μ plus twice σ (i.e., mutation threshold = μ + 2σ). When the DTW cost value of a subpath exceeds this threshold, its corresponding time point is considered to be in an anomalous alignment state, and this time point is marked as a candidate mutation node. This method can effectively identify local anomalous change points from the global evolution trend.

[0097] For time windows marked as candidate mutation nodes, the original remote sensing image sequence is traced back to extract the lake boundary change rate map and the water body centroid offset vector field corresponding to that time period. The boundary change rate is calculated by the overlap between the current and previous frames, while the centroid offset is calculated by the vector difference based on the change in the centroid position of the water body pixels.

[0098] Based on the two spatial characteristics mentioned above, a consistency index is calculated between the intensity of boundary changes and the direction of centroid shift within the abrupt change window (e.g., an angle less than 30 degrees indicates high consistency). If this index exceeds a set threshold (preferably, the condition is met within 70% of the area), the region is identified as a morphological abrupt change region. Finally, a spatial annotation layer is generated on the lake image map for this abrupt change region, which can be used for visualization and integration with dynamic monitoring and early warning systems.

[0099] Based on the spatial distribution of mutation nodes and the amplitude of shoreline curvature change within the morphological transition area, the morphological mutation weight value of the morphological transition area is calculated and compared with the feature deviation of the historical morphological database to form a micro-variation risk probability map.

[0100] A quantitative risk weight assessment is performed on each morphological transition region, which comprehensively considers the local geometric intensity of the current boundary change (amplitude of shoreline curvature change) and the location sensitivity of the abrupt event in the lake's spatial structure (degree of relative center shift). The specific implementation steps are as follows:

[0101] For each identified morphological transition region, all boundary points on its corresponding boundary contour are extracted and arranged into an ordered point set according to the original contour order. A three-point sliding method is used to calculate the local curvature of the point set. Specifically, for each boundary point, a ternary point set is formed with its two adjacent points (one before and one after), a fitted circular arc is constructed, and the reciprocal of the radius of this arc is used as the curvature value of that point. The sliding window size is preferably set to a step size of 3 to 5 pixels to ensure the locality and sensitivity of the calculation. Finally, a continuous curvature change curve for the transition region is generated, which can be used to analyze the intensity of boundary fluctuations.

[0102] Based on the curvature curves described above, the maximum and minimum curvature values ​​are extracted, and the difference between the two is used to define the curvature amplitude ΔK in this region. This index reflects the extreme fluctuation in the complexity of the shoreline morphology within the transition area; a larger value indicates more severe local deformation and greater instability in morphological changes.

[0103] Extract the position coordinates of the current mutation node in the overall space of the lake, denoted as P(xp, yp), and calculate its Euclidean distance D relative to the centroid coordinates G(xg, yg) of the current lake water body. .

[0104] To eliminate the influence of lake-scale differences, the distance D is normalized using the maximum boundary radius Rmax as the normalization benchmark: Normalized Distance Where Rmax is the maximum distance from the water body boundary to the centroid at the current time point, used to represent the maximum possible range of offset. The larger the D' value, the further the abrupt change event deviates from the main structural zone of the lake, which may have a more severe impact on the local boundary but a smaller overall impact; conversely, the smaller the D' value, the more sensitive the abrupt change is to the central structure.

[0105] Based on the curvature amplitude ΔK and the normalized spatial distance D', a comprehensive morphological abrupt change weight W is constructed, and the calculation formula is as follows: Wherein, α and β are weighting coefficients, preferably α = 0.7 and β = 0.3, to emphasize the dominant role of geometrical abrupt changes in the overall weight. The larger the obtained W value, the greater the impact of the abrupt change region on the lake structure at the current moment, serving as the basic input parameter for subsequent risk probability calculations.

[0106] To improve the accuracy of identifying anomalies in abrupt change regions, this invention introduces historical time-series morphological data for comparative analysis. By assessing the degree of difference between the current mutation characteristics and the historical morphological benchmark, the risk of morphological deviation is quantified, and a final risk probability distribution map is formed by combining mutation weights.

[0107] The morphological data of the target lake over past monitoring periods were archived and organized to construct a historical morphological database, including:

[0108] Morphological feature map of each time point (including curvature, compactness, instability factor, and center of gravity position);

[0109] Its corresponding boundary location and region identifier in two-dimensional space;

[0110] Historical abrupt events and their impact tags (such as records of low water levels, siltation, water diversion, etc.).

[0111] The database can be organized using a two-dimensional matrix and geographic index, which facilitates querying the morphological evolution trends and characteristic statistics of a specific spatial region within a specific time window.

[0112] The morphological feature vector of the current abrupt change region is extracted, including the N-dimensional features constituted by the aforementioned indicators (such as ΔK, compactness index, boundary instability factor, etc.). The historical morphological feature vector set corresponding to this spatial location (based on boundary coordinate matching) is then searched in the historical database. The Euclidean distance Dfeature between the current feature and its historical average feature is calculated, which is the square root of the sum of the squares of the differences between the corresponding features. This distance value reflects the degree to which the current abrupt change region deviates from the historical steady state in the morphological feature dimension, serving as an anomaly indicator.

[0113] The mutation weight W and the feature deviation value Dfeature are jointly normalized to unify their values ​​between 0 and 1. Then, a probability mapping function is used to convert them into a risk probability value P, preferably the Sigmoid function, defined as follows: Wherein, γ is a control parameter that regulates the smoothness of the mapping function, preferably γ = 2.0. The probability value P represents the risk of irreversible morphological changes occurring in the transition region at the current moment.

[0114] Based on the micro-variation risk probability map, a lake morphology stability index is output, and a warning marker layer for high-risk areas is generated. This part mainly includes the following steps:

[0115] First, a comprehensive statistical analysis is performed on the constructed micro-variation risk probability map. This risk probability map is represented in two-dimensional grid form, where each abrupt change region is assigned a probability value between 0 and 1 to quantify the risk level of structural variation in that region. To assess the risk distribution characteristics of the entire lake morphology system, the probability values ​​of all abrupt change regions in the map are statistically analyzed, and their mean μ and standard deviation σ are calculated, reflecting the average risk level and risk fluctuation amplitude across the entire lake area, respectively.

[0116] The specific calculation method is as follows: The risk probability values ​​of the n transition regions are sequentially denoted as... The mean μ is the sum of all probability values ​​divided by n; the standard deviation σ is the sum of the squares of the differences between all probability values ​​and the mean, divided by n, and then the square root is taken. This statistical result provides a basis for subsequent stability index modeling and threshold setting.

[0117] To comprehensively summarize the local risk probability values ​​across the entire lake area and form a quantitative index reflecting the overall morphological stability, this invention introduces the area-weighted risk mean as the lake morphological stability index, as shown in the following formula: Morphological Stability Index , where i ranges from 1 to n, and n is the total number of abrupt change regions. The resulting stability index S ranges from 0 to 1. The closer the value is to 0, the more stable the overall morphology of the lake is, and the sparse and low-intensity abrupt change risk is. The closer the value is to 1, the more widespread and strong the trend of abrupt morphological change exists, which requires management attention.

[0118] To accurately identify transitional zones that pose a significant potential threat to ecosystems, it is necessary to establish high-risk assessment thresholds. This invention establishes the following two categories of high-risk identification criteria:

[0119] Overall risk criterion: If the lake morphological stability index S is greater than 0.6, the lake is judged to be in a high-risk state.

[0120] Local risk criterion: If the risk probability value Pi of a single abrupt change region is greater than 0.8, the region is identified as a local high-risk region, and the flag can be triggered even if the overall index is low.

[0121] The aforementioned thresholds were determined based on experimental analysis and ecological management needs, and can be dynamically adjusted according to actual scenarios. Preferably, S = 0.6 and Pi = 0.8 are empirical thresholds, which demonstrate good anomaly sensitivity and risk identification capabilities in multiple typical lake cases.

[0122] Abrupt changes that meet any of the above high-risk criteria are designated as high-risk areas, and spatial layer generation is performed, mainly including:

[0123] Extract the boundary coordinate set of high-risk areas and spatially encode them in vector polygon format;

[0124] High-risk areas are located in the original remote sensing images and displayed in an enhanced image format, such as a heat map or a highlight color overlay.

[0125] Assign a risk level label to each high-risk area in the raster layer, and output the spatial coordinates, risk probability value and risk level information (such as high, medium and low) in combination with the attribute table.

[0126] Export the above layers to a standard spatial format (such as GeoTIFF, Shapefile, etc.) to facilitate integration and use with ecological monitoring platforms, emergency dispatch systems, or map service systems.

[0127] The resulting early warning label layer can be overlaid on the remote sensing base map in real time, clearly showing the spatial distribution of risk hotspots and providing an intuitive and scientific reference for water resource regulation, ecological protection plan formulation, and selection of key inspection areas.

[0128] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for monitoring the shape of a lake based on satellite remote sensing technology, characterized in that: The method comprises the following steps: Obtain multi-source satellite remote sensing image data of the target area at multiple time nodes, including optical images, SAR radar images and thermal infrared images; Perform radiation correction, atmospheric correction and geometric registration processing on the obtained remote sensing images to construct a multi-source image sequence under a unified time reference; Construct a lake initial boundary extraction model based on the fused image sequence, extract the initial shape profile of the lake boundary through multi-scale edge detection and water body spectral feature matching; Convert the initial shape profile into a lake shape feature map, and extract the shoreline curvature distribution index, shape compactness index, boundary instability factor and water body distribution barycenter value; Construct a time series feature matrix based on the shape feature map, establish a lake shape variation path graph model, and use dynamic time warping algorithm to align and match the variation paths of adjacent time periods to identify mutation nodes and shape jump regions; According to the spatial distribution position of the mutation node and the shoreline curvature amplitude in the shape jump region, calculate the shape mutation weight value of the shape jump region, and compare it with the historical shape database to form a micro variation risk probability map; Based on the micro variation risk probability map, output the lake shape stability index and generate a warning identification layer of the high-risk area. 2.The lake shape monitoring method based on satellite remote sensing technology according to claim 1, characterized in that: The method for constructing a lake initial boundary extraction model based on a fused image sequence comprises the following steps: Fuse the preprocessed optical, SAR and thermal infrared images at multiple resolution scales to output multi-scale fused images; Apply the improved Canny operator combined with the Laplace operator to the fused images at each scale for multi-scale edge detection to obtain an edge candidate set, and cross-verify the edge candidate set to generate a boundary confidence map; Based on the boundary confidence map, take the edge with high confidence as a seed, and use adaptive region growing and spline curve fitting to fuse to form the final lake initial shape profile. 3.The lake shape monitoring method based on satellite remote sensing technology according to claim 1, characterized in that: The extraction of the shoreline curvature distribution index, the shape compactness index, the boundary instability factor and the water body distribution barycenter value comprises the following steps: The extraction of the shoreline curvature distribution index comprises the following steps: discretize the lake initial boundary profile into an equidistant sampling point sequence, calculate the curvature value of each point using the three-point sliding method, and segmentally count the curvature values of all points to form a curvature distribution histogram, wherein the curvature is defined as the reciprocal of the circular arc formed by three consecutive points, and the curvature distribution index is the ratio of the maximum curvature segment to the average value of all curvatures; The extraction of the shape compactness index comprises: calculating the area A of the region surrounded by the lake boundary contour and the boundary length P thereof, according to the formula: shape compactness index = A / P ; and performing solving. The extraction of the boundary instability factor comprises the following steps: point-to-point match the boundary profile at the current time node with the boundary profile at the previous time node, count the proportion of boundary segments with a boundary change rate higher than a set threshold, and define it as the boundary instability factor; The extraction of the water body distribution gravity center value includes: converting a lake area into a binary mask image, taking water body pixel coordinates as input, calculating a weighted average position of all water body pixels on a two-dimensional plane, and the water body distribution gravity center value wherein n is the total number of water body pixels, xi and yi are the coordinates of the i th pixel.

4. The lake shape monitoring method based on satellite remote sensing technology according to claim 1, characterized in that: The method for constructing a time series feature matrix based on a shape feature map and establishing a lake shape variation path graph model comprises the following steps: Arrange the shape features extracted at each time node in time sequence to form an original time series feature matrix, and perform zero-mean unit-variance normalization on each feature column; Perform principal component analysis on the normalized time series matrix, retain principal components with a cumulative variance contribution rate not less than 90%, and obtain a low-dimensional time series embedding representation; The low-dimensional representation of each time node is taken as a node of a path graph, and a weighted distance of a dynamic time warping similarity and a cosine similarity is calculated as an edge weight between adjacent time nodes; The constructed weighted path graph is pruned by a threshold, and a shortest path search is used to identify a morphological variation main path and a mutation node.

5. The method for monitoring lake morphology based on satellite remote sensing technology according to claim 4, characterized in that: The alignment and matching of the variation paths of adjacent time periods by using the dynamic time warping algorithm to identify the mutation node and the morphological jump region, comprising: Based on the low-dimensional time sequence characteristics of each node in the lake morphological variation path graph model, a feature sequence pair between adjacent time nodes is extracted, and a local variation sub-path set is constructed; The dynamic time warping algorithm is used for nonlinear alignment of each pair of sub-path sequences, and the optimal paired path and its alignment cost matrix are calculated as a quantitative expression of the influence of the time scale change on the morphological evolution; The DTW cost value of all paired paths is normalized, and an abnormal alignment region is identified by setting a mutation threshold, which is marked as a candidate mutation node; In the time window of the candidate node, the corresponding boundary change rate graph and the center of gravity offset vector field in the original image sequence are traced back, and the morphological jump region is finally determined and the spatial annotation result is generated by combining the morphological jump strength and the spatial position overlap degree. 6.The lake shape monitoring method based on satellite remote sensing technology according to claim 1, characterized in that: The morphological mutation weight value of the morphological jump region is calculated according to the spatial distribution position of the mutation node and the shoreline curvature amplitude in the morphological jump region, comprising: The curvature values of all boundary points in the identified morphological jump region are extracted, a three-point sliding method is used for local curvature calculation, and a curvature change curve is constructed; The difference between the maximum curvature value and the minimum curvature value in the jump region is calculated as a curvature amplitude index of the region; According to the spatial distribution position of the mutation node in the whole lake region, the normalized distance from the lake center of gravity is calculated to reflect the central sensitivity of the mutation; The curvature amplitude index and the normalized spatial distance are combined by weighting to construct the morphological mutation weight value.

7. The method for monitoring lake morphology based on satellite remote sensing technology according to claim 6, characterized in that: The feature deviation comparison with the historical morphological database is performed to form a micro-variation risk probability graph, comprising: A historical morphological database covering the target lake is established to store the morphological feature spectrum and the corresponding spatial distribution information of multiple time nodes in the past; The Euclidean distance between the morphological feature vector of the current mutation region and the average feature vector of the historical region at the same position is calculated to obtain the feature deviation value; The morphological mutation weight value and the feature deviation value are jointly normalized, and a probability mapping function is used to convert the micro-variation risk probability value; The risk probability values of all jump regions are mapped to the spatial coordinate system to form a micro-variation risk probability graph containing the mutation degree and the spatial position. 8.The lake shape monitoring method based on satellite remote sensing technology according to claim 1, characterized in that: Based on the micro-variation risk probability graph, a lake morphological stability index is output, comprising: Statistical analysis is performed on the risk probability values of all jump regions in the micro-variation risk probability graph, and the mean and standard deviation are calculated; Based on the risk probability values of each jump region and the spatial area, an area-weighted risk mean value is constructed, which is defined as the lake morphological stability index; A high-risk judgment threshold is set, and the lake morphological stability index is greater than 0.6 or the risk probability value of a single region exceeds 0.8, which is judged as a high-risk region. The area meeting the high risk condition is highlighted on the remote sensing image in the form of a layer, a warning identification layer is generated, and spatial coordinates and corresponding risk levels are output.

Images