CN113657028B - Online aerosol optical thickness prediction method based on multi-source information - Google Patents

Abstract

The invention relates to an aerosol optical thickness online prediction method based on multi-source information, and belongs to the technical field of atmospheric channel optical signal transmission. The method comprises the steps of establishing a multi-source data acquisition system, carrying out feature extraction on atmospheric environment parameters, constructing a nonlinear mapping model based on a dynamic extreme learning machine, carrying out feature extraction through principal component analysis, and realizing fusion of features such as wind speed, temperature, humidity, sky background brightness and the like, wherein the output parameters of the model are aerosol optical thickness. The method has the advantages that the parameter estimation is carried out by using an artificial bee colony algorithm, so that the optimal selection of the weight and the bias is realized; according to the data fragments which arrive in real time, a nonlinear mapping model is updated, an online aerosol optical thickness prediction function is realized, a single hidden layer neural network based on a dynamic extreme learning machine is formed, the bottleneck of insufficient precision and difficult calculation caused by a physical model driving mode is broken through, and a solution is provided for constructing online estimation of optical parameters under various conditions.

Description

An online prediction method for aerosol optical depth based on multi-source information

Technical Field

The invention belongs to the technical field of atmospheric channel optical signal transmission, and in particular relates to an online prediction method for aerosol optical thickness based on multi-source information.

Background Art

Optical parameter estimation is a key element in optical information transmission. In atmospheric channels, due to the influence of turbulence, the light spot is often offset and the arrival angle fluctuates, which makes it difficult to estimate the optical parameters. Aerosol optical depth is an important support for judging optical parameters, and is of great significance for determining the values and changing laws of optical parameters. Model establishment and inversion calculations using physical methods will cause bottlenecks that are difficult to model or cannot be solved, affecting the accuracy of optical parameter estimation. How to predict aerosol optical depth (AOD) online in a data-driven way is of great significance for accurately estimating optical parameters and designing optical communication systems.

Chinese patent "AOD vertical correction effect evaluation method and system based on aerosol ground-based data" publication number 106446307A proposes a method for calculating aerosol data, including estimating the near-ground aerosol extinction coefficient and inverting the aerosol extinction coefficient based on visibility meter observation data, and using a simulation method to vertically correct the AOD. However, it does not consider the impact of different weather conditions on aerosol thickness, and is unable to predict the changing trend of AOD.

Chinese patent "Method and system for fusion of satellite AOD data based on inverse variance weighted average" Publication No. 106407634A provides a method and system for fusion of satellite AOD data based on inverse variance weighted average. This invention can improve the coverage rate and fusion accuracy requirements of satellite aerosol thickness data. Chinese patent "A method for correcting satellite AOD products by integrating multi-source characteristic geographic parameters" Publication No. 109213964A discloses an aerosol correction method, which corrects the parameters measured by satellite through random forests and correction thresholds to provide a reference for studying the effective concentration of atmospheric PM2.5. However, these two methods are difficult to meet the conditions of online changes and have a strong dependence on the selection of thresholds.

China's patent publication number 106096246A, "Aerosol optical thickness estimation method based on PM2.5 and PM10", provides an aerosol estimation method that estimates through a gravitational neural network model to solve the problem of low inversion accuracy and difficulty in real-time acquisition. However, this method does not consider the relationship between measurement data changes and model updates, and cannot meet the automatic adjustment of the model. China's patent publication number 110186823A, "A method for inverting aerosol optical thickness", establishes a bidirectional reflectivity distribution function parameter database through historical data, and obtains the corresponding aerosol optical thickness value through a table lookup method, but this method cannot calculate the generalization performance of the model, and the calculation results are only valid under specific assumptions.

To meet different needs, extreme learning machine modeling is widely used, such as the Chinese patent "A power fault prediction method based on ELM-CHMM" publication number 109615003A, and the Chinese patent "A wireless sensor network intrusion detection system and method integrating KPCA and ELM" publication number 109818798A. Due to the characteristics of fast convergence speed and clear analytical expression, the extreme learning machine effectively supports the software functional modules of complex systems.

Studying an online prediction method for aerosol optical depth based on multi-source information can not only effectively estimate optical parameters, but also promote the design of atmospheric optical communication systems, which has practical application value.

Summary of the invention

The present invention provides an online prediction method for aerosol optical thickness based on multi-source information, which fully considers the multi-source data information and the change trend of the data, and provides a solution for online estimation of optical parameters under various conditions.

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

(1) Establish a multi-source information acquisition system to collect and store atmospheric environment parameters and optical parameters;

(2) Extract the characteristics of atmospheric environmental parameters to obtain wind speed characteristics WS _D , temperature characteristics TE _D , humidity characteristics HU _D , and sky background brightness characteristics LU _D , which constitute the input parameter x, and the aerosol optical thickness is the output parameter y;

(3) Using data-driven thinking, a nonlinear mapping model based on a dynamic extreme learning machine (ELM) was established, the relationship between input parameters and output parameters was determined, and the artificial bee colony algorithm was used for optimization to determine the optimal solution for the input weight and bias of the nonlinear mapping model;

(4) The nonlinear mapping model is updated online, the error between the real value and the predicted value of the real-time data segment is calculated, the nonlinear mapping model is quantitatively updated according to the error, the output weight of the nonlinear mapping model is analytically calculated, and the predicted value of the aerosol optical thickness of the next data segment is obtained to realize the online prediction of the aerosol optical thickness. Through feature extraction, nonlinear model and model updating strategy, the multi-source features are integrated to realize the online prediction of the aerosol optical thickness.

The multi-source information acquisition system in step (1) of the present invention includes a data acquisition device, an information exchange terminal, and a database. The data acquisition device includes a laser radar, a ground meteorological station, a background radiometer, and a corresponding data interface. The database uses SQL Server to store data. The information exchange terminal exchanges information with the data acquisition device and saves the collected data information in the database. The ground meteorological station, the background radiometer, and the laser radar are used to realize the collection of atmospheric environment parameters and optical parameters. The atmospheric environment parameters include the atmospheric data wind speed WS, temperature TE, humidity HU provided by the ground meteorological station, and the sky background brightness LU collected by the background radiometer. The optical parameter is the aerosol optical thickness, which is collected by the laser radar.

The atmospheric environment parameter feature extraction method in step (2) of the present invention is as follows:

The atmospheric data stored in the data acquisition system are subjected to principal component analysis, including wind speed WS, temperature TE, humidity HU, and sky background brightness LU. Considering the cumulative effect of M days, the wind speed sample fragment WS _K on the Kth day is written in the following matrix form:

Among them, ws _KM represents the wind speed value of the Mth day in the Kth sample segment, and each wind speed value in the wind speed sample matrix segment WS _K is standardized:

Wherein, ws′ _km represents the normalized wind speed value of the mth day in the kth sample segment, μ _m and σ _m represent the mean and variance of the wind speed sample segment of the mth day, respectively. After processing, the wind speed normalization matrix WS′ _K is obtained, and the wind speed eigenvalues λ _WS,1 ,λ _WS,2 ,…λ _WS,M corresponding to the wind speed eigenvalues η _WS,1 ,η _WS,2 ,…η _WS,M are calculated. The wind speed eigenvalues are sorted from large to small, and the contribution of the wind speed eigenvalues is calculated. When the contribution of the selected wind speed eigenvalue is greater than the threshold θ ₁ , the number of principal components Num ₁ is determined;

The corresponding wind speed feature vector set is the wind speed feature

Considering the cumulative effect of M days, the temperature sample fragment TE _K of the Kth day is standardized to obtain the temperature standardized matrix TE′ _K , and the temperature eigenvalues λ _TE,1 ,λ _TE,2 ,…λ _TE,M corresponding to the temperature eigenvalues η _TE,1 ,η _TE,2 ,…η _TE,M are calculated. The temperature eigenvalues are sorted from large to small, and the contribution of the temperature eigenvalues is calculated. When the contribution of the selected temperature eigenvalue is greater than the threshold θ ₂ , the number of principal components Num ₂ is determined;

The corresponding temperature feature vector set is the temperature feature

Considering the cumulative effect of M days, the humidity sample segment HU _K of the Kth day is standardized to obtain the humidity standardized matrix HU′ _K , and the humidity eigenvalues λ _HU,1 ,λ _HU,2 ,…λ _HU,M corresponding to the humidity eigenvalues η _HU,1 ,η _HU,2 ,…η _HU,M are calculated. The humidity eigenvalues are sorted from large to small, and the contribution of the humidity eigenvalues is calculated. When the contribution of the selected humidity eigenvalue is greater than the threshold θ ₃ , the number of principal components Num ₃ is determined.

The corresponding humidity feature vector set is the humidity feature

Considering the cumulative effect of M days, the sky background brightness sample fragment LU _K of the Kth day is standardized to obtain the sky background brightness standardized matrix LU′ _K , the sky background brightness eigenvalues λ _LU,1 ,λ _LU,2 ,…λ _LU,M corresponding to the sky background brightness eigenvalues η _LU,1 ,η _LU,2 ,…η _LU,M are calculated, the sky background brightness eigenvalues are sorted from large to small, the contribution of the sky background brightness eigenvalues is calculated, and when the contribution of the selected sky background brightness eigenvalue is greater than the threshold θ ₄ , the number of principal components Num ₄ is determined;

The corresponding sky background brightness feature vector set is the sky background brightness feature

The input parameter x is a set of wind speed feature TE _D , temperature feature TE _D , humidity feature HU _D , and sky background brightness feature LU _D , that is, x=(WS _D , TE _D , HU _D , LU _D ).

The construction of the nonlinear mapping model in step (3) of the present invention includes a three-layer neural network structure of an input layer, a hidden layer, and an output layer. Considering the newly arrived sample data, the data sampling length is N. In the initial stage, the input layer has N nodes, the hidden layer has L nodes, and the output layer has 1 node. The input parameter x and the output parameter y are sampled and merged to obtain the sample data. _Xi represents the i-th input parameter of x, and _Yi represents the i-th output parameter of y. The input-output relationship of the nonlinear mapping model is as follows:

Among them, a and b are input weight and bias respectively, β is the output layer weight, L is the number of hidden layer nodes of the neural network, is the activation function. Considering the mapping relationship between nodes as a whole, according to the extreme learning machine theory, the weights and biases between the input layer nodes and the hidden layer are randomly selected. The output weight β is the only variable that needs to be solved in the entire network. Through Moore-Penrose generalized inverse calculation, formula (7) can be expressed in the following matrix form:

y＝Hβ······························(8)

Among them, y = [y ₁ ,y ₂ ,...,y _N ] ^T represents the aerosol optical thickness of the current segment. The hidden matrix H is introduced into the nonlinear model for calculation. H is a random matrix with N rows and L columns. The specific expression is:

in, To activate the function, nonlinear mapping is realized, which ensures that the function has derivatives in the interval. In order to further calculate the output weight β, the Moore-Penrose generalized inverse of H is calculated.

H ⁺ =(H ^T H) ^-1 H ^T ··················(10)

By solving the Moore-Penrose generalized inverse, we can obtain the optimal solution β ^* of the output weight β:

β ^* =H ⁺ y·····························(11)

Further optimize the input weights and biases, encode the weights and biases into strings, and obtain the candidate solution vector set E = (a, b), which is in the following form:

The boundaries of the candidate vector set are the maximum and minimum solution vectors, which are expressed here as E _max = {max(e)|e∈E} and E _min = {min(e)|e∈E}, where e represents a set of elements that satisfy the boundary conditions. According to the artificial bee colony algorithm, the solution vector of the input weight and bias under the initial conditions is expressed as:

E _h,j =E _min,j +ω×(E _max,j -E _min,j )···················(13)

E _min,j and E _max,j represent the maximum and minimum values of the jth input solution vector respectively. ω is evenly distributed in the interval [-1,1]. The optimal solution vector search is performed and the following update strategy is obtained:

Among them, E _h,j represents the current solution vector, Represents a new feasible solution vector, h, j and r are all indicators, u _{h, j} changes randomly in the interval [-1,1], calculates the fitness, and selects the optimal solution vector:

Among them, num represents the number of iterations. After the iteration is completed, the optimal solution vector will be updated to obtain the optimal input weight a ^* and bias b ^* . Therefore, the optimal solution of the expression of the nonlinear model is:

The nonlinear mapping model updating method in step (4) of the present invention is:

In the initial stage, the output weights are obtained New sample snippet After arrival, the error between the true value and the predicted value of the real-time arriving data segment is calculated, expressed as:

Among them, N ₀ and N ₁ represent the starting position and end position of the current sample segment, H ₁ and y ₁ represent the hidden layer output and label set of the newly arrived segment, β ₁ represents the current nonlinear model output weight. Based on the current output weight, the expression of β ₁ is:

Divide the above expression into two parts and introduce intermediate variables According to the associative law of matrix multiplication, the expression of the first part of the matrix is:

In the current fragment state, H ₁ and K ₀ are both known conditions. Substituting them into the overall output layer weight expression, we can get the expression of the output weight β ₁ in the model update phase:

Further iterative calculations yield the intermediate variable K _k+1 expression corresponding to the (k+1)th segment:

Further iterative calculations are performed to obtain the optimal solution for the output weight β _k+1 of the nonlinear mapping model corresponding to the (k+1)th segment:

The model expression after online update is:

By substituting the input parameter x into the updated model, the aerosol optical depth y can be calculated, thus realizing the online prediction of aerosol optical depth.

The present invention establishes a multi-source data acquisition system, and realizes information interaction from the data acquisition device to the database through the information interaction terminal. The data acquisition device realizes the collection of data such as wind speed, temperature, humidity, sky background brightness, etc., and transmits them to the information interaction terminal through the data interface, and realizes data storage through data processing operations. Based on the multi-source data acquisition system, a nonlinear mapping model is constructed, a neural network structure is designed, and the optimization of the neural network weights and biases is realized to obtain a mapping model of aerosol optical thickness. The nonlinear model is updated through the real-time collected data to realize an online prediction method for aerosol optical thickness.

The advantages of the present invention are:

(1) The present invention takes into account a multi-source data acquisition system, stores environmental parameters and optical parameters, establishes a nonlinear mapping model through the idea of data-driven, realizes the online prediction of aerosol optical parameters, and breaks through the bottleneck of insufficient accuracy and difficulty in calculation caused by the physical model-driven method.

(2) Feature extraction is performed through principal component analysis, and a neural network is constructed using a dynamic extreme learning machine. The weights and biases are optimized to obtain an input-output relationship that can adapt to the online environment.

(3) Relying on a multi-source information processing solution, the model is quantitatively adjusted according to the current data fragment, avoiding the need to relearn the model and automating the entire prediction process, providing a solution for the intelligent processing of complex systems and having a certain degree of universality in engineering practice.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG1 is a multi-source data acquisition system of the present invention;

FIG2 is a nonlinear mapping model of the present invention;

FIG. 3 is a flow chart of the present invention.

DETAILED DESCRIPTION

The following steps are involved:

(1) Establish a multi-source information acquisition system to collect and store atmospheric environment parameters and optical parameters;

(2) Extract the characteristics of atmospheric environment parameters to obtain wind speed characteristics WS _D , temperature characteristics TE _D , humidity characteristics HU _D , and sky background brightness characteristics LU _D , which constitute the input parameter x, and the aerosol optical thickness is the output parameter y;

(3) Using data-driven thinking, a nonlinear mapping model based on a dynamic extreme learning machine (ELM) was established, the relationship between input parameters and output parameters was determined, and the artificial bee colony algorithm was used for optimization to determine the optimal solution for the input weight and bias of the nonlinear mapping model;

(4) The nonlinear mapping model is updated online, the error between the real value and the predicted value of the real-time data segment is calculated, the nonlinear mapping model is quantitatively updated according to the error, the output weight of the nonlinear mapping model is analytically calculated, and the predicted value of the aerosol optical thickness of the next data segment is obtained to realize the online prediction of aerosol optical thickness. Through feature extraction, nonlinear model and model updating strategy, multi-source features are integrated to realize the online prediction of aerosol optical thickness.

The multi-source information acquisition system in step (1) of the present invention includes a data acquisition device, an information exchange terminal, and a database. The data acquisition device includes a laser radar, a ground meteorological station, a background radiometer, and a corresponding data interface. The database uses SQL Server to store data. The information exchange terminal exchanges information with the data acquisition device and saves the collected data information into the database. The ground meteorological station, the background radiometer, and the laser radar are used to collect atmospheric environment parameters and optical parameters. The atmospheric environment parameters include atmospheric data wind speed WS, temperature TE, humidity HU provided by the ground meteorological station, and sky background brightness LU collected by the background radiometer. The optical parameter is aerosol optical thickness, which is collected by the laser radar.

The atmospheric environment parameter feature extraction method in step (2) of the present invention is as follows:

The atmospheric data stored in the data acquisition system are subjected to principal component analysis, including wind speed WS, temperature TE, humidity HU, and sky background brightness LU. Considering the cumulative effect of M days, the wind speed sample fragment WS _K on the Kth day is written in the following matrix form:

Among them, ws _KM represents the wind speed value of the Mth day in the Kth sample segment, and each wind speed value in the wind speed sample matrix segment WS _K is standardized:

Wherein, ws′ _km represents the normalized wind speed value of the mth day in the kth sample segment, μ _m and σ _m represent the mean and variance of the wind speed sample segment of the mth day, respectively. After processing, the wind speed normalization matrix WS′ _K is obtained, and the wind speed eigenvalues λ _WS,1 ,λ _WS，2 ,…λ _WS,M corresponding to the wind speed eigenvalues η _WS,1 ,η _WS,2 ,…η _WS,M are calculated. The wind speed eigenvalues are sorted from large to small, and the contribution of the wind speed eigenvalues is calculated. When the contribution of the selected wind speed eigenvalue is greater than the threshold θ ₁ , the number of principal components Num ₁ is determined;

The corresponding wind speed feature vector set is the wind speed feature

Considering the cumulative effect of M days, the temperature sample fragment TE _K of the Kth day is standardized to obtain the temperature standardized matrix TE′ _K , and the temperature eigenvalues λ _TE,1 ,λ _TE,2 ,…λ _TE,M corresponding to the temperature eigenvalues η _TE,1 ,η _TE,2 ,…η _TE,M are calculated. The temperature eigenvalues are sorted from large to small, and the contribution of the temperature eigenvalues is calculated. When the contribution of the selected temperature eigenvalue is greater than the threshold θ ₂ , the number of principal components Num ₂ is determined;

The corresponding temperature feature vector set is the temperature feature

Considering the cumulative effect of M days, the humidity sample segment HU _K of the Kth day is standardized to obtain the humidity standardized matrix HU′ _K , and the humidity eigenvalues λ _HU,1 ,λ _HU,2 ,…λ _HU,M corresponding to the humidity eigenvalues η _HU,1 ,η _HU,2 ,…η _HU,M are calculated. The humidity eigenvalues are sorted from large to small, and the contribution of the humidity eigenvalues is calculated. When the contribution of the selected humidity eigenvalue is greater than the threshold θ ₃ , the number of principal components Num ₃ is determined.

The corresponding humidity feature vector set is the humidity feature

Considering the cumulative effect of M days, the sky background brightness sample fragment LU _K of the Kth day is standardized to obtain the sky background brightness standardized matrix LU′ _K , the sky background brightness eigenvalues λ _LU,1 ,λ _LU,2 ,…λ _LU,M corresponding to the sky background brightness eigenvalues η _LU,1 ,η _LU,2 ,…η _LU,M are calculated, the sky background brightness eigenvalues are sorted from large to small, the contribution of the sky background brightness eigenvalues is calculated, and when the contribution of the selected sky background brightness eigenvalue is greater than the threshold θ ₄ , the number of principal components Num ₄ is determined;

The corresponding sky background brightness feature vector set is the sky background brightness feature

The input parameter x is the set of wind speed feature TE _D , temperature feature TE _D , humidity feature HU _D , and sky background brightness feature LU _D , that is, x=(WS _D , TE _D , HU _D , LU _D );

The construction of the nonlinear mapping model in step (3) of the present invention includes a three-layer neural network structure of an input layer, a hidden layer, and an output layer. Considering the newly arrived sample data, the data sampling length is N. In the initial stage, the input layer has N nodes, the hidden layer has L nodes, and the output layer has 1 node. The input parameter x and the output parameter y are sampled and merged to obtain the sample data. _Xi represents the i-th input parameter of x, and _Yi represents the i-th output parameter of y. The input-output relationship of the nonlinear mapping model is as follows:

Among them, a and b are input weight and bias respectively, β is the output layer weight, L is the number of hidden layer nodes of the neural network, is the activation function. Considering the mapping relationship between nodes as a whole, according to the extreme learning machine theory, the weights and biases between the input layer nodes and the hidden layer are randomly selected. The output weight β is the only variable that needs to be solved in the entire network. Through Moore-Penrose generalized inverse calculation, formula (7) can be expressed in the following matrix form:

y＝Hβ·····························(8)

Among them, y = [y ₁ ,y ₂ ,...,y _N ] ^T represents the aerosol optical thickness of the current segment. The hidden matrix H is introduced into the nonlinear model for calculation. H is a random matrix with N rows and L columns. The specific expression is:

in, To activate the function, nonlinear mapping is realized, which ensures that the function has derivatives in the interval. In order to further calculate the output weight β, the Moore-Penrose generalized inverse of H is calculated.

H ⁺ =(H ^T H) ^-1 H ^T ··················(10)

By solving the Moore-Penrose generalized inverse, we can obtain the optimal solution β ^* of the output weight β:

β ^* =H ⁺ y·····························(11)

Further optimize the input weights and biases, encode the weights and biases into strings, and obtain the candidate solution vector set E = (a, b), which is in the following form:

The boundaries of the candidate vector set are the maximum and minimum solution vectors, which are expressed here as E _max = {max(e)|e∈E} and E _min = {min(e)|e∈E}, where e represents a set of elements that satisfy the boundary conditions. According to the artificial bee colony algorithm, the solution vector of the input weight and bias under the initial conditions is expressed as:

E _h,j =E _min,j +ω×(E _max,j -E _min,j )···················(13)

E _min,j and E _max,j represent the maximum and minimum values of the jth input solution vector respectively. ω is evenly distributed in the interval [-1,1]. The optimal solution vector search is performed and the following update strategy is obtained:

Among them, E _h,j represents the current solution vector, represents a new feasible solution vector, h, j and r are all indicators. _{u h,j} changes randomly in the interval [-1,1]. Calculate the fitness and select the optimal solution vector:

Among them, num represents the number of iterations. After the iteration is completed, the optimal solution vector will be updated to obtain the optimal input weight a ^* and bias b ^* . Therefore, the optimal solution of the expression of the nonlinear model is:

The nonlinear mapping model updating method in step (4) of the present invention is:

In the initial stage, the output weights are obtained New sample snippet After arrival, the error between the true value and the predicted value of the real-time arriving data segment is calculated, expressed as:

Among them, N ₀ and N ₁ represent the starting position and end position of the current sample segment, H ₁ and y ₁ represent the hidden layer output and label set of the newly arrived segment, β ₁ represents the current nonlinear model output weight. Based on the current output weight, the expression of β ₁ is:

Divide the above expression into two parts and introduce intermediate variables According to the associative law of matrix multiplication, the expression of the first part of the matrix is:

In the current fragment state, H ₁ and K ₀ are both known conditions. Substituting them into the overall output layer weight expression, we can get the expression of the output weight β ₁ in the model update phase:

Further iterative calculations yield the intermediate variable K _k+1 expression corresponding to the (k+1)th segment:

Further iterative calculations are performed to obtain the optimal solution for the output weight β _k+1 of the nonlinear mapping model corresponding to the (k+1)th segment:

The model expression after online update is:

Substituting the input parameter x into the updated model, the aerosol optical thickness y can be calculated, realizing the online prediction of aerosol optical thickness. During the entire solution process, the nonlinear model does not need to be relearned, and the model weights can be updated according to the received data to realize the online prediction function of aerosol optical thickness.

The present invention will be described in detail below with reference to the accompanying drawings.

(1) Build a multi-source data acquisition system. The specific form is shown in Figure 1. The laser radar collects aerosol optical thickness information, the ground meteorological station collects wind speed, temperature, and humidity information, and the background radiometer collects sky background brightness information. The information exchange terminal is an ordinary computer that can realize information exchange with the data acquisition equipment. SQLServer software is loaded in the database, and the data storage function is realized by maintaining the data table;

The laser radar is connected to the information interaction terminal through the RJ45 network interface. The background radiometer and the ground meteorological station use the RS232 interface to communicate with the information interaction terminal. During the information collection process, the information interaction terminal sends a collection command to the data collection device, and the data collection device transmits the data information to the information interaction terminal through the data interface to realize information interaction.

(2) Feature extraction. The collected data is sorted and the cumulative effects of wind speed, temperature, humidity, and sky background brightness are considered through principal component analysis. The data collection period is one day, and each half-hour interval is one sampling point. There are 48 sample vectors per day. Considering the data of the 7 days before the current time node, a sample segment matrix with 48 rows and 7 columns is formed;

The principal component analysis of the atmospheric data stored in the data acquisition system, including wind speed WS, temperature TE, humidity HU, and sky background brightness LU, is performed. Considering the cumulative effect of M days, the wind speed sample fragment WS _K on the Kth day can be written in the following matrix form:

Among them, ws _KM represents the wind speed value of the Mth day in the Kth sample segment, and each wind speed value in the wind speed sample matrix segment WS _K is standardized:

Wherein, ws′ _km represents the normalized wind speed value of the mth day in the kth sample segment, μ _m and σ _m represent the mean and variance of the wind speed sample segment of the mth day, respectively. After processing, the wind speed normalized matrix WS′ _K is obtained, and the wind speed characteristic vectors η _WS,1 ,η _WS,2 ,…η _WS,M corresponding to the wind speed characteristic values λ _WS,1 ,λ _WS,2 ,…λ _WS,M are calculated. The wind speed characteristic values are sorted from large to small, and the contribution of the wind speed characteristic values is calculated. When the contribution of the selected wind speed characteristic value is greater than the threshold value θ ₁ =0.90, the number of principal components Num ₁ is determined;

The corresponding wind speed feature vector set is the wind speed feature

Considering the cumulative effect of M days, the temperature sample fragment TE _K of the Kth day is standardized to obtain the temperature standardized matrix TE′ _K , and the temperature eigenvalues λ _TE,1 ,λ _TE,2 ,…λ _TE,M corresponding to the temperature eigenvalues η _TE,1 ,η _TE,2 ,…η _TE,M are calculated. The temperature eigenvalues are sorted from large to small, and the contribution of the temperature eigenvalues is calculated. When the contribution of the selected temperature eigenvalue is greater than the threshold value θ ₂ =0.85, the number of principal components Num ₂ is determined.

The corresponding temperature feature vector set is the temperature feature

Considering the cumulative effect of M days, the humidity sample segment HU _K of the Kth day is standardized to obtain the humidity standardized matrix HU′ _K , and the humidity eigenvalues λ _HU,1 ,λ _HU,2 ,…λ _HU,M corresponding to the humidity eigenvalues η _HU,1 ,η _HU,2 ,…η _HU,M are calculated. The humidity eigenvalues are sorted from large to small, and the contribution of the humidity eigenvalues is calculated. When the contribution of the selected humidity eigenvalue is greater than the threshold value θ ₃ =0.85, the number of principal components Num ₃ is determined.

The corresponding humidity feature vector set is the humidity feature

Considering the cumulative effect of M days, the sky background brightness sample fragment LU _K of the Kth day is standardized to obtain the sky background brightness standardized matrix LU′ _K , and the sky background brightness eigenvalues λ _LU,1 ,λ _LU,2 ,…λ _LU,M corresponding to the sky background brightness eigenvalues η _LU,1 ,η _LU,2 ,…η _LU,M are calculated. The sky background brightness eigenvalues are sorted from large to small, and the contribution of the sky background brightness eigenvalues is calculated. When the contribution of the selected sky background brightness eigenvalue is greater than the threshold value θ ₄ =0.80, the number of principal components Num ₄ is determined;

The corresponding sky background brightness feature vector set is the sky background brightness feature

The input parameter x is the set of wind speed feature TE _D , temperature feature TE _D , humidity feature HU _D , and sky background brightness feature LU _D , that is, x=(WS _D , TE _D , HU _D , LU _D );

(3) Using data-driven thinking, a nonlinear mapping model based on a dynamic extreme learning machine (ELM) was established, the relationship between input parameters and output parameters was determined, and the artificial bee colony algorithm was used for optimization to determine the optimal solution for the input weight and bias of the nonlinear mapping model;

The construction of the nonlinear mapping model includes a three-layer neural network structure of input layer, hidden layer, and output layer. Considering the newly arrived sample data, the data sampling length is 48. In the initial stage, the input layer has 48 nodes, the hidden layer has 50 nodes, and the output layer has 1 node. Among them, the input layer nodes and the hidden layer are fully connected, and the hidden layer and the output layer are fully connected. The input layer weight is a matrix of 48 rows and 50 columns, and the bias is a matrix of 1 row and 50 columns. The input parameters x and output parameters y are sampled and merged to obtain sample data. _Xi represents the i-th input parameter of x, and _Yi represents the i-th output parameter of y. The input-output relationship of the nonlinear mapping model is as follows:

Among them, a and b are input weight and bias respectively, β is the output layer weight, L is the number of hidden layer nodes of the neural network, is the activation function. Considering the mapping relationship between nodes as a whole, according to the extreme learning machine theory, the weights and biases between the input layer nodes and the hidden layer are randomly selected. The output weight β is the only variable that needs to be solved in the entire network. Through Moore-Penrose generalized inverse calculation, formula (7) can be expressed in the following matrix form:

y＝Hβ·····················(8)

Among them, y = [y ₁ ,y ₂ ,...,y _N ] ^T represents the aerosol optical thickness of the current segment. The hidden matrix H is introduced into the nonlinear model for calculation. H is a random matrix with 48 rows and 50 columns. The specific expression is:

Among them, the input parameters of the nonlinear model can be expressed as x = [x ₁ , x ₂ , ..., x _N ] ^T , which integrates multiple features;

Calculate the output weights. The weights and biases between the input layer nodes and the hidden layer are randomly selected, and the weights between the hidden layer and the output layer are calculated using the MP generalized inverse. The output weights are obtained using the MP generalized inverse, and the generalized inverse of the hidden layer matrix H is

H ⁺ =(H ^T H) ^-1 H ^T ··········(10)

To meet the minimum error condition, β is the only variable that needs to be solved in the entire network, and y = [y ₁ ,y ₂ ,...,y _N ] ^T represents the aerosol optical thickness value of the current segment. Through the MP generalized inverse, the expression of the output weight is obtained;

β ^* =H ⁺ y·····························(11)

The weights and biases are further optimized, and the weights and biases are string-encoded to obtain a set of candidate vectors, which are in the following form:

The boundaries of the candidate vector set are the solution vectors of the maximum and minimum values, e represents a set of sets whose elements satisfy the boundary conditions, E _max = {max(e)|e∈E} and E _min = {min(e)|e∈E}. The optimization is performed according to the artificial bee colony algorithm. The artificial bee colony algorithm is analogous to bees searching for food to achieve multi-parameter optimization. Three forms are mainly considered, including employed bees, follower bees and scout bees. The three types of bees complete the information search task. In the process of solving the optimal weights and biases, the employed bees first search for information and feedback the results of the search. The follower bees receive the feedback information and continue to search around the optimal vector. If a better solution is found, it is updated; otherwise, the current solution is abandoned and the scout bees continue to search. Finally, the employed bees realize the local search near the optimal solution, the scout bees complete the global search, the entire algorithm converges, and the optimal solution vector of the weights and bias of the nonlinear mapping model is obtained through optimization;

Based on the maximum and minimum values E _min,j and E _max,j of the jth input solution vector, the boundary combination under the initial conditions is calculated, and the weighted solution vector is expressed as:

E _h,j =E _min,j +ω×(E _max,j -E _min,j )···················(13)

Determine the update parameters, and get u _h,j that changes randomly in the interval [-1,1]. Then, through the optimal vector search, get the candidate vector update strategy:

Through the iteration of random variables u _h,j , the fitness calculation under different weights and bias conditions is realized:

Here, the number of iterations num is 50. After the iteration is completed, the optimal solution vector will be updated to obtain the optimal input weight a ^* and bias b ^* . By integrating the weights and biases of each layer of the neural network, the relationship model between the input parameters and the output parameters is obtained:

In the initial stage, the neural network input layer has 48 nodes, the hidden layer has 50 nodes, and the output layer has 1 node. After determining all the parameters of the nonlinear mapping model, the relationship between the input parameters and the output parameters is obtained, that is, the analytical relationship between wind speed, temperature, humidity, sky background brightness and aerosol optical thickness. The output weight obtained in the initial stage is expressed as:

(4) Update the model and add new sample segments After arrival, the approximate error of the sample is calculated based on the collected data fragments, and the objective function is expressed as:

N ₀ and N ₁ represent the starting position and end position of the current sample segment, respectively. H ₁ and Y ₁ represent the hidden layer output and label set of the newly arrived segment, respectively. β ₁ represents the current nonlinear model weight. Based on the current output weight, the output layer weight expression is:

Consider the front and back parts of the expression separately and introduce intermediate variables By using the associative law of matrix multiplication, the expression for calculating the first part of the matrix is:

Consider the latter part of the expression of β ₁ and expand it using the principle of matrix multiplication to obtain:

On the basis of the current output weight, the new input parameters are substituted into the nonlinear model, and the expression of the output weight is obtained through matrix transformation. At the current moment, _H1 and _K0 are both known conditions, and the output weight of the calculation model update stage is:

During the calculation process, the intermediate variable random matrix K is introduced to realize the mapping expression of the input parameters and output parameters of the current segment, and the expression of the (k+1)th segment after each iteration is further obtained:

Calculate the aerosol optical thickness, substitute the input parameter information of the new segment into the nonlinear model, and obtain the output parameter value through neural network calculation;

Finally, through inverse normalization, the final prediction result of aerosol optical thickness is obtained.

Claims (2)

1. An online prediction method for aerosol optical thickness based on multi-source information, characterized in that it comprises the following steps: (1) Establish a multi-source information acquisition system to collect and store atmospheric environment parameters and optical parameters; (2) Extract the characteristics of atmospheric environmental parameters to obtain wind speed characteristics WS _D , temperature characteristics TE _D , humidity characteristics HU _D , and sky background brightness characteristics LU _D , which constitute the input parameter x, and the aerosol optical thickness is the output parameter y; The atmospheric environment parameter feature extraction method is as follows: The atmospheric data stored in the data acquisition system are subjected to principal component analysis, including wind speed WS, temperature TE, humidity HU, and sky background brightness LU. Considering the cumulative effect of M days, the wind speed sample fragment WS _K on the Kth day is written in the following matrix form: Among them, ws _KM represents the wind speed value of the Mth day in the Kth sample segment, and each wind speed value in the wind speed sample matrix segment WS _K is standardized: Wherein, ws′ _km represents the normalized wind speed value of the mth day in the kth sample segment, μ _m and σ _m represent the mean and variance of the wind speed sample segment of the mth day, respectively. After processing, the wind speed normalization matrix WS′ _K is obtained, and the wind speed eigenvalues λ _WS,1 ,λ _WS,2 ,...λ _WS,M corresponding to the wind speed eigenvalues η _WS,1 ,η _WS,2 ,...η _WS,M are calculated. The wind speed eigenvalues are sorted from large to small, and the contribution of the wind speed eigenvalues is calculated. When the contribution of the selected wind speed eigenvalue is greater than the threshold θ ₁ , the number of principal components Num ₁ is determined; The corresponding wind speed feature vector set is the wind speed feature Considering the cumulative effect of M days, the temperature sample fragment TE _K of the Kth day is standardized to obtain the temperature standardized matrix TE′ _K , and the temperature eigenvalues λ _TE,1 ,λ _TE,2 ,...λ _TE,M corresponding to the temperature eigenvalues η _TE,1 ,η _TE,2 ,...η _TE,M are calculated. The temperature eigenvalues are sorted from large to small, and the contribution of the temperature eigenvalues is calculated. When the contribution of the selected temperature eigenvalue is greater than the threshold θ ₂ , the number of principal components Num ₂ is determined; The corresponding temperature feature vector set is the temperature feature Considering the cumulative effect of M days, the humidity sample segment HU _K of the Kth day is standardized to obtain the humidity standardized matrix HU′ _K , and the humidity eigenvalues λ _HU,1 ,λ _HU,2 ,...λ _HU,M corresponding to the humidity eigenvalues η _HU,1 ,η _HU,2 ,...η _HU,M are calculated. The humidity eigenvalues are sorted from large to small, and the contribution of the humidity eigenvalues is calculated. When the contribution of the selected humidity eigenvalue is greater than the threshold θ ₃ , the number of principal components Num ₃ is determined. The corresponding humidity feature vector set is the humidity feature Considering the cumulative effect of M days, the sky background brightness sample fragment LU _K of the Kth day is standardized to obtain the sky background brightness standardized matrix LU′ _K , the sky background brightness eigenvalues λ _LU,1 ,λ _LU,2 ,...λ _LU,M corresponding to the sky background brightness eigenvalues η _LU,1 ,η _LU,2 ,...η _LU,M are calculated, the sky background brightness eigenvalues are sorted from large to small, the contribution of the sky background brightness eigenvalues is calculated, and when the contribution of the selected sky background brightness eigenvalue is greater than the threshold θ ₄ , the number of principal components Num ₄ is determined; The corresponding sky background brightness feature vector set is the sky background brightness feature The input parameter x is the set of wind speed feature TE _D , temperature feature TE _D , humidity feature HU _D , and sky background brightness feature LU _D , that is, x=(WS _D , TE _D , HU _D , LU _D ); (3) Using data-driven thinking, a nonlinear mapping model based on a dynamic extreme learning machine (ELM) was established, the relationship between input parameters and output parameters was determined, and the artificial bee colony algorithm was used for optimization to determine the optimal solution for the input weight and bias of the nonlinear mapping model; The construction of the nonlinear mapping model includes a three-layer neural network structure of input layer, hidden layer and output layer. Considering the newly arrived sample data, the data sampling length is N. In the initial stage, the input layer is N nodes, the hidden layer is L nodes, and the output layer is 1 node. The input parameter x and the output parameter y are sampled and merged to obtain the sample data. _Xi represents the i-th input parameter of x, and _Yi represents the i-th output parameter of y. The input-output relationship of the nonlinear mapping model is as follows: Among them, a and b are input weight and bias respectively, β is the output layer weight, L is the number of hidden layer nodes of the neural network, is the activation function. Considering the mapping relationship between nodes as a whole, according to the extreme learning machine theory, the weights and biases between the input layer nodes and the hidden layer are randomly selected. The output weight β is the only variable that needs to be solved in the entire network. Through Moore-Penrose generalized inverse calculation, formula (7) can be expressed in the following matrix form: y=Hβ·································(8) Among them, y = [y ₁ ,y ₂ ,...,y _N ] ^T represents the aerosol optical thickness of the current segment. The hidden matrix H is introduced into the nonlinear model for calculation. H is a random matrix with N rows and L columns. The specific expression is: in, To activate the function, nonlinear mapping is realized, which ensures that the function has derivatives in the interval. In order to further calculate the output weight β, the Moore-Penrose generalized inverse of H is calculated. H ⁺ =(H ^T H) ^-1 H ^T ···························(10) By solving the Moore-Penrose generalized inverse, we can obtain the optimal solution β ^* of the output weight β: β ^* =H ⁺ y····························(11) Further optimize the input weights and biases, encode the weights and biases into strings, and obtain the candidate solution vector set E = (a, b), which is in the following form: The boundaries of the candidate vector set are the maximum and minimum solution vectors, which are expressed here as E _max = {max(e)|e∈E} and E _min = {min(e)|e∈E}, where e represents a set of elements that satisfy the boundary conditions. According to the artificial bee colony algorithm, the solution vector of the input weight and bias under the initial conditions is expressed as: E _h,j =E _min,j +ω×(E _max,j -E _min,j )····················(13) E _min,j and E _max,j represent the maximum and minimum values of the jth input solution vector respectively. ω is evenly distributed in the interval [-1,1]. The optimal solution vector search is performed and the following update strategy is obtained: Among them, E _h,j represents the current solution vector, Represents a new feasible solution vector, h, j and r are all indicator marks, u _{h, j} changes randomly in the interval [-1,1], calculates the fitness, and selects the optimal solution vector: Among them, num represents the number of iterations. After the iteration is completed, the optimal solution vector will be updated to obtain the optimal input weight a ^* and bias b ^* . Therefore, the optimal solution of the expression of the nonlinear model is: (4) The nonlinear mapping model is updated online, the error between the real value and the predicted value of the real-time data segment is calculated, the nonlinear mapping model is quantitatively updated according to the error, the output weight of the nonlinear mapping model is calculated analytically, and the predicted value of the aerosol optical thickness of the next data segment is obtained to realize the online prediction of the aerosol optical thickness. Through feature extraction, nonlinear model, and model update strategy, multi-source features are integrated to realize the online prediction of the aerosol optical thickness. The nonlinear mapping model updating method is: In the initial stage, the output weights are obtained: New sample snippet After arrival, the error between the true value and the predicted value of the real-time arriving data segment is calculated, expressed as: Among them, N ₀ and N ₁ represent the starting position and end position of the current sample segment, H ₁ and y ₁ represent the hidden layer output and label set of the newly arrived segment, β ₁ represents the current nonlinear model output weight. Based on the current output weight, the expression of β ₁ is: Divide the above expression into two parts and introduce intermediate variables According to the associative law of matrix multiplication, the expression of the first part of the matrix is: In the current fragment state, H ₁ and K ₀ are both known conditions. Substituting them into the overall output layer weight expression, we can get the expression of the output weight β ₁ in the model update phase: Further iterative calculations yield the intermediate variable K _k+1 expression corresponding to the (k+1)th segment: Further iterative calculations are performed to obtain the optimal solution for the output weight β _k+1 of the nonlinear mapping model corresponding to the (k+1)th segment: The model expression after online update is: Substitute the input parameter x into the updated model and calculate the aerosol optical depth y to achieve online prediction of aerosol optical depth. 2. According to claim 1, a method for online prediction of aerosol optical depth based on multi-source information is characterized in that the multi-source information acquisition system in step (1) includes a data acquisition device, an information exchange terminal, and a database, the data acquisition device includes a laser radar, a ground meteorological station, a background radiometer, and a corresponding data interface, the database uses SQL Server to store data, the information exchange terminal exchanges information with the data acquisition device, and saves the collected data information into the database; wherein the ground meteorological station, the background radiometer and the laser radar are used to realize the collection of atmospheric environment parameters and optical parameters, the atmospheric environment parameters include the atmospheric data wind speed WS, temperature TE, humidity HU provided by the ground meteorological station, and the sky background brightness LU collected by the background radiometer, and the optical parameter is the aerosol optical thickness, which is collected by the laser radar.