Multi-time-sequence water eutrophication evaluation method based on satellite images, electronic equipment and storage medium

Abstract

The invention discloses a multi-time-sequence water eutrophication evaluation method based on satellite images, electronic equipment and a storage medium, and belongs to the technical field of satellite remote sensing quantitative inversion. In order to solve the problem that precision cannot meet application requirements, a shuffling module is designed, a DSAT1DShuffleNetV2 module is constructed, channel shuffling, multi-scale feature extraction, deep and shallow feature fusion and multi-attention information sharing are realized, feature information mining is effectively performed, and model feature information collection is guided more accurately; designing a double-weight logarithmic root-mean-square error loss function Wwrmse to guide the optimization of a multi-time-sequence water eutrophication evaluation model; the designed ShuffleGRUnet model can complete model training associated with multi-period data, and realizes data prediction of a single period and data prediction based on multi-period linkage. The method provided by the invention can stably and accurately realize multi-time-sequence satellite remote sensing image water eutrophication evaluation.

Description

Technical Field

[0001] The invention belongs to the technical field of satellite remote sensing quantitative inversion, and in particular relates to a multi-time series water eutrophication assessment method based on satellite images, electronic equipment and storage medium. Background Art

[0002] Eutrophication of lake water has become one of the global environmental problems. Eutrophication affects aquatic ecosystems worldwide. The continuous enrichment of nitrogen and phosphorus in lakes due to the drainage of farmland and the discharge of production and domestic wastewater has attracted widespread attention from countries around the world. However, human infringement on aquatic ecosystems is still increasing at an unprecedented rate. How to timely, accurately and comprehensively evaluate and grasp the status of lake eutrophication and monitor its changing trends, and for relevant departments to formulate corresponding countermeasures and measures to reduce and control lake eutrophication and give full play to the various functions that lake water bodies should have, is of great significance.

[0003] At present, the conventional monitoring method for inland water quality in my country is to first collect water samples at water quality sampling points in lakes, then conduct water quality laboratory analysis, and use single parameter evaluation index method or multi-parameter comprehensive evaluation method to evaluate water quality based on the analysis data. Although this method can make accurate analysis and evaluation of many water quality indicators, it is time-consuming, labor-intensive and uneconomical; and the number of water samples collected and analyzed is very limited. For the entire lake water body, these measurement point data only have local and typical representative significance.

[0004] The emergence and development of remote sensing technology has opened up new ways for the monitoring and research of lake eutrophication. Using satellite or aerial remote sensing information and based on water color remote sensing theory, water quality remote sensing monitoring applies empirical, semi-empirical or theoretical analysis methods, selects appropriate remote sensing band data, and establishes a water quality parameter remote sensing quantitative inversion model to invert the concentration of water quality parameters in the water body. Remote sensing means can be used to conduct large-scale eutrophication surveys and evaluations. The traditional statistical relationship model directly based on the single-scene remote sensing data and the measured water quality results can only fuzzily calculate the water quality results, and the monitoring effect is not good; moreover, the national general lake water quality model used by the existing lake water eutrophication monitoring means cannot accurately adapt to the water quality monitoring of each lake. Furthermore, due to poor adaptability, the monitoring effect is poor and the accuracy cannot meet the requirements.

[0005] With the comprehensive development of technologies combining artificial intelligence with remote sensing satellites, higher requirements have been put forward for the automation, intelligence and precision of water eutrophication assessment. In addition, the existing technologies only consider the correlation between a single period and cannot effectively associate the mutual constraints between data in multiple periods, resulting in poor generalization ability and stability of the model and the accuracy cannot meet application requirements. Summary of the invention

[0006] The present invention aims to solve the problems of poor model generalization ability and stability and inability to meet application requirements caused by water body eutrophication assessment based on single-period satellite image data, and proposes a multi-time series water body eutrophication assessment method based on satellite images, an electronic device and a storage medium.

[0007] To achieve the above object, the present invention is implemented through the following technical solutions:

[0008] A multi-time series water eutrophication assessment method based on satellite images includes the following steps:

[0009] S1. Collect paddy field sample data in the first, second, third and fourth quarters respectively, and use the collected sample data to calculate the nutritional status index of each sample point based on the comprehensive nutritional status index calculation formula as the target parameter of the model;

[0010] S2. Obtain multispectral satellite image data with a similar collection time to the sample data and preprocess them, and use the preprocessed satellite image data to calculate the spectral index of each sample point based on the spectral index calculation formula as the input parameter of the model;

[0011] S3. Design shuffle modules shuffle_block1 and shuffle_block2 respectively, and build DSAT_1D_ShuffleNetV2 module based on the designed shuffle modules for feature data extraction;

[0012] S4. Design a dual-weight logarithmic root mean square error loss function Wwrmse to guide the optimization of multi-time series water eutrophication assessment model;

[0013] S5. Based on the DSAT_1D_ShuffleNetV2 module designed in S3 and the Wwrmse loss function designed in step S4, combined with the gated recurrent neural network, a multi-time series water eutrophication assessment model Shuffle_GRUnet is designed;

[0014] S6. Input the target parameters obtained in step S1 and the input parameters obtained in step S2 into the multi-time series water eutrophication assessment model Shuffle_GRUnet designed in step S5 for model training to obtain the optimal multi-time series water eutrophication assessment model;

[0015] S7. Using the optimal multi-time series water eutrophication assessment model obtained in step S6, perform a multi-time series water eutrophication assessment in the target area.

[0016] Furthermore, the sample data of step S1 includes chlorophyll a data chl_a, total phosphorus data TP, total nitrogen data TN, transparency data SD, and chemical oxygen demand data CODmn.

[0017] Furthermore, the calculation formula of the spectral index in different bands of step S2 is as follows:

[0018] SPI b4,b1 =b4 / b1 (1)

[0019] SPI b4,b2 =b4 / b2 (2)

[0020] SPI b4,b3 =b4 / b3 (3)

[0021] SPI b3,b4,b1 =(b3+b4) / b1 (4)

[0022] SPI b3,b4,b2 =(b3+b4) / b2 (5)

[0023] SPI b1,b2,b4,b3 =(b1+b2+b4) / b3 (6)

[0024] SPI b2,b3,b4,b1 =(b2+b3+b4) / b1 (7)

[0025] SPI b1,b3,b4,b2 =(b1+b3+b4) / b2 (8)

[0026] SPI b3,bb4,b1,b2 =(b3+b4) / (b1+b2) (9)

[0027] Among them, SPI represents the spectral index, b1 represents the blue band, b2 represents the green band, b3 represents the infrared band, and b4 represents the near-infrared band.

[0028] Furthermore, the specific implementation method of step S3 includes the following steps:

[0029] S3.1. Design the shuffle_block1 module, which consists of 1 channel segmentation module Channel_Split, 5 batch normalization modules BN, 2 improved 1D multi-scale convolution attention MSCA_1D, and 3 1D convolution Conv1D 1,1 , 1 depth 1D convolution DWConv1D 3,1 , 3 Relu activation functions and 1 channel shuffle module Channel_Shuffle. The calculation process expression is as follows:

[0030] [in x1 ,in x2 ]=Channel_Split(in x ,0.5) (10)

[0031]

[0032] in x2_1 =Conv1D 1,1 (BN(DWConv1D 3,1 (Relu(BN(Conv1D 1,1 (in x2 )))))) (12)

[0033]

[0034] out y1 =Channel_Shuffle(Concat(in x1_1 ,in x2_2 )) (14)

[0035] Among them, x Indicates the input feature information of the shuffle_block1 module, in x1 ,in x2 Respectively represent the channel information after segmentation according to the first channel and the second channel, in x1_1 Indicates in x1 After MSCA_1D operation and in x1 The result of multiplication, in x2_1 Indicates in x2 After Conv1D 1,1 , BN, Relu, DWConv1D 3,1 The feature information after operation, in x2_2 Indicates in x2_1 After BN and MSCA_1D operations and in x2_1 The result of the multiplication is then passed through Conv1D 1,1 , BN, Relu after the feature information, out y1 Indicates in x1_1 ,in x2_2 The feature information after matrix connection and Channel_Shuffle operation, Channel_Split(tz,bl) means splitting the feature information tz into channels according to the ratio of bl, MSCA_1D(.) means the designed 1D multi-scale convolution attention module, BN(.) means batch normalization operation, Relu(.) means activation function operation, Conv1D 1,1 (.) indicates a 1D convolution operation with a kernel size of 1 and a stride of 1, DWConv1D 3,1(.) indicates a depthwise 1D convolution operation with a kernel size of 3 and a stride of 1. Concat(.) indicates a matrix concatenation operation. Channel_Shuffle(.) indicates a channel shuffle operation. Represents matrix multiplication;

[0036] S3.2. The designed shuffle_block2 module consists of 6 batch normalization modules BN, 2 improved 1D multi-scale convolution attention MSCA_1D, and 4 1D convolution Conv1D 1,1 , 2 depth 1D convolution DWConv1D 3,2 , 3 Relu activation functions and 1 channel shuffle module Channel_Shuffle. The calculation process expression is as follows:

[0037] in 4_1 =DWConv1D 3,2 (in x4 )(15)

[0038]

[0039] in x4_3 =Conv1D 1,1 (BN(DWConv1D 3,2 (Relu(BN(Conv1D 1,1 (in x4 ))))))(17)

[0040]

[0041] out y2 =Channel_Shuffle(Concat(in x4_2 ,in x4_4 ))(19)

[0042] Among them, x4 Indicates the input feature information of the shuffle_block2 module, in x4_1 Indicates that after DWConv1D 3,2 The feature information after operation, in x4_2 Indicates in x4_1 After BN, MSCA_1D, Conv1D 1,1 The result of the operation is then combined with in x4_1 The feature information after multiplication and BN and Relu operation, in x4_3 Indicates in x4 After Conv1D 1,1 , BN, Relu, DWConv1D 3,2 、Conv1D1,1 Feature information after operation, inx 4_4 Indicates in x4_3 After BN, MSCA_1D, Conv1D 1,1 The result of the operation is then in x4_3 After multiplication, the feature information after BN and Relu operation is obtained. y2 Indicates in x4_2 ,in x4_4 The feature information after matrix connection and Channel_Shuffle operation, MSCA_1D(.) represents the 1D multi-scale convolution attention module designed by the present invention, BN(.) represents batch normalization operation, Relu(.) represents activation function operation, Conv1D 1,1 (.) indicates a 1D convolution operation with a kernel size of 1 and a stride of 1, DWConv1D 3,2 (.) represents a depth 1D convolution operation with a kernel size of 3 and a stride of 2. Concat(.) represents a matrix concatenation operation. Channel_Shuffle(.) represents a channel shuffle operation. Represents matrix multiplication;

[0043] S3.3. The designed 1D multi-scale convolutional attention module MSCA_1D consists of 4 1D convolutions of different scales and depths DWConv1D, 1 1D convolution Conv1D 1,1 The calculation process is as follows:

[0044] in x5_1 =DWConv1D5(in x5 )(20)

[0045]

[0046] Among them, x5 Represents the input feature information of the 1D multi-scale convolutional attention module MSCA_1D, in x5_1 Indicates in x5 Feature information after DWConv1D5 operation, in x5_2 ,in x5_3 ,in x5_4 Respectively represent in x5_1 After DWConv1D7, DWConv1D 11 、DWConv1D 21 Feature information after operation, out y3 Indicates in x5_1 ,in x5_2 ,in x5_3 ,in x5_4The result after matrix addition and Conv1D1 operation is the same as in x5 Multiply feature information, DWConv1D q Represents the depth 1D convolution operation of q convolution kernel, q∈[1,5,7,11,21], Conv1D 1,1 1 indicates that the convolution kernel is 1 and the step size is 1, which is a convolution operation. Add indicates matrix addition. Represents matrix multiplication;

[0047] S3.4. The calculation process of channel shuffle Channel_Shuffle is as follows:

[0048] S3.4.1. Determine the number of groups of feature maps x and the number of channels C per group g The calculation formula is as follows:

[0049]

[0050] Among them, C represents the number of channels of the feature map x, G represents the number of groups to be divided, and C g Indicates the number of channels in each group;

[0051] S3.4.2.reshape operation:

[0052] Expand the channel dimension of the feature map x from one dimension to two dimensions, one is the number of groups G, and the other is the number of channels in each group C g ; If the batch size of the feature map is N and the spatial dimension is H×W, the reshape operation can be expressed as:

[0053] x reshaped =x.view(N,G,C g ,H,W) (24)

[0054] S3.4.3.transpose operation:

[0055] The reshaped feature map is transposed in terms of the number of groups and the number of channels in each group, so that the channels in each group can be disrupted. The transposition operation can be expressed as:

[0056] x transposed =x reshaped .transpose(1,2)(25)

[0057] Among them, transpose(1,2) means exchanging the second dimension and the third dimension;

[0058] S3.4.4.flatten operation:

[0059] The transposed feature map xtransposed , flatten back to the original number of channels C, and get the final channel shuffled feature map. The flatten operation is expressed as:

[0060] x shuffled =x transposed .view(N,C,H,W)(26)

[0061] Among them, view(N,C,H,W) restores the feature map to its original four-dimensional shape;

[0062] S3.5. The designed DSAT_1D_ShuffleNetV2 consists of 1 1D convolution Conv1D 3,2 , 1 1D convolution Conv1D 1,1 , 1 1-dimensional maximum pooling MaxPool 3,2 , 1 global average pooling GlobaPool, 3 designed shuffle_block2, 13 designed shuffle_block1, the calculation process expression is as follows:

[0063] Stage 0 = MaxPool 3,2 (Conv1D 3,2 (input x )) (27)

[0064]

[0065] Stage 2_1 =Concat(interpolate(Stage1,Stage2.shape()),Stage2)(29)

[0066] Stage 3_1 =Concat(interpolate(Stage 2_1 ,Stage3.shape()),Stage3)(30)

[0067] output y =GlobaPool(Conv1D 1,1 (Stage 3_1 )) (31)

[0068] Among them, input x Indicates the input feature information of the designed DSAT_1D_ShuffleNetV2, Stage0 represents input x After Conv1D 3,2 、MaxPool 3,2The feature information after the operation. Stage1, Stage2, and Stage3 represent the feature information after Stage0 is operated by shuffle_block2 and shuffle_block1 respectively. 2_1 It indicates the characteristic information obtained by interpolating the result of Stage1 with the size of Stage2 and connecting it with the matrix of Stage2. 3_1 Represents Stage 2_1 After interpolation operation based on the Stage 3 size, the interpolation result is connected with the Stage 3 matrix to obtain the feature information, output y Represents Stage 3_1 After Conv1D 1,1 , feature information of GlobaPool operation, Conv1D 1,1 Represents a 1D convolution operation with a convolution kernel of 1 and a stride of 1, Conv1D 3,2 (.) indicates a 1D convolution operation with a kernel size of 3 and a stride of 2. MaxPool 3,2 represents a 1D maximum pooling operation with a convolution kernel of 3 and a stride of 2. shuffle_block2(,v) represents a shuffle_block2 module operation designed v times, v∈[1,1,1], shuffle_block1(,w) represents a shuffle_block1 module operation designed w times, w∈[3,7,3], interpolate(.) represents an interpolation operation, Concat(.) represents a matrix concatenation operation, and GlobaPool(.) represents a global average pooling operation.

[0069] Furthermore, the calculation formula of the dual-weight logarithmic root mean square error loss function Wwrmse designed in step S4 is as follows:

[0070]

[0071] Among them, i represents the i-th period, i.e., the first quarter, the second quarter, the third quarter, and the fourth quarter, n represents the total number of periods, j represents the j-th sample data, m represents the total number of samples, and W i represents the weight of the i-th period, wrmse represents the logarithmic root mean square error of the weight designed by the present invention, represents the forecast value for the i-th period, represents the true value of the i-th period, represents the average value of the true value in period i, It represents the maximum value of the true value in the i-th period, cov represents the covariance operation, and σ represents the standard deviation operation.

[0072] Furthermore, the specific implementation method of step S5 includes the following steps:

[0073] The Shuffle_GRUnet model designed in S5.1. consists of 1 input layer, 1 DSAT_1D_ShuffleNetV2 layer designed in step S3, 1 gated recurrent neural network layer, and 1 fully connected layer. Each layer contains n components; input layer X = {1st period, 2nd period, 3rd period, 4th period...nth period}, each component input parameter contains 9 steps 2 Calculate the spectral index (SPI b4,b1 、SPI b4,b2 、SPI b4,b3 、SPI b3,b4,b1 、SPI b3,b4,b1 、SPI b1,b2,b4,b3 、SPI b2,b3,b4,b1 、SPI b1,b3,b4,b2 、SPI b3,b4,b1,b2 ), input in the format of a three-dimensional tensor, so that DSAT_1D_ShuffleNetV2 and GRU can be used for feature extraction and model training. The three-dimensional tensor used is (number of samples × time step × number of features), the number of features is 9, and the spectral index stacking method is used to expand it to 224. The insufficient part is supplemented with 0. The calculation process expression is as follows:

[0074]

[0075] Among them, X i represents the input parameters of the i-th period, Indicates the output parameter of the i-th period, i.e., the comprehensive nutritional status index, DSAT_1D_ShuffleNetV2 i (.) indicates the DSAT_1D_ShuffleNetV2 module operation designed in step 3 for the i-th period, GRU i (.) indicates that the GRU module operation is performed in the i-th period, fc i (.) indicates that the i-th period has undergone a full connection operation;

[0076] S5.2. The referenced GRU consists of an update gate and a reset gate. First, the update gate Z is calculated. t and reset gate r t The value of the input variable x t and the hidden layer result h at the previous moment t-1 The concatenated matrix is ​​transformed by sigmoid nonlinearity and then input into the update gate; the reset gate value will act on h t-1 and the input variable x t The concatenation is nonlinearly transformed and then activated by the tanh function to obtain the state of the candidate hidden layer at the current moment. By 1-Z ttimes h t-1 Store the information of the previous moment through Z t times Record the information at the current moment, and add the results of the two as the state output h of the hidden layer at the current moment t Finally, the weight matrix of the output layer is applied to h t On the top, we get the current model output y t , the calculation formula is as follows:

[0077] Z t =sigmoid(W Z [h t-1 ,x t ]) (36)

[0078] r t =sigmoid(W r [h t-1 ,x t ]) (37)

[0079]

[0080] y t =sigmoid(W o h t ) (40)

[0081] Among them, sigmoid represents the sigmoid activation function operation, tanh represents the tanh activation function operation, and Z t Indicates that the door state is updated at time t, r t Indicates that the door state is reset at time t, represents the candidate hidden door state at time t, h t represents the hidden layer state at time t, y t represents the model output at time t, x t represents the model input at time t, h t-1 represents the hidden layer state at the previous moment, x t , W Z represents the update gate weight matrix, W r Represents resetting the gate weight matrix, represents the candidate hidden layer weight matrix, W o Represents the output layer weight matrix.

[0082] Furthermore, the implementation method of step S6 is to input the target parameters prepared in step S1 and the input parameters prepared in step S2 into the multi-time series water eutrophication assessment model Shuffle_GRUnet designed in step S5 for model training, and use the Pearson correlation coefficient (r pcc ), mean root mean square error The model accuracy is evaluated to obtain the optimal multi-time series water eutrophication assessment model.

[0083] An electronic device comprises a memory and a processor, wherein the memory stores a computer program, and the processor implements the steps of a multi-time series water eutrophication assessment method based on satellite images when executing the computer program.

[0084] A computer-readable storage medium stores a computer program, which, when executed by a processor, implements a multi-time series water eutrophication assessment method based on satellite images.

[0085] Beneficial effects of the present invention:

[0086] The multi-time series water eutrophication assessment method based on satellite images described in the present invention designs a DSAT_1D_ShuffleNetV2 feature extraction module, which realizes channel shuffling, multi-scale feature extraction, deep and shallow feature fusion, and multi-scale attention information mining, and realizes multi-dimensional and multi-scale context information acquisition and sharing, effectively retains high-quality information, and more accurately guides the model feature information collection; the Shuffle_GRUnet model is designed to realize multi-period associated model training, realizes single-period data prediction and data prediction based on multi-period linkage, and designs a dual-weight logarithmic root mean square error loss function Wwrmse to guide the optimization of the multi-time series water eutrophication assessment model. The method of the present invention can efficiently and accurately realize multi-period water eutrophication assessment, evaluate and monitor water eutrophication faster, more comprehensively, and more economically in a larger range, comprehensively and timely grasp the eutrophication state of water bodies, and provide a scientific basis for formulating timely control measures for lake eutrophication. BRIEF DESCRIPTION OF THE DRAWINGS

[0087] Figure 1 It is a flow chart of a multi-time series water eutrophication assessment method based on satellite images according to the present invention;

[0088] Figure 2 The sampling point distribution diagram of the present invention, wherein (a) is the sampling point distribution in the first quarter, (b) is the sampling point distribution in the second quarter, (c) is the sampling point distribution in the third quarter, and (d) is the sampling point distribution in the fourth quarter;

[0089] Figure 3 This is the Shuffle_GRUnet model architecture diagram described in the present invention;

[0090] Figure 4 This is the DSAT_1D_ShuffleNetV2 module architecture diagram of the present invention;

[0091] Figure 5 This is the shuffle_block1 module architecture diagram of the present invention;

[0092] Figure 6 This is the shuffle_block2 module architecture diagram of the present invention;

[0093] Figure 7 This is the MSCA_1D module architecture diagram of the present invention;

[0094] Figure 8 This is a diagram of the Channel_Shuffle module architecture of the present invention;

[0095] Fig. 9 These are the water eutrophication results maps described in the present invention, wherein (a) is the water eutrophication results map for the first quarter, (b) is the water eutrophication results map for the second quarter, (c) is the water eutrophication results map for the third quarter, and (d) is the water eutrophication results map for the fourth quarter. DETAILED DESCRIPTION

[0096] In order to make the purpose, technical solution and advantages of the present invention more clear, the present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention, that is, the specific embodiments described are only part of the embodiments of the present invention, rather than all of the specific embodiments. The components of the specific embodiments of the present invention described and shown in the drawings herein can be arranged and designed in various different configurations, and the present invention can also have other embodiments.

[0097] Therefore, the following detailed description of the specific embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents the selected specific embodiments of the present invention. Based on the specific embodiments of the present invention, all other specific embodiments obtained by those skilled in the art without making creative work are within the scope of protection of the present invention.

[0098] In order to further understand the content, features and effects of the present invention, the following specific implementation methods are given as examples, and the attached Figure 1 -Attached Fig. 9 The detailed instructions are as follows:

[0099] Embodiment 1:

[0100] A multi-time series water eutrophication assessment method based on satellite images includes the following steps:

[0101] S1. Collect paddy field sample data in the first, second, third and fourth quarters respectively, and use the collected sample data to calculate the nutritional status index of each sample point based on the comprehensive nutritional status index calculation formula as the target parameter of the model;

[0102] The sample data of step S1 includes chlorophyll a data chl_a, total phosphorus data TP, total nitrogen data TN, transparency data SD, and chemical oxygen demand data CODmn;

[0103] Furthermore, in step S1, samples were collected from a total of 50 lakes, ponds, and reservoirs, and 150 sample points were collected each quarter, totaling 600. The collection, storage, and transportation of water quality samples strictly follow the relevant standards for lake water quality sampling technology, sample storage, and processing specified in the "Technical Guide for Sampling of Water Quality Lakes and Reservoirs" (GB / T 14581-93). Use a water quality sample collection device to collect water quality samples 0.5m below the water surface. Add the corresponding fixative to water quality samples of different indicators, store them in an environment of 2 to 5°C, and deliver them to the laboratory within no more than 48 hours. Except for transparency (SD), which requires field measurement using the Secco plate visual method, the other four indicators are tested in the laboratory using the corresponding methods. Among them, total nitrogen (TN), total phosphorus (TP), and chlorophyll (chl_a) are measured by spectrophotometry, and chemical oxygen demand (CODmn) is measured by potassium dichromate method.

[0104] Furthermore, the calculation formula of the comprehensive nutritional status index in step S1 is as follows:

[0105]

[0106] Among them, TLI (∑) represents the comprehensive nutritional status index, W k represents the relevant weight of the nutritional status index of the kth parameter, TLI(k) represents the nutritional status index of the kth parameter, and K represents the number of evaluation parameters.

[0107] Taking chlorophyll a (chl_a) as the benchmark parameter, the normalized correlation weight calculation formula of the kth parameter is as follows:

[0108]

[0109] Among them, r fk Represents the correlation coefficient between the kth parameter and the reference parameter chl_a.

[0110] The relevant system statistics are as follows:

[0111] Table 1: Correlation coefficient table of benchmark parameter chl_a

[0112] parameter chl_a TP TN SD COD <![CDATA[r fk ]]> 1 0.84 0.82 0.83 0.83 <![CDATA[r fk 2 ]]> 1 0.7056 0.6724 0.6889 0.6889

[0113] The calculation formula of the kth parameter nutritional status index TLI(k) is as follows, where the parameters are chlorophyll a (chl_a), total phosphorus (TP), total nitrogen (TN), transparency (SD), and chemical oxygen demand (CODmn).

[0114] TLI(chl_a)=10×(2.5+1.086lnchl_a); (3)

[0115] TLI(TP)=10×(9.436+1.624lnTP); (4)

[0116] TLI(TN)=10×(5.453+1.694lnTN); (5)

[0117] TLI(SD)=10×(5.118-1.94lnSD); (6)

[0118] TLI(CODmn)=10×(0.109+2.661lnCODmn); (7)

[0119] The unit of chlorophyll a (chl_a) is mg / m 3 , the unit of transparency (SD) is m; the units of total phosphorus (TP), total nitrogen (TN) and chemical oxygen demand (CODmn) are all mg / L.

[0120] S2. Obtain multispectral satellite image data with a similar collection time to the sample data and preprocess them, and use the preprocessed satellite image data to calculate the spectral index of each sample point based on the spectral index calculation formula as the input parameter of the model;

[0121] Furthermore, in step S2, 20 satellite image data of March, June, September and December are obtained respectively, and the image data are preprocessed by radiometric calibration, atmospheric correction, RPC orthorectification, image fusion and geometric precision correction.

[0122] Furthermore, the spectral index calculation formula based on the pre-processed satellite image data in step S2 is as follows:

[0123] SPI b4,b1 =b4 / b1; (8)

[0124] SPI b4,b2 =b4 / b2; (9)

[0125] SPI b4,b3 =b4 / b3; (10)

[0126] SPI b3,b4,b1 =(b3+b4) / b1; (11)

[0127] SPI b3,b4,b2 =(b3+b4) / b2; (12)

[0128] SPI b1,b2,b4,b3 =(b1+b2+b4) / b3; (13)

[0129] SPI b2,b3,b4,b1 =(b2+b3+b4) / b1; (14)

[0130] SPI b1,b3,b4,b2 =(b1+b3+b4) / b2; (15)

[0131] SPI b3,b4,b1,b2 =(b3+b4) / (b1+b2); (16)

[0132] Among them, SPI represents different spectral indices, b1 represents the blue band, b2 represents the green band, b3 represents the infrared band, and b4 represents the near-infrared band.

[0133] S3. Design shuffle modules shuffle_block1 and shuffle_block2 respectively, and build DSAT_1D_ShuffleNetV2 module based on the designed shuffle modules for feature data extraction;

[0134] Furthermore, the specific implementation method of step S3 includes the following steps:

[0135] S3.1. Design the shuffle_block1 module, which consists of 1 channel segmentation module Channel_Split, 5 batch normalization modules BN, 2 improved 1D multi-scale convolution attention MSCA_1D, and 3 1D convolution Conv1D 1,1 , 1 depth 1D convolution DWConv1D 3,1 , 3 Relu activation functions and 1 channel shuffle module Channel_Shuffle. The calculation process expression is as follows:

[0136] [in x1 ,in x2 ]=Channel_Split(in x ,0.5); (17)

[0137]

[0138] in x2_1 =Conv1D 1,1 (BN(DWConv1D 3,1 (Relu(BN(Conv1D 1,1 (inx2 ))))))); (19)

[0139]

[0140] out y1 =Channel_Shuffle(Concat(in x1_1 ,in x2_2 )); (twenty one)

[0141] Among them, x Indicates the input feature information of the shuffle_block1 module, in x1 ,in x2 Respectively represent the channel information after segmentation according to the first channel and the second channel, in x1_1 Indicates in x1 After MSCA_1D operation and in x1 The result of multiplication, in x2_1 Indicates in x2 After Conv1D 1,1 , BN, Relu, DWConv1D 3,1 The feature information after operation, in x2_2 Indicates in x2_1 After BN and MSCA_1D operations and in x2_1 The result of the multiplication is then passed through Conv1D 1,1 , BN, Relu after the feature information, out y1 Indicates in x1_1 ,in x2_2 The feature information after matrix connection and Channel_Shuffle operation, Channel_Split(tz,bl) means splitting the feature information tz into channels according to the ratio of bl, MSCA_1D(.) means the designed 1D multi-scale convolution attention module, BN(.) means batch normalization operation, Relu(.) means activation function operation, Conv1D 1,1 (.) indicates a 1D convolution operation with a kernel size of 1 and a stride of 1, DWConv1D 3,1 (.) indicates a depthwise 1D convolution operation with a kernel size of 3 and a stride of 1. Concat(.) indicates a matrix concatenation operation. Channel_Shuffle(.) indicates a channel shuffle operation. Represents matrix multiplication;

[0142] S3.2. The designed shuffle_block2 module consists of 6 batch normalization modules BN, 2 improved 1D multi-scale convolution attention MSCA_1D, and 4 1D convolution Conv1D 1,1, 2 depth 1D convolution DWConv1D 3,2 , 3 Relu activation functions and 1 channel shuffle module Channel_Shuffle. The calculation process expression is as follows:

[0143] in 4_1 =DWConv1D 3,2 (in x4 );(twenty two)

[0144]

[0145] in x4_3 =Conv1D 1,1 (BN(DWConv1D 3,2 (Relu(BN(Conv1D 1,1 (in x4 )))))); (twenty four)

[0146]

[0147] out y2 =Channel_Shuffle(Concat(in x4_2 ,in x4_4 )); (26)

[0148] Among them, x4 Indicates the input feature information of the shuffle_block2 module, in x4_1 Indicates that after DWConv1D 3,2 The feature information after operation, in x4_2 Indicates in x4_1 After BN, MSCA_1D, Conv1D 1,1 The result of the operation is then combined with in x4_1 The feature information after multiplication and BN and Relu operation, in x4_3 Indicates in x4 After Conv1D 1,1 , BN, Relu, DWConv1D 3,2 、Conv1D 1,1 Feature information after operation, inx 4_4 Indicates in x4_3 After BN, MSCA_1D, Conv1D 1,1 The result of the operation is then in x4_3 After multiplication, the feature information after BN and Relu operation is obtained. y2 Indicates in x4_2 ,in x4_4The feature information after matrix connection and Channel_Shuffle operation, MSCA_1D(.) represents the 1D multi-scale convolution attention module designed by the present invention, BN(.) represents batch normalization operation, Relu(.) represents activation function operation, Conv1D 1,1 (.) indicates a 1D convolution operation with a kernel size of 1 and a stride of 1, DWConv1D 3,2 (.) represents a depth 1D convolution operation with a kernel size of 3 and a stride of 2. Concat(.) represents a matrix concatenation operation. Channel_Shuffle(.) represents a channel shuffle operation. Represents matrix multiplication;

[0149] S3.3. The designed 1D multi-scale convolutional attention module MSCA_1D consists of 4 1D convolutions of different scales and depths DWConv1D, 1 1D convolution Conv1D 1,1 The calculation process is as follows:

[0150] in x5_1 =DWConv1D5(in x5 ); (27)

[0151]

[0152]

[0153] Among them, x5 Represents the input feature information of the 1D multi-scale convolutional attention module MSCA_1D, in x5_1 Indicates in x5 Feature information after DWConv1D5 operation, in x5_2 ,in x5_3 ,in x5_4 Respectively represent in x5_1 After DWConv1D7, DWConv1D 11 、DWConv1D 21 Feature information after operation, out y3 Indicates in x5_1 ,in x5_2 ,in x5_3 ,in x5_4 The result after matrix addition and Conv1D1 operation is the same as in x5 Multiply feature information, DWConv1D q Represents the depth 1D convolution operation of q convolution kernel, q∈[1,5,7,11,21], Conv1D 1,1 1 indicates that the convolution kernel is 1 and the step size is 1, which is a convolution operation. Add indicates matrix addition. Represents matrix multiplication;

[0154] S3.4. The calculation process of channel shuffle Channel_Shuffle is as follows:

[0155] S3.4.1. Determine the number of groups of feature maps x and the number of channels C per group g The calculation formula is as follows:

[0156]

[0157] Among them, C represents the number of channels of the feature map x, G represents the number of groups to be divided, and C g Indicates the number of channels in each group;

[0158] S3.4.2.reshape operation:

[0159] Expand the channel dimension of the feature map x from one dimension to two dimensions, one is the number of groups G, and the other is the number of channels in each group C g ; If the batch size of the feature map is N and the spatial dimension is H×W, the reshape operation can be expressed as:

[0160] x reshaped =x.view(N,G,C g ,H,W); (31)

[0161] S3.4.3.transpose operation:

[0162] The reshaped feature map is transposed in terms of the number of groups and the number of channels in each group, so that the channels in each group can be disrupted. The transposition operation can be expressed as:

[0163] x transposed =x reshaped .transpose(1,2);(32)

[0164] Among them, transpose(1,2) means to transform the second dimension (number of groups G) and the third dimension (number of channels per group C) g ) for exchange;

[0165] S3.4.4.flatten operation:

[0166] The transposed feature map x transposed , flatten back to the original number of channels C, and get the final channel shuffled feature map. The flatten operation can be expressed as:

[0167] x shuffled =x transposed.view(N,C,H,W); (33)

[0168] Among them, view(N,C,H,W) restores the feature map to its original four-dimensional shape;

[0169] S3.5. Based on the designed shuffle_block1 and shuffle_block2 modules, the DSAT_1D_ShuffleNetV2 module is designed for feature data extraction. The specific implementation method includes the following steps:

[0170] The designed DSAT_1D_ShuffleNetV2 consists of 1 1D convolution Conv1D 3,2 , 1 1D convolution Conv1D 1,1 , 1 1-dimensional maximum pooling MaxPool 3,2 , 1 global average pooling GlobaPool, 3 designed shuffle_block2, 13 designed shuffle_block1, the calculation process expression is as follows:

[0171] Stage 0 = MaxPool 3,2 (Conv1D 3,2 (input x )); (34)

[0172]

[0173] Stage 2_1 =Concat(interpolate(Stage1,Stage2.shape()),Stage2); (36)

[0174] Stage 3_1 =Concat(interpolate(Stage 2_1 ,Stage3.shape()),Stage3);(37)

[0175] output y =GlobaPool(Conv1D 1,1 (Stage 3_1 )); (38)

[0176] Among them, inputx represents the input feature information of the designed DSAT_1D_ShuffleNetV2, and Stage0 represents input x After Conv1D 3,2 、MaxPool 3,2The feature information after the operation. Stage1, Stage2, and Stage3 represent the feature information after Stage0 is operated by shuffle_block2 and shuffle_block1 respectively. 2_1 It indicates the characteristic information obtained by interpolating the result of Stage1 with the size of Stage2 and connecting it with the matrix of Stage2. 3_1 Represents Stage 2_1 After interpolation operation based on the Stage 3 size, the interpolation result is connected with the Stage 3 matrix to obtain the feature information, output y Represents Stage 3_1 After Conv1D 1,1 , feature information of GlobaPool operation, Conv1D 1,1 Represents a 1D convolution operation with a convolution kernel of 1 and a stride of 1, Conv1D 3,2 (.) indicates a 1D convolution operation with a kernel size of 3 and a stride of 2. MaxPool 3,2 represents a 1D maximum pooling operation with a convolution kernel of 3 and a stride of 2. shuffle_block2(,v) represents a shuffle_block2 module operation designed v times, v∈[1,1,1], shuffle_block1(,w) represents a shuffle_block1 module operation designed w times, w∈[3,7,3], interpolate(.) represents an interpolation operation, Concat(.) represents a matrix concatenation operation, and GlobaPool(.) represents a global average pooling operation.

[0177] S4. Design a dual-weight logarithmic root mean square error loss function Wwrmse to guide the optimization of multi-time series water eutrophication assessment model;

[0178] Furthermore, the calculation formula of the dual-weight logarithmic root mean square error loss function Wwrmse designed in step S4 is as follows:

[0179]

[0180] Among them, i represents the i-th period, i.e., the first quarter, the second quarter, the third quarter, and the fourth quarter, n represents the total number of periods, j represents the j-th sample data, m represents the total number of samples, and W i represents the weight of the i-th period, wrmse represents the logarithmic root mean square error of the weight designed by the present invention, represents the forecast value for the i-th period, represents the true value of the i-th period, represents the average value of the true value in period i, It represents the maximum value of the true value in the i-th period, cov represents the covariance operation, and σ represents the standard deviation operation.

[0181] S5. Based on the DSAT_1D_ShuffleNetV2 module designed in S3 and the Wwrmse loss function designed in step S4, combined with the gated recurrent neural network, a multi-time series water eutrophication assessment model Shuffle_GRUnet is designed;

[0182] Furthermore, the specific implementation method of step S5 includes the following steps:

[0183] The Shuffle_GRUnet model designed in S5.1. consists of 1 input layer, 1 DSAT_1D_ShuffleNetV2 layer designed in step S3, 1 gated recurrent neural network layer, and 1 fully connected layer. Each layer contains n components; input layer X = {1st period, 2nd period, 3rd period, 4th period...nth period}, each component input parameter contains 9 steps 2 Calculate the spectral index (SPI b4,b1 、SPI b4,b2 、SPI b4,b3 、SPI b3,b4,b1 、SPI b3,b4,b1 、SPI b1,b2,b4,b3 、SPI b2,b3,b4,b1 、SPI b1,b3,b4,b2 、SPI b3,b4,b1,b2 ), input in the format of a three-dimensional tensor, so that DSAT_1D_ShuffleNetV2 and GRU can be used for feature extraction and model training. The three-dimensional tensor used is (number of samples × time step × number of features), the number of features is 9, and the spectral index stacking method is used to expand it to 224. The insufficient part is supplemented with 0. The calculation process expression is as follows:

[0184]

[0185] Among them, X i represents the input parameters of the i-th period, Indicates the output parameter of the i-th period, i.e., the comprehensive nutritional status index, DSAT_1D_ShuffleNetV2 i (.) indicates the DSAT_1D_ShuffleNetV2 module operation designed in step 3 for the i-th period, GRU i (.) indicates that the GRU module operation is performed in the i-th period, fc i (.) indicates that the i-th period has undergone a full connection operation;

[0186] S5.2. The referenced GRU consists of an update gate and a reset gate. First, the update gate Z is calculated. t and reset gate rt The value of the input variable x t and the hidden layer result h at the previous moment t-1 The concatenated matrix is ​​transformed by sigmoid nonlinearity and then input into the update gate; the reset gate value will act on h t-1 and the input variable x t The concatenation is nonlinearly transformed and then activated by the tanh function to obtain the state of the candidate hidden layer at the current moment. By 1-Z t times h t-1 Store the information of the previous moment through Z t times Record the information at the current moment, and add the results of the two as the state output h of the hidden layer at the current moment t Finally, the weight matrix of the output layer is applied to h t On the top, we get the current model output y t , the calculation formula is as follows:

[0187] Z t =sigmoid(W Z [h t-1 ,x t ]); (43)

[0188] r t =sigmoid(W r [h t-1 ,x t ]); (44)

[0189]

[0190]

[0191] y t =sigmoid(W o h t ); (47)

[0192] Among them, sigmoid represents the sigmoid activation function operation, tanh represents the tanh activation function operation, and Z t Indicates that the door state is updated at time t, r t Indicates that the door state is reset at time t, represents the candidate hidden door state at time t, h t represents the hidden layer state at time t, y t represents the model output at time t, x t represents the model input at time t, h t-1 represents the hidden layer state at the previous moment, x t , W Zrepresents the update gate weight matrix, W r Represents resetting the gate weight matrix, represents the candidate hidden layer weight matrix, W o Represents the output layer weight matrix.

[0193] S6. Input the target parameters obtained in step S1 and the input parameters obtained in step S2 into the multi-time series water eutrophication assessment model Shuffle_GRUnet designed in step S5 for model training to obtain the optimal multi-time series water eutrophication assessment model;

[0194] Furthermore, in step S6, the target parameters prepared in step S1 and the input parameters prepared in step S2 are input into the multi-time series water eutrophication assessment model Shuffle_GRUnet designed in step S5 for model training. pcc ), mean root mean square error The model accuracy was evaluated to obtain the optimal multi-time series water eutrophication assessment model. The calculation formula of the evaluation index is as follows:

[0195]

[0196] Among them, r pcc represents the Pearson correlation coefficient, represents the mean root mean square error, i represents the i-th period, i.e. the first quarter, the second quarter, the third quarter, and the fourth quarter, j represents the j-th sample data, and m represents the total number of samples. represents the forecast value for the i-th period, represents the true value of the i-th period, It represents the average value of the true value in the i-th period, cov represents the covariance operation, and σ represents the standard deviation operation.

[0197] S7. Using the optimal multi-time series water eutrophication assessment model obtained in step S6, perform a multi-time series water eutrophication assessment in the target area.

[0198] Furthermore, in step S7, a known nutrient status classification method (such as the evaluation method in the "Surface Water Environmental Quality Evaluation Method (Trial)") is referred to to distinguish the eutrophication degree of the water body. In this way, based on the spectral characteristic information of the eutrophic water body, a eutrophication segmentation classification index is constructed, and the specific classification standards are shown in the following table:

[0199] Table 2: Eutrophication classification standards

[0200]

[0201] To sum up: this implementation case makes full use of the characteristic information of the satellite image spectral index, links the characteristic information of different periods, realizes the multi-time series water body eutrophication assessment, and solves the problem that the model's generalization ability and stability are poor and the accuracy cannot meet the application requirements when conducting water body eutrophication assessment based on single-period satellite image data.

[0202] Embodiment 2:

[0203] An electronic device includes a memory and a processor, wherein the memory stores a computer program, and when the processor executes the computer program, the steps of a multi-time series water eutrophication assessment method based on satellite images described in Example 1 are implemented.

[0204] The computer device of the present invention may be a device including a processor and a memory, such as a single chip microcomputer including a central processing unit, etc. Moreover, the processor is used to implement the steps of the multi-time series water eutrophication assessment method based on satellite images when executing the computer program stored in the memory.

[0205] The processor may be a central processing unit (CPU), other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or any conventional processor, etc.

[0206] The memory may mainly include a program storage area and a data storage area, wherein the program storage area may store an operating system, an application required for at least one function (such as a sound playback function, an image playback function, etc.), etc.; the data storage area may store data created according to the use of the mobile phone (such as audio data, a phone book, etc.), etc. In addition, the memory may include a high-speed random access memory, and may also include a non-volatile memory, such as a hard disk, a memory, a plug-in hard disk, a smart memory card (Smart Media Card, SMC), a secure digital (Secure Digital, SD) card, a flash card (Flash Card), at least one disk storage device, a flash memory device, or other volatile solid-state storage devices.

[0207] Embodiment 3:

[0208] A computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements a multi-time series water eutrophication assessment method based on satellite images as described in Example 1.

[0209] The computer-readable storage medium of the present invention can be any form of storage medium that can be read by a processor of a computer device, including but not limited to non-volatile memory, volatile memory, ferroelectric memory, etc. A computer program is stored on the computer-readable storage medium. When the processor of the computer device reads and executes the computer program stored in the memory, the steps of the above-mentioned multi-time series water eutrophication assessment method based on satellite images can be implemented.

[0210] The computer program includes computer program code, which may be in source code form, object code form, executable file or some intermediate form, etc. The computer readable medium may include: any entity or device capable of carrying the computer program code, recording medium, USB flash drive, mobile hard disk, magnetic disk, optical disk, computer memory, read-only memory (ROM), random access memory (RAM), electric carrier signal, telecommunication signal and software distribution medium, etc. It should be noted that the content contained in the computer readable medium may be appropriately increased or decreased according to the requirements of legislation and patent practice in the jurisdiction. For example, in some jurisdictions, according to legislation and patent practice, computer readable media do not include electric carrier signals and telecommunication signals.

[0211] It should be noted that relational terms such as "first" and "second" are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device. In the absence of further restrictions, the elements defined by the sentence "comprise a ..." do not exclude the existence of other identical elements in the process, method, article or device including the elements.

[0212] Although the present application has been described above with reference to specific embodiments, various modifications may be made thereto and parts thereof may be replaced with equivalents without departing from the scope of the present application. In particular, as long as there is no structural conflict, the various features in the specific embodiments disclosed in the present application may be used in combination with each other in any manner, and the fact that these combinations are not exhaustively described in this specification is only for the sake of omitting space and saving resources. Therefore, the present application is not limited to the specific embodiments disclosed herein, but includes all technical solutions falling within the scope of the claims.

Claims

1. A multi-temporal water eutrophication assessment method based on satellite images, characterized in that: The steps include: S1. Collect paddy field sample data in the first, second, third and fourth quarters respectively, and use the collected sample data to calculate the nutritional status index of each sample point based on the comprehensive nutritional status index calculation formula as the target parameter of the model; S2. Obtain multispectral satellite image data with a similar collection time to the sample data and perform preprocessing, and use the preprocessed satellite image data to calculate the spectral index of each sample point based on the spectral index calculation formula as the input parameter of the model; S3. Design shuffle modules shuffle_block1 and shuffle_block2 respectively, and build DSAT_1D_ShuffleNetV2 module based on the designed shuffle modules for feature data extraction; S4. Design a dual-weight logarithmic root mean square error loss function Wwrmse to guide the optimization of multi-time series water eutrophication assessment model; S5. Based on the DSAT_1D_ShuffleNetV2 module designed in S3 and the Wwrmse loss function designed in step S4, combined with the gated recurrent neural network, a multi-time series water eutrophication assessment model Shuffle_GRUnet is designed; S6. Input the target parameters obtained in step S1 and the input parameters obtained in step S2 into the multi-time series water eutrophication assessment model Shuffle_GRUnet designed in step S5 for model training to obtain the optimal multi-time series water eutrophication assessment model; S7. Using the optimal multi-time series water eutrophication assessment model obtained in step S6, perform a multi-time series water eutrophication assessment in the target area.

2. The multi-time series water eutrophication assessment method based on satellite images according to claim 1 is characterized in that: The sample data of step S1 includes chlorophyll a data chl_a, total phosphorus data TP, total nitrogen data TN, transparency data SD, and chemical oxygen demand data CODmn.

3. The multi-time series water eutrophication assessment method based on satellite images according to claim 2 is characterized in that: The calculation formula of the spectral index in different bands of step S2 is as follows: SPI b4,b1 =b4 / b1 (1) SPI b4,b2 =b4 / b2 (2) SPI b4,b3 =b4 / b3 (3) SPI b3,b4,b1 =(b3+b4) / b1 (4) SPI b3,b4,b2 =(b3+b4) / b2 (5) <h2 style=";text-align:left;direction:ltr">SPI<h2 style=";text-align:left;direction:ltr"> b1,b2,b4,b3 <h2 style=";text-align:left;direction:ltr"> (b1+b2+b4) / b3 (6) <h2 style=";text-align:left;direction:ltr">SPI<h2 style=";text-align:left;direction:ltr"> b2,b3,b4,b1 <h2 style=";text-align:left;direction:ltr"> (b2+b3+b4) / b1 (7) <h2 style=";text-align:left;direction:ltr">SPI<h2 style=";text-align:left;direction:ltr"> b1,b3,b4,b2 <h2 style=";text-align:left;direction:ltr"> (b1+b3+b4) / b2 (8) <h2 style=";text-align:left;direction:ltr">SPI<h2 style=";text-align:left;direction:ltr"> b3,b4,b1,b2 <h2 style=";text-align:left;direction:ltr"> (b3+b4) / (b1+b2) (9) Among them, SPI represents the spectral index, b1 represents the blue band, b2 represents the green band, b3 represents the infrared band, and b4 represents the near-infrared band.

4. The multi-time series water eutrophication assessment method based on satellite images according to claim 3 is characterized in that: The specific implementation method of step S3 includes the following steps: S3.

1. Design the shuffle_block1 module, which consists of 1 channel segmentation module Channel_Split, 5 batch normalization modules BN, 2 improved 1D multi-scale convolution attention MSCA_1D, and 3 1D convolution Conv1D 1,1 , 1 depth 1D convolution DWConv1D 3,1 , 3 Relu activation functions and 1 channel shuffle module Channel_Shuffle. The calculation process expression is as follows: [in x1 ,in x2 ]=Channel_Split(in x ,0.5) (10) in x2_1 =Conv1D 1，1 (BN(DWConv1D 3,1 (Relu(BN(Conv1D 1，1 (in x2 )))))) (12) out v1 =Channel_Shuffle(Concat(in x1_1 ,in x2_2 )) (14) Among them, x Indicates the input feature information of the shuffle_block1 module, in x1 , in x2 Respectively represent the channel information after segmentation according to the first channel and the second channel, in x1_1 Indicates in x1 After MSCA_1D operation and in x1 The result of multiplication, in x2_1 Indicates in x2 After Conv1D 1，1 , BN, Relu, DWConv1D 3，1 The feature information after operation, in x2_2 Indicates in x2_1 After BN and MSCA_1D operations and in x2_1 The result of the multiplication is then passed through Conv1D 1，1 , BN, Relu after the feature information, out y1 Indicates in x1_1 , in x2_2 The feature information after matrix connection and Channel_Shuffle operation, Channel_Split(tz, bl) means splitting the feature information tz into channels according to the ratio of bl, MSCA_1D(.) means the designed 1D multi-scale convolution attention module, BN(.) means batch normalization operation, Relu(.) means activation function operation, Conv1D 1,1 (.) indicates a 1D convolution operation with a kernel size of 1 and a stride of 1, DWConv1D 3，1 (.) indicates a depth 1D convolution operation with a kernel size of 3 and a stride of 1. Concat(.) indicates a matrix concatenation operation. Channel_Shufffle(.) indicates a channel shuffle operation. Represents matrix multiplication; S3.

2. The designed shuffle_block2 module consists of 6 batch normalization modules BN, 2 improved 1D multi-scale convolution attention MSCA_1D, and 4 1D convolution Conv1D 1，1 , 2 depth 1D convolution DWConv1D 3，2 , 3 Relu activation functions and 1 channel shuffle module Channel1_Shuffle, the calculation process expression is as follows: in 4_1 =DWConv1D 3,2 (in x4 ) (15) in x4_3 =Conv1D 1,1 (BN(DWConv1D 3,2 (Relu(BN(Conv1D 1，1 (in x4 )))))) (17) Among them, x4 Indicates the input feature information of the shuffle_block2 module, in x4_1 Indicates that after DWConv1D 3，2 The feature information after operation, in x4_2 Indicates in x4_1 After BN, MSCA_1D, Conv1D 1，1 The result of the operation is then combined with in x4_1 The feature information after multiplication and BN and Relu operation, in x4_3 Indicates in x4 After Conv1D 1，1 , BN, Relu, DWConv1D 3，2 、Conv1D 1,1 Feature information after operation, inx 4_4 Indicates in x4_3 After BN, MSCA_1D, Conv1D 1,1 The result of the operation is then in x4_3 After multiplication, the feature information after BN and Relu operation is obtained. y2 Indicates in x4_2 , in x4_4 The feature information after matrix connection and Channel_Shuffle operation, MSCA_1D(.) represents the 1D multi-scale convolution attention module designed by the present invention, BN(.) represents batch normalization operation, Relu(.) represents activation function operation, Conv1D 1,1 (.) indicates a 1D convolution operation with a kernel size of 1 and a stride of 1, DWConv1D 3,2 (.) represents a depth 1D convolution operation with a kernel size of 3 and a stride of 2. Concat(.) represents a matrix concatenation operation. Channel_Shuffle(.) represents a channel shuffle operation. Represents matrix multiplication; S3.

3. The designed 1D multi-scale convolutional attention module MSCA_1D consists of 4 1D convolutions of different scales and depths DWConv1D, 1 1D convolution Conv1D 1，1 The calculation process is as follows: in x51 =DWConv1D5(in x5 )(20) Among them, x5 Represents the input feature information of the 1D multi-scale convolutional attention module MSCA_1D, in x5_1 Indicates in x5 Feature information after DWConv1D5 operation, in x5_2 , in x5_3 , in x5_4 Respectively represent in x5_1 After DWConv1D7, DWConv1D 11 、DWConv1D 21 The feature information after operation, Indicates in x5_1 , in x5_2 , in x5_3 , in x5_4 The result after matrix addition and Conv1D1 operation is the same as in x5 Multiply feature information, DWConv1D q Represents the depth 1D convolution operation of q convolution kernel, q∈[1,5,7,11,21], Conv1D 1,1 1 indicates that the convolution kernel is 1 and the step size is 1, which is a convolution operation. Add indicates matrix addition. Represents matrix multiplication; S3.

4. The calculation process of channel shuffle Channel_Shuffle is as follows: S3.4.

1. Determine the number of feature map x groups and the number of channels C per group g The calculation formula is as follows: Among them, C represents the number of channels of the feature map x, G represents the number of groups to be divided, and C g Indicates the number of channels in each group; S3.4.2.reshape operation: Expand the channel dimension of the feature map x from one dimension to two dimensions, one is the number of groups G, and the other is the number of channels in each group C g ; If the batch size of the feature map is N and the spatial dimension is H×W, the reshape operation can be expressed as: x reshaped =x.view(N,G,C g ,H,W) (24) S3.4.3.transpose operation: The reshaped feature map is transposed in terms of the number of groups and the number of channels in each group, so that the channels in each group can be disrupted. The transposition operation can be expressed as: x transposed =x reshaped .transpose(1,2) (25) Among them, transpose(1,2) means exchanging the second dimension and the third dimension; S3.4.

4. flatte n operate: The transposed feature map X transposed , flatten back to the original number of channels C, and get the final channel shuffled feature map. The flatten operation is expressed as: x shuffled =x transposed .view(N,C,H,W) (26) Among them, view(N, C, H, W) restores the feature map to its original four-dimensional shape; S3.

5. The designed DSAT_1D_ShuffleNetV2 consists of 1 1D convolution Conv1D 3,2 , 1 1D convolution Conv1D 1，1 , 1 1-dimensional maximum pooling MaxPool 3,2 , 1 global average pooling GlobaPool, 3 designed shuffle_block2, 13 designed shuffl e block1, the calculation process expression is as follows: Stage0=MaxPool 3,2 (Conv1D 3,2 (input x )) (27) Stage 2_1 =Concat(interpolate(Stage1,Stage2.shape()),Stage2) (29) Stage 3_1 =Concat(interpolate(Stage 2_1 ,Stage3.shape()),Stage3) (30) output y =GlobaPool(Conv1D 1,1 (Stage 3_1 )) (31) Among them, input x Indicates the input feature information of the designed DSAT_1D_ShuffleNetV2, Stage0 represents input x After Conv1D 3，2 、MaxPool 3,2 The feature information after the operation. Stage1, Stage2, and Stage3 represent the feature information after Stage0 is operated by shuffle_block2 and shuffle_block1 respectively. 21 It indicates the characteristic information obtained by interpolating the result of Stage1 with the size of Stage2 and connecting it with the matrix of Stage2. 31 Represents Stage 21 After interpolation operation based on the Stage 3 size, the interpolation result is connected with the Stage 3 matrix to obtain the feature information, output y Represents Stage 3_1 After Conv1D 1，1 , feature information of GlobaPool operation, Conv1D 1，1 Represents a 1D convolution operation with a convolution kernel of 1 and a stride of 1, Conv1D 3,2 (.) indicates a 1D convolution operation with a kernel size of 3 and a stride of 2. MaxPool 3,2 represents a 1D maximum pooling operation with a convolution kernel of 3 and a stride of 2. shuffle_block2(, v) represents a shuffle_block2 module operation designed v times, v∈[1,1,1], shuffle_block1(, w) represents a shuffle_block1 module operation designed w times, w∈[3,7,3], interpolate(.) represents an interpolation operation, Concat(.) represents a matrix concatenation operation, and GlobaPool(.) represents a global average pooling operation.

5. The multi-time series water eutrophication assessment method based on satellite images according to claim 4 is characterized in that: The calculation formula of the dual-weight logarithmic root mean square error loss function Wwrmse designed in step S4 is as follows: Among them, i represents the i-th period, i.e., the first quarter, the second quarter, the third quarter, and the fourth quarter, n represents the total number of periods, j represents the j-th sample data, m represents the total number of samples, and W i represents the weight of the i-th period, wrmse represents the logarithmic root mean square error of the weight designed by the present invention, represents the forecast value for the i-th period, represents the true value of the i-th period, represents the average value of the true value in period i, It represents the maximum value of the true value in the i-th period, cov represents the covariance operation, and σ represents the standard deviation operation.

6. The multi-time series water eutrophication assessment method based on satellite images according to claim 5 is characterized in that: The specific implementation method of step S5 includes the following steps: The Shuffle_GRUnet model designed in S5.

1. consists of 1 input layer, 1 DSAT_1D_ShuffleNetV2 layer designed in step S3, 1 gated recurrent neural network layer, and 1 fully connected layer. Each layer contains n components; input layer X = {1st period, 2nd period, 3rd period, 4th period...nth period}, each component input parameter contains 9 steps 2 Calculate the spectral index (SPI b4,b1 、SPI b4,b2 、SPI b4,b3 、SPI b3,b4,b1 、SPI b3,b4,b1 、SPI b1,b2,b4,b3 、SPI b2,b3,b4,b1 、SPI b1,b3,b4,b2 、SPI b3,b4,b1,b2 ), input in the format of a three-dimensional tensor, so that DSAT-1D_ShuffleNetV2 and GRU can be used for feature extraction and model training. The three-dimensional tensor used is (number of samples × time step × number of features), the number of features is 9, and the spectral index stacking method is used to expand it to 224. The insufficient part is supplemented with 0. The calculation process expression is as follows: Among them, X i represents the input parameters of the i-th period, Indicates the output parameter of the i-th period, i.e., the comprehensive nutritional status index, DSAT_1D_ShuffleNetV2 i (.) indicates the DSAT_1D_ShuffleNetV2 module operation designed in step 3 for the i-th period, GRU i (.) indicates that the GRU module operation is performed in the i-th period, fc i (.) indicates that the i-th period has undergone a full connection operation; S5.

2. The referenced GRU consists of an update gate and a reset gate. First, the update gate Z is calculated. t and reset gate r t The value of the input variable x t and the hidden layer result h at the previous moment t-1 The concatenated matrix is ​​transformed by sigmoid nonlinearity and then input into the update gate; the reset gate value will act on h t-1 and the input variable x t The concatenation is nonlinearly transformed and then activated by the tanh function to obtain the state of the candidate hidden layer at the current moment. By 1-Z t times h t-1 Store the information of the previous moment through Z t times Record the information at the current moment, and add the results of the two as the state output h of the hidden layer at the current moment t Finally, the weight matrix of the output layer is applied to h t On the top, we get the current model output y t , the calculation formula is as follows: Z t =sigmoid(W Z [h t-1 ,x t ]) (36) r t =sigmoid(W r [h t-1 ,x t ]) (37) y t =sigmoid(W o h t ) (40) Among them, sigmoid represents the sigmoid activation function operation, tanh represents the tanh activation function operation, and Z t Indicates that the door state is updated at time t, r t Indicates that the door state is reset at time t, represents the candidate hidden door state at time t, h t represents the hidden layer state at time t, y t represents the model output at time t, x t represents the model input at time t, h t-1 represents the hidden layer state at the previous moment, x t , W Z represents the update gate weight matrix, W r Represents resetting the gate weight matrix, represents the candidate hidden layer weight matrix, W o Represents the output layer weight matrix.

7. The multi-time series water eutrophication assessment method based on satellite images according to claim 6 is characterized in that: The implementation method of step S6 is to input the target parameters prepared in step S1 and the input parameters prepared in step S2 into the multi-time series water eutrophication assessment model Shuffle_GRUnet designed in step S5 for model training, and use the Pearson correlation coefficient (r pcc ), mean root mean square error The model accuracy is evaluated to obtain the optimal multi-time series water eutrophication assessment model.

8. An electronic device, characterized in that: The method comprises a memory and a processor, wherein the memory stores a computer program, and the processor implements the steps of a multi-time series water eutrophication assessment method based on satellite images as described in any one of claims 1 to 7 when executing the computer program.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the multi-time series water eutrophication assessment method based on satellite images as described in any one of claims 1 to 7 is implemented.

Images