A laser radar synthetic wind speed and direction processing method based on Kalman filtering

Abstract

The application discloses a kind of laser radar synthetic wind speed and wind direction processing method based on Kalman filtering, through the data detected by laser radar is calculated by Kalman filtering algorithm, specifically for the weight of sample point and sample point set one-step prediction, system state weighted mean and covariance matrix, measurement point set one-step prediction, measurement system weighted mean and covariance matrix, Kalman filtering gain and update system state value, update covariance value calculation, with prior state input and prior covariance matrix selects sample, obtains new system state weighted mean and covariance as the parameter of Kalman filtering by the nonlinear transformation of inversion algorithm, under the premise that the mean of wind speed and wind direction inversion error at hub is zero and unchanged, reduce the standard deviation of inversion error.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of wind measurement signal processing, and particularly relates to a laser radar synthetic wind speed and direction processing method based on Kalman filtering. BACKGROUND

[0002] In the prior art, a nacelle type laser radar can provide atmospheric wind field information at different distances by coherently detecting wind speed in different directions, and calculating the horizontal wind speed and direction at the hub, thereby improving the utilization rate of wind energy of a wind turbine and enhancing the safety of the wind turbine during operation. However, there is still a certain error between the calculation results of the existing hub wind speed and direction inversion algorithm and the theoretical values of the wind speed and direction, because the linear wind shear model used in the existing method cannot accurately describe the nonlinear wind shear phenomenon, resulting in the accumulation effect of measurement outliers and noise.

[0003] Although various frequency domain filtering techniques and estimation filtering techniques can filter out noise at corresponding frequencies, and can eliminate outliers and predict reasonable inversion data, due to the nonlinear effect of vertical wind shear and the difficulty of effective data prediction, the existing linear filtering method cannot achieve good processing effect.

[0004] Therefore, there is a need for a new hub wind speed and direction data processing method to meet the needs in practical applications. SUMMARY

[0005] The application aims to provide a laser radar synthetic wind speed and direction processing method based on Kalman filtering, which can solve the technical problems of the prior art, such as inapplicability to nonlinear systems, large measurement error, and low outlier correction efficiency.

[0006] The technical solution of the application to solve the technical problem is as follows:

[0007] A laser radar synthetic wind speed and direction processing method based on Kalman filtering, specifically comprising the following steps:

[0008] S1: reading the laser radar storage data in a unit of time, extracting the values of the laser radar lower plane swing angle, lower plane roll angle, upper plane swing angle, upper plane roll angle, the corresponding wind direction of the four emitted beams, and the detection distance L as the state input value X of this time;

[0009] S2: determining whether the filtering number of the state input value of this time is initial, if yes, performing initial filtering on the state input value of this time, obtaining the state output Y as the initial filtering result by using the state transition matrix F, and obtaining the measurement value Z by using the state output Y and the measurement matrix H; if no, going to step S3;

[0010] S3: take the previous state input value as the current state input value, input to the Kalman core algorithm; filter the current state input value to obtain the measurement value of the wind speed and direction at the hub, and input to the Kalman core algorithm; input the covariance matrix of the current state to the Kalman core algorithm;

[0011] S4: the Kalman core algorithm calculates the input quantity in step S3, and the specific calculation steps are as follows:

[0012] S4.1: calculate the weight value corresponding to the sample point, and the calculation formula is:

[0013] λ=α 2 (n+κ)-n, wherein m is the mean symbol, c is the covariance symbol, λ is the scaling parameter, κ is the to-be-selected parameter, α is the parameter for controlling the distribution state of the sample point, β is the non-negative weight coefficient, is the weight coefficient corresponding to the mean of the first sample point; is the weight coefficient corresponding to the covariance of the first sample point; is the weight coefficient corresponding to the mean of the i-th sample point; is the weight coefficient corresponding to the covariance of the i-th sample point; n is the dimension of the state input value X;

[0014] S4.2: calculate 2n+1 sets of sample point sets, and the calculation formula is:

[0015] wherein, is the i-th column of the sample point set obtained by adding the i-th column of the square root of the state input value X and the covariance matrix P, X(k-1) is the state input value X at k-1 time, P(k-1) is the covariance matrix at k-1 time, indicates the i-th column of the square root of the covariance matrix, and the sample point set in the vector form is obtained by arranging the sample point formula:

[0016]

[0017] S4.3: calculate one-step prediction of 2n+1 sets of sample point sets, and the calculation formula is:

[0018] X set (k|k-1)=F(k|k-1)·X set (k-1), wherein F is a state transition matrix, X set (k|k-1) is the sample point set prediction result at k time obtained by the sample point set at k-1 time and the state transition matrix F;

[0019] S4.4: calculate the system state weighted mean and the covariance matrix, and the calculation formula is:

[0020]

[0021]

[0022] wherein, Q is the process covariance matrix, X(k|k-1) is the mean of the sample point set at k time obtained by the prediction result of the sample point set at k time and the mean weighting coefficient, P xx (k|k-1) is the autocovariance matrix at k time obtained by the prediction result of the sample point set at k time and the covariance weighting coefficient and the process covariance matrix Q;

[0023] S4.5: the one-step prediction of the sample point set of the 2n+1 groups is brought into the measurement equation to obtain the measurement value, and the calculation formula is:

[0024] Z set (k|k-1) = H(k)·X set (k|k-1); wherein, H is the measurement matrix, Z set (k|k-1) is the measurement value point set at k time obtained by the prediction result of the sample point set at k time and the measurement matrix;

[0025] S4.6: the prediction mean and covariance of the measurement system are obtained by weighting, and the calculation formula is:

[0026]

[0027]

[0028]

[0029] wherein, R is the measurement covariance matrix, Z(k|k-1) is the mean of the measurement value point set at k time obtained by the measurement value point set at k time and the mean weighting coefficient, P xz (k|k-1) is the cross-covariance at k time obtained by the prediction result of the sample point set at k time, the mean of the sample point set at k time, the measurement value point set at k time, the mean of the measurement value point set at k time and the covariance weighting coefficient, P zz (k|k-1) is the autocovariance at k time obtained by the measurement value point set at k time, the mean of the measurement value point set at k time and the covariance weighting coefficient;

[0030] S4.7: the Kalman filter gain is calculated, and the calculation formula is:

[0031] K(k) is the gain at k time,

[0032] S4.8: the system output result is calculated, including the filtering result and the covariance update value, and the calculation formula is:

[0033] Xkf(k) = X(k|k-1) + K(k)[Z(k) - Z(k|k-1)],

[0034] Xkf(k) is the filtering result of the core algorithm at k time, P(k) is the covariance matrix at k time, P(k|k-1) is the covariance prediction matrix predicted from k-1 time to k time;

[0035] S4.9: judge whether the length of the filtered signal reaches the set value, if yes, output the filtering result; if no, add the filtering times, and the system covariance update value is taken as the new state input value, and steps S1-S3 are repeated.

[0036] The specific steps of filtering the current state input value in step S3 to obtain the measurement value of the wind speed and the wind direction at the hub are as follows:

[0037] D1: substituting the current state input value X into the state transition equation F of the Kalman filter to obtain the state output Y of the upper plane wind speed, the lower plane wind speed, the vertical wind shear, the upper plane wind direction, the lower plane wind direction, the vertical wind direction change rate, the wind speed at the hub, the wind direction at the hub and the detection distance;

[0038] D2: substituting the current state output Y in D1 into the measurement equation H of the Kalman filter to obtain the measurement value Z of the wind speed and the wind direction at the hub.

[0039] The state transition equation of the Kalman filter is: Y(k) = F(k) * X(k), wherein Y(k) is the state output at k time; F(k) is the state transition matrix; and X(k) is the state input at k time.

[0040] The measurement equation of the Kalman filter is: Z(k) = H(k) * Y(k), wherein Z(k) is the measurement value at k time; and H(k) is the measurement matrix.

[0041] The value of the alpha is 0.01.

[0042] The beneficial effects of the application are as follows: the data detected by the laser radar is calculated through the Kalman filter algorithm, specifically, the weight of the sample point and the one-step prediction of the sample point set, the weighted mean and the covariance matrix of the system state, the one-step prediction of the measurement point set, the weighted mean and the covariance matrix of the measurement system, the Kalman filter gain and the update of the system state value and the covariance value are calculated, the samples are selected from the prior state input and the prior covariance matrix, the new weighted mean and the covariance of the system state are obtained through the nonlinear transformation of the inversion algorithm as the parameters of the Kalman filter, under the premise that the inversion error mean of the wind speed and the wind direction at the hub is zero and unchanged, the standard deviation of the inversion error is reduced. Moreover, the method is suitable for nonlinear systems, and the processing of abnormal values is more efficient. BRIEF DESCRIPTION OF DRAWINGS

[0043] Figure 1 is a schematic diagram of a wind measuring laser radar beam structure referred to by the present application;

[0044] Figure 2 is a flow chart of a laser radar combined wind speed and wind direction processing method of the present application. DETAILED DESCRIPTION

[0045] The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application.

[0046] As shown in Figure 1 , Figure 2 , the present application comprises a laser radar combined wind speed and wind direction processing method based on Kalman filtering, specifically comprising the following steps:

[0047] S1: reading the laser radar storage data in a unit time, extracting the values of the lower plane swing angle, the lower plane roll angle, the upper plane swing angle, the upper plane roll angle, the four corresponding view direction wind speed and the detection distance L of the laser radar as the state input value X of this time;

[0048] S2: judging whether the filtering number of the state input value of this time is initial, if yes, performing initial filtering on the state input value of this time, obtaining the state output Y by using the state transition matrix F as the initial filtering result, and obtaining the measurement value Z by using the state output Y and the measurement matrix H; if no, entering step S3;

[0049] S3: taking the state input value of the previous time as the current state input value and inputting it into the Kalman core algorithm; filtering the state input value of this time to obtain the measurement value of the hub wind speed and wind direction and inputting it into the Kalman core algorithm; inputting the covariance matrix of the state of this time into the Kalman core algorithm; initializing the state value generally as X(1) and initializing the covariance matrix generally as P(1).

[0050] The specific steps of filtering the state input value of this time to obtain the measurement value of the hub wind speed and wind direction in step S3 are as follows:

[0051] D1: substituting the state input value X of this time into the state transition equation F of the Kalman filtering to obtain the state output Y of the upper plane wind speed, the lower plane wind speed, the vertical wind shear, the upper plane wind direction, the lower plane wind direction, the vertical wind direction change rate, the hub wind speed, the hub wind direction and the detection distance;

[0052] The state transition equation of the Kalman filtering is: Y(k)=F(k)·X(k), wherein Y(k) is the state output at time k; F(k) is the state transition matrix; X(k) is the state input at time k.

[0053] D2: the current state output Y in D1 is substituted into the Kalman filter measurement equation H to obtain the measurement value Z of the wheel hub wind speed and wind direction.

[0054] The Kalman filter measurement equation is: Z(k) = H(k) Y(k), wherein Z(k) is the measurement value at k time; H(k) is a measurement matrix.

[0055] S4: the Kalman core algorithm calculates the input quantity in step S3, and the specific calculation steps are as follows:

[0056] S4.1: the weight of the sample point is calculated, and the calculation formula is:

[0057] λ = α 2 (n + κ) - n, wherein m is a mean symbol, c is a covariance symbol, λ is a scaling parameter used to reduce the total prediction error; κ is a to-be-selected parameter, which has no specific value limit, but generally should ensure that (n + λ) P is a semi-definite matrix, α is a parameter for controlling the distribution state of the sample point, and the value of the α is 0.01; β is a non-negative weight coefficient, which can combine the dynamic range of high-order terms in the equation, so that the influence of the high-order terms can be included; is the weight coefficient corresponding to the mean value of the first sample point; is the weight coefficient corresponding to the covariance of the first sample point; is the weight coefficient corresponding to the mean value of the i-th sample point; is the weight coefficient corresponding to the covariance of the i-th sample point; n is the dimension of the state input value X;

[0058] S4.2: 2n + 1 sets of sample point sets are calculated, and the calculation formula is:

[0059] wherein, is the i-th column of the sample point set obtained by adding the i-th column of the square root of the state input value X and the covariance matrix P, X(k-1) is the state input value X at k-1 time, P(k-1) is the covariance matrix at k-1 time, indicates the i-th column of the square root of the covariance matrix, and the sample point formula is arranged to obtain a vector form of the sample point set:

[0060]

[0061] S4.3: one-step prediction of 2n + 1 sets of sample point sets is calculated, and the calculation formula is:

[0062] X set (k|k-1) = F(k|k-1) X set(k|k-1), wherein F is a state transition matrix, that is, a matrix expression of a wind speed and wind direction inversion algorithm at a hub;

[0063] X set (k|k-1) is a sample point set prediction result at the k moment obtained through the sample point set at the k-1 moment and the state transition matrix F;

[0064] S4.4: Calculate the weighted mean and covariance matrix of the system state, and the calculation formula is:

[0065]

[0066]

[0067] wherein Q is a process covariance matrix, X(k|k-1) is a sample point set mean at the k moment obtained through the sample point set prediction result at the k moment and the mean weighted coefficient, P xx (k|k-1) is a self-covariance matrix at the k moment obtained through the sample point set prediction result at the k moment, the covariance weighted coefficient and the process covariance matrix Q;

[0068] S4.5: Bring the one-step prediction of the sample point set of the 2n+1 groups into the measurement equation to obtain the measurement value, and the calculation formula is:

[0069] Z set (k|k-1) = H(k)·X set (k|k-1); wherein H is a measurement matrix, Z set (k|k-1) is a measurement value point set at the k moment obtained through the sample point set prediction result at the k moment and the measurement matrix;

[0070] S4.6: Obtain the measurement system prediction mean and covariance through weighting, and the calculation formula is:

[0071]

[0072]

[0073]

[0074] wherein R is a measurement covariance matrix, Z(k|k-1) is a measurement value point set mean at the k moment obtained through the measurement value point set at the k moment and the mean weighted coefficient, P xz (k|k-1) is a cross-covariance at the k moment obtained through the sample point set prediction result at the k moment, the sample point set mean at the k moment, the measurement value point set at the k moment, the measurement value point set mean at the k moment and the covariance weighted coefficient, P zz(k|k-1) is the self-covariance at k time obtained by the k time measurement point set, the k time measurement point set mean value and the covariance weighted coefficient;

[0075] S4.7: Calculate the Kalman filter gain, the calculation formula is:

[0076] K(k) is the gain at k time,

[0077] S4.8: Calculate the system output result, including the filtering result and the covariance update value, the calculation formula is:

[0078] Xkf(k)=X(k|k-1)+K(k)[Z(k)-Z(k|k-1)],

[0079] Xkf(k) is the filtering result of the core algorithm at k time, P(k) is the covariance matrix at k time, and P(k|k-1) is the covariance prediction matrix predicted from k-1 time to k time;

[0080] S4.9: Determine whether the length of the filtering signal reaches the set value, if yes, output the filtering result, if no, add the filtering times, the system covariance update value is used as the new state input value, and steps S1-S3 are repeated.

[0081] The hub wind speed and the hub wind direction in the output Xkf(k) are the results processed by the method, which is more accurate to estimate the state of the nonlinear system by a group of sample points based on the Kalman filtering method based on the sample probability density distribution, the Kalman filtering method based on the sample probability density distribution is more flexible and efficient in error processing, can accurately estimate the covariance of the error, and apply the information to state estimation and prediction, so that the subsequent calculation result is more accurate.

[0082] The Kalman filtering method based on the sample probability density distribution in the application transmits the uncertainty in the Gaussian distribution to the sample point set, and the transmission of the uncertainty is more accurate. Compared with the existing filtering method, the standard deviation of the error can be more effectively reduced under the premise that the error mean is zero and unchanged.

[0083] Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor belong to the scope of protection of the application.

Claims

1. A method for processing synthetic wind speed and direction data from a lidar system based on Kalman filtering, characterized in that, Specifically, the following steps are included: S1: Read the data stored in the lidar within a unit of time, and extract the values ​​L of the lidar's lower plane sway angle, lower plane roll angle, upper plane sway angle, upper plane roll angle, line-of-sight wind speed corresponding to the four emitted beams, and detection distance as the current state input value X; S2: Determine whether the current state input value is filtered for the first time. If so, perform the initial filtering on the current state input value, use the state transition matrix F to obtain the state output Y as the initial filtering result, and obtain the measured value Z through the state output Y and the measurement matrix H; otherwise, proceed to step S3. S3: Take the previous state input value as the current state input value and input it into the Kalman core algorithm; filter the current state input value to obtain the measured values ​​of wind speed and wind direction at the hub and input them into the Kalman core algorithm; input the covariance matrix of the current state into the Kalman core algorithm. S4: The Kalman core algorithm calculates the input from step S3. The specific calculation steps are as follows: S4.1: Calculate the corresponding weights of the sample points using the following formula: λ=α 2 (n+κ)-n, where m is the mean sign, c is the covariance sign, λ is the scaling parameter, κ is a candidate parameter, α is a parameter controlling the distribution of sample points, and β is a non-negative weighting coefficient. The weight coefficients are the values ​​corresponding to the mean of the first sample point. These are the weighting coefficients corresponding to the covariance of the first sample point; The weight coefficients corresponding to the mean of the i-th sample point; is the weight coefficient corresponding to the covariance of the i-th sample point; n is the dimension of the state input value X; S4.2: Calculate the 2n+1 sets of sample points using the following formula: in, Let X(k-1) be the i-th column of the sample point set obtained by adding the i-th column of the square root of the state input value X to the square root of the covariance matrix P, where X(k-1) is the state input value X at time k-1, and P(k-1) is the covariance matrix at time k-1. The i-th column represents the root of the covariance matrix. The sample point set is then rearranged using the formula to obtain a vector form. S4.3: Calculate the one-step prediction for the 2n+1 sample point set. The calculation formula is as follows: X set (k|k-1)=F(k|k-1)·X set (k-1), where F is the state transition matrix, X set (k|k-1) represents the prediction result of the sample point set at time k obtained from the sample point set at time k-1 and the state transition matrix F; S4.4: Calculate the system state-weighted mean and covariance matrix. The calculation formula is as follows: Where Q is the process covariance matrix, X(k|k-1) is the mean of the sample point set at time k obtained by using the prediction results of the sample point set at time k and the mean weighting coefficient, and P xx (k|k-1) is the autocovariance matrix at time k obtained by using the prediction results of the sample point set at time k, the covariance weighting coefficients, and the process covariance matrix Q; S4.5: Substitute the one-step prediction of the sample point set of group 2n+1 into the measurement equation to obtain the measured value. The calculation formula is as follows: Z set (k|k-1)=H(k)·X set (k|k-1); where H is the measurement matrix, Z set (k|k-1) represents the set of measurement values ​​at time k obtained by predicting the sample point set at time k and using the measurement matrix; S4.6: The predicted mean and covariance of the measurement system are obtained through weighted average. The calculation formula is as follows: Where R is the measurement covariance matrix, Z(k|k-1) is the mean of the measurement point set at time k obtained by weighting the measurement point set at time k and the mean, and P xz (k|k-1) represents the cross-covariance at time k, obtained by using the prediction results of the sample point set at time k, the mean of the sample point set at time k, the measured value point set at time k, the mean of the measured value point set at time k, and the covariance weighting coefficient. P zz (k|k-1) is the autocovariance at time k obtained by using the set of measurement points at time k, the mean of the set of measurement points at time k, and the weighted coefficient of covariance; S4.7: Calculate the Kalman filter gain using the following formula: K(k) is the gain at time k. S4.8: Calculate the system output, including the filtering results and the updated covariance value. The calculation formula is as follows: Xkf(k)=X(k|k-1)+K(k)[Z(k)-Z(k|k-1)], Xkf(k) is the filtering result of the core algorithm at time k, P(k) is the covariance matrix at time k, and P(k|k-1) is the covariance prediction matrix for time k predicted from time k-1. S4.9: Determine whether the length of the filtered signal has reached the set value. If yes, output the filtering result; if no, increment the number of filtering iterations, use the updated system covariance value as the new state input value, and repeat steps S1-S3.

2. The method for processing synthetic wind speed and direction using lidar based on Kalman filtering according to claim 1, characterized in that: The specific steps in step S3 to filter the current state input value to obtain the measured values ​​of wind speed and wind direction at the hub are as follows: D1: Substitute the current state input value X into the state transition equation F of the Kalman filter to obtain the state output Y of the upper plane wind speed, lower plane wind speed, vertical wind shear, upper plane wind direction, lower plane wind direction, vertical wind direction change rate, wind speed at the hub, wind direction at the hub, and detection distance. D2: Substitute the current state output Y from D1 into the Kalman filter measurement equation H to obtain the measured values ​​Z of wind speed and wind direction at the hub.

3. The method for processing synthetic wind speed and direction using lidar based on Kalman filtering according to claim 2, characterized in that: The state transition equation of the Kalman filter is: Y(k)=F(k)·X(k), where Y(k) is the state output at time k; F(k) is the state transition matrix; and X(k) is the state input at time k.

4. The method for processing synthetic wind speed and direction using lidar based on Kalman filtering according to claim 2, characterized in that: The Kalman filter measurement equation is: Z(k) = H(k)·Y(k), where Z(k) is the measurement value at time k; H(k) is the measurement matrix.

5. The method for processing synthetic wind speed and direction using lidar based on Kalman filtering according to claim 1, characterized in that: The value of α is 0.01.

Images