-
由于激光测风雷达直接探测得到的径向风矢量为沿扫描光束方向的风场分量。因此,需要选取适当的方法对径向风速进行反演。风场反演算法在国内外被普遍应用,在20世纪60年代速度方位显示方法由LHERMITTE提出,之后由BROWNING等人将其完善[3]。到20世纪70年代,在局部地区均匀风的条件下有人提出了速度体积处理技术,但此方法需要进行体扫,运算量大,分辨率也只有10m。到21世纪,由NEWSOM和BANAKH等人提出并改进了DBS四波束合成计算的方法[4-6],此方法至今仍被广泛运用。而在国内风场反演方法发展也十分迅速,JIANG[7]和QIU[8]等人分别提出了涡度-散度方法和简单伴随函数方法,对风场反演技术的发展提供了支持。
现阶段风场反演算法以速度方位显示和DBS四波束合成计算两种算法居多,两种算法均基于假定风场稳定的情况下对风场进行多波束扫描,再通过反演计算得到高精度的风速风向信息,并计算出垂直气流信息。
-
VAD风场反演算法的基本原理是假定在均匀风和线性风的条件下,利用雷达在某一固定仰角上作360°的全方位圆周扫描,获取不同方位上的径向风速值后, 把一系列的径向风速拟合为方位角的余弦函数,最后求得扫描区域内的平均风场。所拟合的3个参量代表 3维风场信息[9-12]。激光光束的扫描方式如图 1所示。
设风矢量为(u, v, w),在雷达直角坐标系中:u沿着x轴,v沿着y轴,w沿着z轴,圆锥扫描半角γ是激光束指向与z正轴的夹角,扫描方位角θ以x坐标正轴为0°起点。
选取正东方向为x轴正向,正北方向为y轴正向,测得的径向风速为vr,水平风速为vh和垂直气流为w,则有:
$ {v_{\rm{r}}} = {v_{\rm{h}}}\sin \gamma + w\cos \gamma $
(1) 式中,γ为圆锥扫描半角。
根据扫描角θ,水平风速vh可由水平分量u和w合成,则得到公式:
$ {v_{\rm{h}}} = u\cos \theta + v\sin \theta $
(2) 将(2)式带入(1)式,可得:
$ {v_{\rm{r}}} = \left( {u\cos \theta + v\sin \theta } \right)\sin \gamma + w\cos \gamma $
(3) 令$a = \sqrt {{u^2} + {v^2}} $,由(3)式得:
$ {v_{\rm{r}}} = b\cos \left( {a - \theta } \right) + c $
(4) 式中,$b = \sqrt {{u^2} + {v^2}} \sin \gamma , c = w\cos \gamma $。
综上推理,VAD扫描方式中,每一个高度层都对应一个余弦曲线,余弦曲线的振幅、相位、偏距等参量均包含了水平风速、风向以及垂直气流等信息。
-
DBS原理即用四波束合成计算方法进行风场反演。DBS四波束扫描过程中,假设风场矢量保持不变,如图 2所示。设风矢量为(u,v,w),在雷达直角坐标系中:u沿着x轴,v沿着y轴,w沿着z轴,圆锥扫描半角γ是激光束指向与z正轴的夹角,扫描方位角θ以x坐标正轴为0°起点。
径向风速与风向矢量的关系可表示为:
$ {v_{\rm{r}}} = u\cos \theta \sin \gamma + v\sin \theta \sin \gamma + w\cos \gamma $
(5) 设4个方向测得的径向风速分别为vr, 1,vr, 2,vr, 3,vr, 4,对应的方位角分别为θ1,θ2,θ3,θ4,则有:
$ \left\{ \begin{array}{l} {v_{{\rm{r,1}}}} = u\cos {\theta _1}\sin \gamma + v\sin {\theta _1}\sin \gamma + w\cos \gamma \\ {v_{{\rm{r,2}}}} = u\cos {\theta _2}\sin \gamma + v\sin {\theta _2}\sin \gamma + w\cos \gamma \\ {v_{{\rm{r,3}}}} = u\cos {\theta _3}\sin \gamma + v\sin {\theta _3}\sin \gamma + w\cos \gamma \\ {v_{{\rm{r,4}}}} = u\cos {\theta _4}\sin \gamma + v\sin {\theta _4}\sin \gamma + w\cos \gamma \end{array} \right. $
(6) 对于激光雷达,4个扫描方位角与N,E,S,W这4个方向对准,即θ1,θ2,θ3,θ4,分别为0°,90°,180°,270°。
则(6)式可简化为:
$ \left\{ \begin{array}{l} {v_{{\rm{r,E}}}} = {v_{{\rm{r,1}}}} = u\sin \gamma + w\cos \gamma \\ {v_{{\rm{r,N}}}} = {v_{{\rm{r,2}}}} = v\sin \gamma + w\cos \gamma \\ {v_{{\rm{r,W}}}} = {v_{{\rm{r,3}}}} = - u\sin \gamma + w\cos \gamma \\ {v_{{\rm{r,S}}}} = {v_{{\rm{r,4}}}} = - v\sin \gamma + w\cos \gamma \end{array} \right. $
(7) 进而推导出:
$ \left\{ \begin{array}{l} u = \frac{{{v_{{\rm{r,E}}}} - {v_{{\rm{r,W}}}}}}{{2\sin \gamma }}\\ v = \frac{{{v_{{\rm{r,N}}}} - {v_{{\rm{r,S}}}}}}{{2\sin \gamma }}\\ w = \frac{{{v_{{\rm{r,E}}}} + {v_{{\rm{r,N}}}} + {v_{{\rm{r,W}}}} + {v_{{\rm{r,S}}}}}}{{4\cos \gamma }} \end{array} \right. $
(8) -
可以看出,VAD和DBS原理均在雷达水平放置的情况下来进行风速测量和风场反演。而当雷达随其放置载体(例如车载雷达)经过各种地形时,会遇到上坡、下坡、转弯等各类状况,使得测风雷达不能保证时刻处于水平状态进行测量,导致风场反演结果与真实情况不符。如图 3所示,其中箭头为假定水平风和垂直气流。
-
假定雷达在车体上放置时与车体底盘平行,则由车体倾斜而引起的雷达相对于水平地面的倾斜方向与车体一致。雷达倾斜示意图如图 4所示。
当车体沿轴向方向倾斜,用β表示,β>0°表示朝前倾斜,β < 0°表示朝后倾斜;车体垂直于轴向方向倾斜(横向倾斜),用α表示,α>0°表示朝右边倾斜,α < 0°表示朝左边倾斜。
由于雷达倾斜的角度与车体倾斜角度相一致,因此雷达倾斜情况同样为:横向倾斜α,轴向倾斜为β,由此相关计算如下。
如图 4所示,设立方体ABCD-A1B1C1D为雷达,AD方向为雷达轴向,ABCD所处平面为水平地面,雷达D处为车头一侧,车体经横向倾斜α,轴向倾斜β以后,使雷达到达GHFC位置,于是:∠ABE=∠DCF=α,∠EFH=∠BCG=β。
设立方体边长为1,由于BRFC为平行四边形,所以EF⊥AA1,又由于平面BRFC绕着轴CF旋转β角到达位置GHFC,易知GHFC为平行四边形。
延长HF与AD的延长线交于点P,连接C和P,过点H做HK⊥CP,连接A和K,则∠AKH为平面GHFC与平面ABCD的二面角。
因为正方体ABCD-A1B1C1D棱长为1,所以有AE=ABtanα=tanα,经计算,可求出:
$ \left\{ \begin{array}{l} AH = \tan \alpha + \tan \beta \\ AK = \frac{{\tan \alpha + \tan \beta }}{{\sqrt {{{\tan }^2}\alpha + {{\tan }^2}\beta } }} \end{array} \right. $
(9) 因此,车体倾斜的角度为:
$ \angle AKH = \arctan \sqrt {{{\tan }^2}\alpha + {{\tan }^2}\beta } $
(10) 同理可得,倾斜方向与车体轴向的夹角为:
$ \angle PAK = \arccos \frac{{\tan \beta }}{{\sqrt {{{\tan }^2}\alpha + {{\tan }^2}\beta } }} $
(11) 由于雷达倾斜具有一定的方向性,通过分析可知:
(1) 当α≥0°, β>0°时:
$ \left\{ \begin{array}{l} \angle AKH = \arctan \sqrt {{{\tan }^2}\alpha + {{\tan }^2}\beta } \\ \angle PAK = \arccos \frac{{\tan \beta }}{{\sqrt {{{\tan }^2}\alpha + {{\tan }^2}\beta } }} \end{array} \right. $
(12) (2) 当α≤0°, β>0°时:
$ \left\{ \begin{array}{l} \angle AKH = \arctan \sqrt {{{\tan }^2}\alpha + {{\tan }^2}\beta } \\ \angle PAK = 2{\rm{ \mathsf{ π} }} - \arccos \frac{{\tan \beta }}{{\sqrt {{{\tan }^2}\alpha + {{\tan }^2}\beta } }} \end{array} \right. $
(13) (3) 当α≤0°, β<0°时:
$ \left\{ \begin{array}{l} \angle AKH = \arctan \sqrt {{{\tan }^2}\alpha + {{\tan }^2}\beta } \\ \angle PAK = 2{\rm{ \mathsf{ π} }} - \arccos \frac{{\tan \beta }}{{\sqrt {{{\tan }^2}\alpha + {{\tan }^2}\beta } }} \end{array} \right. $
(14) (4) 当α≥0°, β<0°时:
$ \left\{ \begin{array}{l} \angle AKH = \arctan \sqrt {{{\tan }^2}\alpha + {{\tan }^2}\beta } \\ \angle PAK = \arccos \frac{{\tan \beta }}{{\sqrt {{{\tan }^2}\alpha + {{\tan }^2}\beta } }} \end{array} \right. $
(15) (5) 当α=0°, β=0°时:
$ \left\{ \begin{array}{l} \angle AKH = 0\\ \angle PAK = 0 \end{array} \right. $
(16) -
当雷达与正北方向成γ角,同时自身倾斜角度为θ时,设原坐标系中的风矢量为(u, v, w),倾斜后新坐标系中的风矢量为(vx, vy, vz),${v_{\rm{h}}} = \sqrt {v_x^2 + v_y^2} $为倾斜以后测得的水平风速值的大小,ξ为倾斜后测得的风向的值,vz为倾斜后测得的垂直气流,如图 5所示。A为倾斜前的水平面,B为倾斜以后的面,OP为倾斜方向,O′P为OP在A平面内的投影。
vx-O-vy和u-O′-v都为地理坐标系(原点与雷达重合,vy轴和v轴指向正北方向,vx轴和u轴指向正东方向,vz轴和w轴指向垂直向下方向)。
则在vx-O-vy平面内,以OP为x轴,OQ为y轴建立新的直角坐标系,将风矢量为[vx, vy, vz]重新分解,如图 6所示。
(1) 当ξ≥γ时,风矢量与倾斜方向的夹角为π-(ξ-γ),水平风速在OP轴和OQ轴上的分量分别为:
$ \left\{ \begin{array}{l} OP = {v_{\rm{h}}}\cos \left( {{\rm{ \mathsf{ π} }} + \gamma - \xi } \right) = - {v_{\rm{h}}}\cos \left( {\gamma - \xi } \right)\\ OQ = {v_{\rm{h}}}\sin \left( {{\rm{ \mathsf{ π} }} + \gamma - \xi } \right) = {v_{\rm{h}}}\sin \left( {\xi - \gamma } \right) \end{array} \right. $
(17) (2) 同理,当ξ < γ时,水平风速在OP轴和OQ轴上的分量分别为:
$ \left\{ \begin{array}{l} OP = {v_{\rm{h}}}\cos \left( {{\rm{ \mathsf{ π} }} + \gamma - \xi } \right) = - {v_{\rm{h}}}\cos \left( {\gamma - \xi } \right)\\ OQ = {v_{\rm{h}}}\sin \left( {{\rm{ \mathsf{ π} }} + \gamma - \xi } \right) = {v_{\rm{h}}}\sin \left( {\xi - \gamma } \right) \end{array} \right. $
(18) 可见,无论ξ与γ的关系如何,水平风速在OP轴和OQ轴上的分量计算结果相同,即:
$ \left\{ \begin{array}{l} OP = - {v_{\rm{h}}}\cos \left( {\xi - \gamma } \right)\\ OQ = {v_{\rm{h}}}\sin \left( {\gamma - \xi } \right) \end{array} \right. $
(19) 同理,在u-O′-v平面内,如图 7所示。以O′P为x轴、O′Q′为y轴建立新的直角坐标系,则由图 4可知,∠OPO′=θ,O′P与v轴正向的夹角为γ。
经计算可得:
$ \left\{ \begin{array}{l} u = O'P\sin \gamma - O'Q'\cos \gamma \\ v = O'P\cos \gamma + O'Q'\sin \gamma \end{array} \right. $
(20) 将图 7中坐标系POQ-vz与坐标系PO′Q′-w相互叠加,使原点重合,得两坐标系的关系如图 8所示。其中,M代表的是原来的O和O′,MS代表O′P,MK代表OP,MD代表O′Q′和OQ,可得到:
$ \left\{ \begin{array}{l} MS = MK\cos \theta - {v_z}\sin \theta \\ MD = MD\\ w = MK\sin \theta + {v_z}\cos \theta \end{array} \right. $
(21) 进而由(20)式和(21)式可推导出:
$ \left\{ \begin{array}{l} u = - {v_{\rm{h}}}\cos \left( {\xi - \gamma } \right)\cos \theta \sin \gamma - {v_z}\sin \theta \cdot \\ \;\;\;\;\sin \gamma - {v_{\rm{h}}}\sin \left( {\xi - \gamma } \right)\cos \gamma \\ v = - {v_{\rm{h}}}\cos \left( {\xi - \gamma } \right)\cos \theta \cos \gamma - {v_z}\sin \theta \cdot \\ \;\;\;\;\cos \gamma + {v_{\rm{h}}}\sin \left( {\xi - \gamma } \right)\sin \gamma \\ w = - {v_{\rm{h}}}\cos \left( {\xi - \gamma } \right)\sin \theta + {v_z}\cos \theta \end{array} \right. $
(22) 对于风向计算,设倾斜前风向为ξx,k为中间变量,则相关计算如下:
$ k = \left\{ \begin{array}{l} \alpha \tan \left( {v{u^{ - 1}}} \right),\left( {u \ge 0,v > 0} \right)\\ \alpha \tan \left( {v{u^{ - 1}}} \right) + {\rm{ \mathsf{ π} }},\left( {u \le 0,v > 0} \right)\\ \alpha \tan \left( {v{u^{ - 1}}} \right) + {\rm{ \mathsf{ π} }},\left( {u \le 0,v < 0} \right)\\ \alpha \tan \left( {v{u^{ - 1}}} \right) + 2{\rm{ \mathsf{ π} }},\left( {u \ge 0,v < 0} \right)\\ 0,\left( {u = 0,v = 0} \right) \end{array} \right. $
(23) 从而推导出:
$ {\xi _x} = \left\{ \begin{array}{l} 3{\rm{ \mathsf{ π} }}/2 - k,\left( {0 \le k \le {\rm{ \mathsf{ π} }}/2} \right)\\ 3{\rm{ \mathsf{ π} }}/2 - k,\left( {{\rm{ \mathsf{ π} }}/2 \le k \le {\rm{ \mathsf{ π} }}} \right)\\ 3{\rm{ \mathsf{ π} }}/2 - k,\left( {{\rm{ \mathsf{ π} }} \le k \le 3{\rm{ \mathsf{ π} /2}}} \right)\\ 7{\rm{ \mathsf{ π} }}/2 - k,\left( {3{\rm{ \mathsf{ π} }}/2 \le k \le 2{\rm{ \mathsf{ π} }}} \right) \end{array} \right. $
(24) 最后通过ξx×180/π将风向由弧度换算为角度便得出风向信息。
在雷达随车体横向倾斜α、而轴向倾斜β后,雷达沿着与车头方向夹角为∠PAK的方向倾斜∠AKH。如果已知雷达所在车体轴向与正北方向的夹角为δ的话,则倾斜方向与正北方向的夹角为γ=δ+∠PAK,(vh, ξ, vz)分别为倾斜后测得的水平风速、风向、垂直气流,则可推出真实的水平风速为:
$ {V_x} = \sqrt {{u^2} + {v^2}} $
(25) 修正后垂直气流(即倾斜以前的垂直气流)为:
$ \begin{array}{*{20}{c}} {w = {v_z}\cos \angle AKH - {v_{\rm{h}}}\cos \left( {\xi - \gamma } \right) \cdot }\\ {\sin \angle AKH = {v_z}\cos \angle AKH - }\\ {{v_{\rm{h}}}\cos \left( {\xi - \delta - \angle PAK} \right)\sin \angle AKH} \end{array} $
(26) 在雷达倾斜的情况下(倾斜角度在10°以内),对于垂直气流测量的影响较大,而对水平风速风向的影响并不大[13](将在下一节的实验例证中运用数据对比曲线来证明),面对此类问题,作者提出运用空间坐标变换原理,将雷达的倾斜状况进行整体修正,计算量小且修正效果明显。
机动型激光测风雷达倾斜风场修正算法研究
Correction method of tilt wind field of mobile wind lidar
-
摘要: 为了解决激光测风雷达处于倾斜放置的情况下测出的风场信息与实际风场之间存在的误差,分析了激光雷达的测量原理、放置情况以及存在的问题,采用空间坐标变换的方法对雷达倾斜状态进行理论计算并给出了修正计算公式。对比了水平放置和横向倾斜7.2°放置的激光雷达的实测数据,发现修正前由于倾斜放置导致雷达所测垂直气流数据间产生的误差较大,加入修正算法进行仿真后,数据曲线吻合良好,修正效果明显。结果表明,空间坐标变换的修正算法可使激光雷达在处于非特定位置时仍能正确地反映实际风场信息,从而验证了该修正算法的正确性和可行性。Abstract: In order to correct the error between wind field information measured by a tilt lidar and the actual wind field, after analyzing the measurement principle, placement and exist problems of laser radar, theoretical calculation was carried out based on space coordinate transformation method and correction formula was put forward. After comparing the measured data of laser radar located at horizontal position and lateral tilt 7.2°, it was found that the error of vertical airflow data measured by tilt lidar was larger. After modification, data curve is in good agreement with the modified algorithm. The results indicate that could get the real wind field information can be obtained correctly based on space coordinate transformation algorithm for tilt lidar. Correctness and feasibility of the algorithm are proved.
-
Key words:
- laser technique /
- lidar /
- correcting algorithm /
- coordinate transformation /
- wind field /
- tilt
-
[1] LAI D, CHEN Y, ZHOU D F, et al. Beam scanning of lidar and the simulation of the improved VAD inversion methods[J]. Laser Technology, 2008, 32(6): 584-586(in Chinese). [2] LHERMITTE R M, ATLAS D. Precipitation motion by pulse Doppler[C]//American Meteorological Society. 9th Weather Radar Confe-rence.New York, USA: American Meteor Society, 1961: 218-223. [3] BROWNING K A, WEXLER R. The determination of kinematic properties of a wind field using Doppler radar[ J]. Journal of Applied Meteorology, 1968, 7(1):105-113. doi: 10.1175/1520-0450(1968)007<0105:TDOKPO>2.0.CO;2 [4] LI L, XIE Y F, DONG G Y, et al. Wind field inversion technique for scanning wind lidar[J]. Chinese Optics, 2013, 6(2):251-258(in Chinese). [5] NEWSOM R K, LIGON D, CALHOON R, et al. Retrieval of micro scale wind and temperature fields from single-and dual-Doppler lidar data[J]. Journal of Applied Meteorology and Climatology, 2005, 44(9): 1324-1344. doi: 10.1175/JAM2280.1 [6] BANAKH V A, WERNER C. Retrieval of the turbulence parameters from spatial statistics of Doppler estimates[J]. Proceedings of the SPIE, 2004, 5575: 55-66. doi: 10.1117/12.566458 [7] JIANG H Y, GE R Sh. A new retrieval technique for single Doppler radar[J]. Quarterly Journal of Applied Meteorology, 1997, 8(2):219-223. [8] QIU C J, XU Q. A simple adjoint method of wind analysis for single-Doppler data[J]. Journal of Atmospheric and Oceanic Technology, 1992, 9(5): 588-598. doi: 10.1175/1520-0426(1992)009<0588:ASAMOW>2.0.CO;2 [9] WAN R, TANG D Z, ZHANG P, et al. VAD eiementary analysis of nonlinear wind field[J]. Scientia Meteorologica Sinica, 2003, 23(3):314-324(in Chinese). [10] BAI J, ZHOU X B, WANG L K, et al. Error analysis of doppler radar expanded VAD technique[J]. Acta Scientiarum Naturalium Universitatis Pekinensis, 2001, 37(1):48-54(in Chinese). [11] REN P, WU Y L, ZHOU D F, et al. Numerical optimization in VAD invertion technique for wind lidar[J]. Laser Technology, 2009, 33(6):664-669(in Chinese). [12] XUE S Q, YANG H P, GAO Y C. Comparison between VAD retrieved wind and radiosonde data[J].Scientia Meteorologica Sinica, 2010, 30(6):856-861(in Chinese). [13] XU M Q, WANG B, XU H Z. 3D-radar altimetry definition analysis and error correction[J].Modern Defence Technology, 2013, 41(2):123-127(in Chinese).