1 引言

大气风场是重要的大气物理参数,精确的大气风场探测对防治大气污染^[1]、提高航空航天的安全性^[2]、预报**环境^[3]、改进气候研究模型、提高长期天气预报的准确性^[4]等具有重大意义。测风激光雷达作为主动大气遥感仪器,具有风场探测精度高、时间和空间分辨率高,以及不受地物杂波影响等优点^[5],非常适用于快速准确的风场探测。

现有测风激光雷达一般包括多普勒测风激光雷达和相关方法测风激光雷达。多普勒测风激光雷达通过测量随风场运动的大气分子和气溶胶粒子对激光造成的多普勒频移实现风场探测,根据探测实现方式的不同分为相干探测和直接探测。多普勒相干探测利用回波信号和本振激光进行光学混频^[6-7],差频信号的频率减去起始的频率偏置就是由运动粒子引起的多普勒频移,相干探测对激光的相干性要求较高,所以光学器件加工难度大、成本高、光路复杂。多普勒直接探测则是利用光学鉴频器或光谱分析仪将光频率移动转化为可直接测量的光功率或光功率空间分布的变化^[8-10],对激光器频率稳定性要求很高,系统复杂。

相关方法测风激光雷达只是基于气溶胶不均匀结构造成的回波信号强度起伏来反演风速^[11],对光源相干性要求不高,对光谱宽度限制也不严格,光路简单,可以在简单、稳定的系统下实现大气风场探测。目前相关方法测风激光雷达采用的测风方法主要有信号直接相关^[12-13]和信号合成图像相关^[14-16]两种。其中信号直接相关方法是直接根据时间序列上的回波信号强度曲线进行相关计算,方法简单,但是只能反演特定方向的风速信息;信号合成图像相关方法是先将激光雷达信号合成为时空演变图,再对时空演变图进行相关处理,此方法可以反演平面任意方向的风速信息,但是涉及到图像的处理,运算复杂,不利于风场的快速探测。

本文提出了基于激光雷达信号相关计算的三光束二维风场反演方法,可以通过三个同一平面发射的激光雷达,反演探测平面上任意方向的风速信息,方法简单。基于三光束二维风场反演方法,使用三台激光雷达搭建了相关二维测风激光雷达系统。系统所采用激光雷达的波长为532 nm,单脉冲能量为10 μJ,重复频率为2.5 kHz,可以实现30 m距离分辨率、1 s时间分辨率的二维风场测量。利用该系统进行了水平风场探测实验,并与探测路径上的风速、风向传感器进行了对比。

2 基于相关计算的激光雷达二维风场探测方法

2.1 基本原理

近地层大气中充满了大小不等的气溶胶微粒,在大气湍流的作用下,悬浮在大气中的气溶胶出现不均匀结构^[17],从而形成大量浓度和尺度不同的气溶胶团块。这些浓度和尺度不等并在大气中随风飘移的气溶胶团块,使激光雷达回波出现幅度和尺度的不均匀起伏^[18-20],从而成为激光雷达相关计算测风的探测信息。

米散射激光雷达的大气后向散射回波信号能量^[21]可表示为

Pr(z)=PtCz-2[βa(z)+βm(z)]exp·-2∫0z[αa(z')+αm(z')]dz',(1)

式中:P_r(z)为激光雷达接收到的高度z处的大气后向散射回波信号能量;P_t为激光雷达的发射能量;C为雷达常数;β_a(z)、β_m(z)分别为气溶胶和大气分子的后向散射系数;α_a(z)、α_m(z)分别为气溶胶和大气分子的消光系数。

气溶胶后向散射系数和消光系数可以分别表示为β_m=n_sγ_m,α_m=n_sσ_m,其中n_s为散射气溶胶粒子浓度,γ_m为单个气溶胶散射粒子对入射激光光束的后向散射截面,σ_m为单个气溶胶散射粒子的消光截面。假设在同一时刻,激光雷达路径上不同高度的气溶胶后向散射截面及消光截面一样,则(1)式可改写为

Pr(z)=PtCz-2[βa(z)+ns(z)γm]exp·-2∫0z[αa(z')+ns(z')σm]dz'。(2)

在低层大气中,气溶胶粒子的散射和吸收远大于气体分子的散射和吸收,因此激光雷达回波信号能量和气溶胶粒子数浓度呈正相关。当激光光束照射到粒子浓度高的区域时,回波信号强度高;反之,回波信号强度弱。

2.2 三光束二维风场反演方法

基于激光雷达回波相关性计算风速的原理如图1所示,大气中气溶胶的不均匀结构在风的带动下^[22],从距离激光雷达1 L₁处的探测点移动到距离激光雷达2 L₂处的探测点,由(2)式可知,不同密度的大气造成激光雷达回波信号强度的起伏不同,所以距离激光雷达1 L₁处和距离激光雷达2 L₂处的回波信号随时间起伏曲线类似,但是有一定的延迟时间t。

延迟时间t可以通过距离激光雷达1 L₁处和距离激光雷达2 L₂处的回波信号时间序列的相关性求得。假设距离激光雷达1 L₁处t₁到t₂时间内的信号为一维矩阵D₁[t₁:t₂],距离激光雷达2 L₂处t₁+Δt到t₂+Δt时间内的信号为一维矩阵D₂[t₁+Δt:t₂+Δt],两个矩阵的相关系数为R₁₂,则R₁₂的曲线如图2所示。R₁₂最大值对应的时间Δt处,即为气溶胶不均匀结构从距离激光雷达1 L₁处移动至距离激光雷达2 L₂处的时间t。距离激光雷达1 L₁处和距离激光雷达2 L₂处两点的距离s可以用两个激光雷达光束夹角、L₁和L₂计算得到。求得t之后,根据公式v=s/t就可以计算出风速。

图 1. 基于激光雷达回波相关性的测风原理示意图。(a)激光雷达探测气溶胶不均匀结构;(b)距离激光雷达1 L1处和距离激光雷达2 L2处的信号曲线

Fig. 1. Schematic of wind measurement based on lidar echo correlation. (a) Measurement of aerosol inhomogeneity by lidar; (b) signal curves of lidar 1 at distance of L1 and lidar 2 at distance of L2

下载图片 查看所有图片

图 2. 相关系数R12与延迟时间Δt的关系曲线

Fig. 2. Correlation coefficient R12 versus delay time Δt

下载图片 查看所有图片

基于此原理,可以用1台激光雷达某个距离处的信号和另1台激光雷达不同距离处的信号进行相关计算并反演二维风场,也可以通过3台呈一定角度的激光雷达相同距离处的信号进行相关计算并反演二维风场。用2台激光雷达成本低,但是由于不同距离处激光雷达信号的信噪比不一样,因此不同距离处激光雷达信号的相关性相比相同距离处激光雷达信号的相关性较差,反演风场的精度较低,且1台激光雷达某个距离和另1台激光雷达不同距离的夹角随着距离的变化而变化,计算较复杂;用3台激光雷达虽然成本较高,但是二维风场的反演精度较高,计算简单,所以采用3台激光雷达来反演二维风场。

三光束二维风场反演方法的原理如图3所示。3台激光雷达在同一平面放置,如图3(a)所示,3台激光雷达发射光束也在同一平面,激光雷达1与激光雷达2的夹角以及激光雷达2与激光雷达3的夹角均为θ。已知激光雷达望远镜的光斑距离为a,则距离3台激光雷达探测L处3点之间的间距s为

s=a+2Lsinθ2。(3)

计算激光雷达1与激光雷达2的信号相关系数R₁₂以及激光雷达2与激光雷达3的信号相关系数R₂₃,然后求出R₁₂最大值对应的时间t₁₂和R₂₃最大值对应的时间t₂₃,进而根据v₁₂=s/t₁₂,v₂₃=s/t₂₃,求出距离激光雷达1 L处到距离激光雷达2 L处方向上的风速v₁₂和距离激光雷达2 L处到距离激光雷达3 L处方向上的风速v₂₃。如图3(b)所示,在v₁₂>0、v₂₃>0情况下,图中虚线l₁垂直于v₂₃方向,虚线l₂垂直于v₁₂方向,v_w为真实风速。假设v_w与l₂的夹角为α,v_w与l₁的夹角为β,则有

|v12|=vw/sin(θ+β)|v23|=vw/sinβα+β=π-θ。(4)

当v₁₂<0、v₂₃<0时,同样满足(4)式。在v₁₂>0、v₂₃<0和v₁₂<0、v₂₃>0情况下,有

|v12|=vw/sin(θ+β)|v23|=vw/sinβα+β=θ。(5)

根据(4)、(5)式可求得

vw=2v12v23tan(θ/2)[tan2(θ/2)+1][(v12-v23)2+(v12+v23)2tan2(θ/2)],(6)α=φ-arcsin2v12tan(θ/2)[tan2(θ/2)+1][(v12-v23)2+(v12+v23)2tan2(θ/2)],(7)β=arcsin2v12tan(θ/2)[tan2(θ/2)+1][(v12-v23)2+(v12+v23)2tan2(θ/2)]。(8)

在v₁₂>0、v₂₃>0和v₁₂<0、v₂₃<0时,φ=π-θ;在v₁₂>0、v₂₃<0和v₁₂<0、v₂₃>0时,φ=θ。

最终的风向和激光雷达2光束方向的夹角γ为

γ=θ2+β(v12>0,v23>0),(9)γ=θ2+π-β(v12<0,v23>0),(10)γ=θ2+π+β(v12<0,v23<0),(11)γ=2π-β+θ2(v12>0,v23<0)。(12)

这样就可以根据三光速二维风场反演方法得到距离激光雷达2光束方向L处的风速和风向数据,按照同样的方法计算不同距离处的风速和风向,就可以得到整个探测路径上的二维风场数据。

图 3. 三光束二维风场反演方法原理。(a)三台激光雷达位置及角度;(b)风速风向反演示意图

Fig. 3. Principle of three-beam two-dimensional wind field inversion method. (a) Position and angle of three lidars; (b) schematic of wind velocity and direction inversion

下载图片 查看所有图片

2.3 相关二维测风激光雷达系统

基于三光束二维风场反演方法,搭建了相关二维测风激光雷达系统,系统所用激光雷达是在课题组微脉冲激光雷达^[23]基础上研制的。雷达采用532 nm Nd:YAG激光器,单脉冲激光能量为10 μJ,重复频率为2.5 kHz,探测器为单光子探测器,信号采集的时间分辨率为1 s。表1为激光雷达的主要参数。

相关二维测风激光雷达系统结构及实物图如图4所示。系统由3台上述激光雷达组成,激光雷达摆放在一个机械调整平台上,两边的激光雷达与中间激光雷达的夹角θ可以调整。激光雷达1与激光雷达2、激光雷达2与激光雷达3的望远镜光斑距离均为22 mm。3台激光雷达的控制线和信号线均连接到1台工控机上,由工控机同时控制和采集。由于3台激光雷达信号采集的距离分辨率均为30 m,时间分辨率均为1 s,可以实现激光雷达2探测路径上30 m距离分辨率、1 s时间分辨率的互相关计算,根据三光束二维风场反演方法,可以实现激光雷达2探测路径上30 m距离分辨率、1 s时间分辨率的风场探测。

图 4. 相关二维测风激光雷达系统。 (a)结构图;(b)实物图

Fig. 4. Correlation-based two-dimensional wind measurement lidar system. (a) Structural diagram; (b) physical map

下载图片 查看所有图片

3 二维风场探测实验与结果分析

为了验证相关二维测风激光雷达系统,对系统和风塔上的风速、风向传感器进行了对比实验。系统实验位置为合肥市科学岛中国科学院安徽光学精密机械研究所1号楼7楼平台,如图5上的红点所示,风塔位于图上的黑点标识处。系统探测方向为水平探测,光束方向为图上3条绿色直线所示。为了避免激光光束直接射在风塔上,激光雷达2的探测路径稍微偏离风塔。风塔高为30 m,每隔5 m高度有一组风速、风向传感器,风速传感器为杯式风速传感器,风向传感器为风向标式传感器, 每10 min采集一组数据。实验点和风塔的直线距离为460 m,风向与正南方向的偏角为34.4°。相关二维测风激光雷达系统所在高度距离地面约为27 m,由于系统水平探测,系统风塔处数据和风塔25 m高度上的风速风向传感器数据最为接近,所以选取风塔25 m高度上的风速风向数据作为对比。

图 5. 相关二维测风激光雷达系统实验位置与风塔位置

Fig. 5. Positions of correlation-based two-dimensional wind measurement lidar system and wind tower

下载图片 查看所有图片

相关二维测风激光雷达系统中激光雷达间的夹角非常重要,对测量结果精度有较大影响。为了确定激光雷达间的夹角,在2018年11月26日19:40—11月27日00:30,相关二维测风激光雷达系统每隔30 min调整激光雷达1与激光雷达2、激光雷达2与激光雷达3之间的夹角,两个夹角相同,依次调整为0.5°,1.0°,1.5°,2.0°,2.5°,3.0°,3.5°,4.0°,4.5°,5.0°,根据反演结果和实际风速、风向传感器数据对比,确定合适夹角。将相关二维测风激光雷达系统中激光雷达夹角调整为合适夹角后,在2018年11月27日4:10—4:55及2018年11月27日8:05—8:50,利用相关二维测风激光雷达系统进行了连续的测风实验。

3.1 气溶胶不均匀结构

图6(a)~(c)分别为系统激光雷达1、激光雷达2、激光雷达3在450 m距离处且在2018年11月26日20:00:00—20:08:20时间段内的信号起伏,信号用光子数表示。图6(d)~(f)分别为激光雷达1、激光雷达2、激光雷达3在相同500 s时间段的信号时空演变图,从图中可以明显看出大气中的气溶胶不均匀结构,且3台雷达探测到的气溶胶不均匀结构变化十分相似,验证了激光雷达信号相关计算的原理。由于3台雷达系统不可能完全相同,所以3台激光雷达的信号也有所区别,激光雷达2的信噪比最高,激光雷达3的信噪比稍低。

图 6. 3台激光雷达在450 m距离处的信号起伏曲线及时空演变图。(a)激光雷达1的信号起伏曲线;(b)激光雷达2的信号起伏曲线;(c)激光雷达3的信号起伏曲线;(d)激光雷达1的信号时空演变图;(e)激光雷达2的信号时空演变图;(f)激光雷达3的信号时空演变图

Fig. 6. Signal curves at distance of 450 m and signal temporal variations for three lidars. (a) Signal curve of lidar 1; (b) signal curve of lidar 2; (c) signal curve of lidar 3; (d) signal temporal variation of lidar 1; (e) signal temporal variation of lidar 2; (f) signal temporal variation of lidar 3

下载图片 查看所有图片

图7(a)~(i)分别为距离激光雷达2在150,300,450,600,750,900,1050,1200,1350 m处的回波信号于2018年11月26日20:00:00—20:08;20时间段内(500 s)的起伏。从图中可以看出:在1350 m处因信噪比低,气溶胶不均匀结构不明显;在其他距离处都能看到明显由于气溶胶不均匀结构造成的信号起伏。不过在150 m和300 m距离处,因为处于几何修正区内,信噪比较低。

图 7. 激光雷达2在不同距离处的信号起伏曲线。(a) 150 m; (b) 300 m; (c) 450 m; (d) 600 m; (e) 750 m; (f) 900 m; (g) 1050 m; (h) 1200 m; (i) 1350 m

Fig. 7. Signal curves of lidar 2 at different distances. (a) 150 m; (b) 300 m; (c) 450 m; (d) 600 m; (e) 750 m; (f) 900 m; (g) 1050 m; (h) 1200 m; (i) 1350 m

下载图片 查看所有图片

3.2 激光雷达夹角及相关计算时间序列长度的确定

为了确定激光雷达间合适的夹角,将相关二维测风激光雷达系统中激光雷达1与激光雷达2、激光雷达2与激光雷达3之间的夹角(两个夹角相同)依次调整为0.5°,1.0°,1.5°,2.0°,2.5°,3.0°,3.5°,4.0°,4.5°,5.0°,在每个角度下测量30 min,并将测量结果与风塔上的风速、风向传感器的3组数据进行对比。从之前气溶胶不均匀结构分析中可以看出,500 s的信号起伏曲线明显表现出不同激光雷达之间的相关性,所以风场计算过程中激光雷达信号矩阵D[t₁:t₂]的所用时间序列长度t₂-t₁可以初步选定为500 s。图8为系统激光雷达在不同夹角情况下的探测结果与风速风向传感器数据的对比。在激光雷达夹角分别为0.5°,1.0°,1.5°,2.0°,2.5°,3.0°,3.5°,4.0°,4.5°,5.0°的情况下,风速平均测量误差分别为88.44%,76.58%,24.06%,12.17%,12.44%,59.55%,79.18%,46.72%,88.31%,68.63%;风向平均测量误差分别为40.84%,56.57%,17.61%,14.03%,14.46%,22.44%,80.61%,88.68%,40.14%,61.81%。从对比结果中可以看出:当系统中激光雷达间的夹角为2.0°时,测量的风速和风向的平均误差最小。因此选定系统中激光雷达的夹角为2.0°进行测风实验。

在相关系数计算中,如果激光雷达信号矩阵D[t₁:t₂]的所用时间序列长度t₂-t₁过大,计算量就会变大,从而影响计算效率;如果t₂-t₁过小,可能无法得到有效的相关系数。图9(a)~(f)分别为激光雷达1和激光雷达2在450 m距离处、t₂-t₁为500,300,200,100,50,20 s时的相关系数曲线,可见在t₂-t₁为500,300,200 s时,根据相关系数曲线可以准确地判断出最大值,但是100 s时曲线已经变差,50 s和20 s时已经无法准确判断相关系数的最大值,所以选择t₂-t₁为200 s。

图 8. 相关二维测风激光雷达系统中激光雷达间不同夹角下的风速及风向探测结果。 (a)风速;(b)风向

Fig. 8. Measurement results of wind speed and direction by correlation-based two-dimensional wind measurement lidar system under different lidar angles. (a) Wind speed; (b) wind direction

下载图片 查看所有图片

图 9. 不同时间序列长度时的相关曲线。(a) 500 s; (b) 300 s; (c) 200 s; (d) 100 s; (e) 50 s; (f) 20 s

Fig. 9. Correlation curves for different lengths of time sequences. (a) 500 s; (b) 300 s; (c) 200 s; (d) 100 s; (e) 50 s; (f) 20 s

下载图片 查看所有图片

3.3 探测结果分析

基于三光束二维风场反演方法,相关二维测风激光雷达系统在2018年11月27日4:10探测的风速和风向数据分别如图10(a)、(b)所示,图中曲线为相关二维测风激光雷达系统探测的风速、风向,方块为风塔上25 m处的风速、风向传感器数据,将其与450 m处的曲线进行对比。与风速、风向传感器相比,系统探测的风速测量误差为12.50%,风向测量误差为0.45%。从图10中可以看出,风速和风向在近距离100 m以内和远距离的相差较大,这可能与近距离路径上的建筑物遮挡有关。而在700 m以后,风速和风向的波动明显变大,这主要是由于远距离信噪比降低,导致反演误差变大。

图 10. 相关二维测风激光雷达系统探测的单个时刻风速和风向曲线。(a)风速;(b)风向

Fig. 10. Curves of wind speed and direction measured by correlation-based two-dimensional wind measurement lidar system at specific time. (a) Wind speed; (b) wind direction

下载图片 查看所有图片

图11为相关二维测风激光雷达系统在2018年11月27日4:10—4:55及2018年11月27日8:05—8:50(共90 min)内探测的450 m距离处的风速风向数据,并将其与风塔25 m高度上的风速风向传感器数据进行对比。从图中可以看出,相关二维测风激光雷达系统的时间和距离不是连续变化的,而是按时间分辨率和距离分辨率的倍数变化的,导致反演的风速、风向结果也不是连续变化的,变化曲线呈阶梯状。系统夜间的探测结果和风速、风向传感器所测数据吻合得较好,风速平均测量误差为9.02%,风向平均测量误差为2.55%。白天激光雷达系统受太阳背景光的影响,信号的信噪比相比夜间低,所以探测结果的误差比夜间的大,风速平均测量误差为27.45%,风向平均测量误差为22.70%,但是总体趋势吻合得较好。

图 11. 相关二维测风激光雷达系统在450 m距离处探测得到的2700 s夜间、白天的风速和风向曲线。 (a)夜间风速;(b)夜间风向;(c)白天风速;(d)白天风向

Fig. 11. Curves of 2700 s wind speed and direction at night and in daytime measured by correlation-based two-dimensional wind lidar measurement system at distance of 450 m. (a) Wind speed at night; (b) wind direction at night; (c) wind speed in daytime; (d) wind direction in daytime

下载图片 查看所有图片

4 结论

提出了基于激光雷达信号相关计算的三光束二维风场反演方法,并基于此方法搭建了相关二维测风激光雷达系统。相关二维测风激光雷达系统所用的激光雷达采用单脉冲能量为10 μJ、重复频率为2.5 kHz的532 nm激光雷达,可以实现距离分辨率为30 m、时间分辨率为1 s的二维风场探测。在相关系数计算中,时间序列长度的选择十分重要,根据实际探测数据计算,最后选择时间序列长度为200 s。相关二维测风激光雷达系统在合肥科学岛进行了水平测风实验,并和探测路径上风塔对应高度处的风速、风向传感器进行了对比,结果表明,相关二维测风激光雷达系统夜间探测的风速、风向数据与风塔上传感器的数据误差分别为9.02%和2.55%,白天的探测数据误差分别为27.45%和22.70%。说明相关二维测风激光雷达系统在白天探测受背景噪声的影响较大,因此需要进行信号去噪的研究;系统目前实现的是二维风场的快速探测,在以后的工作中需进行三维风场探测的方法研究及系统搭建。