文章来源：基于激光雷达测风在距离加权下的风速重构算法 - 中国知网 (cnki.net)

本文针对文献Real-time three dimensional wind field reconstruction from nacelle LiDAR measurements中提出的“在计算出多个测量平面对应的转子有效风速的基础上，利用求和取平均值的方法计算风轮面转子有效风速”这一方法，认为其计算是不准确的。

贝兹(Betz)极限表明风机从自然风中所能索取的能量是有限的，风机叶片不能全部吸收来自流风的能量，由于空气流动具有连续性的特点，风速从风机前方经过风轮面流向风机后方，经过风轮叶片的吸收，风速经历的是逐渐降低的过程，不是一个瞬时衰减的过程。所以雷达所测10个测量距离上的风速对于计算风轮面转子有效风速所贡献的量是不同的，本文在此基础上开发出基于最小二乘法求解10个测量距离上对应的转子有效风速的加权因子进而拟合出风轮面转子有效风速的方法，相比于求和取平均值来讲，提高了计算风轮面转子有效风速的精确性。

第 i i i光束在第 j j j个测量距离上的测点纵向坐标、横向坐标、竖向坐标可以表示为： [ x i , j y i , j z i , j ] = L j [ 1 tan ⁡ ( θ 1 ) cos ⁡ ( θ 2 2 + i ∗ π ) tan ⁡ ( θ 1 ) sin ⁡ ( θ 2 2 + i ∗ π ) ] \left[\begin{array}{l} x_{i, j} \\ y_{i, j} \\ z_{i, j} \end{array}\right]=L_j\left[\begin{array}{c} 1 \\ \tan \left(\theta_1\right) \cos \left(\frac{\theta_2}{2}+i^* \pi\right) \\ \tan \left(\theta_1\right) \sin \left(\frac{\theta_2}{2}+i^* \pi\right) \end{array}\right] ​xi,j​yi,j​zi,j​​ ​=Lj​ ​1tan(θ1​)cos(2θ2​​+i∗π)tan(θ1​)sin(2θ2​​+i∗π)​ ​ 其中 L j L_j Lj​为激光雷达的测量距离。 i = 0 , 1 , 2 , 3 i=0,1,2,3 i=0,1,2,3， j = 0 ⋯ 9 j=0 \cdots9 j=0⋯9，表示10个测量距离。

激光雷达测风数据是以激光雷达坐标系为参考输出的，雷达会跟随机舱振动，导致实测数据产生测量误差，由此将机舱激光雷达的振动分解为3个方向的运动。

设绕Z轴旋转的为偏航角 Ψ L \Psi_L ΨL​，绕X轴旋转的为翻滚角 Φ L \Phi_L ΦL​，绕Y轴旋转的为俯仰角 Θ L \Theta_L ΘL​。建立如下的惯性坐标系。

将激光雷达测风数据由激光雷达坐标系向惯性坐标系变换来进行数值修正，消除机舱振动带来的测量误差。设从激光雷达坐标系到惯性坐标系的转换矩阵为 T I L T_{IL} TIL​，其定义如下： T I L = [ cos ⁡ ( Ψ L ) − sin ⁡ ( Ψ L ) 0 sin ⁡ ( Ψ L ) cos ⁡ ( Ψ L ) 0 0 0 1 ] T yaw  [ cos ⁡ ( Θ L ) 0 sin ⁡ ( Θ L ) 0 1 0 − sin ⁡ ( Θ L ) 0 cos ⁡ ( Θ L ) ] T p i t c h [ 1 0 0 0 cos ⁡ ( Φ L ) sin ⁡ ( Φ L ) 0 − sin ⁡ ( Φ L ) cos ⁡ ( Φ L ) ] T roll  , \begin{aligned} & T_{\mathrm{IL}}=\left[\begin{array}{ccc} \cos \left(\Psi_{\mathrm{L}}\right) & -\sin \left(\Psi_{\mathrm{L}}\right) & 0 \\ \sin \left(\Psi_{\mathrm{L}}\right) & \cos \left(\Psi_{\mathrm{L}}\right) & 0 \\ 0 & 0 & 1 \end{array}\right]_{T_{\text {yaw }}}\left[\begin{array}{ccc} \cos \left(\Theta_{\mathrm{L}}\right) & 0 & \sin \left(\Theta_{\mathrm{L}}\right) \\ 0 & 1 & 0 \\ -\sin \left(\Theta_{\mathrm{L}}\right) & 0 & \cos \left(\Theta_{\mathrm{L}}\right) \end{array}\right]_{T_{pitch}} \\ & {\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos \left(\Phi_{\mathrm{L}}\right) & \sin \left(\Phi_{\mathrm{L}}\right) \\ 0 & -\sin \left(\Phi_{\mathrm{L}}\right) & \cos \left(\Phi_{\mathrm{L}}\right) \end{array}\right]_{T_{\text {roll }}}, } \end{aligned} ​TIL​= ​cos(ΨL​)sin(ΨL​)0​−sin(ΨL​)cos(ΨL​)0​001​ ​Tyaw ​​ ​cos(ΘL​)0−sin(ΘL​)​010​sin(ΘL​)0cos(ΘL​)​ ​Tpitch​​ ​100​0cos(ΦL​)−sin(ΦL​)​0sin(ΦL​)cos(ΦL​)​ ​Troll ​​,​ 这三个矩阵分别为偏航修正矩阵，俯仰修正矩阵，翻滚修正矩阵。

惯性坐标系下的测点坐标为： [ x i , j , I y i , j , I z i , j , I ] = T I L [ x i , j ,   L y i , j ,   L z i , j ,   L ] + [ x L , I y L , I z L , I ] \left[\begin{array}{l} x_{i, j, \mathrm{I}} \\ y_{i, j, \mathrm{I}} \\ z_{i, j, \mathrm{I}} \end{array}\right]=T_{\mathrm{IL}}\left[\begin{array}{l} x_{i, j, \mathrm{~L}} \\ y_{i, j, \mathrm{~L}} \\ z_{i, j, \mathrm{~L}} \end{array}\right]+\left[\begin{array}{l} x_{\mathrm{L}, \mathrm{I}} \\ y_{\mathrm{L}, \mathrm{I}} \\ z_{\mathrm{L}, \mathrm{I}} \end{array}\right] ​xi,j,I​yi,j,I​zi,j,I​​ ​=TIL​ ​xi,j, L​yi,j, L​zi,j, L​​ ​+ ​xL,I​yL,I​zL,I​​ ​ 通过建立激光雷达测量坐标系与惯性坐标系的转换模型，修正了机舱振动带来的数据波动，提高了风速重构的精确性。为解决叶片遮挡的问题，采用以上一秒的有效数据填充下一秒缺失或无效的数据的方法。

基于动态模型的风场重建思想，来得到风场参数，包括在第 j j j个测量距离平面的平均风速 u 0 , j u_{0,j} u0,j​、水平剪切 s h , j s_{h,j} sh,j​、竖直剪切 s v , j s_{v,j} sv,j​，基于这3个风场参数，即可求得在这一测量平面任意一点的速度。 v i , j , l o s = x i , j , I f i , j u i , j , I + y i , j , I f i , j v i , j , I + z i , j , I f i , j w i , j , I v_{i, j, \mathrm{los}}=\frac{x_{i, j, \mathrm{I}}}{f_{i, j}} u_{i, j, \mathrm{I}}+\frac{y_{i, j, \mathrm{I}}}{f_{i, j}} v_{i, j, \mathrm{I}}+\frac{z_{i, j, \mathrm{I}}}{f_{i, j}} w_{i, j, \mathrm{I}} vi,j,los​=fi,j​xi,j,I​​ui,j,I​+fi,j​yi,j,I​​vi,j,I​+fi,j​zi,j,I​​wi,j,I​

f i , j = x i , j , I 2 + y i , j , I 2 + z i , j , I 2 v ⃗ i , j ( y i , j , z i , j , t k ) = ( u 0 , j ( t k ) + s h , j ( t k ) y i , j + s v , j ( t k ) z i , j ) [ 1 0 0 ] \begin{gathered} f_{i, j}=\sqrt{x_{i, j, \mathrm{I}}^2+y_{i, j, \mathrm{I}}^2+z_{i, j, \mathrm{I}}^2} \\ \vec{v}_{i, j}\left(y_{i, j}, z_{i, j}, t_k\right)=\left(u_{0, j}\left(t_k\right)+s_{\mathrm{h}, j}\left(t_k\right) y_{i, j}+s_{\mathrm{v}, j}\left(t_k\right) z_{i, j}\right)\left[\begin{array}{l} 1 \\ 0 \\ 0 \end{array}\right] \end{gathered} fi,j​=xi,j,I2​+yi,j,I2​+zi,j,I2​ ​v i,j​(yi,j​,zi,j​,tk​)=(u0,j​(tk​)+sh,j​(tk​)yi,j​+sv,j​(tk​)zi,j​) ​100​ ​​

s h , j ( t k ) s_{\mathrm{h}, j}\left(t_k\right) sh,j​(tk​)是 t k t_k tk​时刻第 j j j个测量距离上的水平剪切系数， s v , j ( t k ) s_{\mathrm{v},j}\left(t_k\right) sv,j​(tk​)是 t k t_k tk​时刻第 j j j个测量距离上的竖直剪切系数。

将上式联立可得： v ⃗ i , j ( y i , j , z i , j , t k ) = f i , j x i , j , I v i , j , l o s ( t k ) [ 1 0 0 ] \vec{v}_{i, j}\left(y_{i, j}, z_{i, j}, t_k\right)=\frac{f_{i, j}}{x_{i, j, \mathrm{I}}} v_{i, j, \mathrm{los}}\left(t_k\right)\left[\begin{array}{l} 1 \\ 0 \\ 0 \end{array}\right] v i,j​(yi,j​,zi,j​,tk​)=xi,j,I​fi,j​​vi,j,los​(tk​) ​100​ ​

f i , j x i , j , I v i , j ,  los  ( t k ) = u 0 , j ( t k ) + s h , j ( t k ) y i , j + s v , j ( t k ) z i , j \frac{f_{i, j}}{x_{i, j, \mathrm{I}}} v_{i, j, \text { los }}\left(t_k\right)=u_{0, j}\left(t_k\right)+s_{\mathrm{h}, j}\left(t_k\right) y_{i, j}+s_{\mathrm{v}, j}\left(t_k\right) z_{i, j} xi,j,I​fi,j​​vi,j, los ​(tk​)=u0,j​(tk​)+sh,j​(tk​)yi,j​+sv,j​(tk​)zi,j​

然后通过线性最小二乘法就能求出 t k t_k tk​时刻 j j j测量距离的纵向平均风速、竖直剪切、水平剪切系数。 [ f i , j x i , j , I v 1 , j , l o s ( t k ) ⋮ f i , j x i , j , I v i , j , l o s ( t k ) ] ⏟ m = [ 1 y 1 , j , I z 1 , j , I ⋮ ⋮ ⋮ 1 y i , j , I z i , j , I ] ⏟ A [ u i , j ( t k ) s h , j ( t k ) s v , j ( t k ) ] ⏟ s \underbrace{\left[\begin{array}{c} \frac{f_{i, j}}{x_{i, j, \mathrm{I}}} v_{1, j, \mathrm{los}}\left(t_k\right) \\ \vdots \\ \frac{f_{i, j}}{x_{i, j, \mathrm{I}}} v_{i, j, \mathrm{los}}\left(t_k\right) \end{array}\right]}_m=\underbrace{\left[\begin{array}{ccc} 1 & y_{1, j, \mathrm{I}} & z_{1, j, \mathrm{I}} \\ \vdots & \vdots & \vdots \\ 1 & y_{i, j, \mathrm{I}} & z_{i, j, \mathrm{I}} \end{array}\right]}_A \underbrace{\left[\begin{array}{l} u_{i, j}\left(t_k\right) \\ s_{\mathrm{h}, j}\left(t_k\right) \\ s_{\mathrm{v}, j}\left(t_k\right) \end{array}\right]}_s m ​xi,j,I​fi,j​​v1,j,los​(tk​)⋮xi,j,I​fi,j​​vi,j,los​(tk​)​ ​​​=A ​1⋮1​y1,j,I​⋮yi,j,I​​z1,j,I​⋮zi,j,I​​ ​​​s ​ui,j​(tk​)sh,j​(tk​)sv,j​(tk​)​ ​​​ 这样就可以求出该平面上任意一点的风速。

提出基于最小二乘先进先出循环嵌套算法，计算出不同测量平面对应不同的加权因子，进而计算出更加准确的风轮面转子有效风速。

定义10个加权因子 [ μ 0 , ⋯   , μ 9 ] [\mu_0,\cdots,\mu_9] [μ0​,⋯,μ9​]

建立线性方程组： [ v reus  , 0 , t 1 ⋯ v reus  , 9 , t 1 ⋮ ⋮ ⋮ v reus  , 0 , t a ⋯ v reus  , 9 , t a ] [ μ 0 ⋮ μ 9 ] = [ v true  , t 1 ⋮ v true  , t a ] \left[\begin{array}{ccc} v_{\text {reus }, 0, t 1} & \cdots & v_{\text {reus }, 9, t 1} \\ \vdots & \vdots & \vdots \\ v_{\text {reus }, 0, t a} & \cdots & v_{\text {reus }, 9, t a} \end{array}\right]\left[\begin{array}{c} \mu_0 \\ \vdots \\ \mu_9 \end{array}\right]=\left[\begin{array}{c} v_{\text {true }, t 1} \\ \vdots \\ v_{\text {true }, t a} \end{array}\right] ​vreus ,0,t1​⋮vreus ,0,ta​​⋯⋮⋯​vreus ,9,t1​⋮vreus ,9,ta​​ ​ ​μ0​⋮μ9​​ ​= ​vtrue ,t1​⋮vtrue ,ta​​ ​ 利用最小二乘法得到 μ \mu μ的值。

然后对下一个时间戳的10个平面的风速进行加权求和，得到有效风速。

对于下一个时间戳就将窗口向下移动一格，然后重新计算 μ \mu μ，然后计算风速。