Comprehensive assessment method of water ecological environment based on remote sensing and spatial big data technology

Abstract

The present invention provides a comprehensive water ecological environment assessment method based on remote sensing and space big data technology, which belongs to the field of data analysis and solves the problem of low efficiency of water environment assessment. The method is specifically as follows: Step S1: obtaining river parameters; Step S2: obtaining standard rainfall and historical rainfall in a target area, and estimating the rainfall in the target area this month; Step S3: obtaining historical river parameters in the target area, constructing a relationship equation of river flow, vegetation, and sand content and a simulation equation of pollution diffusion; using the relationship equation and the simulation equation, analyzing river parameters, and obtaining a water environment analysis report; Step S4: continuously obtaining river parameters and updating the water environment analysis report; The present invention improves the efficiency of water environment assessment by acquiring, analyzing, and mathematically modeling relevant data of rivers in the target area, so that managers can better understand changes in the water ecological environment.

Description

Technical Field

[0001] The present invention discloses a comprehensive water ecological environment assessment method based on remote sensing and space big data technology, and relates to the field of data analysis. Background Art

[0002] The existing comprehensive assessment methods for water ecological environment have the following shortcomings:

[0003] Large data demand: Existing assessment methods usually require a large amount of monitoring data to support the assessment process, but most of this data has little or no impact on the water ecological environment. In addition, too much data will increase the burden on data processing engineering and analysis processes.

[0004] Lack of assessment methods: Most existing assessment methods use numerical comparison or threshold analysis methods, which can only analyze the water environment of a certain area in the short term and are not suitable for long-term water environment analysis. In addition, when using numerical comparison or threshold analysis methods, professionals are required to calculate the threshold, which affects the efficiency of the analysis.

[0005] Lack of changes in time and space: When conducting water environment analysis, most existing assessment methods assume that the water environment in the target area remains unchanged, and then fit the acquired data; this analysis method lacks changes in the water environment over time or space, is not real-time and has large errors. Summary of the invention

[0006] In view of the shortcomings of the prior art, the purpose of the present invention is to provide a comprehensive assessment method of water ecological environment based on remote sensing and spatial big data technology, aiming to solve the problem of low efficiency of water environment assessment.

[0007] In order to achieve the above object, the present invention is implemented through the following technical solutions: A comprehensive water ecological environment assessment method based on remote sensing and space big data technology includes:

[0008] Step S1: Obtain an aerial picture of the target area, and obtain the river area and vegetation area of ​​the target area based on the aerial picture; obtain the average flow, average sediment content and average pollution of the rivers in the target area as river parameters;

[0009] Step S2: Obtain the standard rainfall and historical rainfall of the target area, and determine whether there is regularity in the rainfall in the target area according to the standard rainfall; if there is, use the historical rainfall as a benchmark and use maximum likelihood estimation to estimate the rainfall in the target area this month; if there is no, use the historical rainfall and standard rainfall as a benchmark and use the state transfer matrix to estimate the rainfall in the target area this month;

[0010] Step S3: Obtain historical river parameters of the target area, and construct a relationship equation of river flow-vegetation-sediment content and a simulation equation of pollution diffusion based on the historical river parameters; use the relationship equation and the simulation equation to analyze the river parameters and obtain a water environment analysis report;

[0011] Step S4: Continuously obtain river parameters and update the water environment analysis report.

[0012] Furthermore, the specific steps of step S2 are as follows:

[0013] Step S21: Get the current month, recorded as m; take the last three months of the target area as standard months, that is, m3, m2 and m1 as standard months;

[0014] Get the daily rainfall in the standard month as the standard rainfall, recorded as rc1, rc2~rc dm ; where dm is the total number of days in a standard month; rc1 is the rainfall on the first day; rc2 is the rainfall on the second day; and so on. dm represents the rainfall on day dm;

[0015] Calculate rc1~rc dm The average value is recorded as arc; calculate rc1~rc dm The variance of 2 ;

[0016] Step S22: taking m3, m2, m1 and m as reference months;

[0017] Get the daily rainfall in the reference month in the past three years as the historical rainfall, and get:

[0018] In the first year, refer to the daily rainfall in the month: rⅠ1, rⅠ2~rⅠ dm ~rⅠ dm＋d ;

[0019] Among them, rⅠ1 represents the rainfall on the first day of the reference month in the first year; rⅠ2 represents the rainfall on the second day of the reference month in the first year; and so on. dm represents the rainfall on the dmth day of the reference month in the first year; rⅠ dm＋d It represents the rainfall on the (dm+d)th day in the reference month of the first year;

[0020] In the second year, refer to the daily rainfall in the month: rⅡ1, rⅡ2~rⅡ dm ~rⅡ dm＋d ;

[0021] Among them, rⅡ1 represents the rainfall on the first day of the reference month in the second year; rⅡ2 represents the rainfall on the second day of the reference month in the second year; and so on. dmrepresents the rainfall on the dmth day of the reference month in the second year; rⅡ dm＋d It represents the rainfall on the (dm+d)th day in the reference month of the second year;

[0022] In the third year, refer to the daily rainfall in the month: rⅢ1, rⅢ2~rⅢ dm ~rⅢ dm＋d ;

[0023] Among them, rⅢ1 represents the rainfall on the first day of the reference month in the third year; rⅢ2 represents the rainfall on the second day of the reference month in the third year; and so on. dm represents the rainfall on the dmth day of the reference month in the third year; rⅢ dm＋d It represents the rainfall on the (dm+d)th day in the reference month of the third year;

[0024] Where d represents the number of days in month m;

[0025] Step S23: define the judgment condition of the change point:

[0026] Assume {rc1~rc dm} in a set of breakpoints: rc r and rc p , i.e. rc r and rc p ∈{rc1~rc dm};

[0027] If rc r and rc p If the following conditions are met, it means that the rainfall in the target area conforms to the normal distribution, (rc r , rc p ) is rc1~rc dm A set of change points, rc r ~rc p It is rc r and rc p Regular section; rc r and rc p The conditions met are as follows:

[0028] rc r ≠rc p ;

[0029] [N (arc （r→p） , S 2 ), rc r ~rc p ];

[0030] [(rc1≤rc r )∩(rc p ≤rc dm )];

[0031] Among them, arc （r→p） Indicates rc1~rc dm Medium r ~rc p The average value of .

[0032] Furthermore, the subsequent steps of step S23 are as follows:

[0033] Step S24: record the rainfall in the target area this month as ee;

[0034] According to the judgment condition of the change point, judge rc1~rc dm Is there a change point in the

[0035] If yes, it means that there is regularity in the rainfall in the target area, and the rainfall in the target area this month is estimated using maximum likelihood estimation, and the process goes to step S25;

[0036] If it does not exist, it means that the rainfall in the target area is irregular; based on the historical rainfall and standard rainfall, the state transfer matrix is ​​used to estimate the rainfall in the target area this month, and the process goes to step S26;

[0037] Step S25: Count rc1 to rc dm The number of change points in is recorded as he; the total number of regular segments corresponding to each group of change points is counted as zo, zo = he / 2;

[0038] Calculate the mean and variance of all data in the first regular segment and record them as ax(1) and S 2 (1);

[0039] Calculate the mean and variance of all data in the second regular segment and record them as ax(2) and S 2 (2);

[0040] Similarly, calculate the mean and variance of all data in the z0th regular segment and record them as ax(z0) and S 2 (zo);

[0041] Use maximum likelihood estimation to estimate the rainfall in the target area and determine the value of ee;

[0042] Step S26: Based on the historical rainfall and standard rainfall, the state transfer matrix is ​​used to estimate the rainfall in the target area this month and determine the value of ee;

[0043] Step S27: Use ee as the input of the relationship equation of river flow-vegetation-sediment content and the simulation equation of pollution diffusion, and enter step S3.

[0044] Furthermore, the specific steps of step S25 are as follows:

[0045] Step S251: Set the number of the regular segment to j, and the value range of j is 1 to zo;

[0046] Construct rc1~rc dm The primary mean periodic function is denoted as ax(j), and the primary variance periodic function is denoted as S 2 (j);

[0047] Primary mean periodic function: ax(j) = Ao×cos(ω×j) + Bo×sin(ω×j) + Co;

[0048] Primary variance periodic function: S 2 (j)=Do×cos(ω×j)+Eo×sin(ω×j)+Fo;

[0049] Among them, Ao, Bo and Co are the parameters to be estimated of the primary mean periodic function; Do, Eo and Fo are the parameters to be estimated of the primary variance periodic function; ω represents the angular frequency, ω = (2π) / dm;

[0050] Step S252: Combine ax(1) to ax(zo) and S 2 (1) ~S 2 Substitute (zo) into the primary mean periodic function and the primary variance periodic function in turn, calculate the corresponding values ​​of Ao, Bo and Co, recorded as A, B and C; calculate the values ​​of Do, Eo and Fo, recorded as D, E and F;

[0051] The mean periodic function is recorded as ax-f(j), ax-f(j)=A×cos(ω×j)+B×sin(ω×j)+C;

[0052] The variance periodic function is denoted as S 2 -f(j),S 2 -f(j)=D×cos(ω×j)+E×sin(ω×j)+F;

[0053] Step S253: let the date be dd, the time variable be t, t=dd+dm;

[0054] Using Fourier series, ax-f(j) and S 2 -f(j) is fitted and the representation range of j is broadened;

[0055] Step S254: read the maximum rainfall value in the 1st to zoth regular sections and record it as rc(max);

[0056] Get the total rainfall in month m1, recorded as ry;

[0057] The function expression of rainfall in the target area and time is recorded as rr(t), and the mathematical expression of rr(t) is:

[0058] ; Wherein, exp represents the exponential function with e as the base;

[0059] Step S255: Get the current date and record it as dt;

[0060] Substitute dt into t in rr(t) and calculate to get rr(dt);

[0061] Let rr(dt) be the value of ee.

[0062] Furthermore, the specific steps of step S253 are as follows:

[0063] Step S2531: replace the variable j in ax-f(j) with the variable t to obtain the mixed function ax-f(t), and fit ax-f(j);

[0064] The general formula of the Fourier expansion of ax-f(t) is denoted as fa(t), and the mathematical expression of fa(t) is:

[0065] ;

[0066] Among them, i represents the time corresponding to the jth regular segment;

[0067] ao represents the coefficient of the zero-order term, an represents the coefficient of the cosine term, and bn represents the coefficient of the sine term; the values ​​of ao, an, and bn are all 1;

[0068] Step S2532: taking ax-f(j) as the integrand, calculate the integrals of ao, an and bn in sequence to obtain aao, aan and abn;

[0069] The calculation formula of aao is as follows:

[0070] ;

[0071] The calculation formula of aan is as follows:

[0072] ;

[0073] The calculation formula of abn is as follows:

[0074] ;

[0075] Use aao, aan and abn to replace ao, an and bn in fa(t) to obtain the time-mean equation of fa(t), denoted as a(t);

[0076] Step S2533: Replace S with variable t 2 - The variable j in f(j) gives the mixed function S2 -f(t), repeat the same steps of fitting xf(t) to obtain the time-variance equation, denoted as S 2 (t).

[0077] Furthermore, the specific steps of step S26 are as follows:

[0078] Step S261: Calculate rⅠ1~rⅠ dm The average value is recorded as rⅠa; calculate rⅠ dm+1 ~rⅠ dm＋d The average value of is denoted as rab;

[0079] Calculate rⅡ1~rⅡ dm The average value is recorded as rⅡa; calculate rⅡ dm+1 ~rⅡ dm＋d The average value is denoted as rbb;

[0080] Calculate rⅢ1~rⅢ dm The average value is recorded as rⅢa; calculate rⅢ dm+1 ~rⅢ dm＋d The average value of is denoted as rcb;

[0081] Step S262: Record the rainfall change rate from the third year to the second year as up （3-2） , up （3-2） The calculation formula is as follows:

[0082] up （3-2） =1-(rⅢa / rⅡa);

[0083] The rate of change of rainfall from the second year to the first year is recorded as up （2-1） ,up （2-1） The calculation formula is as follows:

[0084] up （2-1） =1-(rⅡa / rⅠa);

[0085] The rate of change of rainfall from the first year to this year is recorded as up （0） , up （0） The calculation formula is as follows:

[0086] up （0） =1-(rⅠa / arc);

[0087] Step S263: Up （0） 、up （3-2） and up （2-1） , as well as rab, rbb and rcb; construct a state transfer matrix, denoted as matrix B;

[0088] The mathematical expression of matrix B is as follows:

[0089] ; Where n represents the number of iterations of matrix B, × represents matrix multiplication; n is a natural number, and its initial value is 1; each time matrix B is iterated, n is increased by 1 on the original basis;

[0090] Use the NumPy library function to iterate matrix B and obtain the steady-state transfer matrix, which is recorded as matrix C;

[0091] Define a 2×1 order receiving matrix {1 / 2,1 / 2}, denoted as matrix D;

[0092] Calculate the value of matrix C multiplied by matrix D as the value of ee.

[0093] Furthermore, the specific steps of step S3 are as follows:

[0094] Step S31: Obtain the average river flow, average vegetation coverage area and average river sediment content of the target area in the past five years;

[0095] The average river flow in the past five years is recorded as ri1, ri2, ri3, ri4, and ri5;

[0096] The average vegetation cover area in the past five years is recorded as ve1, ve2, ve3, ve4, and ve5;

[0097] The average river sediment content in the past five years is recorded as tu1, tu2, tu3, tu4, and tu5;

[0098] Among them, ri1, ve1 and tu1 represent the average river flow, average vegetation coverage area and average river sediment content in the past year, respectively; ri2, ve2 and tu2 represent the average river flow, average vegetation coverage area and average river sediment content in the past two years, respectively; and so on, ri5, ve5 and tu5 represent the average river flow, average vegetation coverage area and average river sediment content in the past five years, respectively;

[0099] Step S32: Calculate the average value of ri1-ri5 and record it as Ari; calculate the average value of ve1-ve5 and record it as Ave; calculate the average value of tu1-tu5 and record it as Atu;

[0100] The standard deviation of ri1~ri5 is recorded as Sri; the standard deviation of ve1~ve5 is recorded as Sve; the standard deviation of tu1~tu5 is recorded as Stu;

[0101] Calculate the coefficient of river flow on vegetation coverage, denoted as b1; the calculation formula of b1 is as follows:

[0102] ; Among them, ri i represents the average river flow in the past i years; vei represents the average vegetation coverage area in the past i years;

[0103] Calculate the relationship coefficient between vegetation coverage area and river sediment content, denoted as b2; the calculation formula of b2 is as follows:

[0104] ; Among them, tu i represents the average river sediment content in the past i years;

[0105] Step S33: taking the river flow as the independent variable x;

[0106] Define the relationship between river and vegetation coverage area, recorded as equation p(x); the mathematical expression of p(x) is:

[0107] p(x)=b1×x+sc0; where sc0 represents the parameter to be estimated of equation p(x), and the initial value of sc0 is 0;

[0108] Substitute ri1~ri5 and ve1~ve5 into p(x), calculate the value of sc0, and get sc; the final expression of p(x) is:

[0109] p(x)=b1×x+sc;

[0110] Step S34: Taking the vegetation coverage area as the independent variable z, define the relationship equation between river flow, vegetation and sediment content, which is recorded as equation g(x); the mathematical expression of g(x) is:

[0111] g(x)=p(x)+b2×z+cs0; where cs0 represents the parameter to be estimated of equation g(x), and the initial value of cs0 is 0;

[0112] Substitute ri1~ri5, ve1~ve5 and tu1~tu5 into g(x), calculate the value of cs0, and get cs; the final expression of g(x) is:

[0113] g(x)=p(x)+b2×z+cs.

[0114] Furthermore, the subsequent steps of step S34 are as follows:

[0115] Step S35: Obtain the area of ​​the target region, recorded as SC; record the square root of SC as J, and perform continuous square root extraction on SC until J∈(0,10];

[0116] Based on the aerial photos of the target area, GIS is used to extract and construct a grid model of the river area;

[0117] Taking the northeast end of the target area as the origin and J×J as the division size, the grid model is divided and the river area is mapped to a two-dimensional coordinate system; the points in the two-dimensional coordinate system represent the location of the river;

[0118] According to the two-dimensional coordinate system, the simulation equation of river pollution diffusion is constructed;

[0119] Step S351: Assume that there is pollution at (p, q) in the river, and the pollutant concentration at (p, q) is C (p, q);

[0120] Assume that the pollution diffusion coefficient in the horizontal axis direction is Dx, the pollution diffusion coefficient in the vertical axis direction is Dy, and the attenuation coefficient of the pollutant is kk;

[0121] The relationship between Dx and Dy satisfies:

[0122] Dy=[(Dx×r H +Dx×r B ) / 2]×(g×r H ×H J ) 1 / 2 ; Among them, r H represents the river depth, r B Indicates the width of the river, H J It represents the hydraulic gradient of the river;

[0123] Step S352: define the simulation equation for river pollution diffusion, denoted as t(p,q); the mathematical expression of t(p,q) is:

[0124] t(p,q)={Dx×[C(p,q)] 2 / p 2}+{Dy×[C(p,q)] 2 / q 2}-{kk×C(p,q)};

[0125] Step S36: The river area in the river parameters is recorded as RS, the vegetation area is recorded as VS, the average flow is recorded as rd, the average sediment content is recorded as tb, and the average pollution amount is recorded as ap;

[0126] Calculate the increase in river flow, denoted as Δri, Δri = ee / RS;

[0127] Substitute (rd + Δri) and VS into equations p(x) and g(x), respectively, to calculate the rated sediment content of the target area, denoted as ttu;

[0128] Compare the size of ttu and tb;

[0129] If ttu≤tb, it means that soil erosion in the target area is relatively light;

[0130] If ttu>tb, it means that soil erosion in the target area is serious;

[0131] Step S37: Obtain the location of the sewage outlet in the target area and record it as (m, n), and the concentration of the pollutant and record it as po;

[0132] Replace the parameter p with m, replace the parameter q with n, and replace the parameter C(p,q) with po. Substitute (m,n) and po into the equation t(p,q). Use the finite difference algorithm to iterate the equation t(p,q) to obtain the convergence value of the equation t(p,q), which is recorded as tuu.

[0133] Compare the size of tuu and ttu;

[0134] If tuu<ttu, it means that the water pollution in the target area is relatively light;

[0135] If tuu ≥ ttu, ​​it means that the water pollution in the target area is serious;

[0136] Step S38: Summarize the analysis results of steps S36 to S37 as a water environment analysis report and feed it back to the user.

[0137] Compared with the prior art, the present invention has the following beneficial effects:

[0138] Objectivity and accuracy: The present invention adopts the analysis method of data fitting and mathematical modeling to analyze the water environment of the target area; while ensuring the integrity and validity of the data, it analyzes the impact of various indicators in the rivers in the target area on the water environment from a more objective perspective, avoids subjective arbitrariness, and improves the accuracy of the assessment.

[0139] Comprehensiveness and holisticness: The present invention analyzes the water environment from more than three representative indicators, namely rainfall, soil erosion and water pollution; it enables users to have a more comprehensive understanding of the health status of the water ecological environment, and the calculation method and evaluation process of the present invention are easy to implement in the computer field and easy for relevant technical personnel to operate, and has strong practicality.

[0140] Strong dynamic monitoring capability: The present invention designs mathematical equations for the changes of water environment over time and water pollution over space. It has the ability to analyze time continuity and monitor spatial changes, and can monitor the changing trends of the water ecological environment in real time. Combined with big data technology, the present invention can also realize the ability to conduct long-term dynamic monitoring of the water ecological environment. BRIEF DESCRIPTION OF THE DRAWINGS

[0141] Other features, objects and advantages of the present invention will become more apparent from the detailed description of non-limiting embodiments made with reference to the following drawings:

[0142] Figure 1 It is a schematic diagram of the method of the present invention;

[0143] Figure 2 It is a schematic diagram of the change point division of the present invention;

[0144] Figure 3 Schematic diagram of the river area of ​​the present invention. DETAILED DESCRIPTION

[0145] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0146] See also Figure 1 , the comprehensive assessment methods of water ecological environment based on remote sensing and space big data technology include:

[0147] Step S1: Obtain an aerial picture of the target area, and obtain the river area and vegetation area of ​​the target area based on the aerial picture; obtain the average flow, average sediment content and average pollution of the rivers in the target area as river parameters;

[0148] It should be noted that the “target area” in the present invention refers to an area where a comprehensive assessment of the water ecological environment is conducted using the present invention (a comprehensive assessment method of the water ecological environment based on remote sensing and spatial big data technology).

[0149] Step S2: Obtain the standard rainfall and historical rainfall of the target area, and determine whether there is regularity in the rainfall in the target area according to the standard rainfall; if there is, use the historical rainfall as a benchmark and use maximum likelihood estimation to estimate the rainfall in the target area this month; if there is no, use the historical rainfall and standard rainfall as a benchmark and use the state transfer matrix to estimate the rainfall in the target area this month;

[0150] The specific steps of step S2 are as follows:

[0151] Step S21: Get the current month, recorded as m; take the last three months of the target area as standard months, that is, m3, m2 and m1 as standard months;

[0152] For example, if the current month is October, then m3 is July, then m2 is August, then m1 is September; if the current month is January, then m3 is October of last year, then m2 is November of last year, then m1 is December of last year;

[0153] Get the daily rainfall in the standard month (in the target area) as the standard rainfall, recorded as rc1, rc2~rc dm ; where dm is the total number of days in a standard month (ignoring leap years); rc1 is the rainfall on the first day (in a standard month); rc2 is the rainfall on the second day (in a standard month); and so on. dmIndicates the rainfall on the dmth day (in a standard month);

[0154] Calculate rc1~rc dm The average value is recorded as arc; calculate rc1~rc dm The variance of 2 ;

[0155] Step S22: taking m3, m2, m1 and m as reference months;

[0156] Get the daily rainfall in the reference month in the past three years as the historical rainfall, and get:

[0157] In the first year, refer to the daily rainfall in the month: rⅠ1, rⅠ2~rⅠ dm ~rⅠ dm＋d ; (where rⅠ1 represents the rainfall on the first day of the reference month in the first year; rⅠ2 represents the rainfall on the second day of the reference month in the first year; and so on, rⅠ dm represents the rainfall on the dmth day of the reference month in the first year; rⅠ dm＋d It indicates the rainfall on the (dm+d)th day in the reference month of the first year;)

[0158] In the second year, refer to the daily rainfall in the month: rⅡ1, rⅡ2~rⅡ dm ~rⅡ dm＋d ; (where rⅡ1 represents the rainfall on the first day of the reference month in the second year; rⅡ2 represents the rainfall on the second day of the reference month in the second year; and so on, rⅡ dm represents the rainfall on the dmth day of the reference month in the second year; rⅡ dm＋d It indicates the rainfall on the (dm+d)th day in the reference month of the second year;)

[0159] In the third year, refer to the daily rainfall in the month: rⅢ1, rⅢ2~rⅢ dm ~rⅢ dm＋d ; (where rⅢ1 represents the rainfall on the first day of the reference month in the third year; rⅢ2 represents the rainfall on the second day of the reference month in the third year; and so on, rⅢ dm represents the rainfall on the dmth day of the reference month in the third year; rⅢ dm＋d It represents the rainfall on the (dm+d)th day in the reference month of the third year;)

[0160] Where d represents the number of days in month m (ignoring leap years);

[0161] Step S23: define the judgment condition of the change point:

[0162] Assume {rc1~rc dm} in a set of breakpoints: rc r and rc p , i.e. rcr and rc p ∈{rc1~rc dm};

[0163] If rc r and rc p If the following conditions are met, it means that the rainfall in the target area (partially or entirely) conforms to the normal distribution, (rc r , rc p ) is rc1~rc dm A set of change points, rc r ~rc p It is rc r and rc p Regular section; rc r and rc p The conditions met are as follows:

[0164] rc r ≠rc p ;

[0165] [N (arc （r→p） , S 2 ), rc r ~rc p ];

[0166] [(rc1≤rc r )∩(rc p ≤rc dm )];

[0167] Among them, arc （r→p） Indicates rc1~rc dm Medium r ~rc p The average value of

[0168] Step S24: record the rainfall in the target area this month as ee;

[0169] According to the judgment condition of the change point, judge rc1~rc dm Is there a change point in the

[0170] If yes, it means that there is regularity in the rainfall in the target area, and the rainfall in the target area this month is estimated using maximum likelihood estimation, and the process goes to step S25;

[0171] If it does not exist, it means that the rainfall in the target area is irregular; based on the historical rainfall and standard rainfall, the state transfer matrix is ​​used to estimate the rainfall in the target area this month, and the process goes to step S26;

[0172] Step S25: Count rc1 to rc dmThe number of change points in is recorded as he; the total number of regular segments corresponding to each group of change points is counted as zo, zo = he / 2;

[0173] Calculate the mean and variance of all data in the first regular segment and record them as ax(1) and S 2 (1);

[0174] Calculate the mean and variance of all data in the second regular segment and record them as ax(2) and S 2 (2);

[0175] Similarly, calculate the mean and variance of all data in the z0th regular segment and record them as ax(z0) and S 2 (zo);

[0176] Please note that Figure 2 , the "change points" in the present invention can be counted twice; when the regular segments corresponding to two groups of change points are connected, the middle change point is counted twice; Figure 2 An example illustrates:

[0177] Figure 2 "xa, xb and xc" are the change points of the data "x1~xn", and the number of "xa, xb and xc" is 3;

[0178] But in the regular section （xa-xb） In the regular section, xb is the change point of the group (xa, xb); （xb-xc） In the equation, xb is the change point of the group (xb, xc), so xb is counted twice. Therefore, the number of change points in “x1~xn” is (1+1)+(1+1)=4.

[0179] The changing points xa and xb of the (xa, xb) group + the changing points xb and xc of the (xb, xc) group.

[0180] Use maximum likelihood estimation to estimate rainfall in the target area;

[0181] Step S251: Set the number of the regular segment to j, and the value range of j is 1 to zo;

[0182] Construct rc1~rc dm The primary mean periodic function is denoted as ax(j), and the primary variance periodic function is denoted as S 2 (j);

[0183] Primary mean periodic function: ax(j) = Ao×cos(ω×j) + Bo×sin(ω×j) + Co;

[0184] Primary variance periodic function: S 2 (j)=Do×cos(ω×j)+Eo×sin(ω×j)+Fo;

[0185] Among them, Ao, Bo and Co are the parameters to be estimated of the primary mean periodic function; Do, Eo and Fo are the parameters to be estimated of the primary variance periodic function; ω represents the angular frequency, ω = (2π) / dm;

[0186] Step S252: Combine ax(1) to ax(zo) and S 2 (1) ~S 2 Substitute (zo) into the primary mean periodic function and the primary variance periodic function in turn, calculate the corresponding values ​​of Ao, Bo and Co, recorded as A, B and C; calculate the values ​​of Do, Eo and Fo, recorded as D, E and F;

[0187] The mean periodic function is recorded as ax-f(j), ax-f(j)=A×cos(ω×j)+B×sin(ω×j)+C;

[0188] The variance periodic function is denoted as S 2 -f(j),S 2 -f(j)=D×cos(ω×j)+E×sin(ω×j)+F;

[0189] Step S253: let the date be dd, the time variable be t, t=dd+dm;

[0190] Using Fourier series, ax-f(j) and S 2 -f(j) is fitted and the representation range of j is broadened;

[0191] Step S2531: replace the variable j in ax-f(j) with the variable t to obtain the mixed function ax-f(t), and fit ax-f(j);

[0192] The general formula of the Fourier expansion of ax-f(t) is denoted as fa(t), and the mathematical expression of fa(t) is:

[0193] ;

[0194] Where i represents the time corresponding to the jth regular segment (unit: day);

[0195] ao represents the coefficient of the zero-order term, an represents the coefficient of the cosine term, and bn represents the coefficient of the sine term; the values ​​of ao, an, and bn are all 1;

[0196] Step S2532: taking ax-f(j) as the integrand, calculate the integrals of ao, an and bn in sequence to obtain aao, aan and abn;

[0197] The calculation formula of aao is as follows:

[0198] ;

[0199] The calculation formula of aan is as follows:

[0200] ;

[0201] The calculation formula of abn is as follows:

[0202] ;

[0203] Use aao, aan and abn to replace ao, an and bn in fa(t) to obtain the time-mean equation of fa(t), denoted as a(t);

[0204] Step S2533: Replace S with variable t 2 - The variable j in f(j) gives the mixed function S 2 -f(t), repeat the same steps of fitting xf(t) (i.e., steps S2531 to S2532) to obtain the time-variance equation, denoted as S 2 (t);

[0205] Step S254: read the maximum rainfall value in the 1st to zoth regular sections and record it as rc(max);

[0206] Get the total rainfall in month m1, recorded as ry;

[0207] The function expression of rainfall in the target area and time is recorded as rr(t), and the mathematical expression of rr(t) is:

[0208] ; Where exp represents an exponential function with e as the base, (e represents a natural parameter);

[0209] Step S255: Get the current date and record it as dt;

[0210] Substitute dt into t in rr(t) and calculate to get rr(dt);

[0211] Let rr(dt) be the value of ee;

[0212] Step S26: using the state transfer matrix to estimate the rainfall in the target area this month based on the historical rainfall and the standard rainfall;

[0213] Step S261: Calculate rⅠ1~rⅠ dm The average value is recorded as rⅠa; calculate rⅠ dm+1 ~rⅠ dm＋d The average value of is denoted as rab;

[0214] Calculate rⅡ1~rⅡ dmThe average value is recorded as rⅡa; calculate rⅡ dm+1 ~rⅡ dm＋d The average value is denoted as rbb;

[0215] Calculate rⅢ1~rⅢ dm The average value is recorded as rⅢa; calculate rⅢ dm+1 ~rⅢ dm＋d The average value of is denoted as rcb;

[0216] Step S262: Record the rainfall change rate from the third year to the second year as up （3-2） , up （3-2） The calculation formula is as follows:

[0217] up （3-2） =1-(rⅢa / rⅡa);

[0218] The rate of change of rainfall from the second year to the first year is recorded as up （2-1） , up （2-1） The calculation formula is as follows:

[0219] up （2-1） =1-(rⅡa / rⅠa);

[0220] The rate of change of rainfall from the first year to this year is recorded as up （0） , up （0） The calculation formula is as follows:

[0221] up （0） =1-(rⅠa / arc);

[0222] Step S263: Up （0） 、up （3-2） and up （2-1） , as well as rab, rbb and rcb; construct a state transfer matrix, denoted as matrix B;

[0223] The mathematical expression of matrix B is as follows:

[0224] ; Where n represents the number of iterations of matrix B, × represents matrix multiplication; n is a natural number, and its initial value is 1; each time matrix B is iterated, n is increased by 1 on the original basis;

[0225] Use the NumPy library function to iterate matrix B and obtain the steady-state transfer matrix, which is recorded as matrix C;

[0226] Define a 2×1 order receiving matrix {1 / 2,1 / 2}, denoted as matrix D;

[0227] Calculate the value of matrix C multiplied by matrix D as the value of ee.

[0228] Step S3: Obtain historical river parameters of the target area, and construct a relationship equation of river flow-vegetation-sediment content and a simulation equation of pollution diffusion based on the historical river parameters; use the relationship equation and the simulation equation to analyze the river parameters and obtain a water environment analysis report (for the target area);

[0229] The specific steps of step S3 are as follows:

[0230] Step S31: Obtain the average river flow, average vegetation coverage area and average river sediment content of the target area in the past five years;

[0231] The average river flow in the past five years is recorded as ri1, ri2, ri3, ri4, and ri5;

[0232] The average vegetation cover area in the past five years is recorded as ve1, ve2, ve3, ve4, and ve5;

[0233] The average river sediment content in the past five years is recorded as tu1, tu2, tu3, tu4, and tu5;

[0234] Among them, ri1, ve1 and tu1 respectively represent the average river flow, average vegetation coverage area and average river sediment content in the target area in the past year; ri2, ve2 and tu2 respectively represent the average river flow, average vegetation coverage area and average river sediment content in the target area in the past two years; and so on, ri5, ve5 and tu5 respectively represent the average river flow, average vegetation coverage area and average river sediment content in the target area in the past five years;

[0235] Step S32: Calculate the average value of ri1-ri5 and record it as Ari; calculate the average value of ve1-ve5 and record it as Ave; calculate the average value of tu1-tu5 and record it as Atu;

[0236] The standard deviation of ri1~ri5 is recorded as Sri; the standard deviation of ve1~ve5 is recorded as Sve; the standard deviation of tu1~tu5 is recorded as Stu;

[0237] Calculate the coefficient of river flow on vegetation coverage, denoted as b1; the calculation formula of b1 is as follows:

[0238] ; Among them, ri i represents the average river flow in the past i years; ve i represents the average vegetation coverage area in the past i years;

[0239] Calculate the relationship coefficient between vegetation coverage area and river sediment content, denoted as b2; the calculation formula of b2 is as follows:

[0240] ; Among them, tu irepresents the average river sediment content in the past i years;

[0241] Step S33: taking the river flow as the independent variable x;

[0242] Define the relationship between river and vegetation coverage area, recorded as equation p(x); the mathematical expression of p(x) is:

[0243] p(x)=b1×x+sc0; where sc0 represents the parameter to be estimated of equation p(x), and the initial value of sc0 is 0;

[0244] Substitute ri1~ri5 and ve1~ve5 into p(x), calculate the value of sc0, and get sc; the final expression of p(x) is:

[0245] p(x)=b1×x+sc;

[0246] Step S34: Taking the vegetation coverage area as the independent variable z, define the relationship equation between river flow, vegetation and sediment content, which is recorded as equation g(x); the mathematical expression of g(x) is:

[0247] g(x)=p(x)+b2×z+cs0; where cs0 represents the parameter to be estimated of equation g(x), and the initial value of cs0 is 0;

[0248] Substitute ri1~ri5, ve1~ve5 and tu1~tu5 into g(x), calculate the value of cs0, and get cs; the final expression of g(x) is:

[0249] g(x)=p(x)+b2×z+cs;

[0250] Step S35: Obtain the area of ​​the target region, recorded as SC; record the square root of SC as J, and perform continuous square root extraction on SC until J∈(0,10];

[0251] Based on the aerial photos of the target area, use GIS to extract and construct a grid model of the river area (in the target area);

[0252] Taking the northeast end of the target area as the origin and J×J as the division size, the grid model is divided and the river area (of the target area) is mapped into a two-dimensional coordinate system; the points in the two-dimensional coordinate system represent the locations of the rivers (of the target area);

[0253] See also Figure 3 , based on the two-dimensional coordinate system, construct the simulation equation of river pollution diffusion;

[0254] Step S351: Assume that there is pollution at (p, q) in the river, and the pollutant concentration at (p, q) is C (p, q);

[0255] Assume that the pollution diffusion coefficient in the horizontal axis direction (i.e., west direction) is Dx, the pollution diffusion coefficient in the vertical axis direction (i.e., south direction) is Dy, and the attenuation coefficient of the pollutant is kk;

[0256] The relationship between Dx and Dy satisfies:

[0257] Dy=[(Dx×r H +Dx×r B ) / 2]×(g×r H ×H J ) 1 / 2 ; Among them, r H represents the river depth, r B Indicates the width of the river, H J It represents the hydraulic gradient of the river;

[0258] It should be noted that the "pollution diffusion coefficient in the horizontal direction (i.e., the west direction) is Dx" represents the coefficient of spontaneous diffusion of pollutants in a static water environment, which is related to the nature of the pollutants; while the "pollution diffusion coefficient in the vertical direction (i.e., the south direction) Dy" represents the coefficient of passive diffusion of pollutants with the flow of the river in a dynamic water environment, so the value of Dy is affected by Dx;

[0259] Step S352: define the simulation equation for river pollution diffusion, denoted as t(p,q); the mathematical expression of t(p,q) is:

[0260] t(p,q)={Dx×[C(p,q)] 2 / p 2}+{Dy×[C(p,q)] 2 / q 2}-{kk×C(p,q)};

[0261] Step S36: The river area in the river parameters is recorded as RS, the vegetation area is recorded as VS, the average flow is recorded as rd, the average sediment content is recorded as tb, and the average pollution amount is recorded as ap;

[0262] Calculate the increase in river flow, denoted as Δri, Δri = ee / RS;

[0263] Substitute (rd + Δri) and VS into equations p(x) and g(x), respectively, to calculate the rated sediment content of the target area, denoted as ttu;

[0264] Compare the size of ttu and tb;

[0265] If ttu≤tb, it means that soil erosion in the target area is relatively light;

[0266] If ttu>tb, it means that soil erosion in the target area is serious;

[0267] Step S37: Obtain the location of the sewage outlet in the target area and record it as (m, n), and the concentration of the pollutant and record it as po;

[0268] Replace the parameter p with m, replace the parameter q with n, and replace the parameter C(p,q) with po. Substitute (m,n) and po into the equation t(p,q). Use the finite difference algorithm to iterate the equation t(p,q) to obtain the convergence value of the equation t(p,q), which is recorded as tuu.

[0269] Compare the size of tuu and ttu;

[0270] If tuu<ttu, it means that the water pollution in the target area is relatively light;

[0271] If tuu ≥ ttu, ​​it means that the water pollution in the target area is serious;

[0272] Step S38: Summarize the analysis results of steps S36 to S37 as a water environment analysis report and feed it back to the user.

[0273] Step S4: Continuously obtain river parameters and update the water environment analysis report.

[0274] The above formulas are all dimensionless and numerical calculations. The formula is a formula for the most recent real situation obtained by collecting a large amount of data and performing software simulation. The preset parameters in the formula are set by technical personnel in this field according to actual conditions. For example, if there are weight coefficients and proportional coefficients, their set sizes are to quantify each parameter to obtain a specific value, which is convenient for subsequent comparison. Regarding the size of the weight coefficient and the proportional coefficient, it is sufficient as long as it does not affect the proportional relationship between the parameter and the quantized value.

[0275] Finally, it should be noted that the above-described embodiments are only specific implementations of the present invention, which are used to illustrate the technical solutions of the present invention, rather than to limit them. The protection scope of the present invention is not limited thereto. Although the present invention is described in detail with reference to the above-described embodiments, ordinary technicians in the field should understand that any technician familiar with the technical field can still modify the technical solutions recorded in the above-described embodiments within the technical scope disclosed by the present invention, or can easily think of changes, or make equivalent replacements for some of the technical features therein; and these modifications, changes or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be based on the protection scope of the claims.

Claims

1. A comprehensive water ecological environment assessment method based on remote sensing and spatial big data technology, characterized in that: The method comprises: Step S1: Obtain an aerial picture of the target area, and obtain the river area and vegetation area of ​​the target area based on the aerial picture; obtain the average flow, average sediment content and average pollution of the rivers in the target area as river parameters; Step S2: Obtain the standard rainfall and historical rainfall of the target area, and determine whether there is regularity in the rainfall in the target area according to the standard rainfall; if there is, use the historical rainfall as a benchmark, use the maximum likelihood estimation to estimate the rainfall in the target area this month, and determine the value of ee; if there is no, use the historical rainfall and standard rainfall as a benchmark, use the state transfer matrix to estimate the rainfall in the target area this month, and determine the value of ee; Step S3: Obtain historical river parameters of the target area, and construct a relationship equation of river flow-vegetation-sediment content and a simulation equation of pollution diffusion based on the historical river parameters; use the relationship equation and the simulation equation to analyze the river parameters and obtain a water environment analysis report; Step S4: continuously obtain river parameters and update the water environment analysis report; The specific steps of step S3 are as follows: Step S31: Obtain the average river flow, average vegetation coverage area and average river sediment content of the target area in the past five years; Step S32: Calculate the relationship coefficient between river flow and vegetation coverage area, denoted as b1; the calculation formula of b1 is as follows: ; Among them, ri i represents the average river flow in the past i years; ve i represents the average vegetation coverage area in the past i years, Ari represents the average value of ri1~ri5, Ave represents the average value of ve1~ve5, Sve represents the standard deviation of ve1~ve5, and Sri represents the standard deviation of ri1~ri5; Calculate the relationship coefficient between vegetation coverage area and river sediment content, denoted as b2; the calculation formula of b2 is as follows: ; Among them, tu i represents the average river sediment content in the past i years, Atu represents the average value of tu1 to tu5, and Stu represents the standard deviation of tu1 to tu5; Step S33: taking the river flow as the independent variable x; Define the relationship between river and vegetation coverage area, recorded as equation p(x); the mathematical expression of p(x) is: p(x)=b1×x+sc0; where sc0 represents the parameter to be estimated of equation p(x), and the initial value of sc0 is 0; Substitute ri1~ri5 and ve1~ve5 into p(x), calculate the value of sc0, and get sc; the final expression of p(x) is: p(x)=b1×x+sc; Step S34: Taking the vegetation coverage area as the independent variable z, define the relationship equation between river flow, vegetation and sediment content, which is recorded as equation g(x); the mathematical expression of g(x) is: g(x)=p(x)+b2×z+cs0; where cs0 represents the parameter to be estimated of equation g(x), and the initial value of cs0 is 0; Substitute ri1~ri5, ve1~ve5 and tu1~tu5 into g(x), calculate the value of cs0, and get cs; the final expression of g(x) is: g(x)=p(x)+b2×z+cs; The subsequent steps of step S34 are as follows: Step S35: Obtain the area of ​​the target region, recorded as SC; record the square root of SC as J, and perform continuous square root extraction on SC until J∈(0,10]; Based on the aerial photos of the target area, GIS is used to extract and construct a grid model of the river area; Taking the northeast end of the target area as the origin and J×J as the division size, the grid model is divided and the river area is mapped to a two-dimensional coordinate system; the points in the two-dimensional coordinate system represent the location of the river; According to the two-dimensional coordinate system, the simulation equation of river pollution diffusion is constructed to obtain the equation t (p, q); Step S36: The river area in the river parameters is recorded as RS, the vegetation area is recorded as VS, the average flow is recorded as rd, the average sediment content is recorded as tb, and the average pollution amount is recorded as ap; Calculate the increase in river flow, denoted as Δri, Δri = ee / RS; Substitute (rd + Δri) and VS into equations p(x) and g(x), respectively, to calculate the rated sediment content of the target area, denoted as ttu; Compare the size of ttu and tb; Step S37: obtaining the location of the sewage outlet in the target area and recording it as (m, n), and recording the concentration of the pollutant as po; Substitute (m, n) and po into the equation t(p, q), use the finite difference algorithm to iterate the equation t(p, q), and obtain the convergence value of the equation t(p, q), which is recorded as tuu; Compare the size of tuu and ttu; Step S38: Summarize the analysis results of steps S36 to S37 as a water environment analysis report and feed it back to the user; In step S35, the specific process of constructing the simulation equation for river pollution diffusion is as follows: Step S351: Assume that there is pollution at (p, q) in the river, and the pollutant concentration at (p, q) is C (p, q); Assume that the pollution diffusion coefficient in the horizontal axis direction is Dx, the pollution diffusion coefficient in the vertical axis direction is Dy, and the attenuation coefficient of the pollutant is kk; The relationship between Dx and Dy satisfies: Dy=[(Dx×r H +Dx×r B ) / 2]×(g×r H ×H J ) 1 / 2 ; Among them, r H represents the river depth, r B Indicates the width of the river, H J It represents the hydraulic gradient of the river; Step S352: define the simulation equation for river pollution diffusion, denoted as t(p,q); the mathematical expression of t(p,q) is: t(p,q)={Dx×[C(p,q)] 2 / p 2 }+{Dy×[C(p,q)] 2 / q 2 }-{kk×C(p,q)}。 2. The method for comprehensive assessment of water ecological environment based on remote sensing and spatial big data technology according to claim 1 is characterized in that: The specific steps of step S2 are as follows: Step S21: obtaining the current month, recorded as m; taking the last three months of the target area as standard months, namely m3, m2 and m1 as standard months; Get the daily rainfall in the standard month as the standard rainfall, recorded as rc1, rc2~rc dm ; where dm is the total number of days in a standard month; rc1 is the rainfall on the first day; rc2 is the rainfall on the second day; and so on. dm represents the rainfall on day dm; Calculate rc1~rc dm The average value is recorded as arc; calculate rc1~rc dm The variance of 2 ; Step S22: taking m3, m2, m1 and m as reference months; Get the daily rainfall in the reference month in the past three years as the historical rainfall, and get: In the first year, refer to the daily rainfall in the month: rⅠ1, rⅠ2~rⅠ dm ~rⅠ dm＋d ; Among them, rⅠ1 represents the rainfall on the first day of the reference month in the first year; rⅠ2 represents the rainfall on the second day of the reference month in the first year; and so on. dm represents the rainfall on the dmth day of the reference month in the first year; rⅠ dm＋d It represents the rainfall on the (dm+d)th day in the reference month of the first year; In the second year, refer to the daily rainfall in the month: rⅡ1, rⅡ2~rⅡ dm ~rⅡ dm＋d ; Among them, rⅡ1 represents the rainfall on the first day of the reference month in the second year; rⅡ2 represents the rainfall on the second day of the reference month in the second year; and so on. dm represents the rainfall on the dmth day of the reference month in the second year; rⅡ dm＋d It represents the rainfall on the (dm+d)th day in the reference month of the second year; In the third year, refer to the daily rainfall in the month: rⅢ1, rⅢ2~rⅢ dm ~rⅢ dm＋d ; Among them, rⅢ1 represents the rainfall on the first day of the reference month in the third year; rⅢ2 represents the rainfall on the second day of the reference month in the third year; and so on. dm represents the rainfall on the dmth day of the reference month in the third year; rⅢ dm＋d It represents the rainfall on the (dm+d)th day in the reference month of the third year; Where d represents the number of days in month m; Step S23: define the judgment condition of the change point: Assume {rc1~rc dm } in a set of breakpoints: rc r and rc p , i.e. rc r and rc p ∈{rc1~rc dm }; If rc r and rc p If the following conditions are met, it means that the rainfall in the target area conforms to the normal distribution, (rc r , rc p ) is rc1~rc dm A set of change points, rc r ~rc p It is rc r and rc p Regular section; rc r and rc p The conditions met are as follows: rc r ≠rc p ; [N(arc （r→p） ,S 2 ),rc r ~rc p ]; [(rc1≤rc r )∩(rc p ≤rc dm )]; Among them, arc （r→p） Indicates rc1~rc dm Medium r ~rc p The average value of The subsequent steps of step S23 are as follows: Step S24: record the rainfall in the target area this month as ee; According to the judgment condition of the change point, judge rc1~rc dm Is there a change point in the If yes, it means that there is regularity in the rainfall in the target area, and the rainfall in the target area this month is estimated using maximum likelihood estimation, and the process goes to step S25; If it does not exist, it means that the rainfall in the target area is irregular; based on the historical rainfall and standard rainfall, the state transfer matrix is ​​used to estimate the rainfall in the target area this month, and the process goes to step S26; Step S25: Count rc1 to rc dm The number of change points in is recorded as he; the total number of regular segments corresponding to each group of change points is counted as zo, zo = he / 2; Calculate the mean and variance of all data in the first regular segment and record them as ax(1) and S 2 (1); Calculate the mean and variance of all data in the second regular segment and record them as ax(2) and S 2 (2); Similarly, calculate the mean and variance of all data in the z0th regular segment and record them as ax(z0) and S 2 (zo); Use maximum likelihood estimation to estimate the rainfall in the target area and determine the value of ee; Step S26: Based on the historical rainfall and standard rainfall, the state transfer matrix is ​​used to estimate the rainfall in the target area this month and determine the value of ee; Step S27: Use ee as the input of the relationship equation of river flow-vegetation-sediment content and the simulation equation of pollution diffusion, and enter step S3.

3. The method for comprehensive assessment of water ecological environment based on remote sensing and spatial big data technology according to claim 2 is characterized in that: The specific steps of step S25 are as follows: Step S251: Set the number of the regular segment to j, and the value range of j is 1 to zo; Construct rc1~rc dm The primary mean periodic function is denoted as ax(j), and the primary variance periodic function is denoted as S 2 (j); Primary mean periodic function: ax(j) = Ao×cos(ω×j) + Bo×sin(ω×j) + Co; Primary variance periodic function: S 2 (j)=Do×cos(ω×j)+Eo×sin(ω×j)+Fo; Among them, Ao, Bo and Co are the parameters to be estimated of the primary mean periodic function; Do, Eo and Fo are the parameters to be estimated of the primary variance periodic function; ω represents the angular frequency, ω = (2π) / dm; Step S252: Combine ax(1) to ax(zo) and S 2 (1) ~S 2 Substitute (zo) into the primary mean periodic function and the primary variance periodic function in turn, calculate the corresponding values ​​of Ao, Bo and Co, recorded as A, B and C; calculate the values ​​of Do, Eo and Fo, recorded as D, E and F; The mean periodic function is recorded as ax-f(j), ax-f(j)=A×cos(ω×j)+B×sin(ω×j)+C; The variance periodic function is denoted as S 2 -f(j),S 2 -f(j)=D×cos(ω×j)+E×sin(ω×j)+F; Step S253: let the date be dd, the time variable be t, t=dd+dm; Using Fourier series, ax-f(j) and S 2 -f(j) is fitted and the representation range of j is broadened; Step S254: read the maximum rainfall value in the 1st to zoth regular sections and record it as rc(max); Get the total rainfall in month m1, recorded as ry; The function expression of rainfall in the target area and time is recorded as rr(t), and the mathematical expression of rr(t) is: ; where exp represents an exponential function with base e, a(t) represents the time-mean equation, and S 2 (t) represents the time-variance equation; Step S255: Get the current date and record it as dt; Substitute dt into t in rr(t) and calculate to get rr(dt); Let rr(dt) be the value of ee.

4. The method for comprehensive assessment of water ecological environment based on remote sensing and spatial big data technology according to claim 3 is characterized in that: The specific steps of step S253 are as follows: Step S2531: replace the variable j in ax-f(j) with the variable t to obtain the mixed function ax-f(t), and fit ax-f(j); The general formula of the Fourier expansion of ax-f(t) is denoted as fa(t), and the mathematical expression of fa(t) is: ； Among them, i represents the time corresponding to the jth regular segment; ao represents the coefficient of the zero-order term, an represents the coefficient of the cosine term, and bn represents the coefficient of the sine term; the values ​​of ao, an, and bn are all 1; Step S2532: taking ax-f(j) as the integrand, calculate the integrals of ao, an and bn in sequence to obtain aao, aan and abn; The calculation formula of aao is as follows: ； The calculation formula of aan is as follows: ； The calculation formula of abn is as follows: ； Use aao, aan and abn to replace ao, an and bn in fa(t) to obtain the time-mean equation of fa(t), denoted as a(t); Step S2533: Replace S with variable t 2 - The variable j in f(j) gives the mixed function S 2 -f(t), repeat the same steps of fitting xf(t) to obtain the time-variance equation, denoted as S 2 (t).

5. The method for comprehensive assessment of water ecological environment based on remote sensing and spatial big data technology according to claim 2 is characterized in that: The specific steps of step S26 are as follows: Step S261: Calculate rⅠ1~rⅠ dm The average value is recorded as rⅠa; calculate rⅠ dm+1 ~rⅠ dm＋d The average value of is denoted as rab; Calculate rⅡ1~rⅡ dm The average value is recorded as rⅡa; calculate rⅡ dm+1 ~rⅡ dm＋d The average value is denoted as rbb; Calculate rⅢ1~rⅢ dm The average value is recorded as rⅢa; calculate rⅢ dm+1 ~rⅢ dm＋d The average value of is denoted as rcb; Step S262: Record the rainfall change rate from the third year to the second year as up （3-2） , up （3-2） The calculation formula is as follows: up （3-2） =1-(rⅢa / rⅡa); The rate of change of rainfall from the second year to the first year is recorded as up （2-1） ,up （2-1） The calculation formula is as follows: up （2-1） =1-(rⅡa / rⅠa); The rate of change of rainfall from the first year to this year is recorded as up （0） , up （0） The calculation formula is as follows: up （0） =1-(rⅠa / arc); Step S263: Up （0） 、up （3-2） and up （2-1） , as well as rab, rbb and rcb; construct a state transfer matrix, denoted as matrix B; The mathematical expression of matrix B is as follows: ; Where n represents the number of iterations of matrix B, × represents matrix multiplication; n is a natural number, and its initial value is 1; each time matrix B is iterated, n is increased by 1 on the original basis; Use the NumPy library function to iterate matrix B and obtain the steady-state transfer matrix, which is recorded as matrix C; Define a 2×1 order receiving matrix {1 / 2,1 / 2}, denoted as matrix D; Calculate the value of matrix C multiplied by matrix D as the value of ee.

Images