CN116500648A - 一种地基激光雷达目标区风廓线反演方法 - Google Patents

Abstract

本发明属于反演方法，具体涉及一种地基激光雷达目标区风廓线反演方法。它包括：步骤1：根据雷达与目标位置在地理坐标系下的坐标关系，解算雷达波束需要进行扫描的俯仰、方位范围；步骤2：雷达按照步骤1中计算得到的方位俯仰角范围对目标探测区域进行三维体积扫描，得到各激光波束各距离单元的的径向风速数据；步骤3：根据扫描探测得到的数据进行反演，形成目标区域的风廓线数据并输出。本发明的有益效果是：有效提高了激光雷达的风廓线探测性能，并能够实现远程风廓线探测。

Description

一种地基激光雷达目标区风廓线反演方法

技术领域

本发明属于反演方法，具体涉及一种地基激光雷达目标区风廓线反演方法。

背景技术

众所周知，风廓线的反演是激光测风雷达最重要的应用之一。常用的风廓线的反演方法主要包括以下几种：

(1)距离方位显示法(VAD法)：该方法的核心是假定实际风场在同一高度平面上呈线性变化。在这样的假设下，雷达观测到的径向速度V_h在一个距离圈上(同一仰角和径向距离，但方位不同的数据)是方位角θ的三角函数：

V_h＝ucosθsinγ+v sinθsinγ+w cosγ

通过选取某一距离圈上的数据进行计算，即可得到水平风速，再通过计算不同高度距离圈上的数据，即可得到雷达上方的风廓线数据。

(2)多普勒光束摆动法(Doppler beaming swinging，DBS)：激光雷达围绕扫描中心轴进行东、南、西、北四波束扫描，每个波束的方位间隔为90，相邻波束两两正交，依据如下公式：

可得到雷达上方的风廓线数据，此方法因为运算量小、计算简单等优点，一直被广泛应用。

其后出现了VAP、VVP、NVAD方法等多种风场反演方法，但是这些方法都只是针对雷达上方的风廓线进行探测，如图1)所示，无法对几十公里以内的目标区域进行风廓线探测。

然而，在自然风场资源的测绘、精准目标区域风场探测、炮弹落区风场保障、航天系统降落风场保障、空投空降等人员不可达领域，无法将雷达置于目标空域下方的地面处，现有风廓线反演策略将无法发挥作用。因此，对远程风廓线探测有着迫切的需求。

跟本项目相关的坐标系定义如下：

地理坐标系，指WGS-84世界大地坐标系。其几何意义是，坐标原点为地球质心，其地心空间直角坐标系的Z轴指向地球极(CTP)方向，X轴指向零子午面和CTP赤道的交点，Y轴与Z轴、X轴垂直构成右手坐标系，如图9所示。

雷达坐标系，指以雷达位置为原点，X轴指向东边，Y轴指向北边，Z轴指向天顶，，如图10所示。

跟本发明相关的名词术语的定义如下：

径向风速：指沿着激光波束方向的视向风速。

发明内容

本发明的目的是针对现有技术的缺陷，提供一种目标区激光雷达风廓线反演方法。

为了达到上述目的，本发明提供了一种基于地基激光雷达的目标区风廓线反演算法，本发明所采用的具体技术方案如下：一种地基激光雷达目标区风廓线反演方法，包括：

步骤1：根据雷达与目标位置在地理坐标系下的坐标关系，解算雷达波束需要进行扫描的俯仰、方位范围；

步骤2：雷达按照步骤1中计算得到的方位俯仰角范围对目标探测区域进行三维体积扫描，得到各激光波束各距离单元的的径向风速数据；

步骤3：根据扫描探测得到的数据进行反演，形成目标区域的风廓线数据并输出。

如上所述的一种目标区激光雷达风廓线反演方法，其中，所述的步骤1包括

步骤1.1，以地理坐标系为基准，雷达1所在位置坐标为O(Lg1,La1,H1)，探测目标位置2处的坐标为P(Lg2,La2,H2)，目标区域为探测空域3，其中Lgi为经度，Lai为纬度，Hi为海拔高度(其中，i＝1,2)，空域3的探测区域近似为一个圆柱体S，圆柱体S的中心轴在地面的坐标为(Lg2,La2,H2)，圆柱体的半径为R，高度分辨率为h₀，

步骤1.2，根据目标区域风廓线探测高度分辨率将目标区域划分成相应的网络格点，网格点采用立方体划分法，以目标位置高度为中心的立方体区域都计入该高度层，以任一探测海拔高度h为例，高度层h处的探测区域为圆柱体S上一个高度为h₀的圆柱形子区域，记为S_h，S_h下底面的海拔高度为S_h上底面的海拔高度为/>

步骤1.3，对于圆柱体S_h来说，高度h处的圆柱截面记为g，雷达1与g的圆心连线记为m，利用空间立体几何的理论求解平面g内与直线m的法向量平行的直径l₀，

将该直径l₀四等分，取l₀的5个等分点，分别求出5个等分点在地理坐标系下的坐标O₁(x₁,y₁,z₁)，O₂(x_2,y₂,z₂)，O₃(x₃,y₃,z₃)，O₄(x₄,y₄,z₄)，O₅(x₅,y₅,z₅)，分别连接点O与O₁、O₂、O₃、O₄、O₅，则向量(其中,i＝1,2,…,5)，激光波束分别沿着向量/>(其中，i＝1,2,…,5)的方向进行扫描，

步骤1.4，通过坐标变换，将激光波束的坐标从地理坐标系变换到雷达坐标系，在雷达坐标系下的扫描方位角和俯仰角分别为θ_i和γ_i(i＝1,2,…,5)，

测量时，在上位机软件输入雷达坐标、目标位置坐标及目标空域探测高度、高度分辨率等参数信息，上位机软件自动计算目标区域所有高度层的扫描方位角和俯仰角，以二维数组Coef[n][Attitude]建立所有高度层的扫描信息数据库，其中，n代表高度层，Attitude代表雷达扫描波束的方位角俯仰角信息，该数组在上位机软件运行时写入初始化文件。

如上所述的一种目标区激光雷达风廓线反演方法，其中，所述的步骤2包括

步骤2.1，测量前雷达要先对北，控制雷达按照步骤1计算得到的扫描方位俯仰信息对目标探测区域进行三维体积扫描，扫描时按照逐层扫描的方式进行，并不断重复该过程，直至全部测量结束，

步骤2.2，获取每个激光波束的回波时域数据，回波时域数据按照距离大小将其划分成不同的距离单元，距离单元的长度为m₀，其中m₀应不大于数字接收机通过Hilbert变换将雷达信号从中频变换到基频，通过基于变分模态分解的联合时频分析方法将信号进行频谱变换，多个脉冲回波进行非相干累积得到原始功率谱数据，存储在二维数组R[n][gate_num]中，n代表波束序号，gate_num代表距离单元的序号，

步骤2.3，对每个距离单元的功率谱数据R[n][gate_num]先进行平滑去噪等基本处理，然后将回波信号建模为一个随机过程，通过计算随机过程的期望、方差等统计特征量，以“期望+权重系数×方差”作为自适应阈值检测门限，阈值之上的数据认为是信号包络，检测出的每个距离单元的信号包络，利用一阶矩得到包络质心，解算信号波峰位置和强度，进一步计算出多普勒频率及径向风速数据Vri。

如上所述的一种目标区激光雷达风廓线反演方法，其中，所述的步骤3包括

步骤3.1，综合考虑计算精度和计算规模，网格数量为目标区域探测的高度层总数N，网格的高度不超过目标区域的高度分辨率h₀，

对于特定高度层来说，将雷达扫描探测得到的各个距离单元的径向风速数据划分成两类，一类是落入上述网络格点范围内的点，即径向风所在距离单元在地理坐标系下的坐标落入该高度层所在网格点区域，记为C_rI，另一类是径向风所在距离单元在地理坐标系下的坐标未落入上述网络格点范围内的数据，记为C_rO，其中，r＝1，2，...，N，表示目标区域第r个高度层，

步骤3.2，对于目标区域每一个探测高度层来说，C_rI内的距离单元数据不能少于三个；若不满足该要求，则对C_rO内同一激光波束的各距离单元内的径向风速数据进行递归函数曲线拟合，拟合后的曲线记为f，

步骤3.3，依据上述步骤中的拟合函数f将数据插值到所需要的部分格点上，插值后的C_rI变为C‘_rI，实心圆点表示落入网格内的数据，空心圆点表示未落入网格点的数据，三角形表示插值的数据，

步骤3.4，结合各波束的坐标信息，对落入同一高度层风矢量数据集合C‘_rI，依据下述方程组求解，得到该高度层上的三维风矢量(u,v,w)，求解方程组采用共轭梯度法和有约束多元变量线性优化问题求最优解，可以避免出现采集参数的微小变化引起拟合结果的剧烈变化，

其中，V_i,h(i＝1,2,…,n)表示落入高度层h的第i个风速数据，θ_i表示风速数据V_i,h在雷达坐标系下对应的方位角，β_i表示风速数据V_i,h在雷达自身坐标系下对应的俯仰角。

为了实现地基激光雷达目标区的风廓线探测，目标空域3的探测范围按照圆柱体进行高度层划分，根据雷达1与目标位置2在地理坐标系下的经度、纬度及海拔关系，计算雷达进行扫描探测的波束扫描的方位俯仰信息，对目标空域3进行三维体积扫描，采集每个激光波束的回波时域数据，并按照距离大小将其划分成不同的距离单元，并计算各个距离单元的径向风速数据，然后判断径向风速数据落入目标高度区域的数据个数，若数据个数小于三个，则对该波束的各距离单元数据进行拟合插值，插值到该目标高度区域，对插值后所有落入该高度区域的数据采用共轭梯度法和有约束多元变量线性优化问题求最优解，即可得到该高度区域的风廓线数据。

本发明相对于现有技术具有以下有益效果：本发明所述的地基激光雷达目标区风廓线反演方法，根据雷达和目标位置在地理坐标系下的关系，预先计算雷达进行扫描探测的方位俯仰范围，通过对目标空域按照预先设定的方式进行三维体积扫描，对于每个激光波束的回波数据，按照距离划分成不同距离单元进行处理，采用Hilbert变换将雷达信号从中频变换到基频、基于变分模态分解的联合时频分析进行径向风速的计算，再通过网络格点的方式，对落入相应格点的数据通过拟合插值、有约束多元线性最优化等运算，获取目标空域的风廓线数据，有效提高了激光雷达的风廓线探测性能，并能够实现远程风廓线探测，如图1右侧图所示。

附图说明

图1是传统风廓线探测与本发明实施例中风廓线探测对比示意图；

图2是本发明实施例中风廓线探测示意图；

图3是本发明实施例中三维体积扫描的示意图；

图4是本发明实施例中目标区域网格点划分截面图；

图5是本发明实施例中目标区域网格点划分截面图；

图6是本发明实施例中目标区域网格划分点立体图；

图7是本发明实施例中径向速度拟合的示意图；

图8是本发明实施例中径向速度插值处理的示意图；

图9是地理坐标系示意图；

图10是雷达坐标系示意图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步阐述和说明。对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

一种地基激光雷达目标区风廓线反演方法，包括：

步骤1，根据雷达与目标位置在地理坐标系下的坐标关系，解算雷达波束需要进行扫描的俯仰、方位范围；

步骤1.1，如图2所示，以地理坐标系为基准，雷达1所在位置的坐标为O(Lg1,La1,H1)，探测目标位置2处的坐标为P(Lg2,La2,H2)，目标区域为探测空域3，其中Lgi为经度，Lai为纬度，Hi为海拔高度(其中，i＝1,2)。空域3的探测区域近似为一个圆柱体S，圆柱体S的中心轴在地面的坐标为(Lg2,La2,H2)，圆柱体的半径为R，高度分辨率为h₀。

步骤1.2，根据目标区域风廓线探测高度分辨率将目标区域划分成相应的网络格点，网格点采用立方体划分法，以目标位置高度为中心的立方体区域都计入该高度层，以任一探测海拔高度h为例，高度层h处的探测区域为圆柱体S上一个高度为h₀的圆柱形子区域，记为S_h，S_h下底面的海拔高度为S_h上底面的海拔高度为/>

步骤1.3，对于圆柱体S_h来说，高度h处的圆柱截面记为g，雷达1与g的圆心连线记为m，利用空间立体几何的理论求解平面g内与直线m的法向量平行的直径l₀。

将该直径l₀四等分，取l₀的5个等分点，分别求出5个等分点在地理坐标系下的坐标O₁(x₁,y₁,z₁)，O₂(x₂,y₂,z₂)，O₃(x₃,y₃,z₃)，O₄(x₄,y₄,z₄)，O₅(x₅,y₅,z₅)，分别连接点O与O₁、O₂、O₃、O₄、O₅，则向量(其中,i＝1,2,…,5)。激光波束分别沿着向量/>(其中，i＝1,2,…,5)的方向进行扫描。

步骤1.4，通过坐标变换，将激光波束的坐标从地理坐标系变换到雷达坐标系，在雷达坐标系下的扫描方位角和俯仰角分别为θ_i和γ_i(i＝1,2,…,5)。

测量时，在上位机软件输入雷达坐标、目标位置坐标及目标空域探测高度、高度分辨率等参数信息，上位机软件自动计算目标区域所有高度层的扫描方位角和俯仰角，以二维数组Coef[n][Attitude]建立所有高度层的扫描信息数据库，其中，n代表高度层，Attitude代表雷达扫描波束的方位角俯仰角信息，该数组在上位机软件运行时写入初始化文件。

步骤2，雷达按照步骤1中计算得到的方位俯仰角范围对目标探测区域进行三维体积扫描，得到各激光波束各距离单元的的径向风速数据。

步骤2.1，测量前雷达要先对北，如图3所示，控制雷达按照步骤1计算得到的扫描方位俯仰信息对目标探测区域进行三维体积扫描，扫描时按照逐层扫描的方式进行，并不断重复该过程，直至全部测量结束。

步骤2.2，按照步骤1中计算得到的扫描参数保证了在目标区域的任意一个高度层都有扫描数据，获取每个激光波束的回波时域数据，并按照距离大小将其划分成不同的距离单元，距离单元的长度为m₀，其中m₀应不大于如图4所示。

雷达数字接收机接收通过Hilbert变换将雷达信号从中频变换到基频，通过基于变分模态分解的联合时频分析方法将信号进行频谱变换，多个脉冲回波进行非相干累积得到原始功率谱数据，存储在二维数组R[n][gate_num]中，n代表波束序号，gate_num代表距离单元的序号。

步骤2.3，每个距离单元的功率谱数据R[n][gate_num]可以建模为一个服从正太分布的随机过程，通过计算该随机过程的期望、方差等统计特征量，以“期望+权重系数×方差”作为自适应阈值检测门限，阈值之上的数据认为是信号包络。

检测出的每个距离单元的信号包络，利用一阶矩得到包络质心，解算信号波峰位置和强度，进一步计算出多普勒频率及径向风速数据V_ri；

步骤3，根据扫描探测得到的数据形成目标区域的风廓线数据并输出

步骤3.1，图5所示为目标区域网格点截面图，综合考虑计算精度和计算规模，网格数量为目标区域探测的高度层总数N，网格的高度不超过目标区域的高度分辨率h₀。

对于特定高度层来说，将雷达扫描探测得到的各个距离单元的径向风速数据划分成两类，一类是落入上述网络格点范围内的点，即径向风所在距离单元在地理坐标系下的坐标落入该高度层所在网格点区域，记为C_rI，另一类是径向风所在距离单元在地理坐标系下的坐标未落入上述网络格点范围内的数据，记为C_rO，其中，r＝1，2，...，N，表示目标区域第r个高度层。图6所示为目标区域网格点立体图。

本步骤的综合考虑是指综合考虑计算精度和计算规模，计算规模主要和步骤2中扫描径向光速的数量相关，反应到具体计算，就是步骤3.4中方程的数量。因不是本申请关注的重点，且本领域技术人员根据现有技术可以进行相关计算的优化，因此本申请对于该部分技术问题不细致讨论。

步骤3.2，对于目标区域每一个探测高度层来说，C_rI内的距离单元数据不能少于三个；若不满足该要求，则对C_rO内同一激光波束的各距离单元内的径向风速数据进行递归函数曲线拟合，拟合后的曲线记为f，如图7所示。

步骤3.3，依据上述步骤中的拟合函数f将数据插值到所需要的部分格点上，插值后的C_rI变为C‘_rI，如图8所示，实心圆点表示落入网格内的数据，空心圆点表示未落入网格点的数据，三角形表示插值的数据；

步骤3.4，结合各波束的坐标信息，对落入同一高度层风矢量数据集合C‘_rI，依据下述方程组求解，得到该高度层上的三维风矢量(u,v,w)。求解方程组采用共轭梯度法和有约束多元变量线性优化问题求最优解，可以避免出现采集参数的微小变化引起拟合结果的剧烈变化。

其中，V_i,h(i＝1,2,…,n)表示落入高度层h的第i个风速数据，θ_i表示风速数据V_i,h在雷达自身坐标系下对应的方位角，β_i表示风速数据V_i,h在雷达自身坐标系下对应的俯仰角。

改变高度参数，重复执行步骤3.1～3.4，直到所有高度层均完成计算，输出该区域各高度层的风廓线数据。

本申请的几个步骤的作用可以概括如下：

步骤1进行坐标定义和基本的坐标变换；

步骤2是如何利用回波信号反演径向速度；

步骤3是如何利用反演的径向速度反演目标区风廓线。

本领域技术人员应理解的是，在本发明的揭露中，“四等分”、“5个”等指示的数量关系是基于附图所示的位置关系，其仅是为了便于描述本发明和简化描述，而不是指示或暗示本算法必须具有特定数量的操作，因此，上述术语不能理解为对本发明的限制。

本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书界定。

本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书界定。

Claims (4)

1.一种地基激光雷达目标区风廓线反演方法，其特征在于，包括：

步骤1：根据雷达与目标位置在地理坐标系下的坐标关系，解算雷达波束需要进行扫描的俯仰、方位范围；

步骤2：雷达按照步骤1中计算得到的方位俯仰角范围对目标探测区域进行三维体积扫描，得到各激光波束各距离单元的的径向风速数据；

步骤3：根据扫描探测得到的数据进行反演，形成目标区域的风廓线数据并输出。

2.如权利要求1所述的一种目标区激光雷达风廓线反演方法，其特征在于：所述的步骤1包括

步骤1.1，以地理坐标系为基准，雷达1所在位置坐标为O(Lg1,La1,H1)，探测目标位置2处的坐标为P(Lg2,La2,H2)，目标区域为探测空域3，其中Lgi为经度，Lai为纬度，Hi为海拔高度(其中，i＝1,2)，空域3的探测区域近似为一个圆柱体S，圆柱体S的中心轴在地面的坐标为(Lg2,La2,H2)，圆柱体的半径为R，高度分辨率为h₀，

步骤1.2，根据目标区域风廓线探测高度分辨率将目标区域划分成相应的网络格点，网格点采用立方体划分法，以目标位置高度为中心的立方体区域都计入该高度层，以任一探测海拔高度h为例，高度层h处的探测区域为圆柱体S上一个高度为h₀的圆柱形子区域，记为S_h，S_h下底面的海拔高度为S_h上底面的海拔高度为/>

步骤1.3，对于圆柱体S_h来说，高度h处的圆柱截面记为g，雷达1与g的圆心连线记为m，利用空间立体几何的理论求解平面g内与直线m的法向量平行的直径l₀，

将该直径l₀四等分，取l₀的5个等分点，分别求出5个等分点在地理坐标系下的坐标O₁(x₁,y₁,z₁)，O₂(x₂,y₂,z₂)，O₃(x₃,y₃,z₃)，O₄(x₄,y₄,z₄)，O₅(x₅,y₅,z₅)，分别连接点O与O₁、O₂、O₃、O₄、O₅，则向量(其中,i＝1,2,…,5)，激光波束分别沿着向量/>(其中，i＝1,2,…,5)的方向进行扫描，

步骤1.4，通过坐标变换，将激光波束的坐标从地理坐标系变换到雷达坐标系，在雷达坐标系下的扫描方位角和俯仰角分别为θ_i和γ_i(i＝1,2,…,5)，

测量时，在上位机软件输入雷达坐标、目标位置坐标及目标空域探测高度、高度分辨率等参数信息，上位机软件自动计算目标区域所有高度层的扫描方位角和俯仰角，以二维数组Coef[n][Attitude]建立所有高度层的扫描信息数据库，其中，n代表高度层，Attitude代表雷达扫描波束的方位角俯仰角信息，该数组在上位机软件运行时写入初始化文件。

3.如权利要求2所述的一种地基激光雷达目标区风廓线反演方法，其特征在于：所述的步骤2包括

步骤2.1，测量前雷达要先对北，控制雷达按照步骤1计算得到的扫描方位俯仰信息对目标探测区域进行三维体积扫描，扫描时按照逐层扫描的方式进行，并不断重复该过程，直至全部测量结束，

步骤2.2，获取每个激光波束的回波时域数据，回波时域数据按照距离大小将其划分成不同的距离单元，距离单元的长度为m₀，其中m₀应不大于数字接收机通过Hilbert变换将雷达信号从中频变换到基频，通过基于变分模态分解的联合时频分析方法将信号进行频谱变换，多个脉冲回波进行非相干累积得到原始功率谱数据，存储在二维数组R[n][gate_num]中，n代表波束序号，gate_num代表距离单元的序号，

步骤2.3，对每个距离单元的功率谱数据R[n][gate_num]先进行平滑去噪等基本处理，然后将回波信号建模为一个随机过程，通过计算随机过程的期望、方差等统计特征量，以“期望+权重系数×方差”作为自适应阈值检测门限，阈值之上的数据认为是信号包络，检测出的每个距离单元的信号包络，利用一阶矩得到包络质心，解算信号波峰位置和强度，进一步计算出多普勒频率及径向风速数据V_ri。

4.如权利要求3所述的一种目标区激光雷达风廓线反演方法，其特征在于：所述的步骤3包括

步骤3.1，综合考虑计算精度和计算规模，网格数量为目标区域探测的高度层总数N，网格的高度不超过目标区域的高度分辨率h₀，

对于特定高度层来说，将雷达扫描探测得到的各个距离单元的径向风速数据划分成两类，一类是落入上述网络格点范围内的点，即径向风所在距离单元在地理坐标系下的坐标落入该高度层所在网格点区域，记为C_rI，另一类是径向风所在距离单元在地理坐标系下的坐标未落入上述网络格点范围内的数据，记为C_rO，其中，r＝1，2，...，N，表示目标区域第r个高度层，

步骤3.2，对于目标区域每一个探测高度层来说，C_rI内的距离单元数据不能少于三个；若不满足该要求，则对C_rO内同一激光波束的各距离单元内的径向风速数据进行递归函数曲线拟合，拟合后的曲线记为f，

步骤3.3，依据上述步骤中的拟合函数f将数据插值到所需要的部分格点上，插值后的C_rI变为C‘_rI，实心圆点表示落入网格内的数据，空心圆点表示未落入网格点的数据，三角形表示插值的数据，

步骤3.4，结合各波束的坐标信息，对落入同一高度层风矢量数据集合C‘_rI，依据下述方程组求解，得到该高度层上的三维风矢量(u,v,w)，求解方程组采用共轭梯度法和有约束多元变量线性优化问题求最优解，可以避免出现采集参数的微小变化引起拟合结果的剧烈变化，

其中，V_i,h(i＝1,2,…,n)表示落入高度层h的第i个风速数据，θ_i表示风速数据V_i,h在雷达坐标系下对应的方位角，β_i表示风速数据V_i,h在雷达自身坐标系下对应的俯仰角。