一种基于激光雷达风速数据的湍流参数反演方法与流程

本发明涉及湍流探测领域，尤其涉及一种基于激光雷达风速数据的湍流参数反演方法。

背景技术：

大气湍流是气流在三维空间内随空间位置和时间的不规则涨落，伴随着气流的涨落，温度和湿度等大气物质属性的浓度都呈现无规则的涨落。这种不规则运动主要由大气层中的风切变不稳定性、热力对流和波不稳定性引起。

大气湍流效应作为光学领域的一个分支，受到了广泛关注。激光通过大气传输时，大气湍流随机扰动造成大气折射率随机扰动，引起光强起伏、相位起伏、光束漂移以及光束扩展，成为制约自由大气光通信等激光应用技术发展的重要因素之一。当飞机遇到晴空湍流时，会出现左右摇晃或抖动的现象，严重时会影响航空飞行，危及乘客安全。边界层湍流还是影响大气污染扩散的重要因素，大气湍流强弱能直接影响大气对污染物的扩散能力。因此湍流探测对于自适应技术的发展、保障航空安全以及研究大气污染扩散都有重要意义。

折射率结构常数和湍流能量耗散率ε都是描述湍流强度的重要物理量。测量湍流强度的方法有探空法、以恒星为星源的被动测量法和激光雷达法。探空法利用探空气球搭载一对微温探测器测量路径上的温度脉动量来反演折射率结构常数原理简单但实时性不强，气球在高空会随风漂移；以恒星为星源的被动测量法受制于测量时间和双星的相对位置，空间分辨率低；而激光雷达法采用主动测量，利用回波中丰富的湍流信息，以其时空分辨率高，精度高，抗电磁干扰能力强，在恶劣环境下也能保持优良性能的优点，成为最有潜力的测量工具。

基于激光雷达风速数据反演湍流能量耗散率ε的方法中，速度结构函数法通过连续圆锥扫描获得径向风速，由径向风速计算实测速度结构函数，与大气湍流理论vonkarman模型中的理论速度结构函数直接拟合，以此计算湍流能量耗散率ε。但由于采用连续圆锥扫描方式和大气湍流理论中的vonkarman模型并直接拟合，存在径向风速的测速误差和噪声项的计算误差，造成耗散率的测量误差增大。

技术实现要素：

本发明的目的是提供一种基于激光雷达风速数据的湍流参数反演方法，可以测量耗散率以及其他湍流参数，且湍流参数的测量误差小。

本发明的目的是通过以下技术方案实现的：

一种基于激光雷达风速数据的湍流参数反演方法，包括：

控制激光雷达采用非连续圆锥扫描方式获得径向风速；

利用径向风速以及非连续圆锥扫描时的相关角度参数来计算径向风速脉动量；

利用径向风速脉动量计算实测方位角速度结构函数，并结合kolmogorov模型中的理论横向速度结构函数，来计算耗散率；

由径向风速和耗散率计算其他湍流参数。

所述非连续圆锥扫描方式包括：

激光雷达以固定天顶角α发射激光束，围绕垂直轴扫描360°方位角范围，在测量其中一个方位角对应的径向风速时，暂停扫描，保持方位角不变，直至测速结束，再继续扫描至下一个方位角。

所获得的径向风速表示为vr(rm,θn,l)，其中，rm＝r0+mδr是激光雷达与传感中心的距离，r0是激光雷达最小探测距离，m＝1,2,...,m，m是最大探测距离包含的距离分辨率个数，δr是距离分辨率，rm＝r0+mδr是最大探测距离；θn＝nδθ是方位角，n＝1,2,...,n，δθ是方位角分辨率，n＝360°/δθ；l＝1,2,...,l是扫描圈数。

所述利用径向风速以及非连续圆锥扫描时的相关角度参数来计算计算径向风速脉动量的公式为：

vr'(rm,θn,l)＝vr(rm,θn,l)-{cosα,sinαcosθn,sinαsinθn}·{vx,vy,vz}

式中，vr'(rm,θn,l)为径向风速脉动量，{vx,vy,vz}是平均风速矢量，由径向风速vr(rm,θn,l)采用风速反演方法得到。

利用径向风速脉动量计算实测方位角速度结构函数，并结合kolmogorov模型中的理论横向速度结构函数，来计算耗散率包括：

首先，计算实测方位角速度结构函数dr(kδθ)，其公式为：

式中，k≤ψ/δθ，ψ取决于雷达参数；

采用kolmogorov模型中的理论横向速度结构函数da(kδθ)，公式为：

da(kδθ)＝ε^2/3a(kδθ)

式中，φ(κ1,κ2)＝0.0652(κ1²+κ2²)^-4/3[1+(8/3)κ2²/(κ1²+κ2²)]，h||(κ1)＝[exp{-(πδpκ1)²}sinc(πδrκ1)]²是低频滤波器的纵向传输函数，δp＝cσp/2，c是光速，2σp是发射脉冲脉宽(e^-1功率电平)，δr＝ctw/2，tw是瞬时窗长；h⊥(κ2)＝[sinc(πkδθrmκ2)]²是低频滤波器的横向传输函数；

实测方位角速度结构函数dr(kδθ)与理论速度结构函数da(kδθ)满足以下关系式：

式中，是噪声项；

将k＝k1和k＝k2分别代入上式，并作差以消除噪声项：

dr(k1δθ)-dr(k2δθ)＝da(k1δθ)-da(k2δθ)

式中，k1和k2取值需满足k≤ψ/δθ；

最后，计算得到耗散率ε，公式为：

所述其他湍流参数包括：能量、积分尺度、动量通量和大气边界层高度；其中：

由耗散率ε计算能量e，公式为：

式中，g＝ε^2/3[f(k1δθ)+a(k1δθ)/2]，

由耗散率ε和能量e计算积分尺度lv，公式为：

lv＝0.3796e^3/2/ε

由径向风速计算动量通量＜uw＞和＜vw＞，公式为：

式中，是平均风向，j为虚数符号；

由耗散率ε测量大气边界层高度，选取阈值ε＝εblh，εblh对应的高度即为大气边界层高度。

所述风速反演方法包括：直接正弦波拟合法、过滤正弦波拟合法、传统最小二乘法、改进加权最小二乘法、牛顿-高斯最小二乘法与lm最小二乘法。

所述激光雷达为测风激光雷达，包括：相干多普勒测风激光雷达与直接多普勒测风激光雷达。

由上述本发明提供的技术方案可以看出，通过采用非连续圆锥扫描方式，抑制了多普勒频谱展宽，减小了径向风速的测速误差；采用大气湍流理论中的kolmogorov模型，并将速度结构函数作差消除噪声项后再计算耗散率，避免了噪声项带来的计算误差；与现有方法相比，不仅降低了湍流参数的测量误差，还可以测量除耗散率之外的能量、积分尺度、动量通量以及大气边界层高度。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例提供的一种基于激光雷达风速数据的湍流参数反演方法的流程图；

图2为本发明实施例提供的非连续圆锥扫描方式的示意图。

具体实施方式

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

本发明实施例提供一种基于激光雷达风速数据的湍流参数反演方法，如图1所示，其主要包括如下步骤：

步骤1、控制激光雷达采用非连续圆锥扫描方式获得径向风速。

如图2所示，为非连续圆锥扫描方式的示意图。所述非连续圆锥扫描方式是指，激光雷达以固定天顶角α发射激光束，围绕垂直轴扫描360°方位角范围，在测量其中一个方位角对应的径向风速时暂停扫描，保持方位角不变直至测速结束，再继续扫描至下一个方位角。这种方式使得测速期间激光束指向固定的方位角，抑制了多普勒频谱展宽，因此减小了径向风速的测速误差。

所获得的径向风速表示为vr(rm,θn,l)，其中，rm＝r0+mδr是激光雷达与传感中心的距离，r0是激光雷达最小探测距离，m＝1,2,...,m，m是常数，指的是最大探测距离包含的距离分辨率个数，最大探测距离随雷达参数变化，δr是距离分辨率，rm＝r0+mδr是最大探测距离；θn＝nδθ是方位角，n＝1,2,...,n，δθ是方位角分辨率，n＝360°/δθ；l＝1,2,...,l是扫描圈数。

步骤2、利用径向风速以及非连续圆锥扫描时的相关角度参数来计算径向风速脉动量。

本发明实施例中，计算径向风速脉动量vr'(rm,θn,l)的公式为：

vr'(rm,θn,l)＝vr(rm,θn,l)-{cosα,sinαcosθn,sinαsinθn}·{vz,vx,vy}(1)

式中，{vz,vx,vy}是平均风速矢量，由径向风速vr(rm,θn,l)采用风速反演方法得到。

示例性的，所述风速反演方法包括：直接正弦波拟合法、过滤正弦波拟合法、传统最小二乘法、改进加权最小二乘法、牛顿-高斯最小二乘法与lm最小二乘法。

步骤3、利用径向风速脉动量计算实测方位角速度结构函数，并结合kolmogorov模型中的理论横向速度结构函数，来计算耗散率。

本发明实施例中，利用径向风速脉动量vr'(rm,θn,l)计算实测方位角速度结构函数dr(kδθ)，公式为：

式中，k≤ψ/δθ，ψ取决于雷达参数。

示例性的，所述雷达参数中天顶角α＝30°时ψ＝10°，满足k≤ψ/δθ时实测方位角速度结构函数dr(kδθ)与理论横向速度结构函数da(kδθ)近似相等。

采用kolmogorov模型中的理论横向速度结构函数da(kδθ)，公式为：

da(kδθ)＝ε^2/3a(kδθ)(3)

式中，κ1、κ2均为积分变量，φ(κ1,κ2)＝0.0652(κ1²+κ2²)^-4/3[1+(8/3)κ2²/(κ1²+κ2²)]，h||(κ1)＝[exp{-(πδpκ1)²}sinc(πδrκ1)]²是低频滤波器的纵向传输函数，δp＝cσp/2，c是光速，2σp是发射脉冲脉宽(e^-1功率电平)，δr＝ctw/2，tw是瞬时窗长；h⊥(κ2)＝[sinc(πkrmsinαδθκ2)]²是低频滤波器的横向传输函数。

实测方位角速度结构函数dr(kδθ)与理论速度结构函数da(kδθ)满足以下关系式：

式中，是噪声项。

本发明实施例中，将k＝k1和k＝k2分别代入公式(5)，并作差以消除噪声项：

dr(k1δθ)-dr(k2δθ)＝da(k1δθ)-da(k2δθ)(6)

式中，k1和k2取值需满足k≤ψ/δθ。

将公式(2)(3)(4)代入公式(6)，计算得到耗散率ε，公式为：

步骤4、由径向风速和耗散率计算其他湍流参数。

本发明实施例中，所述其他湍流参数包括：能量、积分尺度、动量通量和大气边界层高度；其中：

由耗散率ε计算能量e，公式为：

式中，g＝ε^2/3[f(k1δθ)+a(k1δθ)/2]，

由耗散率ε和能量e计算积分尺度lv，公式为：

lv＝0.3796e^3/2/ε(9)

由径向风速计算动量通量＜uw＞和＜vw＞，公式为：

式中，是平均风向，j为虚数符号。

由耗散率ε测量大气边界层高度，选取合适的阈值ε＝εblh，εblh对应的高度即为大气边界层高度。

示例的，内陆夏季情况下，εblh＝10^-4m²/s^-3。

本发明实施例上述方案中所涉及的激光雷达可以为测风激光雷达，具体可以为相干多普勒测风激光雷达与直接多普勒测风激光雷达。

示例性的，只有固定天顶角时，可以由径向风速和耗散率计算能量。

本发明实施例上述方案主要具有如下有益效果：

1)本发明采用非连续圆锥扫描方式，减小了径向风速的测速误差；

2)本发明采用kolmogorov模型并将速度结构函数作差以消除噪声项，避免了噪声项带来的计算误差；

3)本发明不仅通过上述有益效果降低了湍流参数的测量误差，还可以测量除耗散率之外的能量、积分尺度、动量通量以及大气边界层高度。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例可以通过软件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是cd-rom，u盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。