Lake ecological meteorological disaster risk dynamic monitoring system based on satellite remote sensing

Abstract

This invention relates to the field of remote sensing image monitoring technology and discloses a dynamic monitoring system for lake ecological meteorological disaster risks based on satellite remote sensing. The system comprises four core processing units: multi-source image preprocessing, algal aggregation morphology analysis, shoreline exposure risk assessment, and dynamic early warning generation. The system first performs spatiotemporal normalization on multi-source satellite data and identifies key monitoring areas. Then, by calculating the directional texture intensity and lateral boundary change rate, it converts the two-dimensional distribution area into an equivalent contact length reflecting the striped aggregation morphology. Based on this length, it quantifies the probability of sudden shoreline contact and generates a single observation risk value by combining it with the local accumulation magnification factor. Finally, it introduces cloud interference suppression and time decay mechanisms to generate a cumulative risk score. This invention can effectively identify the instantaneous extreme exposure risk of narrow algal bloom strips at water intakes.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image monitoring technology, and more specifically, to a dynamic monitoring system for lake ecological meteorological disaster risks based on satellite remote sensing. Background Technology

[0002] Lake water environment safety is a crucial component of water supply security and ecological civilization construction. With the development of satellite remote sensing technology, monitoring algal blooms in lakes using multi-source optical satellite imagery has become the mainstream method. Existing technologies typically employ either chlorophyll a concentration inversion or bloom coverage extraction: the former uses band ratios or machine learning algorithms to invert the average concentration of the water surface; the latter divides the image into bloom and non-bloom areas by setting thresholds, then calculates the total coverage area of ​​the entire lake or specific zones. However, these conventional monitoring methods have significant blind spots when dealing with the risk of sudden water quality deterioration at water intakes or highly sensitive coastal areas. Existing technologies often overlook the fact that pollutant distribution is relatively uniform or smoothly gradual. However, in real shallow eutrophic lakes, under the combined influence of wind and water flow, floating algae often do not spread evenly but exhibit a significant strip-like aggregation pattern. This strip-like aggregation pattern has extremely strong spatial anisotropy. Unlike large-area sheet-like aggregations, narrow strips are more easily drifted rapidly by wind. When these high-density algal streaks are pushed by the wind towards the shore or water intake, they can cause localized accumulation in a very short time, leading to an instantaneous increase in algal density near the shoreline to more than a thousand times the background value. Existing monitoring methods based on total area or average concentration essentially average out this extreme situation. They cannot distinguish the essential difference in hazard risk between large areas of low-density thin layers and small areas of high-density streaks. Specifically, existing technologies cannot effectively quantify the potential contact pressure exerted on the shoreline by this streak-like morphology, nor can they capture the sudden, extremely high risk generated by the interaction between morphology and boundaries. Summary of the Invention

[0003] This invention provides a dynamic monitoring system for lake ecological meteorological disaster risks based on satellite remote sensing, which solves the technical problems mentioned in the background art.

[0004] This invention provides a dynamic monitoring system for lake eco-meteorological disaster risks based on satellite remote sensing, comprising: The multi-source image preprocessing unit is used to receive multi-source satellite remote sensing data and map it to a unified grid, generate water reflectance data that eliminates the geometric influence of sunlight, and delineate key monitoring areas including water intakes based on shoreline distance. The algal aggregation morphology analysis unit is used to generate a continuous probability map of algal distribution based on the water reflectance data, calculate the directional texture intensity index reflecting the degree of algal film striping and the lateral boundary change rate index reflecting the gradient of the strip edge, and calculate the weighted distribution area based on the continuous probability map of algal distribution and the directional texture intensity index, and then combine the lateral boundary change rate index to convert the weighted distribution area into the equivalent total strip length. The shoreline exposure risk assessment unit is used to map the total length of the equivalent strip to the unit shoreline contact intensity for the key monitoring area, calculate the probability of sudden shoreline contact using a random shoreline contact model, and generate a single observation risk value by combining the background water quality benchmark and the preset local accumulation amplification factor. The dynamic early warning generation unit is used to generate a cloud interference suppression coefficient based on cloud cover detection, apply the cloud interference suppression coefficient and time decay factor to the single observation risk value to generate a cumulative risk score, and output an alarm signal when the cumulative risk score exceeds a preset alarm threshold.

[0005] The beneficial effects of this invention are as follows: By analyzing the banded aggregation morphology of cyanobacterial blooms under the influence of wind and waves, it effectively solves the technical problem that traditional large-area averaging monitoring methods cannot capture localized sudden high-value risks. This invention utilizes directional texture analysis technology to transform the two-dimensional coverage area into an effective contact length reflecting the disaster-causing potential energy, and combines the shoreline distance field with a random shoreline contact model to quantify the instantaneous extreme exposure pressure caused by wind-driven elongated aggregates on the water intake. This invention can achieve dynamic early warning of low-probability, high-intensity disaster events in highly sensitive shoreline areas using only satellite imagery. Attached Figure Description

[0006] Figure 1 This is a block diagram of the dynamic monitoring system for lake ecological meteorological disaster risks based on satellite remote sensing, as described in this invention. Figure 2 This is a schematic diagram of a specific implementation scenario of the present invention. Detailed Implementation

[0007] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0008] It should be noted that, unless otherwise defined, the technical or scientific terms used in one or more embodiments of the present invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in one or more embodiments of the present invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "include" or "comprising" indicate that the element or object preceding the term encompasses the element or object listed following the term and its equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "up," "down," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0009] like Figure 1 As shown, the dynamic monitoring system for lake eco-meteorological disaster risks based on satellite remote sensing includes: The multi-source image preprocessing unit is used to receive multi-source satellite remote sensing data and map it to a unified grid, generate water reflectance data that eliminates the geometric influence of sunlight, and delineate key monitoring areas including water intakes based on shoreline distance. The algal aggregation morphology analysis unit is used to generate a continuous probability map of algal distribution based on the water reflectance data, calculate the directional texture intensity index reflecting the degree of algal film striping and the lateral boundary change rate index reflecting the gradient of the strip edge, and calculate the weighted distribution area based on the continuous probability map of algal distribution and the directional texture intensity index, and then combine the lateral boundary change rate index to convert the weighted distribution area into the equivalent total strip length. The shoreline exposure risk assessment unit is used to map the total length of the equivalent strip to the unit shoreline contact intensity for the key monitoring area, calculate the probability of sudden shoreline contact using a random shoreline contact model, and generate a single observation risk value by combining the background water quality benchmark and the preset local accumulation amplification factor. The dynamic early warning generation unit is used to generate a cloud interference suppression coefficient based on cloud cover detection, apply the cloud interference suppression coefficient and time decay factor to the single observation risk value to generate a cumulative risk score, and output an alarm signal when the cumulative risk score exceeds a preset alarm threshold.

[0010] Preferably, when the multi-source image preprocessing unit receives multi-source satellite remote sensing data and maps it to a unified grid to generate water reflectance data that eliminates the geometric influence of sunlight, it specifically performs the following steps: Construct a set of multi-temporal surface reflectance products that meet the lake area coverage threshold within a preset time period before the current time, as the multi-source satellite remote sensing data; The multi-source satellite remote sensing data is projected onto a universal transverse Mercator projection coordinate system with the lake centroid as the reference, and then resampled to the unified grid with a resolution of ten meters using bilinear interpolation. Read the solar zenith angle parameter from the multi-source satellite remote sensing data; The water reflectance data, which eliminates the geometric effects of sunlight, is generated based on the following formula: in, Indicates the band, and t represents time. This represents the atmospheric-corrected surface reflectance on the unified grid. Represents pi (π). Indicates the initial remote sensing reflectance. The solar zenith angle parameter is expressed in radians. Represents the cosine function. This represents the water reflectance data after eliminating the geometric effects of sunlight.

[0011] The preset duration is the length of time preceding the current moment selected when constructing a multi-source satellite remote sensing dataset. A preferred duration is 48 hours, because cyanobacterial blooms in shallow eutrophic lakes exhibit rapid redistribution on an hourly basis; a 48-hour duration can cover sufficient spatiotemporal variation information while avoiding reduced processing efficiency due to excessive data volume.

[0012] The lake area coverage threshold is the minimum effective coverage ratio of the imagery for the lake area required when screening satellite remote sensing data. An optimal threshold is 80%, which ensures that the imagery includes the core monitoring area of ​​the lake while providing a reasonable margin for error in image screening to avoid missing valid data.

[0013] The centroid coordinates of a lake are the coordinates of the center of its geographical extent, representing its geometric center. They can be obtained using the centroid calculation tools of a geographic information system (GIS) based on the lake's vector boundary data.

[0014] Grid resolution is the unified pixel spatial resolution after resampling of multi-source satellite remote sensing data, representing the actual geographic scale of the pixels. A resolution of 10 meters is preferred, as this balances spatial detail of the remote sensing imagery with data processing efficiency, and is suitable for the detailed identification of strip-like aggregations of cyanobacterial blooms.

[0015] The solar zenith angle is the angle between sunlight and the local zenith direction, representing the geometric angular characteristics of sunlight. It can be directly read from the metadata file of satellite remote sensing imagery, or calculated using a solar position algorithm combined with the image's imaging time and geographical location.

[0016] Surface reflectance is the ratio of reflected radiant flux to incident radiant flux at the surface of a lake after atmospheric correction, characterizing the reflective optical properties of the surface. It can be directly obtained from secondary surface reflectance products obtained through satellite remote sensing without the need for additional atmospheric correction.

[0017] A band is a range of different spectral wavelengths used by satellite remote sensing sensors to acquire data; different bands correspond to different spectral information. This information can be directly read from the band identifiers and metadata files of satellite remote sensing imagery, with each band corresponding to a fixed center wavelength.

[0018] The observation time is the specific point in time when satellite remote sensing imagery completes imaging of the lake area, representing the temporal attribute of the image. It can be directly read from the metadata file of the satellite remote sensing imagery, containing specific year, month, day, hour, and minute information.

[0019] The spatiotemporal normalization processing of multi-source satellite remote sensing data includes: firstly, projecting all satellite remote sensing data onto a universal transverse Mercator projection coordinate system based on the lake's centroid. This coordinate system ensures that the projection of the lake area is free from significant distortion. Then, resampling is performed using bilinear interpolation, which calculates the pixel value by taking the weighted average of the four pixels surrounding the point to be interpolated. This method improves processing speed while maintaining accuracy. Finally, the data is unified to a 10-meter resolution grid, achieving complete consistency in spatial reference and pixel resolution for multi-source, multi-temporal data. For example, when processing images from Sentinel-2 and Landsat-8, both are resampled to this 10-meter resolution projection grid to eliminate spatial scale differences.

[0020] Optical correction to eliminate the geometric effects of solar illumination includes: first, reading the solar zenith angle from satellite metadata; dividing the atmospherically corrected surface reflectance by pi to complete the standard conversion from hemispherical reflectance to remotely sensed water reflectance, obtaining the initial remotely sensed reflectance; then, dividing the initial remotely sensed reflectance by the cosine of the solar zenith angle, which is positively correlated with the incident solar radiation flux at the Earth's surface, can eliminate reflectance differences caused by different illumination angles at different observation times. For example, after this correction, the reflectance of the same water area in images from 9 AM and 3 PM can accurately represent the true optical characteristics.

[0021] The selection of multi-source satellite remote sensing data includes: selecting only multi-temporal surface reflectance products within a preset time period prior to the current moment that have a lake area coverage rate reaching a threshold. This ensures that the selected images are continuous in time and effectively cover the lake spatially, avoiding the loss of information on algal bloom changes due to excessively large time spans, or invalid monitoring data due to incomplete coverage. For example, if the current moment is 12:00 noon on a certain day, the preset time period is 48 hours, and the coverage rate threshold is 80%, then only images with an effective lake coverage rate of not less than 80% within the previous 48 hours will be selected.

[0022] The preset duration is 48 hours. If no image meets the lake area coverage threshold within 48 hours, the nearest image can be selected within 7 days to ensure that the dataset always has valid data and avoid empty sets.

[0023] The coverage threshold for lake areas is preferably set at 80%, which is a commonly used effective coverage ratio for water body monitoring in remote sensing images, taking into account both data validity and screening flexibility.

[0024] The general transverse Mercator projection is divided into zones based on the latitude and longitude coordinates of the lake's centroid. The EPSG326XX series is selected for the Northern Hemisphere, and the EPSG327XX series is selected for the Southern Hemisphere. XX represents the corresponding 6-degree zone projection number. For example, if the lake's centroid is located at 120 degrees east longitude and 31 degrees north latitude, it corresponds to EPSG32651, ensuring that the projection is distortion-free.

[0025] Bilinear interpolation resampling uses four-neighbor bilinear interpolation. Taking the center of the pixel to be interpolated as the origin, the coordinates and pixel values ​​of the four surrounding neighboring pixels are taken, and the value of the pixel to be interpolated is calculated by the linear interpolation formula. The weight decreases linearly with the distance. This algorithm is simple to implement and its accuracy meets the monitoring requirements.

[0026] The solar zenith angle is uniformly expressed in radians. If the unit read from the metadata is angle, it needs to be converted using the formula that radians equal angle multiplied by pi divided by 180 to ensure that the subsequent cosine value calculation is consistent and to avoid calculation errors.

[0027] Preferably, when the multi-source image preprocessing unit delineates the key monitoring area including the water intake based on the shoreline distance, it specifically performs the following steps: Extract the green and near-infrared band values ​​from the water reflectance data after eliminating the geometric effects of sunlight, and generate a normalized differential water index based on the following formula: A continuous nonlinear mapping is performed on the normalized differential water index using an S-shaped function that includes a kurtosis parameter and a center threshold parameter, and water body soft mask data is generated based on the following formula: The contour lines with a probability value of 0.5 in the water body soft mask data are extracted as dynamic shorelines, and the Euclidean distance from each pixel in the unified grid to the dynamic shoreline is calculated. The key monitoring area is constructed based on the following formula: in, Represents pixel coordinates, and These represent the reflectance values ​​for the green light band and the near-infrared band, respectively. This represents a small constant to prevent the denominator from being zero. This represents the normalized difference water index. This represents the steepness parameter. This represents the center threshold parameter. Represents an exponential function. This represents the data of the soft mask for the water body. B represents the Euclidean distance, and B represents the preset shoreline buffer distance threshold. This represents the coordinates of the pixel point to the i-th preset water intake point. distance, This represents the preset threshold for the neighborhood radius of the water intake. This represents the logical AND operation. This indicates the key monitoring area for the i-th water intake.

[0028] The steepness parameter of the S-curve function of the water body soft mask is a parameter that controls the steepness of the curve change when the normalized differential water index is mapped to the probability value of the water body soft mask. A value of 25 is preferred, as this value allows for a moderate distinction between the probability mapping of water and non-water areas, avoiding both an excessively wide transition zone leading to blurred boundaries and an excessively narrow transition zone causing abrupt hard segmentation.

[0029] The center threshold parameter of the S-shaped function for the water body soft mask is the critical value of the index corresponding to a normalized differential water body index mapped to a water body soft mask probability value of 0.5. A value of 0.05 is preferred, as this value is suitable for the spectral characteristic differences between lake water bodies and shorelines / land, and can accurately distinguish lake water bodies from surrounding non-water areas.

[0030] The small constant used to prevent the denominator from being zero is an extremely small value added to the calculation of the Normalized Difference Water Index (NDDI) to avoid the calculation anomaly of both the numerator and denominator being zero. It is preferably 10 to the power of -6, a very small value that will not change the original calculation result of the NDDI, thus playing a stabilizing role in the numerical calculation.

[0031] The shoreline buffer distance threshold is a critical distance value defined from the dynamic shoreline into the lake interior, used to delineate the monitoring range near the shoreline. A value of 2000 meters is preferred, as this value is suitable for the shoreline aggregation characteristics of cyanobacterial blooms in shallow eutrophic lakes and can cover high-risk areas of shoreline blooms.

[0032] The coordinates of a water intake point are the specific geographical coordinates of a lake's water intake, representing its geographical location. These coordinates can be obtained through on-site measurements using the Global Positioning System (GPS) or directly retrieved from the geographic information archives of water conservancy projects.

[0033] The intake neighborhood radius threshold is a critical value representing the radius of a circular area drawn outward from the intake coordinates, used to define the core monitoring range of the intake. A value of 3000 meters is preferred, as this value covers the core area around the intake where cyanobacterial blooms may reach the shore and affect water intake safety.

[0034] The numerical processing of the normalized difference water index includes adding a small constant to the denominator to prevent it from reaching zero. This method only makes minor corrections at the numerical calculation level and does not affect the differences in index characteristics between water and non-water areas. It also completely avoids calculation interruptions caused by spectral anomalies, ensuring the stability of the algorithm. For example, when the reflectance values ​​in the green and near-infrared bands are both zero, adding this constant allows the index calculation to be completed normally.

[0035] The construction of the water body soft mask involves: using an S-shaped function with kurtosis and center threshold parameters to perform a continuous nonlinear mapping on the normalized differential water index, generating continuous probability values ​​between 0 and 1 as the water body soft mask data, instead of using the traditional hard threshold segmentation method. This method can preserve the fuzzy features of the transition zone between water and non-water bodies and avoid shoreline boundary distortion caused by hard segmentation. For example, the shallow water area near the shore of a lake will be assigned a probability value between 0 and 1, which is more in line with the actual water body distribution characteristics.

[0036] The extraction of dynamic shorelines includes extracting contour lines with a probability value of 0.5 from a soft mask of water bodies as dynamic shorelines. These contour lines serve as the probability boundary between water bodies and non-water bodies, and can adapt to the actual distribution of water bodies in the image without relying on external vector shoreline data. For example, when the lake water level changes slightly, this method can extract shorelines that match the actual water level, solving the problem of discrepancies between fixed vector shorelines and actual shorelines.

[0037] The geometric definition of key monitoring areas includes: calculating the Euclidean distance from each pixel in the grid to the dynamic shoreline; and simultaneously selecting a pixel set based on a dual-condition approach using a shoreline buffer distance threshold and a water intake neighborhood radius threshold. This constructs a key monitoring area that not only identifies high-risk zones along the shoreline but also focuses on the core area of ​​the water intake, achieving precise, scenario-based delineation of the monitoring area. For example, if a water intake is located near the shoreline, this method will define the overlapping area that simultaneously meets the criteria of being within 2000 meters of the shoreline and within 3000 meters of the water intake as the key monitoring area.

[0038] The steepness parameter of the S-shaped function of the water body soft mask is preferably set to 25, and the center threshold parameter is preferably set to 0.05. The combination of these two parameters can achieve the optimal mapping from the normalized differential water body index to the water body probability value, which is suitable for the water body spectral characteristics of most lakes.

[0039] The preferred small constant to prevent the denominator from being zero is 10 to the power of negative 6. This value is a very small constant commonly used in remote sensing image spectral calculations and can effectively avoid division by zero anomalies in the calculation of various water body indices.

[0040] The shoreline buffer distance threshold is preferentially set at 2000 meters, and the intake neighborhood radius threshold is preferentially set at 3000 meters. Both thresholds are set based on the actual scenario of cyanobacterial blooms accumulating along the shore in shallow eutrophic lakes, and can cover the high-risk monitoring range of the intake and shoreline.

[0041] Contour extraction includes: using the step grid method to extract contour lines with a probability value of 0.5 in the soft mask of the water body. This algorithm is a classic algorithm for contour extraction of raster data. It can accurately extract continuous shoreline contour lines, and the algorithm is simple to implement and adaptable to the raster data format of satellite remote sensing images.

[0042] Euclidean distance is calculated in meters by using the grid pixels of the resampled satellite remote sensing image as the unit. The calculation accuracy is consistent with the grid resolution, that is, a 10-meter resolution grid corresponds to a 10-meter calculation accuracy, ensuring the matching of distance calculation with spatial grid.

[0043] The determination of monitoring areas for multiple water intakes includes: delineating independent key monitoring areas for each water intake according to dual distance thresholds, with no mutual influence between the areas, and conducting subsequent algae analysis and risk assessment separately. This rule enables simultaneous and accurate monitoring of multiple water intakes.

[0044] Preferably, when the algal aggregation morphology analysis unit generates a continuous probability map of algal distribution based on the water reflectance data, and calculates the directional texture intensity index reflecting the degree of algal film banding and the lateral boundary change rate index reflecting the gradient of the band edges, it specifically performs the following steps: Extract the red, near-infrared, and short-wave infrared values ​​from the water reflectance data after eliminating the geometric effects of sunlight, and generate the continuous probability map of algal distribution based on the following formula: The structure tensor matrix is ​​calculated using the Sobel gradient operator and the Gaussian smoothing kernel, and the directional texture intensity index and the lateral boundary change rate index are generated based on the following formulas: in, , , These represent the reflectance values ​​for the red, near-infrared, and short-wave infrared bands, respectively. , , Indicates the center wavelength of the corresponding band. Indicates the phytoplankton index, This represents the probability kurtosis parameter of the algal membrane. This represents the threshold parameter at the center of the algal membrane. This represents the data of the soft mask for the water body. This represents a continuous probability map of the algal distribution. and Represents the eigenvalues ​​of the structure tensor matrix and , This represents the directional texture intensity index. Represents the components of the structure tensor matrix. Indicates the angle of the main direction of the strip. Represents the gradient vector. Represents a lateral unit vector. This represents the lateral boundary change rate index.

[0045] The steepness parameter of the S-shaped function of algal film probability is a parameter that controls the steepness of the curve change when the phytoplankton index is mapped to the algal film probability value. The preferred value is 20. The reason for this value is that it allows for a reasonable distinction between the probability mapping of algal film and non-algal film regions, avoiding misjudgment caused by over-sharpening, and also preventing the algal film boundary recognition from being affected by an excessively wide transition zone.

[0046] The central threshold parameter of the S-shaped function for algal film probability is the critical value corresponding to an algal film probability value of 0.5 when the phytoplankton index is mapped to an algal film probability value of 0.0. The preferred value is 0.0, as this value adapts to the spectral differences between cyanobacteria and water bodies, accurately distinguishing cyanobacterial algal films from ordinary water bodies.

[0047] The center wavelength of the red band is the core wavelength value of the red band in satellite remote sensing, representing the spectral position of the red band. It can be directly read from the metadata file of satellite remote sensing imagery, and the center wavelength of the red band has a fixed standard for different satellite sensors.

[0048] The center wavelength of the near-infrared band is the core wavelength value of the near-infrared band in satellite remote sensing, representing the spectral position of the near-infrared band. It can be directly read from the metadata file of satellite remote sensing imagery, adapting to the needs of cyanobacteria spectral identification.

[0049] The center wavelength of the shortwave infrared band is the core wavelength value of the shortwave infrared band in satellite remote sensing, representing the spectral position of the shortwave infrared band. It can be directly read from the metadata file of satellite remote sensing imagery and used to construct a three-band baseline for calculating the phytoplankton index.

[0050] The three-band baseline height difference method for calculating the phytoplankton index involves: first, constructing a linear baseline from red light to shortwave infrared based on the reflectance values ​​of the red, near-infrared, and shortwave infrared bands; then, subtracting the theoretical near-infrared value corresponding to this baseline from the near-infrared band reflectance to obtain the phytoplankton index. This method effectively highlights the spectral characteristics of cyanobacterial films while suppressing water background and other interfering factors. For example, with a red band reflectance of 0.05, a near-infrared reflectance of 0.3, and a shortwave infrared reflectance of 0.1, the theoretical near-infrared value calculated using the baseline is 0.1, resulting in a final phytoplankton index of 0.2, indicating the presence of an algal film.

[0051] The generation of a continuous probability map of algal distribution involves: using an S-shaped function with algal film probability steepness and a center threshold parameter to perform a continuous nonlinear mapping of the phytoplankton index from 0 to 1; then multiplying it pixel-by-pixel with the water body soft mask data to remove interference from non-water areas. This method preserves the gradient information of algal film thickness, avoiding the morphological distortion of algal films caused by traditional binary classification. For example, thin algal film areas are assigned probability values ​​of 0.3 to 0.7, while thick algal film areas are close to 1, which more closely reflects the actual algal film distribution.

[0052] The generation of the structural tensor matrix includes: first, calculating the horizontal and vertical gradient fields of the continuous probability map of algal distribution using a 3×3 Sobel gradient operator to obtain the gradient vector of each pixel; then, calculating the tensor product of the gradient vectors; and finally, performing a Gaussian smoothing kernel convolution operation on the tensor product to smooth out noise interference. This method can effectively extract the texture structure information of algal membranes. For example, the gradient field of striped algal membranes will show obvious directional consistency, and after tensor product and smoothing, a regular structural tensor matrix is ​​formed.

[0053] The calculation of the directional texture intensity index includes: solving for the maximum and minimum eigenvalues ​​of the structure tensor matrix, and using the ratio of the difference between the two to the sum of the two eigenvalues ​​to obtain the directional texture intensity. The closer this ratio is to 1, the higher the degree of striping of the algal membrane; the closer the ratio is to 0, the more uniform the distribution of the algal membrane. For example, for a striped algal membrane with a maximum eigenvalue of 0.8 and a minimum of 0.1, the directional texture intensity is approximately 0.78, accurately quantifying the degree of striping.

[0054] The calculation of the lateral boundary change rate index includes: calculating the principal direction angle of the strip based on the structure tensor matrix, constructing a lateral unit vector perpendicular to this angle, projecting the gradient field onto this vector, and obtaining the lateral boundary change rate of each pixel. This index can characterize the steepness of the strip edge; for example, the lateral boundary change rate of the strip edge pixels is 0.9, while that of the interior pixels is 0.2, clearly distinguishing the strip edge from the interior region.

[0055] The steepness parameter of the S-shaped function of algal film probability is preferably 20, and the center threshold parameter is preferably 0.0. The two parameters are used together to achieve the optimal mapping from phytoplankton index to algal film probability, which is suitable for the spectral characteristics of most lake cyanobacteria.

[0056] To prevent small constants with zero denominators from being prioritized, 10 to the power of -6 is preferred. This value is a commonly used stabilization constant in remote sensing image texture calculation, which can effectively avoid division by zero anomalies in steps such as feature value calculation.

[0057] The Gaussian smoothing kernel preferentially adopts a smoothing scale of 30 meters, corresponding to three pixels with a resolution of 10 meters. This scale can effectively smooth gradient noise while preserving the structural information of algal membrane strips and avoiding texture distortion caused by over-smoothing.

[0058] The Sobel gradient operator uses a 3×3 standard Sobel template, with the horizontal template being [-1,0,1;-2,0,2;-1,0,1] and the vertical template being [-1,-2,-1;0,0,0;1,2,1]. This template can accurately calculate the gradient direction and magnitude of pixels, making it suitable for texture extraction from raster images.

[0059] The construction of the lateral unit vector includes: when the angle of the main direction of the strip is φ, the horizontal coordinate of the lateral unit vector is the negative sine value of φ, and the vertical coordinate is the cosine value of φ. This construction method can ensure that the vector is strictly perpendicular to the main direction of the strip and accurately capture the lateral boundary features of the strip.

[0060] Gradient field projection includes: multiplying the horizontal gradient value of each pixel by the x-coordinate of the lateral unit vector, multiplying the vertical gradient value by the y-coordinate of the vector, adding the two together and taking the absolute value, and adding a small constant to obtain the lateral boundary change rate. This implementation can accurately quantify the distribution characteristics of the gradient in the lateral direction.

[0061] The solution of eigenvalues ​​of the structure tensor matrix includes: using the Jacobi iteration algorithm to solve for the eigenvalues ​​of the structure tensor matrix. This algorithm is a classic method for solving the eigenvalues ​​of symmetric matrices, with high computational accuracy, strong stability, and suitability for the symmetric properties of the structure tensor matrix.

[0062] Preferably, when the algal aggregation morphology analysis unit calculates the weighted distribution area based on the algal distribution continuous probability map and the directional texture intensity index, and then converts the weighted distribution area into the equivalent total strip length by combining the lateral boundary change rate index, it specifically performs the following steps: The weighted distribution area and initial effective strip width within the key monitoring area are calculated based on the following formula: The initial effective strip width is constrained by a soft enhancement function, and the equivalent total strip length is generated based on the following formula: in, This represents a continuous probability map of the algal distribution. This represents the directional texture intensity index. Represents the structural weight distribution data. This refers to the key monitoring area for the i-th water intake. Represents the physical area of ​​a single pixel. Represents the area of ​​the weighted distribution. This represents the lateral boundary change rate index. This represents the initial effective strip width. This indicates the preset minimum width threshold. This indicates the preset maximum width threshold. Represents the natural logarithm function. Represents an exponential function. Indicates the effective width of the stabilized strip. This represents the total length of the equivalent strip.

[0063] The physical area of ​​a single pixel is the actual geographic spatial area corresponding to that pixel in a unified grid, representing the actual coverage scale of the pixel. A preferred value is 100 square meters, with a grid resolution of 10 meters. 10 meters multiplied by 10 meters equals 100 square meters, perfectly matching the spatial scale and ensuring accurate area calculation.

[0064] The minimum effective strip width threshold is a lower limit critical value constraining the initial effective strip width, used to avoid an abnormal increase in the total equivalent strip length due to excessively small width. A value of 10 meters is preferred, as this value closely matches the actual minimum width characteristics of cyanobacteria strips.

[0065] The maximum effective strip width threshold is an upper limit critical value constraining the initial effective strip width, used to prevent the width from becoming too large and deviating from the essential strip structure. A value of 300 meters is preferred, as this value covers the actual width range of most cyanobacterial strips, ensuring the rationality of the strip morphology.

[0066] The unit length constant is a constant added when calculating the equivalent total strip length to avoid calculation anomalies where the denominator is zero. It is preferably 1 meter, as this value is extremely small and does not affect the calculated strip length.

[0067] The generation of structural weight distribution data includes: multiplying the continuous probability map of algal distribution with the directional texture intensity index pixel by pixel to achieve dual weighting of the probability of algal membrane existence combined with the degree of striping, highlighting the contribution of striped algal membranes and suppressing interference from non-striped regions.

[0068] The calculation of the weighted distribution area includes: accumulating the structural weight distribution data of the key monitoring area pixel by pixel, and then multiplying it by the physical area of ​​a single pixel. The resulting weighted distribution area is not the geometric area in the traditional sense, but an effective area that incorporates the characteristics of the strip structure.

[0069] The derivation method of the initial effective strip width is as follows: multiply the structural weight distribution data with the lateral boundary change rate index pixel by pixel to generate weighted gradient energy distribution data. After accumulation, the total weighted lateral gradient energy value is obtained. Then, the weighted distribution area is divided by this total value to realize the structural dimension conversion from two-dimensional area to one-dimensional width.

[0070] The method for constructing a stable effective strip width is to use a soft enhancement function (logarithmic-exponential combination function) to constrain the initial effective strip width for continuity, rather than hard threshold pruning. For example, if the initial width is 5 meters, which is lower than the minimum threshold of 10 meters, after the lower limit constraint, it is 10 meters plus the natural logarithm (1 plus the exponent (5 minus 10)), resulting in approximately 10.0067 meters, which meets the threshold requirement while maintaining continuity.

[0071] The calculation of the equivalent total strip length includes dividing the weighted distribution area by the sum of the stabilized effective strip width and the unit length constant. The resulting length quantifies the potential contact length between the cyanobacterial strip and the shoreline, providing core structural parameters for subsequent risk assessment. For example, with a weighted distribution area of ​​2,000,000 square meters and a stabilized width of 40 meters, the equivalent total strip length is approximately 2,000,000 divided by 41, or approximately 48,780 meters.

[0072] The minimum effective strip width threshold is preferably set to 10 meters, and the maximum effective strip width threshold is preferably set to 300 meters. These two thresholds, together with the soft enhancement function, can cover the strip width of all real-world scenarios.

[0073] The unit length constant is specifically set to 1 meter.

[0074] The soft enhancement function includes: first, calculating the lower bound constraint value as the minimum width threshold plus the natural logarithm (1 plus the exponent (initial width minus minimum threshold)); then calculating the upper bound constraint value as the maximum width threshold minus the natural logarithm (1 plus the exponent (maximum threshold minus lower bound constraint value)); and finally stabilizing the width as the upper bound constraint value.

[0075] The numerical accumulation only applies to all pixels within the key monitoring area, excluding non-monitored pixels outside the area.

[0076] The calculation of weighted gradient energy distribution data includes: first, taking the absolute value of the gradient field projection result, and then adding a small constant (preferably 10 to the power of -6) to prevent the denominator from being zero, thus avoiding numerical anomalies.

[0077] Preferably, when the shoreline exposure risk assessment unit maps the total length of the equivalent strip to the unit shoreline contact intensity for the key monitoring area and calculates the probability of sudden shoreline contact using a random shoreline contact model, it specifically performs the following steps: The physical length of the dynamic shoreline located within the key monitoring area is calculated to generate the regional shoreline length; The multi-source satellite remote sensing data are sorted by time, the time difference between the current data and the previous data is calculated, and an observation time interval is generated. The unit shoreline contact strength is generated based on the following formula: Using the aforementioned random touchdown model, the probability of sudden touchdown is generated based on the following formula: in, Indicates the current moment. This refers to the key monitoring area for the i-th water intake. This indicates the length of the shoreline in the area. This represents the total length of the equivalent strip. This indicates the observation time interval. This represents a preset time scale constant. This indicates the unit shoreline contact strength. This represents an exponential function with the natural constant as its base. This represents the probability of a sudden landslide.

[0078] The timescale constant is a reference constant that standardizes the observation time interval, used to unify the calculation scale of unit shoreline contact intensity under different observation intervals. A timescale of 2 hours is preferred because cyanobacterial blooms in shallow eutrophic lakes exhibit rapid migration characteristics on an hourly scale, and 2 hours aligns with the conventional observation intervals of satellite imagery, ensuring consistency in the temporal dimension of contact intensity calculations.

[0079] The unit length constant is a constant added when calculating the linear contact density to avoid the anomaly of a zero denominator caused by zero shoreline length. It is preferably 1 meter, a very small value that does not affect the calculated contact density.

[0080] The statistics on regional shoreline length include: only counting the physical length of dynamic shorelines within key monitoring areas, rather than the total shoreline length of the entire lake, to achieve scenario-based adaptation of shoreline length and make contact intensity calculations more consistent with the actual situation in key monitoring areas. For example, if the total length of the broken line segment of the dynamic shoreline within the key monitoring area is 5000 meters, this value is the regional shoreline length.

[0081] The generation of observation time intervals includes: first, sorting the multi-source satellite remote sensing data in ascending order of observation time; the time difference between the current time and the previous time is the observation time interval for the corresponding time, ensuring that the time interval strictly matches the image time sequence. For example, after sorting, the times are 8:00, 10:00, and 12:00, and the corresponding observation time intervals are 2 hours and 2 hours, respectively.

[0082] The mapping of contact intensity per unit shoreline includes: first, dividing the total equivalent strip length by the sum of the regional shoreline length and the unit length constant to obtain the linear contact density; then multiplying by the ratio of the observation time interval to the time scale constant to realize the transformation from strip length to time-weighted contact density, thus quantifying the potential contact pressure per unit shoreline per unit time.

[0083] The calculation of the probability of sudden shore contact includes: using the Poisson random shore contact model, with the natural constant as the base, the negative value of the unit shore contact intensity is used to calculate the result of the power operation, and the result is subtracted from 1 to obtain the probability of sudden shore contact. This model can accurately quantify the probability of low-probability, high-impact strip shore contact events.

[0084] The specific value of the time scale constant is preferably 2 hours. This value is set based on the regular revisit cycle of satellite imagery and the hourly migration characteristics of cyanobacterial blooms, which is suitable for the monitoring needs of most lakes.

[0085] The statistical analysis of the shoreline length of the region includes: using a polyline segment length accumulation algorithm, the contour lines of the dynamic shoreline are first discretized into several polyline segments, the Euclidean distance of each polyline segment is calculated and accumulated, the discretization accuracy of the polyline segment is 10 meters, consistent with the grid resolution, to ensure accurate length calculation.

[0086] Time alignment of multi-source satellite data includes: sorting based on the imaging time (UTC time) of satellite images. If there are time zone differences, the data must be converted to the same time zone (such as East 8th time zone) before sorting to avoid time zone deviations causing errors in time interval calculation.

[0087] The observation time interval for the first frame image is padded, including: when there is no preceding observation data for the first frame image, the median of all subsequent observation time intervals is taken as the observation time interval for the first frame, ensuring that all images have a corresponding time interval and avoiding empty values. For example, if the subsequent observation time intervals are 2 hours, 2 hours, and 3 hours, and the median is 2 hours, then the observation time interval for the first frame is 2 hours.

[0088] Preferably, when the shoreline exposure risk assessment unit generates a single observation risk value by combining the background water quality benchmark with a preset local accumulation amplification factor, it specifically performs the following steps: Extract the coastal band, blue band, and green band values ​​from the water reflectance data after eliminating the geometric effects of solar radiation, and generate the background water quality benchmark based on the following formula: The single observation risk value is generated based on the following formula: in, , , These represent the reflectance values ​​for the coastal band, blue band, and green band, respectively. This represents the function that takes the maximum value. This represents a commonly used logarithmic function. Indicates the background chlorophyll index. Indicates the pixel background density proxy value. This refers to the key monitoring area for the i-th water intake. This represents the data of the soft mask for the water body. This indicates the background water quality standard. This represents the probability of the sudden landslide. This indicates the local accumulation amplification factor. This represents the risk value of a single observation.

[0089] The local accumulation amplification factor is a parameter characterizing the degree of amplification of local concentration or toxins after cyanobacterial streaks reach the shore. A value of 1000 is preferred because in shallow, eutrophic lakes, when light winds push cyanobacterial streaks to the shore, the cell or toxin concentration can be amplified to 1000 times the background value, thus reflecting the amplification effect of the actual monitoring scenario.

[0090] Coastal bands are spectral wavelength ranges corresponding to coastal areas acquired by satellite remote sensing sensors, representing specific short-wavelength spectral information. They can be directly read from the metadata files of satellite remote sensing imagery, and different satellite sensors have fixed center wavelengths for their coastal bands.

[0091] The blue light band refers to the spectral wavelength range corresponding to the blue light region acquired by satellite remote sensing sensors, representing short-to-medium wavelength spectral information. It can be directly read from the metadata files of satellite remote sensing images, adapting to the needs of algae spectral feature identification.

[0092] The green band represents the spectral wavelength range corresponding to green areas acquired by satellite remote sensing sensors, characterizing mid-wavelength spectral information. It can be directly read from the metadata file of satellite remote sensing imagery and used in combination with other bands to calculate algae-related indices.

[0093] The calculation of the background chlorophyll index includes: extracting the coastal, blue, and green light band values ​​from the water reflectance data after eliminating the geometric influence of sunlight; first, taking the maximum value between the coastal and blue light band values; then dividing this maximum value by the green light band value; and finally taking the common logarithm of the result to obtain the background chlorophyll index. This method can effectively highlight the spectral characteristics of background algae in the water body and suppress interference from other non-algae components.

[0094] The generation of pixel background concentration proxy values ​​includes: calculating the result of a power operation with base 10 and 0.8 times the background chlorophyll index as the exponent, realizing non-linear quantization scaling of the background chlorophyll index, and converting the exponent value into a proxy value that can characterize the background algae concentration.

[0095] The calculation of the background water quality baseline includes: within key monitoring areas, using water body soft mask data as weights, calculating a weighted average of pixel background concentration proxy values. The weights are positively correlated with the water body probability, which can eliminate the interference of non-water body pixels in the calculation results. For example, in key monitoring areas, some pixels have a water body soft mask probability of 1, while others have 0.5. The weighted average will highlight the contribution of pixels with high water body probability, making the results more consistent with the actual background water quality.

[0096] The method for generating the risk value of a single observation is as follows: The probability of sudden shore contact, the background water quality benchmark, and the local accumulation amplification factor are multiplied together to construct a risk quantification logic of probability × benchmark × amplification. This logic considers both the likelihood of a shore contact event and the background algal intensity and the amplification effect after shore contact, comprehensively characterizing the risk level of a single observation. For example, if the probability of sudden shore contact is 0.8, the background water quality benchmark is 1.74, and the local accumulation amplification factor is 1000, the risk value of a single observation is 0.8 × 1.74 × 1000 = 1392.

[0097] The local accumulation amplification factor is preferably 1000. This value is determined based on a large amount of monitoring data on cyanobacteria reaching the shore in shallow eutrophic lakes. It is a typical magnitude of concentration amplification after reaching the shore and is applicable to risk assessment of most similar lakes.

[0098] The weighted average uses the arithmetic weighted average method. The weight of each pixel is the probability value of the water body soft mask data. During the calculation, the pixel background concentration proxy value is first multiplied by the corresponding weight, then all products are summed, and finally divided by the total weight to ensure that the calculation logic is clear and reproducible.

[0099] The binding of the center wavelengths of the bands includes: the center wavelengths of the coastal band, blue band, and green band must be bound to the selected satellite sensors. For example, in the Sentinel-2 satellite, the center wavelength of the coastal band is 443 nm, the blue band is 490 nm, and the green band is 560 nm; in the Landsat-8 satellite, the center wavelengths are 443 nm for the coastal band, 482 nm for the blue band, and 561 nm for the green band, to ensure consistency in the use of the bands.

[0100] The pixel background concentration proxy value is a dimensionless quantity used only for relative risk characterization and does not directly correspond to the actual concentration unit. Its value is positively correlated with the background algae content.

[0101] Preferably, when generating the cloud interference suppression coefficient based on cloud cover detection, the dynamic early warning generation unit specifically performs the following steps: Extract the blue, green, red, and cirrus band values ​​from the water reflectance data to eliminate the geometric effects of sunlight, and generate the cloud interference suppression coefficient based on the following formula: in, , , , These represent the reflectance values ​​for the blue, green, red, and cirrus light bands, respectively. Indicates spectral brightness index, Represents an exponential function. and These represent the first cloud detection steepness parameter and the first cloud detection center threshold parameter, respectively. and These represent the second cloud detection steepness parameter and the second cloud detection center threshold parameter, respectively. Indicates the probability of pixel cloud amount. This refers to the key monitoring area for the i-th water intake. This represents the data of the soft mask for the water body. This represents the probability of average cloud cover in the region. This represents the penalty intensity parameter. This represents the cloud interference suppression coefficient.

[0102] The kurtosis parameter of the S-curve function for cloud detection is a parameter that controls the steepness of the curve change when the spectral brightness index is mapped to the cloud cover probability. A value of 20 is preferred, as this value allows the probability mapping distinction between cloud and non-cloud areas to be adapted to the cloud detection requirements of satellite imagery, avoiding an excessively wide or narrow transition zone.

[0103] The first cloud detection S-shaped function center threshold parameter is the critical value corresponding to a cloud cover probability of 0.5 when the spectral brightness index is mapped to the cloud cover probability. A value of 0.08 is preferred, as this value closely matches the spectral brightness differences between clouds, water bodies, and land, and can accurately distinguish cloud-covered areas.

[0104] The second cloud detection S-curve steepness parameter is a parameter that controls the steepness of the curve change when the cirrus cloud band values ​​are mapped to cloud cover probability. A value of 40 is preferred, as this results in weaker cirrus cloud spectral characteristics, and a higher steepness enhances the distinction between cirrus and non-cirrus clouds.

[0105] The center threshold parameter of the S-shaped function for the second cloud detection is the critical value corresponding to a cloud cover probability of 0.5 when the cirrus cloud band value is mapped to a cloud cover probability of 0.5. A value of 0.01 is preferred, as this value is suitable for the low reflectivity characteristics of the cirrus cloud band and can accurately capture cirrus cloud signals.

[0106] The cloud cover penalty intensity parameter is used to adjust the sensitivity of the cloud interference suppression coefficient to cloud cover. A value of 3 is preferred, as this balances the suppression effect of cloud cover interference with the preservation of normal data, avoiding excessive or insufficient penalty.

[0107] Cirrus bands are the spectral wavelength ranges corresponding to cirrus clouds acquired by satellite remote sensing sensors, representing specific long-wavelength spectral information. They can be directly read from the metadata files of satellite remote sensing imagery, adapting to the spectral requirements of cirrus cloud detection.

[0108] The calculation of the spectral brightness index includes: extracting the reflectance values ​​of the blue, green, and red light bands, and calculating the arithmetic mean of the three to obtain the spectral brightness index. This index can comprehensively reflect the brightness of a region and highlight the high brightness characteristics of cloud-covered areas. For example, with a blue band value of 0.15, green band value of 0.18, and red band value of 0.16, the spectral brightness index is 0.163, clearly characterizing the high brightness of clouds.

[0109] The dual S-shaped function cloud cover detection method includes: using two independent S-shaped functions to perform continuous nonlinear mapping on the spectral brightness index and the cirrus band values ​​respectively, and then multiplying the two mapping results to generate the pixel cloud cover probability, thereby realizing a two-dimensional fusion detection combining visible light high brightness and cirrus cloud features.

[0110] The calculation of the regional average cloud cover probability includes: within key monitoring areas, using water body soft mask data as weights, a weighted average of the pixel cloud cover probabilities is performed. The weights are positively correlated with the water body probability, which can eliminate cloud cover interference from non-water body pixels within the area. For example, in key areas, the cloud cover probability for water body pixels is 0.7, and for non-water body pixels it is 0.9. The weighted average weakens the influence of non-water body pixels, resulting in a result that more closely reflects the actual cloud cover in water body areas.

[0111] The generation of the cloud interference suppression coefficient includes: calculating the result of a power operation with the negative penalty intensity parameter as the base and the regional average cloud cover probability as the exponent, using the natural constant as the base. This coefficient is negatively correlated with the cloud cover probability, and the more cloud cover there is, the smaller the coefficient becomes, which can continuously suppress the interference of clouds on risk calculation.

[0112] The parameters for the dual S-shaped function are: first cloud detection steepness 20, center threshold 0.08; second cloud detection steepness 40, center threshold 0.01. The combination of these two sets of parameters enables the collaborative detection of visible light clouds and cirrus clouds, and is suitable for cloud detection scenarios in most satellite imagery.

[0113] The cloud cover penalty intensity parameter is preferably set to 3. This value is a commonly used parameter for suppressing interference in remote sensing data, which can retain effective monitoring information while suppressing cloud cover interference.

[0114] The center wavelength binding of the cirrus band includes: the cirrus band must be bound to the selected satellite sensor. For example, the center wavelength of the cirrus band of Sentinel-2 satellite is 1375 nanometers, and that of Landsat-8 is 1373 nanometers, to ensure consistency in the use of the band.

[0115] The use of water body soft mask weights includes: when the probability of water body soft masking is less than 0.1, the weight of the pixel is set to 0, so as to completely eliminate the interference of non-water body areas on cloud amount calculation.

[0116] Preferably, when the dynamic early warning generation unit applies the cloud interference suppression coefficient and the time decay factor to the single observation risk value to generate a cumulative risk score, and outputs an alarm signal when the cumulative risk score exceeds a preset alarm threshold, it specifically performs the following steps: The cumulative risk score is generated based on the following formula: The alarm signal is output based on the following formula: in, The current system time is represented by t, where t represents the observation time of the multi-source satellite remote sensing data. This represents the preset time decay parameter. This represents the time decay factor. This represents the set of times within a preset time period. This represents the risk value of a single observation. This represents the cloud interference suppression coefficient. This represents the cumulative risk score for the i-th water intake. This indicates the preset alarm threshold. This indicates the alarm signal, when An alarm is output when the value is 1.

[0117] The time decay parameter is a parameter that adjusts the sensitivity of the time decay factor to the difference between the observation time and the current system time. A value of 0.35 per hour is preferred, as this balances the weights of recent and long-term observation data, allowing the risk score to better reflect the hourly dynamic changes of cyanobacterial blooms.

[0118] The current system time is the specific point in time when the monitoring system is running, representing the time base for risk calculation. It can be directly obtained through the monitoring system's built-in clock module, containing precise information including year, month, day, hour, and minute.

[0119] The cumulative risk alarm threshold is the critical value used to determine whether to output an alarm signal, representing the maximum acceptable risk level. A value of 1200 is preferred, as this effectively distinguishes between low-risk and high-risk events, meeting the safety monitoring needs of water intakes.

[0120] The preset time period is the time range of satellite imagery used in the cumulative risk score calculation, representing the time span of risk aggregation. A preferred timeframe is 48 hours, which covers the hourly migration cycle of cyanobacterial blooms while avoiding excessively long periods that could lead to delayed risk information.

[0121] The calculation of the time decay factor includes: multiplying a negative time decay parameter by the time difference between the observation time and the current system time, using the natural constant as the base, and then performing a power operation to obtain the time decay factor. This factor decreases as the time difference increases, implementing a time-weighted logic where recent data has a higher risk weight and distant data has a lower risk weight. For example, if the difference between the observation time and the current system time is 2 hours, the time decay parameter is 0.35 per hour, and the decay factor is e to the power of -0.7, approximately 0.497, reflecting the high weight given to recent data.

[0122] The construction of weighted risk units involves multiplying the single observation risk value, the cloud interference suppression coefficient, and the time decay factor to generate a weighted risk unit that integrates spatial risk, cloud interference, and time decay, thus achieving the organic integration of multi-dimensional risk information. For example, if the single observation risk value is 1392, the cloud interference suppression coefficient is 0.223, and the time decay factor is 0.497, the weighted risk unit is approximately 1392 × 0.223 × 0.497 ≈ 155.

[0123] The calculation of the cumulative risk score includes: summing all weighted risk units within a preset time period to obtain the weighted risk sum, summing all time decay factors to obtain the weighted sum, and dividing the two to obtain the cumulative risk score, thus achieving spatiotemporal aggregation of risk values ​​rather than a simple arithmetic average. For example, if the weighted risk sum is 775 and the weighted sum is 2.01, the cumulative risk score is approximately 385.6.

[0124] The time decay parameter is preferably set to 0.35 per hour, in units of one hour, to match the regular observation interval of satellite imagery and ensure the rationality of the time weighting logic.

[0125] The cumulative risk alarm threshold is preferentially set at 1200. This value is determined based on simulation data of cyanobacteria contact risk in a large number of shallow eutrophic lakes, which can effectively balance the false alarm rate and the missed alarm rate.

[0126] The preset time period is initially set to 48 hours. If there are no valid images within 48 hours, it can be extended to 7 days to ensure that there is enough time series data to participate in risk aggregation.

[0127] The cumulative risk assessment of multiple water intakes includes: calculating the cumulative risk score for each water intake independently, taking the maximum value among all scores and comparing it with the alarm threshold; if the score of any water intake exceeds the threshold, an alarm signal is output for the entire system to ensure that no high-risk areas are missed.

[0128] The alarm signal output is a standard digital signal, with 1 indicating an alarm and 0 indicating no alarm. It can be adapted to the existing interface of the water conservancy monitoring system and supports multiple linkage methods such as software prompts, hardware triggers, or SMS notifications.

[0129] like Figure 2 As shown, Figure 2 This is a schematic diagram illustrating the specific implementation scenario of a dynamic monitoring system for lake ecological meteorological disaster risks based on satellite remote sensing. It shows that satellite remote sensing is used as the data source. After being received and processed on the ground, the system identifies the dynamic shoreline of the lake, delineates the shoreline buffer zone, accurately locates the key monitoring area including the water intake, and focuses on the algal bands that need to be tracked. It fully outlines the process from data acquisition to the focus of monitoring objects, and intuitively reflects the targeted monitoring logic for the risk of cyanobacterial blooms touching the shoreline of lake shorelines and water intakes.

[0130] It should be noted that the interval and threshold sizes are set for ease of comparison. The size of the threshold depends on the amount of sample data and the base number set by those skilled in the art for each set of sample data, as long as it does not affect the proportional relationship between the parameter and the quantized value. Furthermore, the above formulas are all dimensionless calculations, and the formulas are derived from software simulations using a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0131] The embodiments of this example have been described above. However, this example is not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms based on the guidance of this example, and all of them are within the protection scope of this example.

Claims

1. A satellite remote sensing-based dynamic monitoring system for lake ecological meteorological disaster risks, characterized in that, include: The multi-source image preprocessing unit is used to receive multi-source satellite remote sensing data and map it to a unified grid, generate water reflectance data that eliminates the geometric influence of sunlight, and delineate key monitoring areas including water intakes based on shoreline distance. The algal aggregation morphology analysis unit is used to generate a continuous probability map of algal distribution based on the water reflectance data, calculate the directional texture intensity index reflecting the degree of algal film striping and the lateral boundary change rate index reflecting the gradient of the strip edge, and calculate the weighted distribution area based on the continuous probability map of algal distribution and the directional texture intensity index, and then combine the lateral boundary change rate index to convert the weighted distribution area into the equivalent total strip length. The shoreline exposure risk assessment unit is used to map the total length of the equivalent strip to the unit shoreline contact intensity for the key monitoring area, calculate the probability of sudden shoreline contact using a random shoreline contact model, and generate a single observation risk value by combining the background water quality benchmark and the preset local accumulation amplification factor. The dynamic early warning generation unit is used to generate a cloud interference suppression coefficient based on cloud cover detection, apply the cloud interference suppression coefficient and time decay factor to the single observation risk value to generate a cumulative risk score, and output an alarm signal when the cumulative risk score exceeds a preset alarm threshold.

2. The satellite remote sensing based lake eco-meteorological disaster risk dynamic monitoring system according to claim 1, characterized in that, Receive multi-source satellite remote sensing data and map it to a unified grid to generate water reflectance data that eliminates the geometric effects of solar illumination, including: Construct a set of multi-temporal surface reflectance products that meet the lake area coverage threshold within a preset time period before the current moment, and use it as the multi-source satellite remote sensing data; The multi-source satellite remote sensing data is projected onto a universal transverse Mercator projection coordinate system based on the lake's centroid, and then resampled to the unified grid with a resolution of ten meters using bilinear interpolation. Read the solar zenith angle parameter from the multi-source satellite remote sensing data; The initial remote sensing reflectance is obtained by dividing the surface reflectance value on the unified grid by pi, and the initial remote sensing reflectance is then divided by the cosine of the solar zenith angle to generate the water reflectance data that eliminates the geometric influence of sunlight.

3. The dynamic monitoring system for lake ecological meteorological disaster risk based on satellite remote sensing according to claim 2, characterized in that, Based on shoreline distance, key monitoring areas including water intakes were delineated, including: Extract the green band and near-infrared band values ​​from the water reflectance data after eliminating the geometric influence of sunlight, calculate the ratio of the difference between the green band and near-infrared band values ​​to the sum of the green band and near-infrared band values, and generate the normalized difference water index. A continuous nonlinear mapping of the normalized differential water index is performed using an S-shaped function containing a steepness parameter and a center threshold parameter to generate water body soft mask data. The contour lines with a probability value of 0.5 in the water body soft mask data are extracted as dynamic shorelines, and the Euclidean distance from each pixel in the unified grid to the dynamic shoreline is calculated. The key monitoring area is constructed by selecting a set of pixels whose Euclidean distance is less than a preset shoreline buffer distance threshold and whose distance from the preset water intake coordinate point is less than a preset water intake neighborhood radius threshold.

4. The dynamic monitoring system for lake ecological meteorological disaster risk based on satellite remote sensing according to claim 3, characterized in that, Based on the water reflectance data, a continuous probability map of algal distribution is generated. An directional texture intensity index reflecting the degree of algal film banding and a lateral boundary change rate index reflecting the gradient of the band edges are calculated, including: The red light band, near-infrared band, and short-wave infrared band values ​​were extracted from the water reflectance data after eliminating the geometric influence of sunlight, so as to calculate the phytoplankton index using the three-band baseline height difference method. The sigmoid index is continuously nonlinearly mapped using an S-shaped function that includes an algal film probability steepness parameter and an algal film center threshold parameter. The continuous nonlinear mapping result is then multiplied with the water body soft mask data to generate the continuous probability map of algal distribution. The gradient fields of the continuous probability map of algal distribution in the horizontal and vertical directions are calculated using the Sobel gradient operator, and the tensor product of the gradient fields is convolved using a Gaussian smoothing kernel to generate a structure tensor matrix. Calculate the maximum and minimum eigenvalues ​​of the structure tensor matrix, and generate the directional texture intensity index using the ratio of the difference between the maximum and minimum eigenvalues ​​to the sum of the maximum and minimum eigenvalues. The principal direction angle of the strip is calculated based on the structural tensor matrix, a lateral unit vector perpendicular to the principal direction angle of the strip is constructed, and the gradient field is projected onto the lateral unit vector to generate the lateral boundary change rate index.

5. The dynamic monitoring system for lake ecological meteorological disaster risk based on satellite remote sensing according to claim 4, characterized in that, The weighted distribution area is calculated based on the continuous probability map of algal distribution and the directional texture intensity index. Then, the weighted distribution area is converted into the equivalent total strip length by combining the lateral boundary change rate index, including: The continuous probability map of algal distribution is multiplied pixel by pixel with the directional texture intensity index to generate structural weight distribution data; The weighted distribution area is generated by numerically summing the structural weight distribution data within the key monitoring area and multiplying it by the physical area of ​​a single pixel. The structural weight distribution data is multiplied pixel by pixel with the lateral boundary change rate index to generate weighted gradient energy distribution data. The weighted gradient energy distribution data in the key monitoring area is numerically accumulated and multiplied by the physical area of ​​a single pixel to generate the total weighted lateral gradient energy value. Divide the weighted distribution area by the total weighted lateral gradient energy to generate the initial effective strip width; The initial effective strip width is constrained between a preset minimum width threshold and a preset maximum width threshold to generate a stable effective strip width; The equivalent total strip length is generated by dividing the weighted distribution area by the sum of the stabilized effective strip width and the unit length constant.

6. The dynamic monitoring system for lake ecological meteorological disaster risk based on satellite remote sensing according to claim 5, characterized in that, The equivalent total length of the strip is mapped to the unit shoreline contact intensity for the key monitoring area, and the probability of sudden shore contact is calculated using a random shore contact model, including: The physical length of the dynamic shoreline located within the key monitoring area is calculated to generate the regional shoreline length; The multi-source satellite remote sensing data are sorted by time, the time difference between the current data and the previous data is calculated, and an observation time interval is generated. The linear contact density is obtained by dividing the total length of the equivalent strip by the sum of the shoreline length of the region and the unit length constant. The linear contact density is then multiplied by the ratio of the observation time interval to the preset time scale constant to generate the unit shoreline contact intensity. Calculate the result of a power operation with the natural constant as the base and the negative value of the unit shoreline contact intensity as the exponent, and take the value of subtracting the result of the power operation as the probability of sudden shore contact.

7. The dynamic monitoring system for lake ecological meteorological disaster risk based on satellite remote sensing according to claim 6, characterized in that, The risk value for a single observation is generated by combining the probability of sudden shore contact with background water quality standards and a preset local accumulation amplification factor, including: Extract the coastal band, blue light band, and green light band values ​​from the water reflectance data after eliminating the geometric influence of sunlight. Calculate the maximum value between the coastal band value and the blue light band value, and use the common logarithm of the ratio of the maximum value to the green light band value as the background chlorophyll index. Calculate the result of the power operation with base 10 and exponent 0.8 times the background chlorophyll index, and generate a pixel background density proxy value; Within the key monitoring area, the background concentration proxy value of the pixel is calculated by weighted averaging using the water body soft mask data to generate the background water quality benchmark. The product of the sudden shore contact probability, the background water quality benchmark, and the local accumulation amplification factor is used as the single observation risk value.

8. The dynamic monitoring system for lake ecological meteorological disaster risk based on satellite remote sensing according to claim 7, characterized in that, Cloud interference suppression coefficients are generated based on cloud cover detection, including: Extract the blue light band values, green light band values, red light band values, and cirrus band values ​​from the water reflectance data after eliminating the geometric influence of sunlight. Calculate the arithmetic mean of the blue light band values, the green light band values, and the red light band values ​​to generate a spectral brightness index. The spectral brightness index is continuously nonlinearly mapped using an S-shaped function that includes the first cloud detection steepness parameter and the first cloud detection center threshold parameter to obtain the first mapping result; The cirrus band values ​​are continuously nonlinearly mapped using an S-shaped function that includes the second cloud detection steepness parameter and the second cloud detection center threshold parameter to obtain the second mapping result; Multiply the first mapping result and the second mapping result to generate the pixel cloud amount probability; Within the key monitoring area, the cloud cover probability of the pixels is calculated by weighted averaging using the water body soft mask data to generate the regional average cloud cover probability. The cloud interference suppression coefficient is generated by calculating the result of a power operation with the natural constant as the base and the negative penalty intensity parameter multiplied by the average cloud cover probability of the region as the exponent.

9. The dynamic monitoring system for lake ecological meteorological disaster risk based on satellite remote sensing according to claim 8, characterized in that, The cloud interference suppression coefficient and time decay factor are applied to the single observation risk value to generate a cumulative risk score, and an alarm signal is output when the cumulative risk score exceeds a preset alarm threshold, including: Obtain the current system time and calculate the time difference between the observation time of the multi-source satellite remote sensing data and the current system time; The time decay factor is generated by calculating the result of a power operation with the natural constant as the base and the negative preset time decay parameter as the exponent, multiplied by the time difference value. The product of the single observation risk value, the cloud interference suppression coefficient, and the time decay factor is used as a weighted risk unit, and all the weighted risk units within a preset time period are accumulated to generate a weighted risk sum. All time decay factors within the preset time period are summed to generate a weighted sum; The cumulative risk score is generated by dividing the weighted sum of risks by the sum of weights. Determine whether the cumulative risk score is greater than a preset alarm threshold. If so, output the alarm signal.

Images