Method for screening water quality risk factors of wetland in cold region

Abstract

The invention relates to a cold region wetland water quality risk factor screening method, and belongs to the technical field of water ecological water quality screening. Establishing a water quality monitoring data file; the water quality screening calculation comprises the steps of data file selection and identification, seasonal screening calculation, data file establishment, principal component scoring, main factor direct calculation score, non-rotation factor sample direct scoring, variance maximum orthogonal rotation, factor post-rotation scoring, Bartlett inspection, risk factor probability screening and conjecturing of over-standard elements in each month of the next year. The method has the advantages that the problem that risk factors cannot be quantitatively analyzed in a factor analysis method can be solved, and identification and analysis can be independently carried out in different seasons. According to a method for performing factor analysis and probability calculation according to seasons, the water quality risk factors are screened by combining factor analysis scores with probable probabilities.

Description

Technical Field

[0001] The present invention belongs to the technical field of water ecological water quality screening, and in particular is a method of quantitatively screening the possible probability of wetland water quality risk factors using a probability statistical analysis model, which is different from factor analysis which can only qualitatively analyze and screen risk factors. Background Art

[0002] Since environmental pollution and water ecological damage directly affect human survival and development, my country pays more and more attention to the protection of water ecology. All kinds of polluted water must be properly treated to meet the water quality standards of recycled water for different purposes. As a result, recycled water quality monitoring and evaluation as a branch of environmental science has also emerged and developed.

[0003] Since water quality monitoring data has a large amount of information and is highly random, as the year progresses, the amount of information increases dramatically, and the amount of repetition is also huge, and the difficulty of screening water quality index information affected by various random factors also increases. Principal component analysis and factor analysis can effectively reduce the dimension of data information and are effective methods to solve this problem. Through water quality monitoring of rivers, lakes, wetlands, reservoirs, urban sewage or industrial wastewater, properly treated recycled water, and analysis of data characteristics and statistical distribution problems in my country, scholars proposed the use of principal component analysis and principal factor analysis to identify, screen, analyze and evaluate water pollution components, and explore the causes of water quality changes and their laws of change.

[0004] The screening of water quality risk factors in cold-region wetlands usually adopts factor analysis methods, using scores to comprehensively evaluate the causes and impacts of risk factors, without quantitatively giving specific risk indicators. Since this type of method only conducts qualitative analysis on the random information of various pollution indicators, it cannot give quantitative indicators. It can only point out the problem, but cannot clarify the possibility and time of occurrence of key polluting elements. Such results cannot enable water quality management departments to clarify which pollution risk factors are most likely to occur, or when it is most economical and effective to deal with certain risk factors in a timely manner. Summary of the invention

[0005] The present invention provides a method for screening water quality risk factors in cold region wetlands, which can not only solve the problem that the factor analysis method cannot quantitatively analyze risk factors, but also can identify and analyze different seasons separately.

[0006] The technical solution adopted by the present invention comprises the following steps:

[0007] (I) Establish water quality monitoring data files

[0008] (1) Establish a water quality monitoring sample content data table

[0009] 1) The number of monitoring samples is the row, and the number of sample detection indicators is the column. There are n water quality monitoring samples, which is called sample capacity. Each sample has m monitoring indicators, which is the dimension of the random variable.

[0010] The monitoring sample matrix is ​​represented as an n×m-order matrix:

[0011] X=(X1,X2,…X j , …X m ) T

[0012] X j is a column vector, j = 1, 2, ..., m, or

[0013] X=Ax(1)

[0014] A is a constant matrix of n×m order consisting of monitoring data A=[a1,a2…,a m ], a j (j=1, 2, ..., m) is the column vector expressing the monitoring data, x=(x1, x2, ..., x m ) T is an m-dimensional random variable;

[0015] 2) To facilitate the understanding of the meaning and use of the monitoring sample data matrix, it is expressed in Excel format and the monitoring sample matrix is ​​supplemented with the following content:

[0016] A. Number the samples on the left side of the monitoring data matrix;

[0017] B. Add three rows at the top of the matrix, the first row is the name of the monitoring indicator, the second row is the number of the monitoring indicator, and the third row is the standard value of the monitoring indicator;

[0018] C. Add a column on the right side of the matrix for the monitoring date of the sample;

[0019] D. Add a row at the bottom of the matrix to determine whether it is greater than or less than the specified standard value;

[0020] The Excel table finally formed is the water quality monitoring sample data table;

[0021] (2) Conditions for sample data capacity required for screening calculations

[0022] Since this algorithm is a factor analysis + probability calculation for overall screening or seasonal screening, the sample size requirements are as follows:

[0023] (1) The sample size of factor analysis is greater than the number of monitoring indicators, which is a condition that must be met if only factor analysis is performed;

[0024] (2) The sample size for probability screening is greater than the number of monitoring indicators, and is generally not less than 30. When performing seasonal screening calculations, the sample size may be no less than 24.

[0025] (3) Overall screening of water quality risks in cold-region wetlands regardless of time period

[0026] For the overall screening of water quality risks in cold-region wetlands regardless of time period, you can skip the "seasonal screening calculation" item and directly conduct screening calculations.

[0027] (4) Three steps to establish the data file for screening calculations

[0028] 1) Create calculation data file

[0029] The data file name is named by the user. For easy searching, it is recommended to start with 0. The file determines the content calculated by the program. The file name cannot start with h or n. The files starting with h or n are classified files after the program recognizes them. The files starting with h are the data of the selected season or time period, and the files starting with n are the remaining data that have not been confirmed after recognition.

[0030] If you do not want to divide the calculation into seasons or time periods, you can directly perform sequential calculations after starting the calculation. It is recommended to use text files, that is, files with the extension .txt;

[0031] The first line of the data file contains the basic information and calculation requirements of the data. Spaces or commas can be used to separate the data. The details are as follows:

[0032] The first data is the number of samples + 1, which is used for factor analysis and probability screening;

[0033] The second data is the number of indicators that the sample participates in the calculation + 1, which is used for factor analysis and probability screening;

[0034] The third data is to require the cumulative variance contribution rate to be taken to decimal places for factor analysis;

[0035] The fourth data is the required cumulative variance contribution rate, which is used for factor analysis;

[0036] The fifth data is the significance level α value, which is used for probability screening;

[0037] The sixth data is the absolute value of the hidden factor loading is less than the specified number, which is used for factor analysis;

[0038] The seventh data is the probability threshold for screening risk factors. If the data is less than or equal to 0, only factor analysis is performed without possible probability analysis, which is used for probability screening.

[0039] The eighth data is the number of probability series calculation groups to determine the availability of parameter estimates, which is used for probability screening;

[0040] The ninth data is the starting month of the selected period, which is used for factor analysis and probability screening;

[0041] The tenth data is the number of months of the duration of the selected period, which is used for factor analysis and probability screening.

[0042] When performing overall screening calculations, the computer automatically performs overall screening without considering the effects of the ninth and tenth data;

[0043] The first column is the row number, which is used to easily check the sample data number; the other columns are the content specified for the monitoring sample, and each sample occupies one row and is entered in the order specified for the record;

[0044] Simplify the data below the second row: fill in the standard value row of the water quality monitoring sample data table to the nth row, and the standard value column to the mth column, and rewrite the "standard value" in the first row to the number "0". This row of data is the national standard value, which is the reference indicator for whether the water quality meets the standard;

[0045] The last line of the data file is an auxiliary supplementary line. The line number of the supplementary line is the same as the first data of the first line. The numbers of the supplementary line correspond to the data of each column. 0 means that the value of the column is less than or equal to the standard value and meets the standard. 1 means that the value of the column is greater than the standard value and meets the standard.

[0046] 2) Create a sample name file

[0047] The sample name file is a special file name that cannot be changed, and its name is: sample name.txt.

[0048] Paste the date column of the water quality monitoring sample data table into the sample name.txt file, and rewrite the date in the first row to the standard value, because the analysis results show the calculation results corresponding to the standard value in this position;

[0049] 3) Create an element name file

[0050] The element name file is a special file name that cannot be changed, and its name is: element name.txt;

[0051] Paste the first row of the water quality monitoring sample data table into element name.txt, delete the sample number and date items, arrange each element in rows, and the calculation results correspond to the input element name;

[0052] (5) After creating three data files, you can enter the screening calculation process.

[0053] The “selection and identification of data files” in the calculation process is to establish preliminary files for the screening process, including the process of “identifying data by season”;

[0054] After this step is completed, if you want to perform the "Screen and Calculate by Season" project, click "Screen and Calculate by Season", otherwise skip the "Screen and Calculate by Season" project and click "Create Data File" to proceed to the next step;

[0055] (II) The specific steps of water quality screening calculation are as follows

[0056] There are 17 steps in the water quality screening calculation, which are expressed in the form of a menu. The calculation process is carried out in the order of the menu. Step 17 "Verify the validity of the probability screening model" is independent and can be run separately as long as the data file is completed. The specific order of the menu content is as follows:

[0057] (1). “Data file selection and identification”

[0058] Click the "Selection and Identification of Data Files" item, and the computer will identify the data file.txt and sample name.txt. This is the starting step of the calculation, including the generation of the initial files for factor analysis and probability screening during the overall screening process;

[0059] 1) “Overall screening calculation”

[0060] During the overall screening calculation, the computer automatically ignores the effects of the ninth and tenth data and directly enters the screening calculation;

[0061] 2) “Data identification by season”

[0062] The ninth data in the first line of the data file is the starting month of the selected period, and the tenth data is the number of months of the selected period. For example, if the calculation data file is named 01.txt, the ninth data is filled with 10, and the tenth data is filled with 6, then the computer selects and identifies the cold season of each year as starting from October and lasting for 6 consecutive months, that is, it identifies the data of October, November, December, January, February, and March, and generates a data file starting with h.txt. At this time, the calculation data file screened by season is generated as h01.txt, the sample name file is generated as sample name_h01.txt, and the element name file is generated as element name_h01.txt;

[0063] (2) “Screening and calculation by season”

[0064] Click "Calculate by Season Screening" to calculate. The computer uses the data in the h01.txt file to generate a file ending with _h01.txt to create the initial file for factor analysis and probability screening.

[0065] (3). "Create data file"

[0066] Click "Create Data File" and the computer will automatically generate the standardized monitoring data matrix, the data format files corresponding to the principal component analysis, factor analysis and probability calculation, including the covariance matrix V of the factor analysis generated by formula (3):

[0067] D(X)=D(Ax)=AD(x)A T =AA T =V (3)

[0068] The correlation matrix R generated by equation (4) has the following elements:

[0069]

[0070] Find the eigenvalue λ of R i (i=1, 2, ..., m) and sorted according to formula (5):

[0071] λ1≥λ2...λ m ≥0 (5)

[0072] The cumulative variance contribution rate is calculated by formula (9):

[0073]

[0074] When the cumulative variance contribution rate is greater than or equal to the cumulative variance contribution rate given by the fourth data in the first line of the data file, it indicates that the first k factors can approximately express the data content of all factors, that is, the information of m factors can be expressed by k factors, achieving the purpose of dimensionality reduction;

[0075] (4). “Principal component score”

[0076] Click "Principal Component Score" to reduce the dimension of the observed data according to the cumulative variance contribution rate given by the fourth data in the first line of the data file. According to formula (7), the first principal component is the largest, and the other k-th principal components are ranked in order:

[0077]

[0078] Use the first k principal components to form a score matrix:

[0079] X k =(u1, u2, ..., u k ) T x or

[0080] X k =(X1, X2, ..., X k ) T

[0081] The score coefficients of the first k principal components are:

[0082]

[0083] Using the factor score formula (21), we can ji Change to X ji

[0084]

[0085] F ji Change to X ji ,Right now Calculate the principal component scores, and then sort the principal component scores by X. Due to space limitations, the results are not given here;

[0086] (5). "Directly calculate the score of the main factor"

[0087] Click "Calculate the score directly for the main factor" and let λ in equation (15) j =1, (j=1, 2, ..., m) to obtain the initial factor loading matrix A. The same as the principal component score, use formula (10) to calculate the score coefficient, and form the score matrix with the first k columns of the initial factor loading matrix A. Use the factor score formula (21) to convert F ji Change to A ji , that is, find the main factors and directly calculate the scores, and then sort the scores;

[0088] (6) “Directly score the sample without rotating the factors”

[0089] Click "Directly score the factor sample without rotating it" and use formula (15) to calculate the factor loading matrix A. The same as the principal component score, use formula (10) to calculate the score coefficient, and form the score matrix with the first k columns of the loading matrix A. Use the factor score formula (21) to convert F ji Change to A ji , that is, to obtain the direct scores of the unrotated factor samples, and then sort the scores;

[0090] Note: Without rotating the factor sample, the direct score should apply the first k eigenvalues ​​λ of the correlation matrix j (j=1, 2, ..., m), and the main factor directly calculates the score j =1, (j=1,2,...,m);

[0091] (7). “Variance Maximization Orthogonal Rotation”

[0092] Click “Variance Maximum Orthogonal Rotation”, build an orthogonal matrix according to formula (16), and calculate the rotation angle according to formula (18):

[0093]

[0094] Will Substituting into equation (17), the factor principal axis completes a rotation;

[0095]

[0096] By changing the values ​​of subscripts j and k in pairs, the next rotation can be performed. When the common factor is greater than 2, one cycle needs to be rotated Generally, one cycle rotation cannot achieve the best result. It takes multiple cycles of rotation until the total variance of the previous and next two cycles does not change much before the rotation can be completed.

[0097] After rotation, the contribution is increased as much as possible, and the score matrix of the maximum variance orthogonal theory is obtained;

[0098] (8). “Score after factor rotation”

[0099] Click "Score after factor rotation" to build the score matrix after factor rotation, and calculate the sorted factor analysis score and then use formula (20) to calculate the score function:

[0100] F=Bx=A T R -1 x (20)

[0101] Use formula (10) to calculate the score coefficient, and use factor score formula (21) to calculate the score, that is, to obtain the score after factor rotation, and then sort the scores;

[0102] (9) “B ARTLETT test"

[0103] Click "Bartlett test" and perform the calculation for m-dimensional independent normal distribution according to the conditions of formula (22). If the calculation result shows "no rejection of homogeneity of variance" on the screen, it means that the data passes the Bartlett test calculation;

[0104] Let n represent the capacity of the monitoring sample and m represent the dimension of the random variable. The Bartlett test is performed on the m-dimensional mutually independent normally distributed random variables using formula (22):

[0105]

[0106] The screen displays "The homogeneity of variance is not rejected";

[0107] The conclusion is that the test passes, otherwise it means that the standard deviation does not meet the homogeneity condition and the data needs to be adjusted;

[0108] (10). “MSA inspection”

[0109] It is not the main procedure used for water quality screening, but is only used as a reference. Due to space limitations, it will not be explained here.

[0110] (11). “KMO test”

[0111] It is not the main procedure used for water quality screening, but is only used as a reference. Due to space limitations, it will not be explained here.

[0112] (12). “Risk factor probability screening”

[0113] Click “Risk Factor Probability Screening”, and the computer will calculate the water quality risk factor threshold β according to the probability threshold specified by the seventh data in the first line of the data file using formula (32):

[0114] β=Φ -1 (Φ β ) (32)

[0115] The monitoring index X exceeding the probability threshold is screened in the matrix X by formula (33):

[0116]

[0117] Obtain the results of the probability analysis of the sample monitoring data, including the probability that several detection indicators may exceed the probability threshold, and give the sample time when the indicator exceeds the probability threshold;

[0118] (13). “Estimation of elements exceeding the standard in each month of next year”

[0119] Click "Predict the elements exceeding the standard in each month of next year" to calculate the possible risk probability and date of the sample monitoring indicator by formula (31) and sample name.txt and element name.txt. Predict the frequency of the possible months exceeding the specified probability threshold in the next year by formula (36):

[0120]

[0121] Estimate the possible months according to the frequency order of pollution indicators;

[0122] (14). “Draw a curve”

[0123] Click "Draw Curve" to draw the scree plot of principal component analysis and factor analysis, and various factor score curves;

[0124] (15) “Data aggregation”

[0125] 1) Data summary of overall screening calculations

[0126] Click "Data Summary" to directly summarize all calculation results in the "Data Summary_Data File.txt" file for easy viewing. For example, if the data file is 01.txt, the calculation results are summarized in the "Data Summary_01.txt" file;

[0127] 2) Summary of data by season screening

[0128] Click "Data Summary" to directly summarize all calculation results in the "Data Summary_h Data File.txt" file for easy viewing. For example, if the data file is 01.txt, the calculation results are summarized in the "Data Summary_h01.txt" file;

[0129] (16). “Reader”

[0130] Click "Reader" to display all files in the displayed directory. Click the file name you want to view to view all data files generated during the screening process;

[0131] (17). “Validation of the effectiveness of the probabilistic screening model”

[0132] Click “Verify the validity of the probability screening model” and use formula (29) to calculate the parameter estimate of the empirical distribution function according to the q group sequence:

[0133]

[0134] The parameter estimate of the empirical distribution function calculated using q groups of sequences is given by equation (31):

[0135]

[0136]

[0137] The validity of the probability model is verified using formula (30) and Table a:

[0138]

[0139] If the conditions in Table a are met, the screen displays:

[0140] The probability model meets the convergence conditions!

[0141] Otherwise, it means that the sample size is insufficient and should be increased;

[0142] For seasonal screening, if the sample size is not less than 24, it can be calculated directly without verification, and the verification process has no impact on the entire analysis process.

[0143] The advantages of the present invention are that the water quality monitoring data has a large amount of information, a large amount of information repetition, strong randomness, and strong seasonality in cold areas. By classifying and analyzing the characteristics of the information and comparing the common characteristics, a method of factor analysis and probability calculation by season is proposed, and the factor analysis score combined with the probability of occurrence is used to screen water quality risk factors. This method can not only improve the accuracy of the factor analysis score, but also clarify the probability and time of occurrence of various water quality indicators screened, so as to solve the problem that the factor analysis method cannot quantitatively analyze risk factors. Therefore, the present invention expands the content that the principal component analysis method and the principal factor analysis method cannot complete, and not only qualitatively screens with the factor analysis method, but also quantitatively screens the risk factors of wetland water quality with the probability method. BRIEF DESCRIPTION OF THE DRAWINGS

[0144] Figure 1 It is a smooth curve graph that approximates the distribution function using the step curve of the empirical distribution function. DETAILED DESCRIPTION

[0145] 1. Establish water quality monitoring data files

[0146] Taking the monitoring data of Yixi Wetland for three years as an example, according to the provisions of "1. Establish a water quality monitoring sample content data summary table", an Excel table based on the summary table is established as shown in Table 1

[0147] Table 1 Wetland water quality index monitoring data table

[0148]

[0149]

[0150] The sample numbers 1 to 48 in the table and the corresponding columns 2 to 13 are the monitoring data matrix. The rightmost column is the sample monitoring time, the leftmost column is the data sequence number, and the top row columns 2 to 13 are the monitoring indicators.

[0151] Since the water quality monitoring of cold-region wetlands is carried out throughout the year or for many years, when calculating the sample data input, the monitoring sample data and the monitoring time must be accurately marked, and the calculation content must be classified into time periods or seasons using the identification time.

[0152] ⑴ Establish calculation data file

[0153] The data file name is named by the user. For easy searching, it is recommended to start with 0. The file determines the content of the program calculation. The file name cannot start with h or n. The files starting with h or n are classified files after the program recognizes them. The files starting with h are the data of the selected season or time period, and the files starting with n are the rest of the data. If you do not want to divide the data into seasons or time periods, you can directly perform sequential calculations after starting the calculation. It is recommended to use text files, that is, files with the extension .txt, for calculation data files.

[0154] The first line of the data file contains the basic information and calculation requirements of the data. Spaces or commas can be used to separate data. The first data is the number of samples + 1; the second data is the number of indicators that the sample participates in the calculation + 1; the third data is the number of decimal places for the cumulative variance contribution rate required by the factor analysis; the fourth data is the cumulative variance contribution rate required by the factor analysis; the fifth data is the significance level α value required for the probability calculation; the sixth data is the absolute value of the hidden factor load that is less than the specified number for the factor analysis calculation to be clear; the seventh data is the threshold for the probability of screening risk factors. If the data is less than or equal to 0, only factor analysis is performed, and no possible probability analysis is performed; the eighth data is the number of probability sequence calculation groups to determine the availability of parameter estimates; the ninth data is the starting month of the selected period; the tenth data is the number of months for the duration of the selected period;

[0155] The first column is the row number, which is used to facilitate the viewing of sample data numbers; the other columns are the contents specified for the monitoring samples, each sample occupies one row and is entered in the order specified for the records.

[0156] For example, create a data file by referring to Table 1, and name the file: 01.txt;

[0157] According to the input rules of the first line of the data file, the first line of the 01.txt file is entered in sequence: 49 13 3 0.64 0.01 0.5 0.75 5 10 6

[0159] The data sequence of the first row is described as follows:

[0160] The first number is the number of rows in the monitoring data matrix, which is 49-1=48;

[0161] The first number is the number of columns in the monitoring data matrix, which is 13-1=12;

[0162] The third number is the cumulative variance contribution rate with 3 decimal places;

[0163] The fourth number requires that the cumulative variance contribution rate of the analysis indicator be greater than 0.64;

[0164] The fifth number specifies the significance index α value of Bartlett's test to be 0.01;

[0165] The sixth number is the hidden factor loading less than 0.5;

[0166] The seventh number specifies the desired probability threshold as 75%;

[0167] The eighth number is to select 5 groups to perform sequence calculations to determine the sample size and verify whether it passes the verification;

[0168] The ninth number is the cold season, which starts from October;

[0169] The tenth number is the length of the cold season, which is 6 months.

[0170] The data below the second line can be simplified: paste the standard value row of Table 1 to row 48 and the standard value column to column number 12 in the second line of the 01.txt file, and rewrite the "standard value" in the first line to the number "0". This line of data is the national standard value, which is the reference data for whether the water quality meets the standard.

[0171] The last line of the data file needs to be supplemented. The line number of the supplemented line is the same as the first data of the first line. The numbers of the supplemented line correspond to the data in each column. 0 means the direction is opposite to the standard value, and 1 means the direction is the same as the standard value. For example, the national standard value in the 7th column requires that it is greater than a certain value to meet the standard, so it is represented by 1. The rest are all national standard values ​​that are less than a certain value to meet the standard, so they are represented by 0.

[0172] The contents of the 01.txt file are shown in Table 2

[0173] Table 2 01.txt file format

[0174]

[0175]

[0176] ⑵Create a sample name file

[0177] The sample name file is a fixed file name that cannot be changed, and its name is: sample name.txt.

[0178] Paste the date column in Table 1 into the sample name.txt file and rewrite the date in the first row to the standard value, because the analysis results show the calculation results corresponding to the standard value in this position.

[0179] The content of the sample name.txt file is shown in Table 3

[0180] Table 3 Sample name.txt file format

[0181] Standard value

[0182] 2021-01-07C

[0183] 2021-01-25C

[0184] 2021-03-01C

[0185] 2021-03-16C

[0186] 2021-04-19C

[0187] 2021-05-12C

[0188] 2021-05-28C

[0189] 2021-06-17C

[0190] 2021-07-19C

[0191] 2021-08-10C

[0192] 2021-08-26C

[0193] 2021-09-10C

[0194] 2021-10-22C

[0195] 2021-11-25C

[0196] 2021-12-07C

[0197] 2021-12-23C

[0198] 2022-01-18C

[0199] 2022-02-16C

[0200] 2022-05-10C

[0201] 2022-06-02C

[0202] 2022-06-26C

[0203] 2022-07-04C

[0204] 2022-07-15C

[0205] 2022-08-02C

[0206] 2022-08-21C

[0207] 2022-09-06C

[0208] 2022-09-26C

[0209] 2022-10-24C

[0210] 2022-11-27C

[0211] 2022-12-12C

[0212] 2023-01-08C

[0213] 2023-02-02C

[0214] 2023-03-10C

[0215] 2023-03-20C

[0216] 2023-04-04C

[0217] 2023-04-18C

[0218] 2023-05-18C

[0219] 2023-06-16C

[0220] 2023-07-12C

[0221] 2023-08-11C

[0222] 2023-09-18C

[0223] 2023-10-18C

[0224] 2023-10-27C

[0225] 2023-11-14C

[0226] 2023-11-21C

[0227] 2023-11-28C

[0228] 2023-12-04C

[0229] 2023-12-20C

[0230] ⑶Create element name file

[0231] The element name file is a fixed file name that cannot be changed, and its name is: element name.txt.

[0232] Paste the first row of Table 1 into element name.txt and delete the sample number and date items. Arrange each element in rows and the calculation results correspond to the input element names.

[0233] The content of the element name.txt file is shown in Table 4

[0234] Table 4 Element name.txt file format

[0235] PH

[0236] COD

[0237] Ammonia nitrogen

[0238] CadmiumCD

[0239] Total Phosphorus

[0240] Total Nitrogen

[0241] Dissolved oxygen

[0242] Nitrite Nitrogen NO2-N

[0243] Iron

[0244] Manganese Mn

[0245] Fluoride F ˉ

[0246] Nitrate Nitrogen NO3-N

[0247] 2. Water quality screening calculation

[0248] Start the "Cold Region Recycled Water Wetland Water Quality Risk Factor Screening" program and begin the water quality risk factor screening calculation. All the following steps correspond to "II. The specific steps of water quality screening calculation are as follows":

[0249] (1) “Data file selection and identification”

[0250] Click "Selection and Identification of Data Files" and enter the data file 01.txt as prompted. Data files h01.txt and n01.txt starting with h and n respectively are automatically generated. The original data file 01.txt remains unchanged. The overall screening content file is established, including the initial files for factor analysis and probability screening during the overall screening process.

[0251] 1) “Overall screening calculation”

[0252] During the overall screening calculation, the computer automatically ignores the effects of the ninth and tenth data and directly enters the screening calculation.

[0253] 2) “Data identification by season”

[0254] The ninth data in the first line of the data file is the starting month of the selected period, and the tenth data is the number of months of the selected period. After identification, the data file h01.txt (Table 5) is generated. The h01.txt file is the cold season calculation data file generated after identification.

[0255] In this example, the data file is 01.txt (Table 2). In the first row of h01.txt (Table 5), only the first data changes. This number minus 1 is the number of samples identified this time. The data in the first column corresponds to the row number identified in the 01.txt file, and also corresponds to the row number of the generated sample name_h01.txt (Table 6).

[0256] In this example, the ninth digit is 10 and the tenth digit is 6. The computer recognizes that the cold season of each year starts from October and lasts for 6 consecutive months. The calculation data file generated after recognition is shown in Table 5, and the sample name file is shown in Table 6:

[0257] Table 5 Calculation data file h01.txt generated after recognition

[0258]

[0259] As can be seen from the table, the first value in the first line of the file is 25, that is, the recognized cold season data is 25-1=24, that is, 24 sample data are recognized.

[0260] Table 6 Sample name file generated after recognition Sample name_h01.txt

[0261] Standard value

[0262] 2021-01-07C

[0263] 2021-01-25C

[0264] 2021-03-01C

[0265] 2021-03-16C

[0266] 2021-10-22C

[0267] 2021-11-25C

[0268] 2021-12-07C

[0269] 2021-12-23C

[0270] 2022-01-18C

[0271] 2022-02-16C

[0272] 2022-10-24C

[0273] 2022-11-27C

[0274] 2022-12-12C

[0275] 2023-01-08C

[0276] 2023-02-02C

[0277] 2023-03-10C

[0278] 2023-03-20C

[0279] 2023-10-18C

[0280] 2023-10-27C

[0281] 2023-11-14C

[0282] 2023-11-21C

[0283] 2023-11-28C

[0284] 2023-12-04C

[0285] 2023-12-20C

[0286] As can be seen from the table, the dates in the file are all dates that specify the cold season, and the file data is 24+1=25 because a row of standard values ​​is added.

[0287] Because the monitoring indicators have not changed, the content of the element name_h01.txt file is the same as the element name.txt content and is not expressed here.

[0288] During the overall screening calculation, the computer automatically ignores the effects of the ninth and tenth data and directly enters the screening calculation;

[0289] (2) “Screening and calculation by season”

[0290] Click "Screen and calculate by season", and the computer will create an initial file for factor analysis and probability screening for the h01.txt file, that is, factor analysis and probability calculation are performed on files ending with h01.txt and _h01.txt. All data files related to h01.txt generated during the calculation process end with _h01.txt. For example, the sample name file obtained after identification is sample name_h01.txt, and so on; the element name content expressed by the monitoring indicator has not changed;

[0291] (3) "Create data file"

[0292] Click "Create Data File" and the computer will automatically generate a standardized monitoring data matrix, data format files corresponding to principal component analysis, factor analysis and probability calculation, including the covariance matrix V of the factor analysis generated by formula (3), the correlation matrix R generated by formula (4), and the eigenvalue λ of R i (i=1,2,…,m) and sorted according to formula (5), and the cumulative variance contribution rate is calculated by formula (9). When the cumulative variance contribution rate is greater than or equal to the cumulative variance contribution rate given by the fourth data in the first line of the data file, it means that the first k factors can approximately express the data content of all factors, that is, the information of m factors can be expressed by k factors, and the dimension reduction is completed. Due to space limitations, the results are not given here.

[0293] (4) “Principal component score”

[0294] Click "Principal Component Score" to reduce the dimension of the observed data according to the cumulative variance contribution rate given by the fourth data in the first line of the data file. Use the first k principal components to form a score matrix, calculate the score using the factor score formula, calculate the principal component score, and then sort the principal component scores by X. Due to space limitations, the results are not given here;

[0295] (5) "Direct calculation of scores for main factors"

[0296] Click "Calculate the score directly for the main factor" and let λ in equation (15) j =1, (j=1,2,…,m) to calculate the initial factor loading matrix A. The same as the principal component score, use formula (10) to calculate the score coefficient, the first k columns of the initial factor loading matrix A form a score matrix, use the factor score formula (21) to calculate the main factor directly to calculate the score, and then sort the scores. Due to space limitations, the results are not given here;

[0297] (6) “Direct scoring of unrotated factor samples”

[0298] Click "Direct score of unrotated factor sample", and use formula (15) to calculate the factor loading matrix A. Same as the principal component score, use formula (10) to calculate the score coefficient, and form the score matrix with the first k columns of the loading matrix A. Use the factor score formula (21) to calculate the direct score of the unrotated factor sample, and then sort the scores. Due to space limitations, the results are not given here;

[0299] (7) “Variance Maximization Orthogonal Rotation”

[0300] Click "Variance Maximum Orthogonal Rotation", build an orthogonal matrix according to formula (16), and calculate the rotation angle according to formula (18)

[0301]

[0302] Will Substituting into formula (17), the principal axis of the factors completes a rotation.

[0303]

[0304] By changing the values ​​of the subscripts j and k in pairs, the next rotation can be performed. For example, if the common factor is 4, the rotations in one cycle are Then select 4 cycles.

[0305] After rotation, the contribution is increased as much as possible, and the score matrix of the maximum variance orthogonal theory is obtained as shown in Table 7.

[0306] Table 7 Common factor score matrix (the last column is the sample comprehensive score)

[0307]

[0308]

[0309] (8) “Score after factor rotation”

[0310] Click "Score after factor rotation" to build the score matrix after factor rotation, and calculate the sorted factor analysis score and then use formula (20) to calculate the score function:

[0311] F=Bx=A T R -1 x (20)

[0312] The score coefficient is calculated using formula (10), and the score is calculated using factor score formula (21). That is, the score after factor rotation is obtained, and then the score is sorted. The sorted scores are shown in Table 8.

[0313] Table 8 After the common factor comprehensive scores are sorted, the third column is the original position number of the second column

[0314]

[0315] As can be seen in the table, the standard value has the lowest score, which means the water quality is the best, which is consistent with the actual situation;

[0316] (9) “B ARTLETT test"

[0317] Click "Bartlett test" to perform a homogeneity test of the normal distribution variance on the monitoring data matrix according to formula (22). In this example, n = 25 represents the capacity of the monitoring sample, m = 12 represents the dimension of the random variable, and the significance level α = 0.01:

[0318]

[0319] in, It is called the pooled variance;

[0320] Given a confidence level of 1-α

[0321] Bartlett's test results

[0322] χ 2 =2064

[0323] α = 0.01 degrees of freedom = 11χ 2 Table value = 14.63

[0324] The screen displays "Do not reject homogeneity of variances"

[0325] The conclusion is that the test has passed;

[0326] (10) “MSA test”

[0327] It is not the main procedure used for water quality screening, but is only used as a reference. Due to space limitations, it will not be explained here.

[0328] (11) “KMO test”

[0329] It is not the main procedure used for water quality screening, but is only used as a reference. Due to space limitations, it will not be explained here.

[0330] (12) “Risk Factor Probability Screening”

[0331] Click "Risk Factor Probability Screening" and the computer will calculate the probability threshold Φ specified by the seventh data in the first line of the data file. β =75%, the water quality risk factor threshold β is calculated by formula (32)

[0332] β=Φ -1 (Φ β )=0.674 (32)

[0333] According to formula (33), the monitoring index X exceeding the probability threshold is screened in the matrix X

[0334]

[0335] The results of the probability analysis of the sample monitoring data are obtained. The results include the probability that several detection indicators may exceed the probability threshold, and the calculation results are shown in Tables 9, 10, and 11.

[0336] Table 9 Ranking of the probability of possible risks in monitoring indicators

[0337]

[0338]

[0339] Table 10 Indicators of samples that may have a risk of exceeding the probability threshold of 75%

[0340] 2 94.02% COD

[0341] 4 77.14% Cadmium CD

[0342] 6 99.84% Total Nitrogen

[0343] Table 11: Probability and date of risk of exceeding the probability threshold of 75% that may occur for sample elements

[0344]

[0345] Note: The overall probability of occurrence of the three elements may exceed the probability threshold, but the probability of their occurrence in each time period may not exceed the probability threshold.

[0346] (13) “Estimation of elements exceeding the standard in each month of next year”

[0347] Click "Predict the elements exceeding the standard in each month of next year" to calculate the possible risk probability and date of the sample monitoring indicator by formula (31) and sample name.txt and element name.txt. Predict the frequency of the possible months exceeding the specified probability threshold in the next year by formula (36)

[0348]

[0349] The inference results are shown in Table 12

[0350] Table 12: The frequency of occurrence of months that may exceed the specified probability threshold of 75% in the next year

[0351]

[0352]

[0353] The possible order of months is inferred according to the frequency order of pollution indicators, and the calculation results are shown in Table 13.

[0354] Table 13: Possible month rankings based on the frequency order of pollution indicators

[0355]

[0356]

[0357] (14) “Draw a curve”

[0358] Click "Draw Curve" to draw the scree plot of principal component analysis and factor analysis, and various factor score curves.

[0359] (15) “Data aggregation”

[0360] This example only calculates the data summary by season screening

[0361] Click "Data Summary" to directly summarize all calculation results in the "Data Summary_h01.txt" file, and all calculation results can be viewed in the data summary file.

[0362] (16) “Reader”

[0363] Click "Reader" to display all files in the displayed directory. Click the file name you want to view to view all data files generated by the screening process.

[0364] (17) “Validation of the effectiveness of the probabilistic screening model”

[0365] Click "Verify the validity of the probability screening model" and use formula (29) to calculate the parameter estimate of the empirical distribution function according to the q group sequence

[0366]

[0367] The comparison process of the calculation results of parameter estimates is shown in Tables 14 and 15.

[0368] Table 14 Calculation and result comparison of mathematical expectation estimator sequence

[0369]

[0370] Table 15 Calculation and result comparison of standard deviation estimator series

[0371]

[0372] The parameter estimate of the empirical distribution function calculated using q-group sequences is given by equation (31)

[0373]

[0374] The calculation results are shown in Table 16;

[0375] Table 16 Calculation and comparison of probability series of samples in each group

[0376]

[0377] The calculation using formula (30) is to substitute each data in the probability average row and each data corresponding to the 5th row into formula (30) to obtain the last row, and use this row of results to check against Table a for verification. If the verification passes, it means that the probability model has sufficient accuracy.

[0378]

[0379] The calculation results show that the last row of Table 17 meets the requirements of Table a.

[0380] If the conditions in table a are met, the screen displays

[0381] According to the principle of practical impossibility of low-probability events, that is, the idea that low probability has no effect on the screening conclusion, the model that meets the following conditions converges:

[0382] 1. The model convergence condition is that the relative error is less than 50% when the probability is less than 10%.

[0383] 2. The relative error of the probability greater than 10% and less than or equal to 25% is less than 30%

[0384] 3. The relative error of the probability greater than 25% and less than or equal to 50% is less than 20%

[0385] 4. For probability greater than 65%, the relative error is less than 5%

[0386] The probability model meets the convergence conditions!

[0387] This example is seasonal screening, and the sample size is 24. It is proved that the sample size is 24, so it can be calculated directly without verification. Of course, the verification process has no impact on the entire analysis process.

[0388] The present invention is further described below by theoretical basis.

[0389] 1. Factor analysis calculation model

[0390] 1. Build a data model

[0391] For simplicity, bold letters are used to represent vectors or matrices.

[0392] The number of monitoring samples is taken as rows, and the number of sample detection indicators is taken as columns. The monitoring sample data is represented by a matrix. There are n water quality monitoring samples, which is called the sample capacity. Each sample has m monitoring indicators, which is the dimension of the random variable. The monitoring sample matrix is ​​represented as an n×m-order matrix:

[0393] X=(X1,X2,…X j , ...X m ) T , whose corresponding matrix element is X ij .

[0394] X j is a column vector, j = 1, 2, ..., m, or

[0395] X=Ax (1)

[0396] A is a constant matrix of n×m order consisting of monitoring data A=[a1,a2…,a m ], a j (j=1, 2, ..., m) is the column vector expressing the monitoring data, x=(x1, x2, ..., x m ) T is an m-dimensional random variable. Obviously, X=(X1,X2,…,X m ) T is the matrix of observable data.

[0397] In order to make the data comparable, the elements of X are considered to be dimensionless data after standardization. Then the mathematical expectation E(X j )=0,variance D(X j )=I,A n×m is a constant matrix and I is the identity matrix.

[0398] This model uses the standard deviation to process the original data items into a standard normal distribution, that is, subtracting the mean of the same variable and then dividing it by the standard deviation.

[0399] The standardized matrix is ​​recorded as X', and the element X is ij Convert to

[0400] X′ ij =(X ij -X j ) / S j (2)

[0401] Where:

[0402]

[0403] S j is the standard deviation of the sample

[0404]

[0405] For simplicity of expression, the standardized matrix X' is still represented by X, and X in factor analysis is always a standardized matrix.

[0406] 1. Principal component analysis calculation model

[0407] The principal component analysis method considers that the sample is a linear combination of monitoring data. By identifying the variance relationship in the data and using the method of obtaining the maximum value of the variance, the information carried by the variables is lost as little as possible, thereby achieving the purpose of summarizing the data information.

[0408] Let X = (X1, X2, ..., X m ) T is the observable data matrix. To make the data comparable, the elements of X are considered dimensionless data that have been standardized by formula (2). Then the mathematical expectation E(X) = 0, the variance D(X) = I, and A n×m is a constant matrix.

[0409] According to formula (1), the covariance of X is

[0410] D(X)=D(Ax)=AD(x)A T =AA T =V (3)

[0411] V is the non-negative covariance matrix of X, whose elements are

[0412]

[0413] Then its correlation coefficient matrix R xx

[0414]

[0415] By D(X)=I, that is have

[0416] r ij = cov(X i , X j )=v ij

[0417] That is to say, in the case of standard normal distribution, the covariance matrix V is numerically equal to the correlation coefficient matrix R between X and x. Xx , but the dimensions are different. In the following, R Xx Denoted as R. According to the symmetry of the correlation matrix, R T = R, so there exists an orthogonal matrix U such that

[0418] U T RU=Λ

[0419] Λ is the characteristic value λ i The diagonal matrix formed by the principal component is to obtain the maximum value of the variance, that is, to find MaxR.

[0420] From the extreme value and non-negative definite properties of the eigenvalue, we can see that the eigenvalue of R is λ i ≥0(i=1,2,...,m), sort them from large to small, then we have

[0421] λ1≥λ2...λ m ≥0 (5)

[0422] On the other hand, the maximum variance of R is the trace of R, that is,

[0423]

[0424] Compared with formula (5), λ i It is a linear function The maximum value of the variance of R is λ i The diagonal matrix Λ (i=1, 2, ..., m) is obtained from this. The calculation method of the principal component of X is:

[0425] The correlation coefficient matrix R is obtained by formula (4) and its eigenvalue is obtained. The sorted eigenvalue λ is used to i The corresponding eigenvector u i Construct an orthogonal matrix U, u i is the eigenvector that constitutes the i-th principal component. From equations (1), (4), and (5), we can get the principal component algorithm. The first principal component is the largest, and the other i-th principal components are arranged in order.

[0426]

[0427] In general, the variance contribution rate of the first few largest principal factors can represent most of the information. The larger its variance, the greater its contribution. In the calculation of principal components, it is very convenient to use eigenvalues ​​to calculate the contribution rate of principal components.

[0428] Definition: The variance contribution rate of the jth principal component is

[0429]

[0430] Assume that the number of principal components is k (k<m), and the cumulative contribution rate of the first k principal components is

[0431]

[0432] The score coefficients of the first k principal components are

[0433]

[0434] 2 Factor analysis calculation model

[0435] Factor analysis statistical model

[0436] Factor analysis is to group the size of the correlation between variables using the internal connection of the correlation coefficient matrix, which is called a common factor, so that the relationship between each factor is clear and the reliability and validity analysis of the acquired data can be performed. Thus, a method of factor analysis using the R-type orthogonal factor model with the help of the principal component method was obtained.

[0437] Assume a random vector x, F, ε, x = (x1, x2, ..., x n ) T is a standardized observable random vector, E(x) = 0, cov(x) = V, F = (F1, F2, ..., F m ) T , (m<n) is a common factor, which is an unobservable random vector, F~N(0,I), cov(F)=I, E(F)=0, cov(F,ε)=0, ε=(ε1,ε2,…,ε n ) T is the unobservable random error vector, ε~N(0,∑), E(ε)=0, The covariance is a diagonal matrix. A n×m It is called the factor loading matrix, which is a constant matrix. The factor analysis model can be expressed as

[0438] x=AF+ε (11)

[0439] Remember the expression

[0440]

[0441] Factor analysis is to find the common factor load with the largest contribution possible while keeping the commonality of the load factors unchanged, so as to explain the meaning and role of the factor load.

[0442] According to formula (11),

[0443] D(x)=AA T +∑

[0444] It can be seen that AA T +∑ is the correlation coefficient matrix R between x and F in formula (10) xF .Right now

[0445]

[0446] According to formula (11), It is still a symmetric non-negative definite matrix with rank ≤ m, and its diagonal elements are transformed into common degrees Will Contribution The degrees are arranged in order starting from 1, and k < m contribution degrees are taken. Then there is

[0447]

[0448] Using the Lagrange multiplier method, we can prove The maximum value of In practical applications, the least squares regression combined with the matrix singular decomposition method can prove that when only the first few contributions are taken, Alternative R xF The error is minimal. Let λ i The corresponding standard is orthogonal eigenvector u i , then R xF Available Express

[0449]

[0450]

[0451] Then the expression of A is

[0452]

[0453] 3-factor rotation calculation model

[0454] In order to increase the contribution of common factors as much as possible and to better explain the meaning of common factors, the maximum variance rotation method (Varimax) can be used to construct the standard orthogonal matrix Γ, ΓΓ T =I, ignoring the influence of ∑, formula (10) can be expressed as:

[0455] x=AF=Arr T F

[0456] Let A * =AΓ=B,F * =Γ T F, then

[0457] x=A * F *

[0458] It is also a factor analysis model, calculating the mathematical expectation and variance

[0459] E(F * )=Γ T E(F)=0

[0460] D(F * )=ΓD(F)Γ T =I

[0461] D(x)=A * A *T =AΓ(AΓ) T =AA T =BB T

[0462] Comparing with formula (11), it can be seen that after linear transformation, B is still a factor loading matrix, and

[0463]

[0464] Indicates the commonality of loads after orthogonal transformation No change. Since the contribution is the sum of squares, it means that the orthogonal transformation can be used to adjust the factor contribution without changing the commonality of the factors, so that the contribution can be increased as much as possible. This transformation load matrix is ​​called the orthogonal factor axis rotation with the largest variance. After several rotations, the common factor contribution can be close to the maximum. F before and after rotation j Contribution to x There are changes during the rotation The factors begin to polarize towards 0 and 1. This is also the purpose of factor rotation, which makes the factors have clear practical significance when explaining the original sample information, which is conducive to the interpretation of common factors. Consider two factor loads a j , a k Orthogonal rotation of .

[0465] A=(a j , a k )but i=1,2,...,n

[0466] make

[0467]

[0468] Take the orthogonal matrix

[0469]

[0470] Let B = AΓ, then

[0471]

[0472] Find the sum of the variances of two columns

[0473]

[0474] make The maximum value of G can be obtained

[0475] remember

[0476]

[0477] but

[0478]

[0479] That is When G reaches its maximum value, the factor principal axis completes a rotation. According to the inverse trigonometric function, the principal value interval of this rotation is By changing the values ​​of subscripts j and k in pairs, the next rotation can be performed. When the common factor is greater than 2, one cycle needs to be rotated Generally, one cycle rotation cannot achieve the best result. It takes multiple cycles of rotation until the total variance of the two cycles does not change much, and then the rotation can be completed.

[0480] 4-factor score calculation model

[0481] The ultimate goal of factor analysis is to calculate the comprehensive score of the factors and then conduct a comprehensive evaluation. In the case of orthogonal factors, the main factors are independent of each other, and each main factor represents a certain causal concept. Therefore, the main factor score can be used to describe the spatial performance of this environmental cause. Since the observation data matrix is ​​generally not a square matrix, it can only be reduced by the matrix singular decomposition method combined with the least squares regression method based on the statistical model principle, so that the error of replacing the original matrix with the first m factors is minimized, so special factors can be ignored. By rotating, a large amount of redundancy is eliminated from the m factors, which greatly reduces the subsequent data processing pressure and brings convenience to data analysis. According to formula (10), ignoring the influence of special factors, we have

[0482] x=AF

[0483] Construct the linear combination of the main factor F represented by the random variable x as:

[0484] F=Bx (19)

[0485] Formula (18) is called the main factor score function. Since the number of equations in the formula is greater than the number of variables, the matrix singular decomposition is used and combined with the least squares method to find the optimal solution. Let the generalized minus inverse of A be B, that is

[0486] B=A -

[0487] Using matrix singular decomposition method to reduce dimension combined with least squares regression, it can be proved that Right now

[0488]

[0489] The expression of the common factor matrix F calculated by the standardized observation data matrix x is obtained. The contribution is calculated using formula (13). If the contribution of the first k (k < m) common factors has reached the requirement, the cumulative variance contribution rate is calculated using formula (9), and the score coefficient w is calculated using formula (10). i , calculate the score as follows

[0490]

[0491] A can be the matrix before rotation or after rotation. Since rotation is achieved through orthogonal linear transformation, the factor loading is not unique, the matrix after rotation is not unique, and therefore, the score result is not unique.

[0492] 5. Random variable test of principal component and factor analysis + probability calculation method

[0493] The water quality risk factors of cold-region wetlands are screened by principal component and factor analysis + probability calculation method, and the random variables are generally tested;

[0494] (1) Random variable test of principal component and factor analysis algorithms

[0495] The principal component and factor analysis methods are dimensionality reduction calculations, using test random variables to assist in determining the accuracy of the calculations;

[0496] There are two ways to test random variables using principal component and factor analysis:

[0497] MSA test;

[0498] KM0 test;

[0499] Because the MSA test and KM0 test are routine tests in general principal component analysis and factor analysis, they only work on random variables with small data changes. They can only be used as reference indicators for the test values ​​of cold-region wetland water quality risk factor screening monitoring data. The test results have no effect on the screening results, so they will not be elaborated here.

[0500] (2) Bartlett ARTLETT ) method to test the homogeneity of variance of normal distribution

[0501] n is used to represent the capacity of the monitoring sample, and m is used to represent the dimension of the random variable or the total number. The specific method is as follows:

[0502] Suppose that m samples are randomly selected independently from m normal populations, and the size of each sample is n. i , the sample variance is The hypothesis test is then:

[0503] The variances of the populations are not all equal;

[0504] Under the condition that H0 is established, the Bartlett test statistic is:

[0505]

[0506] in, It is called the pooled variance;

[0507] Given a confidence level of 1-α

[0508] like Then H0 is not rejected;

[0509] like Then reject H0 and accept H1.

[0510] If the conclusion is "do not reject the homogeneity of variance", the test passes;

[0511] If the homogeneity test fails, it means that a certain test indicator has a singular value, and the homogeneity test can be passed by adjusting the sample data;

[0512] 2. Probabilistic Screening Calculation Model

[0513] 1. Monitoring data random variables can be converted into one-dimensional random variables

[0514] The probability screening algorithm is based on the law of large numbers and the central limit theorem, which requires random variables to obey the normal distribution. The water quality observation samples of cold-region wetlands are a set of data composed of multiple monitoring indicators, that is, the water quality observation samples are random variables composed of a set of random vectors, and their distribution is a multivariate function.

[0515] As is well known, the relevant multivariate distribution function analysis is relatively complicated, the difficulty of calculation and solution is relatively large, and the practical application is limited. If it can be proved that the normal distribution is independent of each other, it does not affect the overall sample to obey the normal distribution, then the difficulty of calculation and solution will be greatly reduced, and the problem will be greatly simplified. Therefore, if it is proved that the random vector component is independent and normal and does not affect the overall normal distribution, it is proved that the random variation law of the water quality monitoring data component of the present invention obeys the law of large numbers and the central limit theorem, and the law and theorem can be used for probability analysis. Only when the multiple monitoring indicators of the water quality monitoring sample are considered according to the one-dimensional normal distribution respectively, can the screening calculation of water quality risk factors be realized using one-dimensional probability statistics theory.

[0516] The present invention applies the uniqueness relation of characteristic functions, that is, the necessary and sufficient condition for the identity of two distribution functions F1(x) and F2(x) is their characteristic functions and The identity theorem proves that the n-dimensional random vector X = (X1, X2, ..., X n ) T It is composed of monitoring indicators of observed samples, each monitoring indicator X i , i=1,2,...n are mutually independent normal distributions, which does not affect the sample population X obeying the normal distribution, and vice versa. Therefore, for m-dimensional mutually independent normal distributions, it can be decomposed into m one-dimensional normal distribution functions to solve.

[0517] Since the monitoring data of cold region wetlands are all independent random variables, the multidimensional subvector sequence {X i , i=1,2,…} are all independent and have the same normal distribution. Their mathematical expectation E(X i ) exists and variance 0<D(X i )<+∞, according to the Lindburg-Levy theorem, the multidimensional subvector sequence {X i , i=1,2,…} obey the central limit theorem;

[0518] Because the cold-region wetland monitoring data is observable, the eigenvalues ​​of its variance matrix exist and are positive, positive definite and finite. Therefore, the sub-vectors of the random vector matrix established by the cold-region wetland water quality observation data are independent of each other, obey the law of large numbers and satisfy the central limit theorem conditions, that is, the random variable sequence of the cold-region wetland water quality monitoring data is independent and normally distributed, positive definite and finite, and meets the positive definiteness requirements of the characteristic function. According to the properties of the multidimensional normal distribution, if the m-dimensional normal distribution has correlation, that is, X~N(μ,∑), where E(X)=μ, D(X)=∑, according to the properties of the multidimensional normal distribution, there exists an orthogonal matrix U, Y=U T(X-μ), so that the components Y1, Y2, ..., Y m are independent of each other, and Y i ~N(0,λ i ), i=1, 2,..., m, where λ1, λ2,..., λ m is the eigenvalue of the matrix ∑. That is to say, if there is correlation, the overall distribution can be decomposed into m independent normal distributions by solving the eigenvalue of the covariance matrix. That is to say, if there is correlation in the m-dimensional normal distribution, the overall distribution can be decomposed into m independent normal distributions.

[0519] The present invention applies the uniqueness relation of characteristic functions, that is, the necessary and sufficient condition for the identity of two distribution functions F1(x) and F2(x) is their characteristic functions and The identity theorem proves that the m-dimensional random vector X = (X1, X2, ..., X m ) T It is composed of monitoring indicators of observed samples, each monitoring indicator X i , i=1,2,...m are mutually independent normal distributions, which do not affect the sample population X's compliance with normal distribution, and vice versa. Therefore, for m-dimensional mutually independent normal distributions, it can be decomposed into m one-dimensional normal distribution functions to solve;

[0520] Considering multiple monitoring indicators of water quality monitoring samples according to one-dimensional normal distribution, when the sample size is large enough, it is theoretically possible to use the law of large numbers and the central limit theorem to screen water quality risk factors;

[0521] For the sake of simplicity, the m random variables of the monitoring indicators are considered to be independent random variables;

[0522] According to the law of large numbers, if the random variable sequence of cold region wetland water quality monitoring data {X i , i = 1, 2, ...} are independent and normally distributed. Obviously, its mathematical expectation E(X i ) exists, variance 0<D(X i )<+∞, construct the empirical normal distribution function from the sample Its second moment exists, then for any real number ε>0, when n→∞:

[0523]

[0524] This is called mean square convergence. Equation (4) indicates that as the sample size increases, the second-order moment of the sample converges to the second-order moment of the population, that is, the possibility of a difference between the second-order moment of the sample and the second-order moment of the population will become smaller and smaller. According to the theorem of probability theory, the mathematical expectation of its first-order moment is also convergent;

[0525] According to Glivenko's theorem, the population distribution function of X is F(x), and the empirical distribution function is F n (x), for any real number x, when n→+∞, we have

[0526]

[0527] Grevenko's theorem shows that when the sample size n is large enough, the empirical distribution function F n (x) converges to F(x) with probability 1 about x, which is the theoretical basis for using samples to infer the population. Formula (24) shows that when n is large enough, for all x values, the empirical distribution function F n (x) is very close to the overall distribution function F(x), such as Figure 1 shown.

[0528] That is to say, for a random variable, construct an empirical normal distribution function. When the sample size n is large enough, the empirical normal distribution function F n (x) converges to the normal distribution function F(x). At this time, it can be considered that when the sample size is greater than a certain value n, in the interval where the sample size is greater than n, its empirical distribution function F n (x) is very close to the distribution function F(x), and the parameters of its distribution function are also very close. According to the uniqueness theorem of characteristic function, the normal empirical characteristic function is constructed and the normal distribution characteristic function The first and second order moments of should also be very close.

[0529] 2. Lindberg-Feller Conditions for Satisfying the Central Limit Theorem

[0530] Assume {X n , n=1,2,...} is a sequence of independent random variables with finite mathematical expectation and variance. If the Lindberg condition is satisfied, that is, if for any ε>0,

[0531]

[0532] in, For X i The density function, μ i =E(X i ), Then {X n , n=1, 2, ...} obey the central limit theorem. Obviously, So there is the Feller condition

[0533]

[0534] Formula (25) and (26) constitute the necessary and sufficient conditions for the central limit theorem to hold. If formula (25) and (26) are satisfied, then Then Lyapunov's central limit theorem holds.

[0535] lim n→∞ P(Y n ≤y)}=Φ(y) (27)

[0536] Formula (26) shows that when n is large enough, Y n Each added term is uniformly small in probability, and the limiting distribution of the sum of these random variables is normal distribution, that is, it obeys the central limit theorem under large samples, and the probability can be calculated using the standard normal distribution function.

[0537] Because the mean and variance of water quality monitoring data can be directly obtained and are bounded, the Lindberg-Feller condition must be met. Equation (27) shows that if the Lindberg-Feller condition is met, the normal distribution can be used for probability calculation.

[0538] 3. Parameter estimators of the empirical distribution function

[0539] Let F n (x) is the empirical distribution function, and its characteristic function is and F(x) is the normal distribution function, and its characteristic function is According to the uniqueness of the characteristic function, when the characteristic function of the empirical distribution function Characteristic function that closely approximates the distribution function When F n (x) is very close to F(x). In other words, the parameters of the characteristic function of the empirical distribution function are also close to the characteristic function of the distribution function.

[0540] The present invention uses the moment method or the maximum likelihood method to estimate the parameters of the empirical distribution function. The results of the two algorithms are consistent.

[0541]

[0542] because Revised Estimates is 2 The minimum mean square error estimator and the optimal unbiased estimator are usually used in practical applications. As σ 2 The estimated amount.

[0543]

[0544] In order to unify the dimensions of the calculation,

[0545]

[0546] The mean and variance of the water quality monitoring data are obtained directly. The monitoring samples are regarded as random variables. When the capacity is large enough, the Lindberg-Feller conditions expressed in equations (25) and 26 are satisfied, that is, the necessary and sufficient conditions for applying the central limit theorem are satisfied. The central limit theorem holds and the random variable obeys the normal distribution.

[0547] 4. Determination of the convergence of the empirical distribution function approximation distribution function

[0548] According to Glivenko's theorem, when the sample size is greater than a certain value n, in the interval where the sample size is greater than n, its empirical distribution function F n (x) is very close to the distribution function F(x), and the parameters of its distribution function are also very close. According to the uniqueness theorem of characteristic function, the normal empirical characteristic function is constructed and the normal distribution characteristic function The first and second moments of should also be very close. This is the theoretical basis for using the empirical distribution function to approximate the distribution function parameters.

[0549] From formula (29) is the optimal unbiased, consistent, effective, minimum mean square error estimator, so we use As σ 2 The estimator of the sample capacity is used to calculate the sample capacity value. Taking n as the sample capacity and q as the number of sequence calculation groups, the mean estimator is In order to make the error estimation results comparable, the dimensions should be the same, so we verify σ 2 The error is estimated using the standard deviation The distribution function of the nth group is F n (x) starts with Construct the sequence empirical distribution function F for the parameters (i = n, n + 1, n + 2, ..., n + q-1) i (x)(i=n,n+1,n+2,…,n+q-1), the absolute value of the relative error between the average value of the sequence calculation results and the calculation result when i=n+q-1 is taken as the convergence condition,

[0550]

[0551] Based on the principle of practical impossibility of low-probability events and the idea that low probability has no effect on the screening conclusion, it is considered that the model converges if the relative error w satisfies the sequence calculation of the conditions specified in Table a.

[0552] Table a Convergence conditions calculated according to probability sequence

[0553]

[0554] This convergence condition is essentially: using the empirical distribution function parameter estimate of the i=n+q-1 screening index of the water quality sample and The standardized distribution function F is obtained i The probability value of (x), that is, the approximate probability, is the average value of the empirical distribution function of q samples with a capacity greater than or equal to n. The absolute value w of the relative error between them satisfies formula (30).

[0555] After the conditions set in Table a are met, it is considered that the Grevenko condition (24) and the Lindberg-Feller condition (25) and (26) are established, that is, the sample size fully meets the conditions of the screening model. According to the Lyapunov theorem in formula (27), the probability can be calculated using the standard normal distribution function Φ(x).

[0556] Through a large number of calculations, the inventors believe that as long as the above conditions are met, the sequence calculation with q=5 can satisfy the Grevenko condition (24) and the Lindberg-Feller condition (25) and (26), that is, the distribution function parameters and sample capacity required for probability screening calculations.

[0557] 6. Algorithm for risk factor thresholds in the probabilistic screening model

[0558] When the sample size n meets the conditions for building a probability model, formula (23) is established, and the values ​​of each monitoring indicator are obtained by formula (29): The q-group sequence is calculated by the following formula:

[0559]

[0560]

[0561] When screening risk factors, set the risk probability threshold to Φ β , the threshold of water quality risk factor is

[0562] β=Φ -1 (Φ β ) (32)

[0563] Then the data exceeding the threshold β of the water quality risk factor in the monitoring data matrix X is screened by the following formula:

[0564]

[0565] Thus, the water quality risk factor data exceeding the probability threshold limit can be obtained.

[0566] 7. Standard normal distribution function Φ(x) and its inverse function Φ -1 Algorithm for (x)

[0567] The algorithm of the normal distribution function uses the Pade approximation to construct the confluent hyperratio function, and the continued fraction recursive format method is used to establish the algorithm of the normal distribution function Φ(x), and it is obtained

[0568]

[0569] The accuracy of formula (34) can reach more than 6 decimal places, which is sufficient;

[0570] Inverse function algorithm of standard normal distribution function:

[0571] When screening risk factors, the water quality risk factor threshold β=Φ -1 (x), the present invention adopts the non-iterative algorithm proposed by Zheng Duo, that is, directly solves Φ using the following formula: -1 (x), the method is as follows:

[0572] Assume that the water quality risk factor threshold is β, Φ(β)=Φ0

[0573] Constructor g(x) = Φ(x) - Φ0

[0574] Obviously, g(β)=0

[0575] remember structure

[0576]

[0577] Formula (35) can directly obtain β = Φ -1 (x), the calculation accuracy can reach 6 decimal places, which is sufficient;

[0578] 8. Season and date of screening for probability of exceeding threshold

[0579] Since the four seasons rotate in a cycle every year, the inventor sets a risk probability threshold, uses formula (32) to screen out monitoring data exceeding the risk probability threshold, and uses the month of the sample corresponding to the monitoring date to estimate the frequency of the monitoring data exceeding the probability threshold near the date in the next stage.

[0580]

[0581] Formula (33) is to infer the frequency of sample monitoring indicators exceeding the threshold in a certain month in the next stage.

[0582] Screen the seasons and dates with probabilities exceeding the probability threshold, and use the random variables of the screening monitoring data exceeding the probability threshold to screen out the corresponding seasons and dates. The probability threshold is set by the user according to needs and entered in the data file.

[0583] The following method is used to solve the water quality risk factor threshold from the probability threshold:

[0584] Set the risk probability threshold as Φ β , then the water quality risk factor threshold β=Φ -1 (Φ β ) is calculated by formula (35), and then formula (33) screens the measured value of the random variable X = βσ + μ to obtain the risk factor data that exceeds the probability threshold.

[0585] The inventors used three years of continuous monitoring sample data to conduct analysis, and the estimated results may have satisfactory accuracy. As the number of years increases and the amount of monitoring data increases, the estimated accuracy will significantly increase.

Claims

1. A method for screening water quality risk factors in cold region wetlands, characterized by: The following steps are involved: (I) Establish water quality monitoring data files (1) Establish a water quality monitoring sample content data table 1) The number of monitoring samples is the row, and the number of sample detection indicators is the column. There are n water quality monitoring samples, which is called sample capacity. Each sample has m monitoring indicators, which is the dimension of the random variable. The monitoring sample matrix is ​​represented as an n×m-order matrix: X=(X1,X2,...X j ,…X m ) T X j is a column vector, j = 1, 2, ..., m, or X=Ax (1) A is a constant matrix of n×m order consisting of monitoring data A=[a1,a2…,a m ], a j (j=1,2,…,m) is the column vector expressing the monitoring data, x=(x1,x2,…,x m ) T is an m-dimensional random variable; 2) To facilitate the understanding of the meaning and use of the monitoring sample data matrix, it is expressed in Excel format, and the final Excel table is the water quality monitoring sample data table; (2) Conditions for sample data capacity required for screening calculations Since this algorithm is a factor analysis + probability calculation for overall screening or seasonal screening, the sample size requirements are as follows: (1) The sample size of factor analysis is greater than the number of monitoring indicators, which is a condition that must be met if only factor analysis is performed; (2) The sample size of probability screening is greater than the number of monitoring indicators, generally not less than 30. When calculating seasonal screening, the sample size may be no less than 24; (3) Overall screening of water quality risks in cold-region wetlands regardless of time period For the overall screening of water quality risks in cold-region wetlands regardless of time period, you can skip the "Screening and Calculation by Season" item and directly conduct the screening calculation; (4) Three steps to establish the data file for screening calculations 1) Create calculation data file The data file name is named by the user. For easy searching, it is recommended to start with 0. The file determines the content calculated by the program. The file name cannot start with h or n. The files starting with h or n are classified files after the program recognizes them. The files starting with h are the data of the selected season or time period, and the files starting with n are the remaining data that have not been confirmed after recognition. If you do not want to divide the calculation into seasons or time periods, you can directly perform sequential calculations after starting the calculation. It is recommended to use text files, that is, files with the extension .txt; The first line of the data file contains the basic information and calculation requirements of the data. Spaces or commas can be used to separate the data. When performing overall screening calculations, the computer automatically performs overall screening without considering the effects of the ninth and tenth data; The first column is the row number, which is used to easily check the sample data number; the other columns are the content specified for the monitoring sample, and each sample occupies one row and is entered in the order specified for the record; Simplify the data below the second row: fill in the standard value row of the water quality monitoring sample data table to the nth row, and the standard value column to the mth column, and rewrite the "standard value" in the first row to the number "0". This row of data is the national standard value, which is the reference indicator for whether the water quality meets the standard; The last line of the data file is an auxiliary supplementary line. The line number of the supplementary line is the same as the first data of the first line. The numbers of the supplementary line correspond to the data of each column. 0 means that the value of the column is less than or equal to the standard value and meets the standard. 1 means that the value of the column is greater than the standard value and meets the standard. 2) Create a sample name file The sample name file is a special file name that cannot be changed, and its name is: sample name.txt; Paste the date column of the water quality monitoring sample data table into the sample name.txt file, and rewrite the date in the first row to the standard value, because the analysis results show the calculation results corresponding to the standard value in this position; 3) Create an element name file The element name file is a special file name that cannot be changed, and its name is: element name.txt; Paste the first row of the water quality monitoring sample data table into element name.txt, delete the sample number and date items, arrange each element in rows, and the calculation results correspond to the input element name; (5) After creating three data files, you can enter the screening calculation process. The "selection and identification of data files" in the calculation process is to establish preliminary files for the screening process, including the process of "identifying data by season"; After this step is completed, if you want to perform the "Screen and Calculate by Season" project, click "Screen and Calculate by Season", otherwise skip the "Screen and Calculate by Season" project and click "Create Data File" to proceed to the next step; (II) Specific steps of water quality screening calculation There are 17 steps in the water quality screening calculation, which are expressed in the form of a menu. The calculation process is carried out in the order of the menu. Step 17 "Verify the validity of the probability screening model" is independent and can be run separately as long as the data file is completed. The specific order of the menu content is as follows: (1) "Data file selection and identification" Click the "Selection and Identification of Data Files" item, and the computer will identify the data file.txt and sample name.txt. This is the starting step of the calculation, including the generation of the initial files for factor analysis and probability screening during the overall screening process; 1) "Overall screening calculation" During the overall screening calculation, the computer automatically ignores the effects of the ninth and tenth data and directly enters the screening calculation; 2) "Data identification by season" The ninth data in the first line of the data file is the starting month of the selected period, and the tenth data is the number of months of the selected period. At this time, the calculation data file for seasonal screening is generated as h01.txt, the sample name file is generated as sample name_h01.txt, and the element name file is generated as element name_h01.txt; (2) "Screening and calculation by season" Click "Calculate by Season Screening" to calculate. The computer uses the data in the h01.txt file to generate a file ending with _h01.txt to create the initial file for factor analysis and probability screening. (3). "Create data file" Click "Create Data File" and the computer will automatically generate a standardized monitoring data matrix; (4) "Principal component score" Click "Principal Component Score" to reduce the dimension of the observed data according to the cumulative variance contribution rate given by the fourth data in the first line of the data file; (5) "Direct calculation of scores by main factors" Click "Calculate score directly based on main factors" and let λ j =1, (j=1,2,…,m) to obtain the initial factor loading matrix A. The same as the principal component score, use formula (10) to calculate the score coefficient, and form the score matrix with the first k columns of the initial factor loading matrix A. Use the factor score formula (21) to convert F ji Change to A ji , that is, find the main factors and directly calculate the scores, and then sort the scores; (6) "Directly score the sample without rotating the factors" (7) "Variance maximal orthogonal rotation" (8). "Score after factor rotation" Click "Factor Rotation Score" to build the score matrix after factor rotation, and calculate the sorted factor analysis score and then use formula (20) to calculate the score function: F=Bx=A T R -1 x (20) Use formula (10) to calculate the score coefficient, and use factor score formula (21) to calculate the score, that is, to obtain the score after factor rotation, and then sort the scores; (9) "BARTLETT test" (10) "MSA test" (11) "KMO test" (12) "Risk factor probability screening" (13) "Predict the elements exceeding the standard in each month of next year" Click "Predict the elements exceeding the standard in each month of next year" to calculate the possible risk probability and date of the sample monitoring indicator by formula (31) and sample name.txt and element name.txt. Predict the frequency of the possible months exceeding the specified probability threshold in the next year by formula (36): Estimate the possible months according to the frequency order of pollution indicators; (14). "Draw a curve" Click "Draw Curve" to draw the scree plot of principal component analysis and factor analysis, and various factor score curves; (15). "Data summary" 1) Data summary of overall screening calculations Click "Data Summary" to directly summarize all calculation results in the "Data Summary_Data File.txt" file for easy viewing; 2) Summary of data by season screening Click "Data Summary" to directly summarize all calculation results in the "Data Summary_h Data File.txt" file for easy viewing; (16). "Reader" Click "Reader" to display all files in the displayed directory. Click the file name you want to view to view all data files generated during the screening process; (17). "Verifying the effectiveness of the probabilistic screening model".

2. The method for screening water quality risk factors in cold region wetlands according to claim 1, characterized in that: The step (i) establishes a water quality monitoring data file: (1) establishing a water quality monitoring sample content data summary table; and 2) supplementing the monitoring sample matrix with the following content: A. Number the samples on the left side of the monitoring data matrix; B. Add three rows at the top of the matrix, the first row is the name of the monitoring indicator, the second row is the number of the monitoring indicator, and the third row is the standard value of the monitoring indicator; C. Add a column on the right side of the matrix for the monitoring date of the sample; D. Add a row at the bottom of the matrix to determine whether it is greater than or less than the specified standard value.

3. The method for screening water quality risk factors in cold region wetlands according to claim 1, characterized in that: In the step (1) establishing a water quality monitoring data file and (4) establishing three data files for screening calculation, 1) establishing a calculation data file, the data files are filled in with the following data: The first data is the number of samples + 1, which is used for factor analysis and probability screening; The second data is the number of indicators that the sample participates in the calculation + 1, which is used for factor analysis and probability screening; The third data is to require the cumulative variance contribution rate to be taken to decimal places for factor analysis; The fourth data is the required cumulative variance contribution rate, which is used for factor analysis; The fifth data is the significance level α value, which is used for probability screening; The sixth data is the absolute value of the hidden factor loading is less than the specified number, which is used for factor analysis; The seventh data is the probability threshold for screening risk factors. If the data is less than or equal to 0, only factor analysis is performed without possible probability analysis, which is used for probability screening. The eighth data is the number of probability series calculation groups to determine the availability of parameter estimates, which is used for probability screening; The ninth data is the starting month of the selected period, which is used for factor analysis and probability screening; The tenth data is the number of months of the duration of the selected period, which is used for factor analysis and probability screening.

4. The method for screening water quality risk factors in cold region wetlands according to claim 1, characterized in that: In the specific steps of the water quality screening calculation in step (ii), (3) "establishing data files" includes data format files corresponding to principal component analysis, factor analysis and probability calculation, including the covariance matrix V of the factor analysis generated by formula (3): D(X)=D(Ax)=AD(x)A T =AA T =V (3) The correlation matrix R generated by equation (4) has the following elements: Find the eigenvalue λ of R i (i=1,2,…,m) and sorted according to formula (5): λ1≥λ2…λ m ≥0 (5) The cumulative variance contribution rate is calculated by formula (9): When the cumulative variance contribution rate is greater than or equal to the cumulative variance contribution rate given by the fourth data in the first line of the data file, it indicates that the first k factors can approximately express the data content of all factors, that is, the information of m factors can be expressed by k factors, achieving the purpose of dimensionality reduction.

5. The method for screening water quality risk factors in cold region wetlands according to claim 1, characterized in that: In the specific steps of the water quality screening calculation in step (ii), (4) "principal component score" includes: According to formula (7), the first principal component is the largest, and the other k-th principal components are arranged in order: Use the first k principal components to form a score matrix: X k = (u1, u2, …, u k ) Y xor X k =(X1,X2,…,X k ) T The score coefficients of the first k principal components are: Using the factor score formula (21), we can ji Change to X ji F ji Change to X ji ,Right now The principal component scores are calculated and then X is used to sort the principal component scores.

6. The method for screening water quality risk factors in cold region wetlands according to claim 1, characterized in that: The specific steps of the water quality screening calculation in step (ii) (6) "direct score of unrotated factor samples" include: Click "Directly score the factor sample without rotation" to obtain the factor loading matrix A by formula (15). The same as the principal component score, use formula (10) to calculate the score coefficient, and form the score matrix by the first k columns of the loading matrix A. Use the factor score formula (21) to convert F ji Change to A ji , that is, to obtain the direct scores of the unrotated factor samples, and then sort the scores; Note: Without rotating the factor sample, the direct score should apply the first k eigenvalues ​​λ of the correlation matrix j (j=1,2,…,m), and the main factor directly calculates the score j =1,(j=1,2,…,m).

7. The method for screening water quality risk factors in cold region wetlands according to claim 1, characterized in that: In the specific steps of the water quality screening calculation in step (ii), (7) "variance maximum orthogonal rotation" includes: Click "Variance Maximum Orthogonal Rotation", build an orthogonal matrix according to formula (16), and calculate the rotation angle according to formula (18): Will Substituting into equation (17), the factor principal axis completes a rotation; By changing the values ​​of the subscripts j and k in pairs, the next rotation can be performed. When the common factor is greater than 2, one cycle needs to rotate Generally, one cycle rotation cannot achieve the best result. It takes multiple cycles of rotation until the total variance of the previous and next two cycles does not change much before the rotation can be completed. After rotation, the contribution is increased as much as possible, and the score matrix of the maximum variance orthogonal theory is obtained.

8. The method for screening water quality risk factors in cold region wetlands according to claim 1, characterized in that: The specific steps of the water quality screening calculation in step (ii) (9) "BARTLETT test" include: Click "Bartlett test" and perform the Bartlett test on the m-dimensional independent normal distribution according to the condition of formula (22). If the calculation result shows "no rejection of homogeneity of variance" on the screen, it means that the data passes the Bartlett test. Let n represent the capacity of the monitoring sample and m represent the dimension of the random variable. The Bartlett test is performed on the m-dimensional mutually independent normally distributed random variables using formula (22): The screen displays "Do not reject the requirement of homogeneity of variance"; The conclusion is that the test passes, otherwise it means that the standard deviation does not meet the homogeneity condition and the data needs to be adjusted.

9. The method for screening water quality risk factors in cold region wetlands according to claim 1, characterized in that: The specific steps of the water quality screening calculation in step (ii) (12) "risk factor probability screening" include: Click "Risk Factor Probability Screening", and the computer will calculate the water quality risk factor threshold β according to the probability threshold specified by the seventh data in the first line of the data file using formula (32): β=Φ -1 (F β ) (32) The monitoring index X exceeding the probability threshold is screened in the matrix X by formula (33): The results of the probability analysis of the sample monitoring data are obtained, including the probability that several detection indicators may exceed the probability threshold, and the sample time when the indicator exceeds the probability threshold is given.

10. The method for screening water quality risk factors in cold region wetlands according to claim 1, characterized in that: The specific steps of the water quality screening calculation in step (ii) (17) "verifying the validity of the probability screening model" include: Click "Verify the validity of the probability screening model" to calculate the parameter estimate of the empirical distribution function according to the q group sequence using formula (29): The parameter estimate of the empirical distribution function calculated using q groups of sequences is given by equation (31): The validity of the probability model is verified using formula (30) and Table a: If the conditions in Table a are met, the screen displays: The probability model meets the convergence conditions! Otherwise, it means that the sample size is insufficient and should be increased; For seasonal screening, if the sample size is not less than 24, it can be calculated directly without verification, and the verification process has no impact on the entire analysis process.

Images