Blue-green algae bloom grading early warning system combining vegetation index and water quality parameters

Abstract

The present application relates to the technical field of cyanobacterial bloom grading, and particularly relates to a cyanobacterial bloom grading and early warning system combining vegetation index and water quality parameters. The present application first deduces the drift trajectory and cumulative turbulent energy of a drift water mass based on wind conditions, thereby calculating a surface algae retention index; a flow direction constraint drift model of the drift water mass is constructed, and the model is solved to evaluate the pairing relationship between the drift water mass and neighboring pixels, further combining the surface algae retention index of the drift water mass and the monitored water quality, decoupling the algae spectral response coefficient and the water turbidity background parameter, and then determining the grading and early warning threshold to perform cyanobacterial bloom grading and early warning. The present application reconstructs the wind-driven flow field and introduces turbulent mixing constraints to solve the inconsistency problem of point monitoring and satellite observation in the time and space scales and physical space, and decouples the algae spectral response coefficient and the water turbidity background parameter to adaptively adjust the grading and early warning threshold, thereby improving the accuracy of cyanobacterial bloom grading and early warning.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of cyanobacterial bloom grading, and particularly relates to a cyanobacterial bloom grading and early warning system combining vegetation index and water quality parameters. BACKGROUND

[0002] In the monitoring technology of cyanobacterial bloom in large shallow lakes, fixed-point water quality monitoring and satellite remote sensing monitoring are two main data acquisition methods; ground stations can provide water quality parameters with high temporal resolution, and satellite remote sensing can provide instantaneous planar spectral images that can cover the entire lake area. The existing cyanobacterial bloom grading and early warning usually directly establishes a regression model based on the pixel values in the planar spectral image at the satellite transit time and the water quality parameters sampled by the ground station at that moment, and grades according to a fixed statistical threshold.

[0003] However, at the satellite imaging time, the water mass collected by the ground station will drift with the flow, resulting in a mismatch between the point water quality parameters and the planar spectral image in geographical space; the strong wind and wave in the shallow lake will cause vertical turbulence, which will also cause the blue-green algae to be rolled into the water, resulting in that the surface optical signal observed by the satellite is lower than the total biomass of the water column; the horizontal displacement and vertical distribution difference caused by fluid motion will make the point water quality data and the planar spectral image lack physical consistency; at the same time, the content of suspended sediment and other background substances in the water body fluctuates in real time with the hydrodynamic force, which will change the optical response baseline between the chlorophyll concentration and the vegetation index. The fixed threshold cannot adapt to this dynamic environment background, and the accuracy of the cyanobacterial bloom grading and early warning is low. SUMMARY

[0004] In order to solve the technical problem of low accuracy of cyanobacterial bloom grading and early warning, the purpose of the present application is to provide a cyanobacterial bloom grading and early warning system combining vegetation index and water quality parameters, and the technical solution adopted is as follows:

[0005] The cyanobacterial bloom grading and early warning system combining vegetation index and water quality parameters comprises:

[0006] The data acquisition module is used to acquire the vegetation index of each pixel in the satellite spectral image of the lake surface, and determine a time window based on the satellite imaging time; and acquire the water quality data and wind condition data of each acquisition time within the time window at the ground station;

[0007] The drift trajectory deduction module is used to determine the drift trajectory coordinates of the drift water mass at the ground station at each acquisition time at the satellite imaging time based on the wind condition data, and determine the cumulative turbulence energy of each drift water mass at the satellite imaging time;

[0008] The physical constraint inversion module is configured to obtain a surface algae retention index of each drift water mass according to the accumulated turbulent energy, and construct a flow direction constraint drift model of the drift water mass, solve the model to obtain a pairing coefficient between each drift water mass and each neighborhood pixel at a coordinate of the drift trajectory in the satellite spectral image; and obtain an algae spectral response coefficient and a water turbidity background parameter according to the surface algae retention index of each drift water mass, corresponding water quality data, a vegetation index of each neighborhood pixel, and a pairing weight between each drift water mass and each neighborhood pixel.

[0009] The hierarchical early warning module is configured to determine hierarchical early warning thresholds according to the algae spectral response coefficient and the water turbidity background parameter, and perform hierarchical early warning of water bloom based on a vegetation index of each pixel.

[0010] Further, the method for obtaining the time window comprises:

[0011] Based on a lake flow rate and a coverage radius of a ground station, a pushback tracing duration is obtained; and a time period of the pushback tracing duration is obtained on both sides of a time sequence centering on a satellite imaging time, and a time range corresponding to the two time periods is taken as the time window.

[0012] Further, the method for obtaining the drift trajectory coordinate comprises:

[0013] At each acquisition time within the time window, a lake surface flow rate vector is determined according to the wind condition data, a preset wind flow conversion coefficient, and a preset flow direction deflection matrix; a drift displacement is determined according to the lake surface flow rate vector and a drift duration between adjacent acquisition times; and a drift trajectory coordinate at a next acquisition time is determined based on a superposition result of a position coordinate of a ground station and the drift displacement at a first acquisition time.

[0014] For each drift water mass, within the time window, a drift trajectory coordinate at a next acquisition time is determined based on a superposition result of a position coordinate of a ground station and the drift displacement at a first acquisition time, and the drift trajectory coordinates at each acquisition time are obtained by sequentially superimposing.

[0015] Further, the method for obtaining the accumulated turbulent energy comprises:

[0016] At each acquisition time within the time window, a wind wave stirring power coefficient is determined according to the wind condition data, a stirring duration between adjacent acquisition times is taken, and the turbulent energy is obtained by fusing the wind wave stirring power coefficient and the stirring duration.

[0017] For each drift water mass, within the time window, the accumulated turbulent energy of the drift water mass at the satellite imaging time is determined based on an initial turbulent energy and a cumulative result of the turbulent energy at all acquisition times before the satellite imaging time.

[0018] Further, the method for obtaining the surface algae retention index comprises:

[0019] Obtaining a relative turbulent mixing intensity coefficient of each drift water mass by comparing the cumulative turbulent energy of different drift water masses; and taking a negative correlation normalization result of the relative turbulent mixing intensity coefficient as a surface algae retention index of each drift water mass.

[0020] Further, the obtaining method of the flow-constrained drift model comprises:

[0021] Taking a pair coefficient between each drift water mass and each neighboring pixel as an optimization variable, determining a spatial matching cost term in combination with a distribution relationship between each drift water mass and each neighboring pixel at a corresponding drift trajectory coordinate, determining a ground distribution constraint term in combination with the water quality data of each drift water mass and the surface algae retention index, and determining a satellite observation constraint term in combination with a position distribution of the neighboring pixels; and fusing the spatial matching cost term, the ground distribution constraint term, and the satellite observation constraint term to construct the flow-constrained drift model.

[0022] Further, the obtaining method of the spatial matching cost term comprises:

[0023] For each drift water mass, a displacement vector of each neighboring pixel is pointed to by the drift trajectory coordinate of the drift water mass, and the displacement vector is decomposed into an along-flow drift component and a cross-flow drift component;

[0024] The along-flow drift component and the cross-flow drift component are weighted and summed to obtain a drift cost parameter between the drift water mass and each neighboring pixel; and a weighting weight of the cross-flow drift component is greater than a weighting weight of the along-flow drift component.

[0025] The optimization variable of the flow-constrained drift model is weighted by using the drift cost parameter to obtain the spatial matching cost term.

[0026] Further, the obtaining method of the algae spectral response coefficient and the water turbidity background parameter comprises:

[0027] The water quality data at least comprises a chlorophyll concentration; for each drift water mass, the chlorophyll concentration is weighted by using the surface algae retention index, and a weighted result is taken as a surface chlorophyll concentration;

[0028] A linear regression model of a vegetation index of a pixel in a satellite spectral image is constructed; wherein a variable in the linear regression model is the surface chlorophyll concentration, a regression coefficient is an algae spectral response coefficient, and a bias is a water turbidity background parameter.

[0029] The regression coefficient and bias in a linear regression model are solved based on a weighted least square method, taking the matching coefficient between each drift water mass and each neighborhood pixel as a weight, to obtain an algae spectral response coefficient and a water turbidity background parameter.

[0030] Further, the method for obtaining the hierarchical early warning threshold comprises:

[0031] A preset number of risk levels are determined, the algae spectral response coefficient is weighted by using a preset early warning concentration of each risk level, and a fusion result of the weighted result and the water turbidity background parameter is taken as a hierarchical early warning threshold of the corresponding risk level.

[0032] Further, the water bloom hierarchical early warning based on the vegetation index of each pixel comprises:

[0033] In the satellite spectral image, for each pixel, a warning attribute of the pixel is marked based on a spatial distance between the pixel and the nearest drift trajectory coordinate, wherein the warning attribute of the spatial distance greater than a preset threshold is determined warning, and the warning attribute of the spatial distance less than or equal to the preset threshold is suspected warning;

[0034] The hierarchical early warning interval is determined based on the hierarchical early warning threshold of the adjacent risk level, and for each pixel, when the vegetation index is in the corresponding hierarchical early warning interval, the pixel is marked as the corresponding risk level;

[0035] The water bloom hierarchical early warning is performed based on the warning attribute and the risk level of the pixel.

[0036] The present application has the following beneficial effects:

[0037] The application first acquires the vegetation index of each pixel in the satellite spectral image of the lake surface, and determines a time window based on the satellite imaging time, to provide a time limit for subsequent backtracking; then acquires the water quality data and wind condition data of each ground station at each collection time within the time window, to provide a data basis for subsequent graded early warning analysis; further, based on the wind condition data, the drift trajectory coordinates of each drift water mass at the satellite imaging time are backtracked, and the cumulative turbulent energy of each drift water mass at the satellite imaging time is determined, and then the surface algae retention index of each drift water mass is acquired; then a flow direction constraint drift model of the drift water mass is constructed, and the pairing coefficient between each drift water mass and each adjacent pixel is solved, and further combined with the surface algae retention index of each drift water mass and the corresponding water quality data, the algae spectral response coefficient and the water turbidity background parameter in the satellite spectral signal are decoupled, so as to determine the graded early warning threshold, and based on the vegetation index of each pixel, the algal bloom graded early warning is carried out. The application reconstructs the wind-driven flow field and introduces the turbulent mixing constraint, solves the inconsistency problem of the point monitoring and the satellite observation in the time and space scales and the physical space, decouples the algae spectral response coefficient and the water turbidity background parameter, and adaptively adjusts the graded early warning threshold, so as to improve the accuracy of the graded early warning of the blue-green algal bloom. BRIEF DESCRIPTION OF DRAWINGS

[0038] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, and the advantages thereof, below, a brief introduction will be given to the drawings needed to be used in the embodiments or the prior art description. Obviously, the drawings in the following description only show some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained without creative labor based on these drawings.

[0039] Figure 1 A module diagram of a blue-green algal bloom graded early warning system combining vegetation index and water quality parameters provided by an embodiment of the present application. DETAILED DESCRIPTION

[0040] In order to further illustrate the technical means and effects adopted by the present application to achieve the predetermined application purpose, below, the specific embodiments, structure, features and effects of a blue-green algal bloom graded early warning system combining vegetation index and water quality parameters according to the present application are described in detail in combination with the drawings and preferred embodiments. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. In addition, the specific features, structures or characteristics in one or more embodiments can be combined in any suitable form.

[0041] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the present application belongs.

[0042] Specifically, the application provides a specific scheme of a blue-green algae bloom grading early warning system combining a vegetation index and water quality parameters.

[0043] Please refer to Figure 1 Fig. 1 shows a module diagram of a blue-green algae bloom grading early warning system combining a vegetation index and water quality parameters according to an embodiment of the application, which comprises a data acquisition module 101, a drift trajectory deduction module 102, a physical constraint inversion module 103 and a grading early warning module 104.

[0044] The application establishes an inversion model based on fluid dynamics constraints aiming at the horizontal drift of monitoring objects of ground stations caused by wind-driven current in shallow lakes and the algae vertical mixing phenomenon caused by strong wind waves, and then solves the model to decouple the algae spectral response coefficient representing the color developing ability of algae and the water turbidity background value representing the interference of non-biological suspended solids, and further establishes grading standards for blue-green algae bloom grading early warning.

[0045] The data acquisition module 101 is used to acquire the vegetation index of each pixel in the satellite spectral image of the lake surface, and determine a time window based on the imaging time of the satellite; and acquire the water quality data and wind condition data of the ground station at each acquisition time within the time window.

[0046] In an embodiment of the application, the satellite spectral image of the lake surface is acquired first, and the vegetation index (NDVI) of each pixel in the satellite spectral image is acquired. The vegetation index is a numerical index for quantifying the ground vegetation coverage, which can provide a data basis for subsequent grading early warning, and is a known technology and will not be described here. At the same time, the water quality data of the lake and the wind condition data of the lake surface are collected by the ground station deployed on the ground. The water quality data can help to evaluate the water bloom information subsequently, and the wind condition data can help to evaluate the drift influence of the water mass and the underwater entrainment of blue algae by wind.

[0047] The water quality data at least includes the chlorophyll concentration in the collected water sample, and the wind condition data includes the wind speed and the wind direction, which can also be directly represented by a wind speed vector. The acquisition frequency of the water quality data and the wind condition data is once every 10-30 minutes, and in this embodiment, it is once every 10 minutes. The implementer can also adjust it, but needs to ensure that the water quality data and the wind condition data are synchronized in time sequence.

[0048] It should be noted that the acquisition of the satellite spectral image, the vegetation index, the water quality data and the wind condition data is a prior art and will not be described here. In the grading early warning process, there can be multiple ground stations, and this embodiment takes one ground station as an example for analysis and description.

[0049] Then the satellite imaging time of the satellite spectral image is taken as the system anchor point or the hierarchical early warning time; since the blue-green algae in the lake water body can drift with the flow, a certain pixel in the satellite spectral image can correspond to a water mass that drifts from a certain position upstream at a historical time, and therefore a time window is further determined based on the satellite imaging time, so as to backtrack or deduce the drifting condition of the water mass.

[0050] It should be noted that the hierarchical early warning method at each satellite imaging time is consistent, and only any satellite imaging time is taken as an example for analysis and description.

[0051] Preferably, in an embodiment of the present application, considering that a too long time window not only increases the deduction burden but also can introduce errors, for example, theoretically the long-term drifting trajectory of the water mass can be deduced, but the water mass can have undergone biochemical reaction or metamorphosis in the long-term drifting process, thereby causing a large error in the deduced trajectory; and considering that the lake flow rate will affect the drifting speed of the water mass to a certain extent, and the monitoring coverage radius of the ground station can be regarded as a representative or reliable deduction boundary, and then the deduction or backtracking length can be determined; therefore, the method for obtaining the time window comprises:

[0052] Based on the lake flow rate and the coverage radius of the ground station, the water mass drifting backtracking length is obtained; and the time window is obtained by taking the satellite imaging time as the center and obtaining the time period of the water mass drifting backtracking length on both sides of the time sequence respectively.

[0053] As an example, first, the coverage radius of the ground station is divided by the lake flow rate to obtain the water mass drifting backtracking length; wherein the lake flow rate is an average flow rate calculated based on historical statistical information, for example, 0.1-0.3 m / s, and the middle value is taken in this embodiment; the coverage radius of the ground station needs to be determined based on its design parameters, for example, 5000 m; then the time window is determined by taking the satellite imaging time as the center, that is, ; wherein, is the satellite imaging time, is the water mass drifting backtracking length.

[0054] The water quality data and wind condition data of the ground station at each collection time within the time window can be further determined; wherein the monitoring objects (or collected water samples) of the ground station at different collection times correspond to different drifting water masses.

[0055] In order to facilitate subsequent drifting deduction, the latitude and longitude of the ground measuring point are converted into plane rectangular coordinates by using Gauss-Kruger projection or universal transverse Mercator projection, which is denoted as , wherein i is the code of the drifting water mass, so as to correspond to the geographical projection coordinates of the satellite spectral image. The above projection means are all known technologies and will not be described in detail.

[0056] The drift trajectory deduction module 102 is configured to determine the drift trajectory coordinates of each drift water mass at the satellite imaging time based on the wind condition data at each collection time at the ground site, and determine the cumulative turbulent energy of each drift water mass at the satellite imaging time.

[0057] Since the water body is continuously moving under the wind stress, and the wind condition of the lake surface will directly affect the drift direction and speed of the water mass, the embodiment of the present application determines the drift trajectory coordinates of each drift water mass at the satellite imaging time based on the wind condition data at each collection time at the ground site.

[0058] Preferably, in an embodiment of the present application, considering that wind energy is not directly converted into water flow power, and since the earth is rotating, the moving water flow is affected by the Coriolis force (geostrophic deflection force), the wind condition data can be analyzed based on a preset wind flow conversion coefficient and a preset flow direction deflection matrix, so as to estimate the lake surface flow velocity vector, and further estimate the drift displacement of each drift water mass at each collection time, thereby deducing the drift trajectory. Therefore, the method for obtaining the drift trajectory coordinates comprises:

[0059] At each collection time within the time window, the lake surface flow velocity vector is determined according to the wind condition data, the preset wind flow conversion coefficient and the preset flow direction deflection matrix; the time length between adjacent collection times is taken as the drift time length, and the drift displacement is determined according to the lake surface flow velocity vector and the drift time length.

[0060] For each drift water mass, within the time window, the drift trajectory coordinates of the next adjacent collection time are determined based on the superposition result of the position coordinates of the ground site and the drift displacement at the first collection time, and the drift trajectory coordinates at each collection time are obtained by sequentially superimposing.

[0061] As an example, taking any drift water mass as an example, and taking any collection time within the time window as an example, the analysis and description are as follows:

[0062] First, the lake surface flow velocity vector is determined ; wherein k is the collection time index, is the lake surface flow velocity vector at the kth collection time, is the wind flow conversion coefficient, usually taking a value of 0.02-0.04, and in the present embodiment, the middle value is taken; is the modulus length symbol; is the wind speed vector, which is expressed as , is the horizontal component of the wind speed, is the vertical component of the wind speed; is the preset flow direction deflection matrix, which is expressed as , is the deflection angle (suggested value: 20°-45°, used to describe the flow direction deflection caused by the Coriolis force; Used to characterize the wind speed after the transformation of airflow. Used to characterize the wind speed and direction after the airflow is deflected.

[0063] Due to the satellite imaging time , No. The actual location of the drifting water mass is not at the ground measuring point. If the collection time <Satellite Imaging Time> Explanation of the first The drifting water mass had already drifted downstream of the ground measurement point at the time of satellite imaging; if the acquisition time... >Satellite Imaging Time Explanation of the first If a drifting water mass is still upstream of the ground measurement point at the time of satellite imaging, then the drift inference within the time window also includes forward inference and backward backtracking.

[0064] like (Forward deduction), water mass moves with the current: ;

[0065] like (Reverse flow), water mass flowing backwards: ;

[0066] Where k is the data acquisition time number; Used to characterize the drift trajectory coordinates at the (K+1)th acquisition time; Used to characterize the drift trajectory coordinates at the Kth acquisition time; This represents the lake surface velocity vector at the k-th acquisition time. This refers to the drift duration, which is the duration between adjacent data collection times. .

[0067] Then, the drift trajectory coordinates of each drifting water mass at the time of satellite imaging are determined from the drift simulation trajectory. .

[0068] When the lake surface is calm, cyanobacteria usually float on the surface and can be clearly observed by satellites. However, wind and waves can cause turbulence, which can sweep the cyanobacteria to the bottom, making it possible for satellites to only observe a portion of the cyanobacteria, thus affecting the effectiveness of subsequent graded early warnings. Based on this, this embodiment of the invention further determines the cumulative turbulent energy of each drifting water mass at the time of satellite imaging based on wind condition data, in order to prepare for subsequent assessment of the surface algae retention index of the drifting water mass for accurate graded early warnings.

[0069] It should be noted that in one embodiment of the present application, the wind field of the lake surface is assumed to be uniform, so that at each collection time, even if the drifting water mass drifts to another place, its wind condition is consistent with the wind condition monitored by the ground station; in other embodiments, the implementer can also deploy multiple ground stations, and by default, the wind condition within the coverage radius of each station is consistent, and then determine the wind condition monitored by the adjacent ground station during the drifting process of the drifting water mass, to be applied to the analysis and calculation of the cumulative turbulent energy.

[0070] Preferably, in one embodiment of the present application, based on the classical empirical relationship in fluid dynamics and ocean physics, it can be deduced that the stirring power of the wind input water body is proportional to the cube of the wind speed, and then the turbulent energy caused by the wind wave at each collection time can be obtained by combining the stirring time, and then the cumulative turbulent energy of the drifting water mass at the satellite imaging time is determined; the method for obtaining the cumulative turbulent energy includes:

[0071] At each collection time within the time window, the wind wave stirring power coefficient is determined according to the wind condition data, the time length between adjacent collection times is taken as the stirring time, and the turbulent energy is obtained by fusing the wind wave stirring power coefficient and the stirring time;

[0072] For each drifting water mass, within the time window, based on the initial turbulent energy and the cumulative result of the turbulent energy at all collection times before the satellite imaging time, the cumulative turbulent energy of the drifting water mass at the satellite imaging time is determined.

[0073] As an example, taking any drifting water mass as an example, the analysis is described as follows:

[0074] Within the time window, the cube of the wind speed vector module length of each collection time is taken , and then multiplied by the stirring time, and the product is taken as the turbulent energy at the collection time; further, the cumulative sum of the turbulent energy at all collection times before the satellite imaging time is taken as the cumulative turbulent energy. Wherein, the initial turbulent energy of the drifting water mass is 0.

[0075] The physical constraint inversion module 103 is used to obtain the surface algal retention index of each drifting water mass according to the cumulative turbulent energy, and construct a flow direction constraint drifting model of the drifting water mass, to solve the model to obtain the pairing coefficient between each drifting water mass and each neighborhood pixel at the drifting trajectory coordinates in the satellite spectral image; according to the surface algal retention index of each drifting water mass and the corresponding water quality data, combined with the vegetation index of each neighborhood pixel and the pairing weight between each drifting water mass and each neighborhood pixel, the algal spectral response coefficient and the water body turbidity background parameter are obtained.

[0076] Because some algae are swept underwater by wind, waves, and turbulence, satellite observation signals show nonlinear signal loss compared to ground-based measurements. Therefore, this embodiment of the invention further obtains the surface algae retention index of each drifting water mass based on the accumulated turbulence energy. The surface algae retention index is used to characterize the proportion of cyanobacteria in the drifting water mass that can still be monitored by satellites under the influence of wind, waves, and turbulence. This prepares for the subsequent construction of a flow-direction-constrained drift model of the drifting water mass, in order to calculate the true algae spectral response coefficient from the complex flow field environment, and thus prepare for graded early warning.

[0077] Preferably, in one embodiment of the present invention, based on the Rouse Profile theory in fluid mechanics and the competition mechanism between Archimedes' buoyancy and turbulent diffusion, it can be deduced that cyanobacteria will actively rise to the surface under the influence of Archimedes' buoyancy, forming algal blooms, while turbulent eddies will mix the water, thus dispersing the cyanobacteria in the vertical direction; and the greater the turbulent energy, the more uniform the vertical distribution of cyanobacteria in the turbulence, and the lower the proportion remaining on the lake surface; therefore, the method for obtaining the surface algae retention index includes:

[0078] By comparing the cumulative turbulent energy of different drifting water masses, the relative turbulent mixing intensity coefficient of each drifting water mass is obtained; the negative correlation normalization result of the relative turbulent mixing intensity coefficient is used as the surface algae retention index of each drifting water mass.

[0079] As an example, the cumulative turbulent energy of each drifting water mass is first normalized to its maximum and minimum values, and the result of the normalization is used as the relative turbulent mixing intensity coefficient (dimensionless parameter) of the drifting water mass. Alternatively, the implementer can divide the cumulative turbulent energy of each drifting water mass by the maximum cumulative turbulent energy to obtain the relative turbulent mixing intensity coefficient. Furthermore, a variant of the Sigmoid function is used for negative correlation normalization to determine the surface algae retention index of the drifting water mass.

[0080] Specifically, the surface algae retention index ; , where is the surface algae retention index of the i-th drifting water mass; It is an exponential function with the natural constant e as its base; A preset sensitivity coefficient (within the range of 5-15, and 10 in this embodiment) is used to control the steepness of the curve; Let be the relative turbulent mixing intensity coefficient of the i-th drifting water mass; The preset phase transition midpoint (with a value range of 0.4-0.6, and 0.5 in this embodiment) is used to characterize the critical point at which turbulent mixing begins to significantly suppress buoyancy.

[0081] Since the drifting deduction of the drifting water mass is only a theoretical estimated value, it can be affected by other local turbulent flows during the drifting process, so that the drifting trajectory coordinates can have a certain deviation, that is, the pixels corresponding to the drifting trajectory coordinates in the satellite spectral image and the drifting water mass can not be completely consistent;

[0082] Based on this, the embodiment of the present application further constructs a flow direction constraint drifting model of the drifting water mass, and obtains the pairing coefficient between each drifting water mass and each neighborhood pixel at the drifting trajectory coordinate in the satellite spectral image by solving the model; the pairing coefficient represents the possibility of the pixel corresponding to the drifting water mass, and prepares for subsequent decoupling of the algal spectral response coefficient and the water body turbidity background parameter.

[0083] It should be noted that for each drifting water mass, the neighborhood pixel refers to the pixel in the preset neighborhood with the drifting trajectory coordinate as the center; in this embodiment, the size of the preset neighborhood is 5x5, and the implementer can also adjust it.

[0084] Considering that the drifting water mass usually drifts along the flow direction, but if the drifting water mass deviates or jumps out of the flow direction, it can violate the natural law, so first, the spatial matching cost term can be determined according to the distribution relationship between each drifting water mass and each neighborhood pixel at the corresponding drifting trajectory coordinate; the spatial matching cost term allows the drifting water mass to have a small error in the flow direction, that is, it is hoped that the pixels are matched along the flow direction, and when the pixels are matched vertically to the flow direction, there will be a larger matching cost;

[0085] In addition, considering that wind waves can roll blue-green algae underwater, the actual observed amount of blue-green algae or the drifting water mass observed by the satellite can be less than the sampling amount of the ground measurement point, so the ground distribution constraint term can be determined according to the water quality data and the surface algal retention index of each drifting water mass; the ground distribution constraint term allows data gaps, that is, part of the pixels can not be hard-matched with the drifting water mass, that is, it does not need to be completely matched one by one, so as to avoid matching distortion;

[0086] In addition, in the matching process, in order to pursue the closest distance between the pixels and the drifting water source, multiple drifting water masses can correspond to the same nearest pixel, and each pixel is theoretically best matched with only one drifting water mass, and if the pixels are forcibly matched with the adjacent drifting water mass, the matching of the drifting water mass and the pixel in the entire lake flow field can be directly distorted; therefore, the satellite observation constraint term can be determined according to the position distribution of the neighborhood pixels; the satellite observation constraint term is used to constrain the matching aggregation degree of the pixels.

[0087] Based on this, in one preferred embodiment of the present application, the method for obtaining the flow direction constraint drifting model comprises:

[0088] The pair coefficient between each drift water mass and each neighborhood pixel is taken as an optimization variable, a spatial matching cost term is determined in combination with a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term is determined in combination with water quality data and a surface algal retention index of each drift water mass, and a satellite observation constraint term is determined in combination with a position distribution of the neighborhood pixels; a flow-constrained drift model is constructed by fusing the spatial matching cost term, the ground distribution constraint term and the satellite observation constraint term.

[0089] As an example, an expression form of the flow-constrained drift model is as follows:

[0090] ; wherein, is the flow-constrained drift model; is a minimum value function; j is the serial number of the neighborhood pixel corresponding to the i th drift water mass; is the total number of the neighborhood pixels corresponding to the i th drift water mass; i is the serial number of the drift water mass; is the total number of the drift water masses; is the spatial matching cost term; is a preset regularization parameter (usually 0.01-0.1, and 0.1 is taken in the embodiment), used for controlling the relaxation degree of matching; is a KL divergence function, used for asymmetric measurement, that is, allowing the total mass of two distributions to be unequal; is a coefficient matrix of the pair coefficient between each drift water mass and each neighborhood pixel, that is, an optimization variable, the size of the coefficient matrix is I×J, and the matrix element is a (variable) pair coefficient, which is prepared for subsequent optimization solution; is a transposed matrix of the coefficient matrix; is a full 1 column vector, with a length of I; is a demand distribution vector of the pixel, also a full 1 column vector, with a length of I, used for constraining the aggregated matching of the pixel; is a ground supply distribution vector of the drift water mass, with a column vector length of I, used for hard matching of the pixel. For example, the matching of the drift water mass and the neighborhood pixel can be regarded as filling the pixel with the water mass;

[0091] In the ground distribution constraint term

[0092] , each vector element in the is used to represent the total amount of each drift water mass allocated to each neighborhood pixel, each vector element in the is used to represent the surface algal residue of each drift water mass after the action of turbulent stirring, and then the allocation closeness is evaluated through the KL divergence function, and the smaller the KL divergence function is, the closer the allocation is;

[0093] In the satellite observation constraint term , Each vector element in the vector is used to represent the total amount of all drift water masses allocated to a certain field pixel, Each vector element in the vector is used to represent the total amount of each neighborhood pixel that can be received, and then the allocation closeness is evaluated by the KL divergence function, the smaller the KL divergence function, the closer the allocation;

[0094] Wherein, in a preferred embodiment of the present application, the method for obtaining the spatial matching cost term comprises:

[0095] For each drift water mass, a displacement vector of the drift trajectory coordinates of the drift water mass is directed to each neighborhood pixel, and the displacement vector is decomposed into an along-flow drift component and a cross-flow drift component;

[0096] The along-flow drift component and the cross-flow drift component are weighted and summed to obtain the drift cost parameter between the drift water mass and each neighborhood pixel; the weighting weight of the cross-flow drift component is greater than that of the along-flow drift component;

[0097] The drift cost parameter is used to weight the optimization variable of the flow-constrained drift model to obtain the spatial matching cost term.

[0098] Specifically, for each drift water mass, first determine the displacement vector of the drift trajectory coordinates of the drift water mass directed to each neighborhood pixel, and then decompose the displacement vector into an along-flow drift component and a cross-flow drift component; wherein the direction of the along-flow drift component is consistent with the direction of the water flow, and the cross-flow drift component is perpendicular to the direction of the water flow;

[0099] Then the along-flow drift component and the cross-flow drift component are weighted and summed to obtain the drift cost parameter between the drift water mass and each neighborhood pixel; wherein the weighting weight of the along-flow drift component is set to 1; the length of the cross-flow drift component is set to a value greater than 10, and in this embodiment, 10 is taken, to severely punish the matching that violates the laws of fluid mechanics;

[0100] Further, the spatial matching cost term ; is the drift cost parameter between the i-th drift water mass and the j-th neighborhood pixel; is the (variable) pairing coefficient between the i-th drift water mass and the j-th neighborhood pixel, that is, the matrix element in the coefficient matrix.

[0101] Then solve the flow-constrained drift model so that S is minimized, thereby obtaining the pairing coefficient between each drift water mass and each neighborhood pixel at the drift trajectory coordinates in the satellite spectral image.

[0102] Since the vegetation index of each pixel in the satellite spectral image is partly derived from the contribution of algae and partly from the contribution of the background water body, and the matching coefficient between the pixel and the drift water mass can help to evaluate the correspondence confidence between them, the vegetation index can be further decomposed into and thereby the algae contribution proportional to the surface algal concentration of the water mass and the interference of the background water body can be separated from the mixed satellite signal (pixel vegetation index), so as to obtain the algal spectral response coefficient and the water body turbidity background parameter.

[0103] The algal spectral response coefficient is used to characterize the vegetation index increment generated by unit effective chlorophyll concentration on the satellite spectral image, and the water body turbidity background parameter is used to characterize the vegetation index base generated by background substances such as silt and colored dissolved organic matter in addition to algae.

[0104] In a preferred embodiment of the present application, the method for obtaining the algal spectral response coefficient and the water body turbidity background parameter comprises:

[0105] The water quality data at least includes the chlorophyll concentration; for each drift water mass, the chlorophyll concentration is weighted by the surface algal retention index, and the weighted result is taken as the surface chlorophyll concentration;

[0106] A linear regression model of the vegetation index of the pixel in the satellite spectral image is constructed; wherein the variable in the linear regression model is the surface chlorophyll concentration, the regression coefficient is the algal spectral response coefficient, and the bias is the water body turbidity background parameter;

[0107] The matching coefficient between each drift water mass and each neighborhood pixel is taken as the weight, and the regression coefficient and the bias in the linear regression model are solved based on the weighted least squares method, so as to obtain the algal spectral response coefficient and the water body turbidity background parameter.

[0108] As an example, first, the surface chlorophyll concentration of each drift water mass is calculated by weighting the chlorophyll concentration by the surface algal retention index; then assuming that the jth pixel corresponds to the ith drift water mass, a linear regression model like is constructed; is the vegetation index of the jth pixel; is the algal spectral response coefficient, that is, the regression coefficient; is the surface chlorophyll concentration of the ith drift water mass, that is, the variable; is the water body turbidity background parameter, that is, the bias.

[0109] Further, the matching coefficient between the drift water mass and the neighborhood pixel is taken as the weight, and the regression coefficient and the bias in the linear regression model are solved based on the weighted least squares method, so that the weighted residual sum of squares is minimized, that is, , the solution is obtained and .

[0110] The hierarchical early warning module 104 is configured to determine the hierarchical early warning threshold according to the algal spectral response coefficient and the water turbidity background parameter, and perform the hierarchical early warning of the water bloom based on the vegetation index of each pixel.

[0111] After the algal spectral response coefficient and the water turbidity background parameter are obtained, the static concentration value specified in the environmental monitoring standard can be converted into the vegetation index under the satellite spectrum, so as to determine the hierarchical early warning threshold, and perform the hierarchical early warning of the water bloom based on the vegetation index of each pixel.

[0112] Preferably, in an embodiment of the present application, the method for obtaining the hierarchical early warning threshold comprises:

[0113] The preset number of risk levels is determined, the algal spectral response coefficient is weighted by the preset early warning concentration of each risk level, and the weighted result and the fusion result of the water turbidity background parameter are taken as the hierarchical early warning threshold of the corresponding risk level.

[0114] As an example, the preset number is set to 3, and the risk levels are divided into a mild risk level, a moderate risk level and a severe risk level; then the preset early warning concentration of each risk level is determined based on the local environmental protection emergency standard, such as the upper limit of the chlorophyll concentration corresponding to the risk level, which will not be described again; then the preset early warning concentration of each risk level is multiplied by the algal spectral response coefficient, and the product is added to the water turbidity background parameter, and the sum is taken as the hierarchical early warning threshold of the corresponding risk level.

[0115] Preferably, in an embodiment of the present application, considering that the water dynamics and optical properties of different regions of the lake may have spatial heterogeneity, the farther the region is from the actual drift trajectory, the greater the inversion uncertainty is, and the lower the accuracy of subsequent hierarchical early warning may be, so the warning attribute needs to be determined; then the risk level is determined in combination with the vegetation index; then the hierarchical early warning of the water bloom based on the vegetation index of each pixel comprises:

[0116] In the satellite spectral image, for each pixel, the warning attribute of the pixel is marked based on the spatial distance between the pixel and the nearest drift trajectory coordinate, wherein the warning attribute of the spatial distance greater than the preset threshold is determined warning, and the warning attribute of the spatial distance less than or equal to the preset threshold is suspected warning; the hierarchical early warning threshold of the adjacent risk level is determined to determine the hierarchical early warning interval, for each pixel, when the vegetation index is in the corresponding hierarchical early warning interval, the pixel is marked as the corresponding risk level; the hierarchical early warning of the water bloom is performed based on the warning attribute and the risk level of the pixel.

[0117] As an example, in the satellite spectral image, for each pixel, first measure the spatial distance in the Euclidean distance, and the preset threshold value is 5-10km, and 5km is taken in this embodiment; when the distance between the coordinates of the drift trajectory of the nearest drift water mass is less than 5km, it can be considered that the matching deduction accuracy of the pixel is relatively large, and the hierarchical warning accuracy is also relatively high, and the warning attribute is set to the determined warning, and the rest is set to the suspected warning;

[0118] Then, the hierarchical warning interval is determined, for example, [0, A1], (A1, B1], (B1, C1], (C1, ∞); wherein, A1 is the hierarchical warning threshold value corresponding to the mild risk level, B1 is the hierarchical warning threshold value corresponding to the moderate risk level, and C1 is the hierarchical warning threshold value corresponding to the severe risk level; then, the risk level is marked based on the vegetation index, and the water bloom hierarchical warning is performed in combination with the warning attribute.

[0119] For example, when the vegetation index of the pixel is in [0, A1] and the warning attribute is the determined warning, the water bloom hierarchical warning result of the pixel is the determined no risk; when the vegetation index of the pixel is in (A1, B1] and the warning attribute is the suspected warning, the water bloom hierarchical warning result of the pixel is the suspected mild risk; and the like, which will not be repeated. When the determined no risk is determined, no marking is needed; the implementer can mark different risk levels by using different colors, and can also distinguish the warning attribute by using solid color patches and semi-transparent color patches.

[0120] Finally, the cyanobacterial bloom hierarchical warning map containing the explicit risk level and the warning attribute (confidence identification) is output.

[0121] In summary, the present application first obtains the vegetation index of each pixel in the satellite spectral image of the lake surface, and determines the time window based on the satellite imaging time; obtains the water quality data and wind condition data of each ground station at each collection time within the time window; then, the drift trajectory coordinates of each drift water mass at the satellite imaging time are deduced based on the wind condition, and the cumulative turbulent energy is determined, and then the surface algal retention index of each drift water mass is obtained; further, the flow direction constraint drift model of the drift water mass is constructed, and the model is solved to evaluate the pairing relationship between the drift water mass and the adjacent pixel; further, the algal spectral response coefficient and the water turbidity background parameter are obtained by combining the surface algal retention index of each drift water mass and the corresponding water quality data, and then the hierarchical warning threshold value is determined, and the water bloom hierarchical warning is performed based on the vegetation index of each pixel. The present application reconstructs the wind-driven flow field and introduces the turbulent mixing constraint, solves the inconsistency problem of point monitoring and satellite observation in the time and space scales and the physical space, decouples the algal spectral response coefficient and the water turbidity background parameter, and self-adaptively adjusts the hierarchical warning threshold value, thereby improving the accuracy of the cyanobacterial bloom hierarchical warning.

[0122] It is to be noted that the sequential order of the above-described embodiments of the present application only for the purpose of description, but not the advantages and disadvantages of the embodiments. The processes depicted in the drawings do not necessarily require the specific order or continuous order shown to achieve the desired results. In some embodiments, multi-task processing and parallel processing are also possible or can be advantageous.

[0123] Each of the embodiments in the specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other. Each embodiment focuses on the difference from other embodiments.

Claims

1. A cyanobacterial bloom grading early warning system that combines vegetation index and water quality parameters, characterized in that, The system comprises: a data acquisition module: for obtaining the vegetation index of each pixel in the satellite spectral image of the lake surface, and determining a time window based on the satellite imaging time; obtaining water quality data and wind condition data of the ground station at each collection time within the time window; a drift trajectory derivation module: for determining the drift trajectory coordinates of each drift water mass at the satellite imaging time based on the wind condition data at each collection time, and determining the cumulative turbulent energy of each drift water mass at the satellite imaging time; a physical constraint inversion module: for obtaining the surface algal retention index of each drift water mass according to the cumulative turbulent energy, and constructing a flow direction constrained drift model of the drift water mass to solve the model and obtain the pairing coefficient between each drift water mass and each neighborhood pixel at the drift trajectory coordinates in the satellite spectral image; and obtaining the algal spectral response coefficient and the water turbidity background parameter based on the surface algal retention index of each drift water mass and the corresponding water quality data, the vegetation index of each neighborhood pixel, and the pairing weight between each drift water mass and each neighborhood pixel. a hierarchical early warning module: for determining hierarchical early warning thresholds based on the algal spectral response coefficient and the water turbidity background parameter, and performing water bloom hierarchical early warning based on the vegetation index of each pixel.

2. The system according to claim 1, wherein the system is characterized in that, The method for obtaining the time window comprises: based on the flow velocity of the lake and the coverage radius of the ground station, obtaining a backtracking time length; taking the satellite imaging time as the center, obtaining time periods of the backtracking time length on both sides of the time sequence, and taking the corresponding time ranges of the two time periods as the time window.

3. The system according to claim 1, wherein the system is characterized in that, The method for obtaining the drift trajectory coordinates comprises: at each collection time within the time window, determining the lake surface flow velocity vector based on the wind condition data, a preset wind flow conversion coefficient, and a preset flow direction deflection matrix; taking the time length between adjacent collection times as a drift time length, and determining the drift displacement based on the lake surface flow velocity vector and the drift time length; for each drift water mass, within the time window, determining the drift trajectory coordinates of the next adjacent collection time based on the superposition result of the position coordinates of the ground station and the drift displacement at the first collection time, and sequentially superimposing to obtain the drift trajectory coordinates at each collection time.

4. The system for early warning of cyanobacterial bloom classification according to claim 1, wherein, The method for obtaining the cumulative turbulent energy comprises: at each collection time within the time window, determining a wind wave stirring power coefficient based on the wind condition data, taking the time length between adjacent collection times as a stirring time length, and fusing the wind wave stirring power coefficient and the stirring time length to obtain the turbulent energy; for each drift water mass, within the time window, determining the cumulative turbulent energy of the drift water mass at the satellite imaging time based on the initial turbulent energy and the cumulative result of the turbulent energy at all collection times before the satellite imaging time.

5. The system for early warning of cyanobacterial bloom classification according to claim 1, wherein, The method for obtaining the surface algal retention index comprises: comparing the cumulative turbulent energy of different drift water masses to obtain the relative turbulent mixing intensity coefficient of each drift water mass; and taking the negative correlation normalization result of the relative turbulent mixing intensity coefficient as the surface algal retention index of each drift water mass.

6. The system for early warning of cyanobacterial bloom classification according to the combination of vegetation index and water quality parameters of claim 1, characterized in that, The method for obtaining the flow direction constrained drift model comprises: The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels.

7. The system according to claim 6, wherein the system is characterized in that, The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels.

8. The system for early warning of cyanobacterial bloom classification according to claim 1, wherein, The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels.

9. The system for early warning of cyanobacterial bloom classification according to claim 1, wherein, The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels.

10. The system for early warning of cyanobacterial bloom classification according to claim 2, wherein, The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood pixels. The spatial matching cost term is determined by combining a matching coefficient between each drift water mass and each neighborhood pixel, a distribution relationship between each drift water mass and each neighborhood pixel at a corresponding drift trajectory coordinate, a ground distribution constraint term determined by combining the water quality data of each drift water mass and the surface algae retention index, and a satellite observation constraint term determined by combining a position distribution of the neighborhood

Images