Multi-index soft measurement method based on slow feature over-complete width learning system

Abstract

The invention provides a multi-index soft measurement method based on a slow feature over-complete width learning system. The multi-index soft measurement method is used for simultaneously measuring five-day biochemical oxygen demand (BOD5) and chemical oxygen demand (COD) in water quality indexes. The existing BOD5 and COD measuring method has high requirements on equipment and operation, is difficult to meet the requirement of real-time monitoring and is mostly used for measuring a single index. According to the SFA-OBLS model, a plurality of sensors are installed, easy-to-measure parameters are collected and input into an online monitoring system, and meanwhile, BOD5 and COD values are predicted. The model is essentially used for predicting indexes, and the precision is slightly inferior to that of traditional direct measurement. According to the method, a slow feature window and an over-complete window jointly form a feature window, slowly-varying dynamic feature components and independent components are mapped to an enhancement window, and nonlinear expansion is carried out to enhance features. The above three windows enable the model to effectively process several characteristics of water quality data, which easily influence monitoring precision, and effectively improve monitoring efficiency and precision.

Description

Technical Field

[0001] The invention belongs to the technical field of soft measurement. Background Art

[0002] The current five-day biochemical oxygen demand (BOD 5 ) measurement methods have problems such as long detection cycle, complex operation, and high equipment requirements, which are difficult to meet the needs of real-time monitoring. The chemical oxygen demand (COD) requires the use of strong oxidants for complex operations under high temperature and high pressure conditions, and also requires the treatment of toxic waste liquids. The laboratory safety requirements are high and it is difficult to measure. Traditional soft measurement technologies mostly measure a single indicator and it is difficult to solve the problem of simultaneous measurement of multiple indicators. Technical difficulties such as the accuracy and synchronization of data acquisition, high-dimensional data processing, nonlinear relationship modeling, and data noise and missing value processing increase the complexity of multi-indicator soft measurement. In response to these difficulties, the present invention proposes a multi-indicator soft measurement method based on a slow feature overcomplete width learning system (SFA-OBLS), which estimates BOD in real time through an efficient and accurate intelligent model. 5 Under the premise of ensuring high precision, this model uses advanced algorithm architecture to realize simultaneous measurement of multiple indicators, which significantly improves the measurement speed and accuracy.

[0003] Biochemical oxygen demand (BOD) 5 ) is an important indicator for evaluating the content of organic pollutants in water and is widely used in water quality monitoring in sewage treatment plants. 5 It indicates the amount of dissolved oxygen consumed by microorganisms to decompose organic matter in water within five days at a specific temperature (usually 20°C). This indicator can reflect the concentration of organic pollutants in water and is directly related to the health of the water body. 5 The measurement method is time-consuming and complex, and the experimental steps are cumbersome, making it difficult to monitor in real time. Since it takes five days to complete the cultivation and reaction process of microorganisms, there is a significant time lag in the measurement results, making it difficult to achieve real-time monitoring. This lag makes it impossible for sewage treatment plants to adjust the treatment process in a timely manner when dealing with sudden pollution incidents, affecting the treatment effect. In addition, the measurement process requires strict control of multiple factors such as temperature, initial dissolved oxygen concentration, and microbial inoculum size, which is prone to operational errors and leads to inaccurate measurement results. At the same time, the experimental equipment has high requirements and high maintenance costs, and requires professional technicians to operate and manage it, which increases operating costs and human resource investment.

[0004] Chemical oxygen demand (COD) is an important indicator for assessing the content of organic and some inorganic pollutants in water bodies, and is widely used in the fields of sewage treatment and water quality monitoring. COD indicates the amount of oxygen required by organic matter in a water sample under the action of a strong oxidant. This indicator can quickly reflect the degree of organic pollution in a water body and provide an effective monitoring basis for sewage treatment plants. However, the traditional COD measurement method is complicated to operate, needs to be carried out under high temperature and high pressure, and involves the use of strong oxidants. Toxic waste liquids must be handled with caution during the operation, which poses a safety hazard. The measurement results are often also affected by factors such as reaction conditions and reagent quality, making it difficult to guarantee the accuracy and repeatability of the results. At the same time, the detection cycle of traditional methods is long, and it is difficult to meet the needs of real-time monitoring, especially when responding to sudden pollution incidents, and it is impossible to respond quickly, affecting water quality management and treatment effects.

[0005] In order to overcome the above difficulties, soft sensor technology has gradually become a research hotspot in recent years. Soft sensor technology predicts BOD by establishing a mathematical model and using other parameters that are easy to measure. 5 Indicators that are difficult to measure directly, such as COD, can be used to achieve real-time prediction of difficult-to-measure indicators. Compared with traditional methods, soft measurement technology has the advantages of fast detection speed, simple operation, and low equipment requirements, which can significantly improve monitoring efficiency and response speed. However, most traditional soft measurement technologies can only measure a single indicator, and it is difficult to monitor multiple key indicators at the same time. When measuring multiple indicators at the same time, it will increase the complexity of the model, significantly increase the demand for computing resources, and increase the complexity of model training and optimization, as well as difficulties in real-time and response speed. If multiple indicators cannot be measured at the same time, the monitoring data may be incomplete and the response to water quality changes may be delayed, thereby affecting the optimization of the sewage treatment process and the overall treatment effect. In addition, when processing complex water quality indicator data, the existing soft measurement methods are difficult to handle their slow-changing characteristics, non-Gaussian characteristics, and nonlinear characteristics that will significantly affect the prediction accuracy, resulting in low prediction accuracy and stability, and there are defects in accuracy compared with traditional monitoring methods.

[0006] In view of these difficulties, the present invention proposes a multi-index soft measurement method based on a slow feature overcomplete width learning system (SFA-OBLS). The model uses an overcomplete slow feature width learning system (SFA-OBLS) to process the slow features, non-Gaussian features and nonlinear features of water quality data, effectively solving the defects of traditional soft measurement models in data feature processing. The model can estimate multiple key indicators (BOD5 and COD are measured in the invention) in real time through an efficient and accurate intelligent algorithm architecture. Under the premise of ensuring high accuracy, it can easily achieve simultaneous measurement of multiple indicators and achieve efficiency that satisfies operators, meet the efficient and real-time requirements of modern sewage treatment and water quality monitoring, and provide strong technical support for real-time monitoring and management of multiple indicators. Summary of the invention

[0007] This paper proposes a multi-index soft sensing method based on slow feature overcomplete width learning system (SFA-OBLS) to solve the BOD 5 The measurement method is time-consuming, complex to operate and requires high equipment. The COD measurement method has high laboratory safety requirements and is difficult to measure. In addition, traditional soft measurement methods are not easy to measure multiple indicators at the same time. The SFA-OBLS model uses soft measurement technology to collect easy-to-measure parameters in water samples in real time to achieve real-time monitoring and efficient management of BOD5 and COD. 5 The simultaneous real-time prediction of BOD and COD significantly shortens the measurement time, eliminates the need for manual intervention, and reduces operational complexity and human errors. The prediction model constructed by SFA-OBLS can effectively handle the inevitable slow-changing, non-Gaussian and nonlinear characteristics of water quality indicators that affect the accuracy of the model, and provide high-precision prediction results. The SFA-OBLS model can also achieve real-time and rapid prediction of actual BOD on the basis of accuracy comparable to traditional measurement methods. 5 and COD value.

[0008] The following scheme was used to complete the BOD 5 And the prediction of COD:

[0009] The present invention is to construct BOD 5 The prediction model of BOD and COD is the research object, aiming to improve the prediction of BOD and COD through fast and high-precision data processing methods. 5 and COD detection accuracy and efficiency. 5 The values ​​of COD and BOD are preprocessed and then input into the SFA-OBLS model for training. The SFA-OBLS model uses its fast learning and powerful data processing capabilities to quickly establish a water quality prediction model based on feature nodes and enhanced nodes. Then, the soft sensor is used to measure the real-time and easy-to-measure parameters for BOD prediction. 5 The present invention can predict BOD in real time. 5 The real-time detection of COD is mainly divided into three parts: preprocessing stage, model training stage and online application stage.

[0010] Preprocessing stage:

[0011] Step 1: First, obtain the actual data set of the sewage plant and divide the data set into a training data set and a test data set; then determine the input and output process variables of the detection model; the input variables are represented by array X; the input variables are: influent suspended solids mass concentration (g / m 3 ), mass concentration of biodegradable organic matter (g / m 3), mass concentration of particulate inert organic carbon (g / m 3 ), active heterotrophic bacteria mass concentration (g / m 3 ), active autotrophic bacteria mass concentration (g / m 3 ), particle mass concentration produced by biological attenuation (g / m 3 ), dissolved oxygen mass concentration (g / m 3 ), nitrate nitrogen mass concentration (g / m 3 ), ammonia nitrogen mass concentration (g / m 3 ), soluble biodegradable organic nitrogen mass concentration (g / m 3 ), alkalinity molar concentration (mol / m 3 ), total nitrogen concentration (g / m 3 ), mass concentration of chronically biodegradable organic matter (g / m 3 ); the output variable is represented by array Y; the output variable is the five-day biochemical oxygen demand (BOD 5 )(g / m 3 ) and chemical oxygen demand (COD) (g / m 3 ); the input and output variable names of the test set are the same as those of the training set; the output indicators are predicted by inputting 13 measurable water quality data;

[0012] Step 2: Preprocess the training data set and represent the training set with the symbol X, that is, X∈R NxM ; R is a real number set, N is the number of samples, M is the number of process variables; M-dimensional input signal X = [x 1 (t),x 2 (t),...,x M (t)], where t is the sample;

[0013] Network training phase:

[0014] Step 3: The goal of the slow feature extraction of the SFA-OBLS model is to find a mapping function Where m represents the number of slow features, by inputting the input signal X into The mapping function makes the characteristic component sf=g(X) change as little as possible over time. The slow characteristic information after mapping is:

[0015]

[0016] Where X is the 13 input variables of the training data set: mass concentration of suspended solids in water (g / m 3 ), mass concentration of biodegradable organic matter (g / m 3 ), mass concentration of particulate inert organic carbon (g / m 3 ), active heterotrophic bacteria mass concentration (g / m 3), active autotrophic bacteria mass concentration (g / m 3 ), particle mass concentration produced by biological attenuation (g / m 3 ), dissolved oxygen mass concentration (g / m 3 ), nitrate nitrogen mass concentration (g / m 3 ), ammonia nitrogen mass concentration (g / m 3 ), soluble biodegradable organic nitrogen mass concentration (g / m 3 ), alkalinity molar concentration (mol / m 3 ), total nitrogen concentration (g / m 3 ), mass concentration of chronically biodegradable organic matter (g / m 3 );where sf is the slow feature set extracted by the SFA-OBLS model through the slow features, and sf is obtained according to the change speed from slow to fast 1 ,sf 2 ,...,sf m m slow features; g(X) is The set of mapping functions of ; [·] T The formula represents the transposition; g(x 1 (t)),g(x 2 (t)), and g(x m (t)) is X = [x 1 (t),x 2 (t),...,x M (t)] through the mapping function The first, second and mth slow feature information is obtained;

[0017] Step 4: The slow feature window needs to find the first few features with the slowest change rate, so that the variance Δ(sf i ) is minimized, thus obtaining the following optimization problem:

[0018]

[0019] Among them <·> t is the slow feature sample at N time points The calculated average value; Represents the first-order difference information of the i-th slow feature, that is, Δ(·) is the value for sf i A natural measure of the rate of change; represent The square of

[0020] Step 5: Add the following constraints to the optimization problem:

[0021]

[0022] These three constraints are the zero mean constraint. <sf i > t =0, unit variance constraint and decorrelation constraints <sf i sf j > t =0, Among them, the zero mean constraint requires that the mean of the output signal is 0 in the time dimension, ensuring that the optimization problem has only one solution; the unit variance constraint requires that the mean of the square of the output signal is equal to 1 in the time series; the decorrelation constraint requires that different components are independent of each other, showing that they carry different information; among them <·> t is a slow feature sample at N time points The calculated average value;

[0023] Step 6: The first-order difference of discrete data indicates the speed of change:

[0024]

[0025] For each slow feature sf i and sf i-1 are represented by linear combinations of all input variables, through The mapping function is used to implement the mapping, namely:

[0026]

[0027] sf=WX (6)

[0028] Where W = [w 1 ,w 2 ,...,w m ] T is the load matrix that needs to be optimized through slow features, w i represents the load vector of the load matrix W, w i is the transpose of ; g(X) is The set of mapping functions, sf is the set of slow features extracted by the SFA-OBLS model through slow features; it is easy to prove that if the input variables are to satisfy the zero mean constraint, it is necessary to force the input variables to automatically remove the mean;

[0029] Step 7: Substitute the formula in step 6 into the formulas in steps 4 and 5 to obtain the slow feature optimization target in the linear case:

[0030]

[0031] The formula It is represented by the symbol A, which is the covariance matrix of the first-order difference of the input X, that is, is the first-order difference of the input X; the formula <X X T >Symbol B represents the covariance matrix of input X, that is, B = <X X T > w i The load vector representing the load matrix W; Yes i The transpose of ; where <·> t is a slow feature sample at N time points The calculated average value; Represents the first-order difference information of the i-th slow feature, that is, Δ(·) is the value for sf i A natural measure of the rate of change; yes The square of , from the constraint condition in step 4, we can get the unit variance constraint <sf i 2 > t =1, that is

[0032] Step 8: Spheroidize X to eliminate the correlation between variables and obtain interpretable data; the spheroidization process involves singular value decomposition of the covariance matrix B:

[0033] B=UΛU T (8)

[0034] Where U is an orthogonal matrix composed of singular vectors; Λ is a rectangular diagonal singular value matrix, and the elements on the diagonal of Λ are the singular values ​​of B; T represents transpose;

[0035] Step 9: The original input X is sphericalized by singular value decomposition, and the correlation is removed, that is:

[0036] z=Λ -1 / 2 U T X (9)

[0037] Where U is an orthogonal matrix composed of singular vectors; Λ is a rectangular diagonal singular value matrix, and the elements on the diagonal of Λ are the singular values ​​of B; T represents transposition, and z is the result of spherification of X. <zz T > t =I m , I m is the m-order identity matrix; through the singular value decomposition operation, the original optimization problem is transformed into finding a matrix P, according to the unit variance and decorrelation constraints in step 5, so that sf = Pz and <sfsf T > t =I m ,Will <zz T > t =Im Substituting sf=Pz, we obtain:

[0038] PP T =I m (10)

[0039] T represents transpose; P is an orthogonal matrix; we get Take the minimum solution, the covariance matrix of the first-order difference of z Perform singular value decomposition and solve, is the first-order difference of z;

[0040]

[0041] Where P is an orthogonal matrix composed of eigenvectors; Ω is a rectangular diagonal singular value matrix, and the diagonal elements are the change values ​​of each slow feature; T represents transpose;

[0042] Step 10: The final parameter coefficient matrix W is shown as follows:

[0043] W=PΛ -1 / 2 U T (12)

[0044] Where U is an orthogonal matrix composed of singular vectors; Λ is a rectangular diagonal singular value matrix, and the elements on the diagonal of Λ are the singular values ​​of B; T represents transpose; P is an orthogonal matrix composed of eigenvectors;

[0045] Step 11: Find the slow feature sf:

[0046] sf=WX=PΛ -1 / 2 U T X=Pz (13)

[0047] Step 12: Input the SF obtained from the slow feature analysis into the feature layer window of the SFA-OBLS network for further feature extraction, and form the i-th group of feature nodes SF through the feature mapping function i :

[0048]

[0049] in is the SFA-OBLS network feature mapping function, the activation function of the feature layer is tansig, and the mapping function is used to extract the centralized variable information of the sewage data and map it to the enhancement window; p is the number of feature nodes in the feature layer, sf is the slow feature set extracted by the SFA-OBLS model through slow features, The weight matrix mapped to the feature node of the i-th window, The bias vector mapped to the feature node of the i-th window, and is randomly generated; to prevent weight Affected by randomness, a sparse autoencoder is used to fine-tune it, where the sparse matrix W s Obtained by the following formula:

[0050]

[0051] Where argmin means finding W that minimizes the objective function s ,λ 1 is the L1 regularization parameter, λ 1 =2-8, after sparse feature node SF i =sfW s ;

[0052] Step 13: Output SF to the feature layer p Repeat the previous step p times to obtain;

[0053] SF p =(SF 1 ,SF 2 ,…,SF p ) (16)

[0054] Among them SF 1 , S.F. 2 and SF p Represent the 1st, 2nd and pth slow feature windows respectively;

[0055] Step 14: The purpose of the overcomplete window of the SFA-OBLS model is to obtain S with an overcomplete feature representation:

[0056] S=XL (17)

[0057] The key to obtaining S is to obtain the mixed component L; X is the input variable, L is the mixed component matrix;

[0058] Step 15: Calculate the cumulative generating function of x:

[0059]

[0060] where φ x (·) is the cumulative generating function, σ is a real variable used to adjust and analyze the characteristics of the variable x, is the moment generating function MGF of the input variable x, which means The expected value of Used to describe the characteristics of x, σ is a continuous real number variable used to construct the function x, σ T represents the transpose of σ;

[0061] Step 16: By calculating The coefficient of the Taylor series expansion of the cumulant generating function at is used to find the second-order cumulant of x. The Hessian estimate of the second-order cumulant of x is as follows:

[0062]

[0063] in Represents the zero vector, all components are 0, → represents the direction of the vector, used to distinguish between vectors and scalars, the special one is There is no direction; cov(x) is the covariance matrix, which represents the covariance between different components of the input variable x. represents the cumulant generating function φ x (σ) The second-order partial derivative matrix at ;

[0064] Step 17: Expand σ to non-zero vectors when When , the generalized covariance matrix C is obtained by accumulating the Hessian matrix of the generating function x (σ);

[0065]

[0066] in, represents the cumulant generating function φ x (σ) The second-order partial derivative matrix at , σ is a real variable used to adjust and analyze the characteristics of the variable x, which obeys the Gaussian distribution; Represents the input variable x and its transpose x T Multiply Expected value; express Expected value; is used to describe the characteristics of x, σ is a real number variable used to adjust and analyze the characteristics of variable x, σ T represents the transpose of σ, β x (σ) represents the first-order derivative of the cumulant generating function of x under σ, that is cov(x) is the covariance matrix;

[0067] Step 18: C x (σ) is further expressed as:

[0068] C x (σ)=LC S (d)L T (twenty one)

[0069]

[0070] Among them C x(σ) is the generalized covariance matrix, σ is a real variable used to adjust and analyze the characteristics of variable x, L is the mixture component matrix, L T is the transpose of L, C s (d) represents the generalized covariance matrix of the independent source components, represents the cumulant generating function φ S (d) The second-order partial derivative matrix at , S is an overcomplete feature representation, d = L T σ; β x (σ) represents the first-order derivative of the cumulant generating function of x under σ, that is is an overcomplete feature representation S and The expected value of the product, yes The expected value of Used to describe the characteristics of S, S T represents the transpose of S, where d = L T σ, β S (d) represents the first derivative of the cumulant generating function of S under d;

[0071] Step 19: Derive the generalized covariance matrix C of the independent sources s (d):

[0072]

[0073]

[0074] Among them, C s (d) represents the generalized covariance matrix of the independent source components, represents the cumulant generating function φ S (d) The second-order partial derivative matrix at , S is an overcomplete feature representation, S T represents the transpose of S, d = L T σ, L is the mixture component matrix, L T is the transpose of L, σ is a real variable used to adjust and analyze the characteristics of variable x, is an overcomplete feature representation S and its transpose and The expected value of the product, yes The expected value of Used to describe the characteristics of S, β S (d) represents the first-order derivative of the cumulant generating function of S under d, β S (d) T Represents β S (d) is the transpose of, cov(S) is the covariance matrix of S, is the first derivative of the cumulant generating function of x, Indicates S and The expected value of the product, express The expected value of T represents the transpose of d;

[0075] Step 20: Due to the independence of S, the generalized covariance C of the independent source components s (d) is a diagonal matrix; therefore, C x (σ) is further expressed as:

[0076]

[0077] Where k is the number of overcomplete features S; ω i (σ)=[C s (L T σ)] ii represents the i-th independent source component S i The covariance of i represents the i-th mixed component, Indicates l i The transpose of is a single atom of the mixing matrix L;

[0078] Step 21: The generalized covariance matrix belongs to the subspace W; the subspace spanned by any number of covariance matrices is either equal to W or a subset of W; select a sufficiently large number s of generalized covariance matrices to span the subspace W = [H 1 ,H 2 ,...,H s ]; Set the parameter s to 5 times the parameter k; Based on the given non-zero vector σ 1 ,σ 2 ,...,σ s , construct an equivalent matrix H j =C x (σ j ),j∈{1,2,...,s};

[0079] Step 22: Based on the subspace W, the overcomplete window is semi-positive regularized to complete the single atom estimation Then the mixed component l is obtained i , the semi-positive definite problem of single atom estimation can be specifically expressed as:

[0080]

[0081] in It is an atom Based on the optimal solution obtained by semidefinite programming, argmax is used to find the matrix B that maximizes <G,B> in the space Span, B∈W=Span{H 1 ,H 2 ,...,H s} indicates that the matrix B is in the linear space Span{H 1 ,H 2 ,...,H s}; the goal is to find a matrix B that maximizes the inner product <G,B>; Span is defined by the matrix H 1 ,H 2 ,...,H s The generated linear space, the matrix G determines The reduction direction of atoms; (In the process of independent component analysis, the objective function is used to measure the independence of signals; In order to maximize the independence, the gradient of the objective function is used to update the estimation matrix; The matrix G is determined by the gradient of the objective function; For a solution vector w in the optimization problem, assuming that we use negative entropy J as the objective function, the matrix G is expressed as the outer product of the gradient vector, or directly as the gradient matrix: The gradient here Represents the objective function J relative to the signal w T The rate of change of X, where X is the input signal matrix;)

[0082] Step 23: Constrain the semi-definite problem:

[0083]

[0084] Where B∈W=Span{H 1 ,H 2 ,...,H s} indicates that the matrix B is in the linear space Span{H 1 ,H 2 ,...,H s}; the goal is to find a matrix B that maximizes the inner product <G,B>, H 1 ,H 2 ,...,H s is the matrix constructed in step 21. Span is defined by the matrix H 1 ,H 2 ,...,H s The generated linear space, Tr(·) represents the sum of the elements on the diagonal of the matrix B; when the optimal solution exists and is unique, It is an atom The optimal solution based on semidefinite programming; G in step 22 determines The reduction direction of atoms, that is, the selection of the best atoms, is determined by the matrix G; (In the process of independent component analysis, the objective function is used to measure the independence of signals; In order to maximize independence, the gradient of the objective function is used to update the estimation matrix; The matrix G is determined by the gradient of the objective function; For a solution vector w in the optimization problem, assuming that we use negative entropy J as the objective function, the matrix G is expressed as the outer product of the gradient vector, or directly as the gradient matrix: The gradient here Represents the objective function J relative to the signal w T The rate of change of X, where X is the input signal matrix;)

[0085] Step 24: Matrix G is W (t) The orthogonal complement of When (t) Yes H (t) The first left singular vector of is u (t) The transpose of , then

[0086] Step 25: Express the constraint set in the form of a set K = {B∈W, Tr(B) = 1, B≥0}:

[0087] K={ <B,F j >=0,j∈[mk],Tr(B)=1,B≥0} (28)

[0088] Where B∈W=Span{H 1 ,H 2 ,...,H s} indicates that the matrix B is in the linear space Span{H 1 ,H 2 ,...,H s}; the goal is to find a matrix B that maximizes the inner product <G,B>, with the constraint that B∈W=Span{H 1 ,H 2 ,...,H s} is equivalent to B being orthogonal to any matrix in the complementary basis of the subspace W; N(W) = {F 1 ,F 2 ,...,F m-k} represents the complementary basis of the orthogonal basis of subspace W, where m = M 2 / 2, Tr(·) represents the sum of the elements on the diagonal of the matrix B;

[0089] Step 26: In actual industrial processes, since noise is inevitable, atomic estimation can often only be performed in an approximate form. To solve this problem, relax the original hard constraints based on the original constraints:

[0090]

[0091] K={Tr(B)=1,B≥0} (30)

[0092] Among them B * is an atomic estimate, argmax is used to find the matrix B that maximizes <G,B> in the space Span, B∈W=Span{H 1 ,H 2 ,...,H s} indicates that the matrix B is in the linear space Span{H 1 ,H 2 ,...,H s}; the goal is to find a matrix B that maximizes the inner product <G,B>, and Span is defined by the matrix H 1 ,H 2 ,...,H s The generated linear space, matrix G determines B * The reduction direction of the atom; Use the gradient of the objective function to update the estimation matrix; The matrix G is determined by the gradient of the objective function; For a solution vector w in the optimization problem, assuming that we use negative entropy J as the objective function, the matrix G is expressed as the outer product of the gradient vector, or directly as the gradient matrix: The gradient here Represents the objective function J relative to the signal w T The rate of change of X, X is the input signal matrix; μ>0 is a regularization parameter; B∈W=Span{H 1 ,H 2 ,...,H s} means that B is orthogonal to any matrix in the complementary basis of the subspace W, N(W) = {F 1 ,F 2 ,...,F m-k} represents the complementary basis of the orthogonal basis of the subspace W, K = {Tr(B) = 1, B ≥ 0} is the set form representation of the constraint set, and Tr(·) represents the sum of the elements on the diagonal of the matrix B;

[0093] Step 27: For a solution vector w in the optimization problem, use negative entropy J as the objective function, and the matrix G is expressed as the outer product of the gradient vector, or directly as the gradient matrix: Here the gradient ▽J represents the objective function J relative to the signal w T The rate of change of X, X is the input signal matrix; the specified k-fold number of atoms is selected; then, the selected specified k-fold number of atoms are clustered into k clusters, denoted as D clust ; In each cluster, select a column as the representative atom. Here, the first column is selected by default. Finally, for Dclust The k atoms in are screened, and the specific steps of screening are as follows:

[0094] G clust =|D clust T D clust |-I k (31)

[0095] N um = argmin{max(G clust )} (32)

[0096] Where N um Indicates the minimum index among all maximum values, D clust T Yes D clust The transpose of G is to find the value such that max(G clust )The smallest N um , that is, D clust (N um ) represents the most independent atom, which can be used as an atom screened by the clustering algorithm; if D clust If the maximum values ​​corresponding to the remaining atoms in do not reach the preset threshold standard, they will be retained as new estimated atoms, and the threshold does not exceed the maximum value of the regularization parameter;

[0097] Step 28: If the total number of atoms a screened by the clustering algorithm is less than the preset number k, all the excellent atoms found so far are excluded from the search range to update the constraint set; secondly, the subspace W is updated to the space spanned by the atoms remaining after the estimated atoms are removed, and the atoms found by the clustering algorithm are added to the complement of the orthogonal basis of the subspace W; finally, the fast iterative shrinkage threshold algorithm FISTA and the optimization-maximization principle MMP are used to complete the estimation of the remaining atoms; the negative objective function can be further expressed as follows:

[0098]

[0099] The matrix G determines the reduction direction of atoms; The gradient here Represents the objective function J relative to the signal w T The rate of change of X, X is the input signal matrix; B∈W=Span{H 1 ,H 2 ,...,H s} indicates that the matrix B is in the linear space Span{H 1 ,H 2 ,...,H s}; the goal is to find a matrix B that maximizes the inner product <G,B>, μ>0 is a regularization parameter, N(W)={F 1 ,F 2 ,...,F m-k} indicates that the complementary basis B of the orthogonal basis of subspace W is an atomic estimate, m = M 2 / 2, t∈[a,k] represents the number of atoms currently obtained;

[0100] Step: 29: Estimate the remaining ka atoms; combined with the given constraints K = {Tr(B) = 1, B ≥ 0}, the gradient of the differentiable part of the objective function is:

[0101]

[0102] The matrix G determines the reduction direction of atoms; The gradient here Represents the objective function J relative to the signal w T The rate of change of X, X is the input signal matrix;), B∈W=Span{H 1 ,H 2 ,...,H s} indicates that the matrix B is in the linear space Span{H 1 ,H 2 ,...,H s}; the goal is to find a matrix B that maximizes the inner product <G,B>, μ>0 is a regularization parameter, N(W)={F 1 ,F 2 ,...,F m-k} indicates that the complementary basis B of the orthogonal basis of subspace W is an atomic estimate, m = M 2 / 2, t∈[a,k] represents the number of atoms currently obtained; the optimization-maximization principle will terminate the contraction process in advance when a suitable atom is found;

[0103] Step 30: Complete the atom estimation and obtain k atoms; then perform eigenvalue decomposition after each atom is obtained by matrixing, and select the eigenvector corresponding to the maximum eigenvalue as a mixing component l in the mixing component L i ; The dimension of the potential independent component S obtained based on the mixed component L exceeds the dimension of the number of input features, which makes the neurons of the SFA-OBLS network have more freedom when setting weights or biases;

[0104] Step 31: Get the cumulative output OF of m overcomplete windows m :

[0105] OF m =[OF1 ,OF 2 ,…,OF m ] (35)

[0106] OF 1 ,OF 2 and OF m Represent the 1st, 2nd and mth over-complete windows respectively;

[0107] Step 32: SFA-OBLS maps the slow-varying dynamic feature components reflecting the overall trend of process changes and the potential independent components with overcomplete feature representation to the enhancement window module through the slow feature window and the overcomplete window to achieve nonlinear expansion and obtain enhanced features; the output of the j-th enhancement node of the enhancement window can be expressed as E j , E j The expression is as follows:

[0108]

[0109] Where n represents the total number of n enhancement nodes in the enhancement layer, ξ j As a nonlinear activation function, sigmoid is selected as the activation function of the enhancement layer, which is similar to the weight and bias generation in the feature layer. and It is also a randomly generated connection weight and bias uniformly distributed in the interval (0,1); F m =[SF m |OF m ] represents the total output of the feature window, SF m and OF m They represent the cumulative output of the overcomplete feature window and the slow feature window, and the total output E of the enhanced window module. n For E n =[E 1 ,E 2 ,...,E n ], E 1 , E 2 and E n Represent the 1st, 2nd and nth enhancement windows respectively;

[0110] Step 33: Integrate the output of the slow feature window, the output of the overcomplete window, and the output of the enhancement window into a transformation matrix A 1 Presented to the output layer of the SFA-OBLS network; at this point, the output Y∈R of the SFA-OBLS network NXC , N is the number of samples, C represents the output dimension, and R represents the set of real numbers, which can be expressed as:

[0111]

[0112] The actual output Y of the output layer of the SFA-OBLS network is the output result predicted based on the 13 sets of input data; [SF 1 ,SF 2 ,…,SF p ∣OF 1 ,OF 2 ,…,OF m |E 1 ,E 2 ,…,E n ]Use the expression A 1 express,

[0113] A 1 =[SF m |OF m |E n ], SF 1 , S.F. 2 and SF p Respectively represent the 1st, 2nd and pth slow feature windows, OF 1 ,OF 2 and OF m denote the 1st, 2nd and mth overcomplete windows respectively, E 1 , E 2 and E n Respectively represent the 1st, 2nd and nth enhancement windows, SF m and OF m They represent the cumulative output of the overcomplete feature window and the slow feature window, E n is the total output of the enhanced window module; the network connection weight of the output layer of the first layer is W 1 Indicates that and The weight parameters of the hidden layer and the output layer composed of the slow feature window module of the first layer, the overcomplete window module and the enhanced window module are respectively, and the output of the SFA-OBLS network is Y = (y 1 ,y 2 ,…,y N ) T ∈R N×C represents the label, N is the number of samples, C represents the output dimension, and R represents the real number set. and are obtained through ridge regression technology, ξ 1 and q Both are sigmoid activation functions, that is, and All are tansig feature layer activation functions;

[0114] Step 34: Use the generalized mapping function P 1(·) to generate m groups of slow feature nodes; using the generalized mapping function F 1 (·) to generate m groups of overcomplete feature nodes; using the generalized mapping function Q 1 (·) to generate n groups of enhanced nodes; therefore, Y is further expressed as:

[0115]

[0116] Where sf is the slow feature set extracted by the SFA-OBLS model through the slow features, and XL is the overcomplete feature set extracted by the SFA-OBLS model through the overcomplete features. YesP 1 (·) and F 1 (·) and Q 1 (·) composite function; weight parameter and They are a series of weights generated by the slow feature window, the over-complete window and the enhanced window; is the slow feature window module of the first layer, [A 1 ] + Indicates A 1 The pseudo-inverse matrix of , the weight parameters of the hidden layer and output layer composed of the overcomplete window module and the enhanced window module;

[0117] Step 35: W 1 Use the following optimization problem to solve:

[0118]

[0119] Among them A 1 =[SF m |OF m |E n ], SF m and OF m They represent the cumulative output of the overcomplete feature window and the slow feature window, E n To enhance the total output of the window module, the network connection weight of the output layer of the first layer is denoted by W 1 express, is the weight parameter of the hidden layer and output layer composed of the slow feature window module, the overcomplete window module and the enhanced window module in the first layer, λ 3 is used to constrain W 1 The adjustment parameters are 10 -4 to 10 -2 between but

[0120] μ(W 1 )=(W 1 ) T A T A1 W 1 -W T (A 1 ) T YY T A 1 W 1 +Y T Y+λ 3 (W 1 ) T W 1 ;

[0121] Step 36: Order get:

[0122]

[0123] Among them, λ 3 is used to constrain W 1 The adjustment parameters are 10 -4 to 10 -2 Between, I is the identity matrix, A 1 =[SF m |OF m |E n ], Indicates A 1 The transpose of

[0124] Step 37: The inverse problem degenerates into a least squares problem:

[0125]

[0126]

[0127] Among them A 1 + Indicates A 1 The pseudo-inverse matrix, λ 3 is used to constrain W 1 The adjustment parameters are 10 -4 to 10 -2 Between, I is the identity matrix, A 1 =[SF m |OF m |E n ], SF m and OF m They represent the cumulative output of the overcomplete feature window and the slow feature window, E n To enhance the overall output of the window module, Indicates A 1 The transpose of

[0128] Online application stage:

[0129] Online test dataset X obtained in real time new The network parameters obtained in the SFA-OBLS network training phase were used to calculate the SFA-OBLS network based on 13 input variables: the mass concentration of suspended solids in the water (g / m 3 ), mass concentration of biodegradable organic matter (g / m 3 ), mass concentration of particulate inert organic carbon (g / m 3 ), active heterotrophic bacteria mass concentration (g / m 3 ), active autotrophic bacteria mass concentration (g / m 3 ), particle mass concentration produced by biological attenuation (g / m 3 ), dissolved oxygen mass concentration (g / m 3 ), nitrate nitrogen mass concentration (g / m 3 ), ammonia nitrogen mass concentration (g / m 3 ), soluble biodegradable organic nitrogen mass concentration (g / m 3 ), alkalinity molar concentration (mol / m 3 ), total nitrogen concentration (g / m 3 ), mass concentration of chronically biodegradable organic matter (g / m 3 ), to predict the output variable five-day biochemical oxygen demand (BOD 5 )(g / m 3 ) and chemical oxygen demand (COD) (g / m 3 ), the input and output variable names of the test set are the same as those of the training set;

[0130] Step 38: Use the SFA-OBLS framework obtained in the network training phase to design and construct a sewage treatment water quality index detection model;

[0131] Step 39: Get new online sampling data X new ,Then the test data is preprocessed, and the test set data uses the same input variables as the training set as the input of the model;

[0132] Step 40: Use the SFA-OBLS model obtained in the network training phase to perform online detection to obtain a new detection result y new ;

[0133] Step 44: Regression prediction is used in the online application detection stage. BRIEF DESCRIPTION OF THE DRAWINGS

[0134] Figure 1 It is a structural diagram of the SFA-OBLS network proposed by the present invention;

[0135] Figure 2 Flow chart of modeling of sewage treatment process and water quality index detection based on SFA-OBLS network;

[0136] Figure 3 Schematic diagram of the sewage treatment process for a sewage treatment plant; DETAILED DESCRIPTION

[0137] The data of the present invention comes from the process wastewater treatment benchmark platform (BenchmarkSimulationModelNo.1, BSM1) developed by the International Water Quality Association and the European Union Science and Technology Cooperation. The specific flow chart of wastewater treatment is as follows: Figure 1 As shown. BSM1 has been widely recognized by the industry as an evaluation standard for sewage treatment process detection methods. The platform simulates a complete sewage treatment plant system, from the entry of sewage to the completion of purification. The treatment plant includes five activated sludge biochemical reaction tanks with a total volume of 5999 cubic meters (composed of two 1000 cubic meters of anoxic tanks and three 1333 cubic meters of aerobic tanks), and a secondary sedimentation tank with a total volume of 6000 cubic meters, divided into 10 layers, each layer 4 meters deep. In the sewage treatment process, two reaction modes of internal reflux and external reflux are adopted. The sewage plant treats 20,000 cubic meters of sewage per day, and the chemical oxygen demand for biodegradable is set to 300 mg per liter. The input and output variable parameters are shown in Table 1. By predicting 2 output variables based on the values ​​of 13 input variables, and comparing the predicted values ​​with the actual values, the detection accuracy of the water quality index detection model is evaluated. The present invention uses the actual influent water quality data provided by the BSM1 platform as the experimental data set. The data were sampled every 15 minutes, with 96 samples per day, and a total of 1,344 sample sets were collected over 14 days.

[0138] The SFA-OBLS model proposed in this invention is applied to the actual sewage treatment process and is specifically stated as follows:

[0139] SFA-OBLS network preprocessing and training phases in the sewage treatment process:

[0140] Step 1: First, preprocess the sampling data collected from the sewage treatment plant.

[0141] Step 2: Slow feature extraction is performed on 13 input process variables through the slow feature analysis algorithm. In the actual network training stage, there are 26 groups of slow features, including slowly changing components and rapidly changing components. The slowly changing components contain the key characteristic information of the sewage treatment data, while the rapidly changing components mainly include noise and internal and external disturbances. These disturbances and noise may cause false positives and false negatives in the modeling stage of the sewage treatment abnormal condition detection network. Based on the experience of sewage treatment process experts, the number of slow features is set to 12.

[0142] Step 3: Overcomplete feature extraction is performed on 13 input process variables through the overcomplete feature extraction algorithm. Affected by the complex biochemical reactions of the sewage treatment process, the extracted data inevitably has non-Gaussian characteristics, which will affect the detection performance of the sewage treatment process, resulting in a decrease in the accuracy of the prediction and control models and affecting the treatment efficiency. The non-Gaussian characteristics of sewage data can be well handled through overcomplete feature extraction, and the number of independent sources is set to 24.

[0143] Step 4: Input the slow feature information after slow feature processing and the overcomplete feature information after overcomplete feature processing into the 10 feature windows of the SFA-OBLS network respectively. After further feature extraction in each feature window, the SF p It is input into the enhancement layer for feature enhancement mapping. The number of layers of SFA-OBLS stacking is determined to be 1.

[0144] Step 5: After repeating the experiment 50 times, the network parameters required to establish the SFA-OBLS model are determined and applied to the SFA-OBLS model for the online sewage treatment process application stage.

[0145] B. Online application training phase of sewage treatment process based on SFA-OBLS network:

[0146] Step 6: The online data obtained from the sewage treatment process is first preprocessed.

[0147] Step 7: In the online detection stage, the water quality index detection results are evaluated by regression prediction. The evaluation indicators selected by the present invention are: MSE, RMSE, MAE and R 2 , in order to evaluate the model detection performance.

[0148] The training phase and online application training phase of the SFA-OBLS model are the specific practical applications of the SFA-OBLS model on the sewage treatment platform. In order to average the performance of the sewage treatment water quality index detection model, the mean square error (MSE), root mean square error (RMSE), mean absolute error (MAE) and determination coefficient (R 2 ) as the evaluation index. Table 2 shows the evaluation index and its average value of the SFA-OBLS model when changing the number of windows of feature nodes in each layer, the number of feature nodes and enhanced nodes in each window.

[0149]

[0150]

[0151]

[0152]

[0153] where y i is the true value, is the predicted value, is the true mean and n is the number of samples.

[0154] Table 1

[0155]

[0156]

[0157] With reference to the Pollutant Discharge Standard for Urban Wastewater Treatment Plants (GB 18918-2002), the International Water Association (IWA), the guidelines of the U.S. Environmental Protection Agency (EPA), the process design manual and the standards mentioned in some engineering case studies, the above variables are set with corresponding standards: Influent suspended solids mass concentration: 100-350g / m 3 ; Mass concentration of biodegradable organic matter: 30-80g / m 3 ; Mass concentration of granular inert organic carbon: 10-30g / m 3 ; Active heterotrophic bacteria mass concentration: 1,000-4,000g / m 3 ; Active autotrophic bacteria mass concentration: 100-500g / m 3 ; Mass concentration of particles produced by biological attenuation: 2-10g / m 3 ; Dissolved oxygen mass concentration: 2-4g / m 3 ; Nitrate nitrogen mass concentration: ≤10g / m 3 ; Ammonia nitrogen mass concentration: ≤5g / m 3 ; Soluble biodegradable organic nitrogen mass concentration: 1-5g / m 3 ; Alkalinity molar concentration: 1-4 mol / m 3 ;Total nitrogen concentration: ≤15g / m 3 ; Mass concentration of chronically biodegradable organic matter: 30-80g / m 3 ; Five-day biochemical oxygen demand: ≤20mg / L; Chemical oxygen demand: ≤50mg / L. Table 2

[0158]

[0159]

[0160] By observing Table 2, we can see that R 2 The values ​​are very close to 1, and the values ​​of MSE, RMSE and MAE are very small, indicating that the model fits the data very well. The experimental results verify that the SFA-OBLS model is suitable for BOD 5The prediction effect of COD is very good, which reflects the advantages of the model in terms of high precision and high efficiency.

Claims

1. A multi-index soft sensing method based on a slow feature overcomplete width learning system, characterized by It is divided into three parts: preprocessing stage, model training stage and online application stage; Preprocessing stage: Step 1: First, obtain the actual data set of the sewage plant and divide the data set into a training data set and a test data set; Then determine the input and output process variables of the detection model; Input variables are represented by array X; The input variables are: mass concentration of suspended solids in water (g / m 3 ), mass concentration of biodegradable organic matter (g / m 3 ), mass concentration of particulate inert organic carbon (g / m 3 ), active heterotrophic bacteria mass concentration (g / m 3 ), active autotrophic bacteria mass concentration (g / m 3 ), particle mass concentration produced by biological attenuation (g / m 3 ), dissolved oxygen mass concentration (g / m 3 ), nitrate nitrogen mass concentration (g / m 3 ), ammonia nitrogen mass concentration (g / m 3 ), soluble biodegradable organic nitrogen mass concentration (g / m 3 ), alkalinity molar concentration (mol / m 3 ), total nitrogen concentration (g / m 3 ), mass concentration of chronically biodegradable organic matter (g / m 3 ); the output variable is represented by array Y; The output variable is the five-day biochemical oxygen demand (BOD5) (g / m 3 ) and chemical oxygen demand (COD) (g / m 3 ); the input and output variable names of the test set are the same as those of the training set; the output indicators are predicted by inputting 13 measurable water quality data; Step 2: Preprocess the training data set and represent the training set with the symbol X, that is, X∈R NxM ; R is a real number set, N is the number of samples, and M is the number of process variables; M-dimensional input signal X = [x1(t), x2(t), ..., x M (t)], where t is the sample; Network training phase: Step 3: The goal of the slow feature extraction of the SFA-OBLS model is to find a mapping function Where m represents the number of slow features, by inputting the input signal X into The mapping function makes the characteristic component sf=g(X) change as little as possible over time. The slow characteristic information after mapping is: Where X is the 13 input variables of the training data set: mass concentration of suspended solids in water (g / m 3 ), mass concentration of biodegradable organic matter (g / m 3 ), mass concentration of particulate inert organic carbon (g / m 3 ), active heterotrophic bacteria mass concentration (g / m 3 ), active autotrophic bacteria mass concentration (g / m 3 ), particle mass concentration produced by biological attenuation (g / m 3 ), dissolved oxygen mass concentration (g / m 3 ), nitrate nitrogen mass concentration (g / m 3 ), ammonia nitrogen mass concentration (g / m 3 ), soluble biodegradable organic nitrogen mass concentration (g / m 3 ), alkalinity molar concentration (mol / m 3 ), total nitrogen concentration (g / m 3 ), mass concentration of chronically biodegradable organic matter (g / m 3 ), where sf is the slow feature set extracted by the SFA-OBLS model through slow features, and sf1, sf2, ..., sf are obtained from the slow to fast speed. m m slow features; g(X) is The set of mapping functions of ; [·] T The formula represents the transpose; g(x1(t)), g(x2(t)), and g(x m (t)) is X=[x1(t),x2(t),...,x M (t)] through the mapping function The first, second and mth slow feature information is obtained; Step 4: The slow feature window needs to find the first few features with the slowest change rate, so that the variance Δ(sf i ), thus the following optimization problem is obtained: Among them <·> t is the slow feature sample at N time points The calculated average value; Represents the first-order difference information of the i-th slow feature, that is, Δ(·) is the value for sf i A natural measure of the rate of change; represent The square of Step 5: Add the following constraints to the optimization problem: These three constraints are the zero mean constraint. <sf i > t =0, unit variance constraint and decorrelation constraints <sf i sf j > t =0, Among them, the zero mean constraint requires that the mean of the output signal is 0 in the time dimension, ensuring that the optimization problem has only one solution; the unit variance constraint requires that the mean of the square of the output signal is equal to 1 in the time series; the decorrelation constraint requires that different components are independent of each other, showing that they carry different information; among them <·> t is a slow feature sample at N time points The calculated average value; Step 6: The first-order difference of discrete data indicates the speed of change: For each slow feature sf i and sf i-1 are represented by linear combinations of all input variables, through The mapping function is used to implement the mapping, namely: sf=WX (6) Where W = [w1,w2,...,w m ] T is the load matrix that needs to be optimized through slow features, w i represents the load vector of the load matrix W, w i is the transpose of ; g(X) is The set of mapping functions, sf is the set of slow features extracted by the SFA-OBLS model through slow features; it is easy to prove that if the input variables are to satisfy the zero mean constraint, it is necessary to force the input variables to automatically remove the mean; Step 7: Substitute the formula in step 6 into the formulas in steps 4 and 5 to obtain the slow feature optimization target in the linear case: The formula It is represented by the symbol A, which is the covariance matrix of the first-order difference of the input X, that is, is the first-order difference of the input X; the formula It is represented by symbol B, where B is the covariance matrix of input X, that is, w i The load vector representing the load matrix W; Yes i The transpose of ; where <·> t is a slow feature sample at N time points The calculated average value; Represents the first-order difference information of the i-th slow feature, that is, Δ(·) is the value for sf i A natural measure of the rate of change; yes The square of , from the constraint condition in step 4, we can get the unit variance constraint Right now Step 8: Spheroidize X to eliminate the correlation between variables and obtain interpretable data; the spheroidization process involves singular value decomposition of the covariance matrix B: B=UΛU T (8) Where U is an orthogonal matrix composed of singular vectors; Λ is a rectangular diagonal singular value matrix, and the elements on the diagonal of Λ are the singular values ​​of B; T represents transpose; Step 9: The original input X is sphericalized by singular value decomposition, and the correlation is removed, that is: z=Λ -1 / 2 U T X (9) Where U is an orthogonal matrix composed of singular vectors; Λ is a rectangular diagonal singular value matrix, and the elements on the diagonal of Λ are the singular values ​​of B; T represents transposition, and z is the result of spherification of X. <zz T > t =I m , I m is the m-order identity matrix; through the singular value decomposition operation, the original optimization problem is transformed into finding a matrix P, according to the unit variance and decorrelation constraints in step 5, so that sf = Pz and <sfsf T > t =I m ,Will <zz T > t =I m Substituting sf=Pz, we obtain: PP T =I m (10) T represents transpose; P is an orthogonal matrix; we get Take the minimum solution, the covariance matrix of the first-order difference of z Perform singular value decomposition and solve, is the first-order difference of z; Where P is an orthogonal matrix composed of eigenvectors; Ω is a rectangular diagonal singular value matrix, and the diagonal elements are the change values ​​of each slow feature; T represents transpose; Step 10: The final parameter coefficient matrix W is shown as follows: W=PΛ -1 / 2 U T (12) Where U is an orthogonal matrix composed of singular vectors; Λ is a rectangular diagonal singular value matrix, and the elements on the diagonal of Λ are the singular values ​​of B; T represents transpose; P is an orthogonal matrix composed of eigenvectors; Step 11: Find the slow feature sf: sf=WX=PΛ -12 U T X=Pz (13) Step 12: Input the SF obtained from the slow feature analysis into the feature layer window of the SFA-OBLS network for further feature extraction, and form the i-th group of feature nodes SF through the feature mapping function i : in is the SFA-OBLS network feature mapping function, the activation function of the feature layer is tansig, and the mapping function is used to extract the centralized variable information of the sewage data and map it to the enhancement window; p is the number of feature nodes in the feature layer, sf is the slow feature set extracted by the SFA-OBLS model through slow features, The weight matrix of the feature node mapped to the i-th window, β fi The bias vector mapped to the feature node of the i-th window, and is randomly generated; to prevent weight Affected by randomness, a sparse autoencoder is used to fine-tune it, where the sparse matrix W s Obtained by the following formula: Where argmin means finding W that minimizes the objective function s , λ1 is the L1 regularization parameter, λ1=2-8, after sparsification, the feature node SF i =sfW s ; Step 13: Output SF to the feature layer p Repeat the previous step p times to obtain; SF p =(SF1,SF2,…,SF p ) (16) Among them, SF1, SF2 and SF p Represent the 1st, 2nd and pth slow feature windows respectively; Step 14: The purpose of the overcomplete window of the SFA-OBLS model is to obtain S with an overcomplete feature representation: S=XL (17) The key to obtaining S is to obtain the mixed component L; X is the input variable, L is the mixed component matrix; Step 15: Calculate the cumulative generating function of x: where φ x (·) is the cumulative generating function, σ is a real variable used to adjust and analyze the characteristics of the variable x, is the moment generating function MGF of the input variable x, which means The expected value of Used to describe the characteristics of x, σ is a continuous real number variable used to construct the function x, σ T represents the transpose of σ; Step 16: By calculating The coefficient of the Taylor series expansion of the cumulant generating function at is used to find the second-order cumulant of x. The Hessian estimate of the second-order cumulant of x is as follows: in Represents the zero vector, all components are 0, → represents the direction of the vector, used to distinguish between vectors and scalars, the special one is There is no direction; cov(x) is the covariance matrix, which represents the covariance between different components of the input variable x. represents the cumulant generating function φ x (σ) The second-order partial derivative matrix at ; Step 17: Expand σ to non-zero vectors when When , the generalized covariance matrix C is obtained by accumulating the Hessian matrix of the generating function x (σ); in, represents the cumulant generating function φ x (σ) The second-order partial derivative matrix at , σ is a real variable used to adjust and analyze the characteristics of the variable x, which obeys the Gaussian distribution; Represents the input variable x and its transpose x T Multiply Expected value; express Expected value; is used to describe the characteristics of x, σ is a real number variable used to adjust and analyze the characteristics of variable x, σ T represents the transpose of σ, β x (σ) represents the first-order derivative of the cumulant generating function of x under σ, that is cov(x) is the covariance matrix; Step 18: C x (σ) is further expressed as: C x (σ)=LC S (d)L T (21) Among them C x (σ) is the generalized covariance matrix, σ is a real variable used to adjust and analyze the characteristics of variable x, L is the mixture component matrix, L T is the transpose of L, C s (d) represents the generalized covariance matrix of the independent source components, represents the cumulant generating function φ S (d) The second-order partial derivative matrix at , S is an overcomplete feature representation, d = L T σ; β x (σ) represents the first-order derivative of the cumulant generating function of x under σ, that is is an overcomplete feature representation S and The expected value of the product, yes The expected value of Used to describe the characteristics of S, S T represents the transpose of S, where d = L T σ, β S (d) represents the first derivative of the cumulant generating function of S under d; Step 19: Derive the generalized covariance matrix C of the independent sources s (d): Among them, C s (d) represents the generalized covariance matrix of the independent source components, represents the cumulant generating function φ S (d) The second-order partial derivative matrix at , S is an overcomplete feature representation, S T represents the transpose of S, d = L T σ, L is the mixture component matrix, L T is the transpose of L, σ is a real variable used to adjust and analyze the characteristics of variable x, is an overcomplete feature representation S and its transpose and The expected value of the product, yes The expected value of Used to describe the characteristics of S, β S (d) represents the first-order derivative of the cumulant generating function of S under d, β S (d) T Represents β S (d) is the transpose of, cov(S) is the covariance matrix of S, is the first derivative of the cumulant generating function of x, Indicates S and The expected value of the product, express The expected value of T represents the transpose of d; Step 20: Due to the independence of S, the generalized covariance C of the independent source components s (d) is a diagonal matrix; therefore, C x (σ) is further expressed as: Where k is the number of overcomplete features S; ω i (σ)=[C s (L T σ)] ii represents the i-th independent source component S i The covariance of i represents the i-th mixed component, Indicates l i The transpose of is a single atom of the mixing matrix L; Step 21: The generalized covariance matrix belongs to the subspace W; the subspace spanned by any number of covariance matrices is either equal to W or a subset of W; select a sufficiently large number s of generalized covariance matrices to span the subspace W = [H1, H2, ..., H s ]; Set the parameter s to 5 times the parameter k; Based on the given non-zero vectors σ1,σ2,...,σ s , construct an equivalent matrix H j =C x (σ j ),j∈{1,2,...,s}; Step 22: Based on the subspace W, the overcomplete window is semi-positive regularized to complete the single atom estimation Then the mixed component l is obtained i , the semi-positive definite problem of single atom estimation can be specifically expressed as: in It's an atom Based on the optimal solution obtained by semidefinite programming, argmax is used to find the matrix B that maximizes <G,B> in the space Span, B∈W=Span{H1,H2,...,H s } indicates that the matrix B is in the linear space Span{H1,H2,...,H s }; the goal is to find a matrix B that maximizes the inner product <G,B>; Span is defined by the matrices H1,H2,...,H s The generated linear space, the matrix G determines The reduction direction of atoms; (In the process of independent component analysis, the objective function is used to measure the independence of signals; In order to maximize the independence, the gradient of the objective function is used to update the estimation matrix; The matrix G is determined by the gradient of the objective function; For a solution vector w in the optimization problem, assuming that we use negative entropy J as the objective function, the matrix G is expressed as the outer product of the gradient vector, or directly as the gradient matrix: The gradient here Represents the objective function J relative to the signal w T The rate of change of X, where X is the input signal matrix;) Step 23: Constrain the semi-definite problem: where B∈W=Span{H1,H2,...,H s } indicates that the matrix B is in the linear space Span{H1,H2,...,H s }; the goal is to find a matrix B that maximizes the inner product <G,B>, H1,H2,...,H s is the matrix constructed in step 21. Span is defined by the matrices H1, H2, ..., H s The generated linear space, Tr(·) represents the sum of the elements on the diagonal of the matrix B; when the optimal solution exists and is unique, It's an atom The optimal solution based on semidefinite programming; G in step 22 determines The reduction direction of atoms, that is, the selection of the best atoms, is determined by the matrix G; (In the process of independent component analysis, the objective function is used to measure the independence of signals; In order to maximize independence, the gradient of the objective function is used to update the estimation matrix; The matrix G is determined by the gradient of the objective function; For a solution vector w in the optimization problem, assuming that we use negative entropy J as the objective function, the matrix G is expressed as the outer product of the gradient vector, or directly as the gradient matrix: The gradient here Represents the objective function J relative to the signal w T The rate of change of X, where X is the input signal matrix;) Step 24: Matrix G is W (t) The orthogonal complement of When (t) Yes H (t) The first left singular vector of is u (t) The transpose of , then Step 25: Express the constraint set in the form of a set K = {B∈W, Tr(B) = 1, B≥0}: K={ <B,F j >=0,j∈[mk],Tr(B)=1,B≥0} (28) where B∈W=Span{H1,H2,...,H s } indicates that the matrix B is in the linear space Span{H1,H2,...,H s }; the goal is to find a matrix B that maximizes the inner product <G,B>, with the constraint that B∈W=Span{H1,H2,...,H s } is equivalent to B being orthogonal to any matrix in the complementary basis of the subspace W; N(W) = {F1, F2, ..., F m-k } represents the complementary basis of the orthogonal basis of subspace W, where m = M 2 / 2, Tr(·) represents the sum of the elements on the diagonal of the matrix B; Step 26: In actual industrial processes, since noise is inevitable, atomic estimation can often only be performed in an approximate form. To solve this problem, relax the original hard constraints based on the original constraints: K={Tr(B)=1,B≥0} (30) Among them B * is an atomic estimate, argmax is used to find the matrix B that maximizes <G,B> in the space Span, B∈W=Span{H1,H2,...,H s } indicates that the matrix B is in the linear space Span{H1,H2,...,H s }; the goal is to find a matrix B that maximizes the inner product <G,B>, and Span is defined by the matrices H1,H2,...,H s The generated linear space, matrix G determines B * The reduction direction of the atom; Use the gradient of the objective function to update the estimation matrix; The matrix G is determined by the gradient of the objective function; For a solution vector w in the optimization problem, assuming that we use negative entropy J as the objective function, the matrix G is expressed as the outer product of the gradient vector, or directly as the gradient matrix: The gradient here Represents the objective function J relative to the signal w T The rate of change of X, X is the input signal matrix; μ>0 is a regularization parameter; B∈W=Span{H1,H2,...,H s } means that B is orthogonal to any matrix in the complementary basis of the subspace W, N(W) = {F1, F2, ..., F m-k } represents the complementary basis of the orthogonal basis of the subspace W, K = {Tr(B) = 1, B ≥ 0} is the set form representation of the constraint set, and Tr(·) represents the sum of the elements on the diagonal of the matrix B; Step 27: For a solution vector w in the optimization problem, use negative entropy J as the objective function, and the matrix G is expressed as the outer product of the gradient vector, or directly as the gradient matrix: The gradient here Represents the objective function J relative to the signal w T The rate of change of X, X is the input signal matrix; the specified k-fold number of atoms is selected; then, the selected specified k-fold number of atoms are clustered into k clusters, denoted as D clust ; In each cluster, select a column as the representative atom. Here, the first column is selected by default. Finally, for D clust The k atoms in are screened, and the specific steps of screening are as follows: G clust =|D clust T D clust |-I k (31) N um =argmin{max(G clust )} (32) Where N um Indicates the minimum index among all maximum values, D clust T Yes D clust The transpose of G is to find the value such that max(G clust )The smallest N um , that is, D clust (N um ) represents the most independent atom, which can be used as an atom screened by the clustering algorithm; if D clust If the maximum values ​​corresponding to the remaining atoms in do not reach the preset threshold standard, they will be retained as new estimated atoms, and the threshold does not exceed the maximum value of the regularization parameter; Step 28: If the total number of atoms a screened by the clustering algorithm is less than the preset number k, all the excellent atoms found so far are excluded from the search range to update the constraint set; secondly, the subspace W is updated to the space spanned by the atoms remaining after the estimated atoms are removed, and the atoms found by the clustering algorithm are added to the complement of the orthogonal basis of the subspace W; finally, the estimation of the remaining atoms is completed using the fast iterative shrinkage threshold algorithm FISTA and the optimization-maximization principle MMP; the negative objective function can be further expressed as follows: The matrix G determines the reduction direction of atoms; The gradient here Represents the objective function J relative to the signal w T The rate of change of X, X is the input signal matrix; B∈W=Span{H1,H2,...,H s } indicates that the matrix B is in the linear space Span{H1,H2,...,H s }; the goal is to find a matrix B that maximizes the inner product <G,B>, μ>0 is a regularization parameter, N(W)={F1,F2,...,F m-k } indicates that the complementary basis B of the orthogonal basis of subspace W is an atomic estimate, m = M 2 / 2, t∈[a,k] represents the number of atoms currently obtained; Step: 29: Estimate the remaining ka atoms; combined with the given constraints K = {Tr(B) = 1, B ≥ 0}, the gradient of the differentiable part of the objective function is: The matrix G determines the reduction direction of atoms; The gradient here Represents the objective function J relative to the signal w T The rate of change of X, X is the input signal matrix;), B∈W=Span{H1,H2,...,H s } indicates that the matrix B is in the linear space Span{H1,H2,...,H s }; the goal is to find a matrix B that maximizes the inner product <G,B>, μ>0 is a regularization parameter, N(W)={F1,F2,...,F m-k } indicates that the complementary basis B of the orthogonal basis of subspace W is an atomic estimate, m = M 2 / 2, t∈[a,k] represents the number of atoms currently obtained; the optimization-maximization principle will terminate the contraction process in advance when a suitable atom is found; Step 30: Complete the atom estimation and obtain k atoms; then perform eigenvalue decomposition after each atom is obtained by matrixing, and select the eigenvector corresponding to the maximum eigenvalue as a mixing component l in the mixing component L i ; The dimension of the potential independent component S obtained based on the mixed component L exceeds the dimension of the number of input features, which makes the neurons of the SFA-OBLS network have more freedom when setting weights or biases; Step 31: Get the cumulative output OF of m overcomplete windows m : OF m =[OF1,OF2,…,OF m ] (35) Among them, OF1, OF2 and OF m Represent the 1st, 2nd and mth over-complete windows respectively; Step 32: SFA-OBLS maps the slow-varying dynamic feature components reflecting the overall trend of process changes and the potential independent components with overcomplete feature representation to the enhancement window module through the slow feature window and the overcomplete window to achieve nonlinear expansion and obtain enhanced features; the output of the j-th enhancement node of the enhancement window can be expressed as E j , E j The expression is as follows: Where n represents the total number of n enhancement nodes in the enhancement layer, ξ j As a nonlinear activation function, sigmoid is selected as the activation function of the enhancement layer, which is similar to the weight and bias generation in the feature layer. and It is also a randomly generated connection weight and bias uniformly distributed in the interval (0,1); F m =[SF m |OF m ] represents the total output of the feature window, SF m and OF m They represent the cumulative output of the overcomplete feature window and the slow feature window, and the total output E of the enhanced window module. n For E n =[E1,E2,...,E n ], E1, E2 and E n Represent the 1st, 2nd and nth enhancement windows respectively; Step 33: Integrate the output of the slow feature window, the output of the overcomplete window, and the output of the enhancement window into a transformation matrix A1 and present it to the output layer of the SFA-OBLS network; so far, the output Y∈R NXC , N is the number of samples, C represents the output dimension, and R represents the set of real numbers, which can be expressed as: The actual output Y of the output layer of the SFA-OBLS network is the output result predicted based on 13 sets of input data; SF1, SF2, …, SF p |OF1,OF2,…,OF m |E1,E2,…,E n ] is expressed by expression A1, A1=[SF m |OF m |E n ], SF1, SF2 and SF p Represent the 1st, 2nd and pth slow feature windows, OF1, OF2 and OF m denote the 1st, 2nd and mth overcomplete windows, E1, E2 and E n Respectively represent the 1st, 2nd and nth enhancement windows, SF m and OF m They represent the cumulative output of the overcomplete feature window and the slow feature window, E n is the total output of the enhanced window module; the network connection weight of the output layer of the first layer is W 1 Indicates that and are the weight parameters of the hidden layer and output layer composed of the slow feature window module of the first layer, the overcomplete window module and the enhanced window module, respectively. The SFA-OBLS network output Y = (y1, y2, ..., y N ) T ∈R N×C represents the label, N is the number of samples, C represents the output dimension, and R represents the real number set. and Both are obtained through ridge regression technology, ξ1 and ξ q Both are sigmoid activation functions, that is, and All are tansig feature layer activation functions; Step 34: Use the generalized mapping function P1(·) to generate m groups of slow feature nodes; use the generalized mapping function F1(·) to generate m groups of overcomplete feature nodes; use the generalized mapping function Q1(·) to generate n groups of enhanced nodes; therefore, Y is further expressed as: Where sf is the slow feature set extracted by the SFA-OBLS model through the slow features, and XL is the overcomplete feature set extracted by the SFA-OBLS model through the overcomplete features. is a composite function of P1(·), F1(·) and Q1(·); the weight parameter and They are a series of weights generated by the slow feature window, the over-complete window and the enhanced window; is the first layer slow feature window module, [A1] + represents the pseudo-inverse matrix of A1, the weight parameters of the hidden layer and output layer composed of the overcomplete window module and the enhanced window module; Step 35: W 1 Use the following optimization problem to solve: Where A1 = [SF m |OF m |E n ], SF m and OF m They represent the cumulative output of the overcomplete feature window and the slow feature window, E n To enhance the total output of the window module, the network connection weight of the output layer of the first layer is denoted by W 1 express, is the weight parameter of the hidden layer and output layer composed of the slow feature window module, the overcomplete window module and the enhanced window module in the first layer, and λ3 is used to constrain W 1 The adjustment parameters are 10 -4 to 10 -2 between but μ(W 1 )=(In 1 ) T AND T A1W 1 -IN T (A1) T YY T A1W 1 +Y T Y+λ3(W 1 ) T IN 1 ; Step 36: Order get: Among them, λ3 is used to constrain W 1 The adjustment parameters are 10 -4 to 10 -2 Between, I is the unit matrix, A1 = [SF m |OF m |E n ], represents the transpose of A1; Step 37: The inverse problem degenerates into a least squares problem: A1 + represents the pseudo-inverse matrix of A1, and λ3 is used to constrain W 1 The adjustment parameters are 10 -4 to 10 -2 Between, I is the unit matrix, A1 = [SF m |OF m |E n ], SF m and OF m They represent the cumulative output of the overcomplete feature window and the slow feature window, E n To enhance the overall output of the window module, represents the transpose of A1; Online application stage: Online test dataset X obtained in real time new The network parameters obtained in the SFA-OBLS network training phase were used to calculate the SFA-OBLS network based on 13 input variables: the mass concentration of suspended solids in the water (g / m 3 ), mass concentration of biodegradable organic matter (g / m 3 ), mass concentration of particulate inert organic carbon (g / m 3 ), active heterotrophic bacteria mass concentration (g / m 3 ), active autotrophic bacteria mass concentration (g / m 3 ), particle mass concentration produced by biological attenuation (g / m 3 ), dissolved oxygen mass concentration (g / m 3 ), nitrate nitrogen mass concentration (g / m 3 ), ammonia nitrogen mass concentration (g / m 3 ), soluble biodegradable organic nitrogen mass concentration (g / m 3 ), alkalinity molar concentration (mol / m 3 ), total nitrogen concentration (g / m 3 ), mass concentration of chronically biodegradable organic matter (g / m 3 ), to predict the output variable five-day biochemical oxygen demand (BOD5) (g / m 3 ) and chemical oxygen demand (COD) (g / m 3 ), the input and output variable names of the test set are the same as those of the training set; Step 38: Use the SFA-OBLS framework obtained in the network training phase to design and construct a sewage treatment water quality index detection model; Step 39: Get new online sampling data X new ,Then the test data is preprocessed, and the test set data uses the same input variables as the training set as the input of the model; Step 40: Use the SFA-OBLS model obtained in the network training phase to perform online detection to obtain a new detection result y new ; Step 44: Regression prediction is used in the online application detection stage.

Images