基于风廓线雷达的沙尘暴沙尘质量浓度定量反演估算方法

摘要

本发明公开了一种基于风廓线雷达的沙尘暴沙尘质量浓度定量反演估算方法。该方法包括：第一步，计算沙尘暴过程雷达反射率因子Z；第二步，计算不同高度沙尘粒子的概率p(Dj)；第三步，计算沙尘数浓度N0和沙尘谱N(Dj)；第四步，计算沙尘质量浓度M；第五步，由于反射率因子Z和沙尘质量浓度M都与沙尘谱N(Dj)存在关系，建立雷达反射率因子Z与沙尘质量浓度M之间的关系式，得到Z＝AMb。本发明有效的解决了卫星遥感和微脉冲激光雷达定量探测沙尘暴的不足，实现了对沙尘暴过程高空沙尘质量浓度的定量反演估算。

说明书

技术领域

本发明涉及灾害天气监测领域，尤其涉及一种基于风廓线雷达的沙尘暴监测及其质量浓度定量反演估算方法。

背景技术

卫星遥感和微脉冲激光雷达是目前监测沙尘暴采用的技术手段。

卫星遥感的方法虽然可以有效的监测到沙尘暴的发生范围和移动路径，但无法获取沙尘暴的精细垂直结构、大气含沙量和沙尘质量浓度等定量信息。

微脉冲激光雷达的方法可以全天候实时观测沙尘气溶胶的垂直结构、成分和光学特性；但由于发射功率较小，在强沙尘暴期间激光不能穿透整个沙尘暴剖面，无法对强沙尘暴进行准确的定量监测。

故，亟需一种方法来实现实时对沙尘暴的准确定量监测。

发明内容

本发明的目的在于提供一种基于风廓线雷达的沙尘暴监测及其质量浓度定量反演估算方法，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明采用所述基于风廓线雷达的沙尘暴沙尘质量浓度定量反演估算方法，所述方法包括：

第一步，按照公式(1)计算沙尘暴过程雷达反射率因子Z；

$P_r=\frac{{\pi }^3}{1024*ln2}*\frac{P_t{\mathrm{h\tau cG}}^2\theta \psi }{2{\lambda }^2{R}^{2}L}*|\frac{m^2-1}{m^2+1}{|}^2*Z---\left(1\right);$

在式(1)中，P_r是接收回波功率，P_t是发射功率，c为电磁波传输速度，τ是脉冲宽度，h是雷达的有效探测高度，G是天线增益，L是双程损耗，θ是水平波速宽度，Ψ是垂直波速宽度，λ是发射的电磁波波长，R是雷达到探测目标物的距离，m是沙尘的复折射指数；

第二步，利用公式(2)估算不同高度沙尘粒子的概率p(D_j)；

$p\left(D_j\right)=\underset{j=1}{\overset{N}{\Sigma }}\frac{a}{\sqrt{2\pi }{\mathrm{\sigma D}}_j}\exp (-\frac{{({\mathrm{lnD}}_j-E)}^2}{2{\sigma }^2})\Delta D---\left(2\right);$

在式(2)中，p(D_j)表示直径为D_j的粒子的概率值，N表示离散化的总个数，a是调整因子，σ表示标准方差，D_j表示第j个沙尘粒子的直径，j表示积累变量，E表示期望值，ΔD表示离散化粒子间距。

第三步，根据公式(3)计算沙尘数浓度N₀和沙尘谱N(D)；

$Z=\underset{j=0}{\overset{Max}{\Sigma }}N\left(D_j\right)D_j^6\Delta D=\underset{D_j=0}{\overset{D_{max}}{\Sigma }}{N}_{0}p\left(D_j\right){D_j}^6\Delta D---\left(3\right);$

在式(3)中，Z表示雷达反射率因子，N(D_j)表示直径为D_j的沙尘粒子的沙尘粒子谱，j表示积累变量，D_j表示第j个沙尘粒子的直径，p(D_j)表示直径为D_j的沙尘粒子的概率值，Max表示变量j的最大值；

第四步，根据公式(4)计算沙尘质量浓度M；

$M=\frac{1}{6}\pi \rho \underset{j=1}{\overset{N}{\Sigma }}N\left(D_j\right){D_j}^3\Delta D---\left(4\right);$

在式(4)中，ρ表示沙尘的密度，N表示粒子总个数，N(D_j)表示直径为D_j的沙尘粒子的沙尘粒子谱，ΔD表示表示离散化粒子间距。

第五步，由于反射率因子Z和沙尘质量浓度M都与沙尘谱N(D_j)存在关系，采用最小二乘法方法建立雷达反射率因子Z与沙尘质量浓度M的关系，得到公式(5)；

Z＝AM^b>

式(5)中，A和b都是常数系数，完成基于风廓线雷达的沙尘暴沙尘质量浓度定量反演。

优选地，所述期望值E按照下述公式(6)计算：

E＝0.434e^(-1.827R)>

式(6)中，R是雷达探测目标物距离。

优选地，在第二步骤中，标准方差σ按照公式(7)计算：

σ＝0.081e^(0.847R)>

式(7)中，R是雷达探测目标物距离。

优选地，在浮尘时期，A＝20713.5，b＝0.995；在扬沙时期，A＝22988.3，b＝1.006；在沙尘暴时期，A＝24584.2，b＝1.013。

本发明的反演方法原理：利用风廓线雷达开展了对沙尘暴过程的探测试验，计算了沙尘暴的雷达反射率因子Z(即回波强度)，发现沙尘暴天气现象在雷达反射率因子图(即回波强度图)上有着清晰的印痕，可判识沙尘输送的高度和厚度，在此基础上，基于沙尘暴过程实测沙尘粒子谱资料，将多普勒天气雷达定量估测降水量的思路移植应用于风廓线雷达定量探测沙尘浓度，定量反演估算了沙尘暴期间高空沙尘粒子数浓度和质量浓度，建立了基于风廓线雷达的沙尘暴定量探测反演新方法。

本发明的有益效果是：

本发明利用风廓线雷达探测的功率谱数据计算得到了沙尘数浓度、沙尘谱和沙尘质量浓度，有效的解决了卫星遥感和微脉冲激光雷达探测沙尘暴天气的不足，实现了实时定量监测和分析沙尘暴天气的精细垂直结构和时空演变特点。本发明建立的反射率因子和沙尘质量浓度的关系式丰富了雷达气象学理论，可以进行扩展，具有很大的应用价值。

附图说明

图1是2014.05.22 14:00～05.23 10:00塔克拉玛干沙漠地区的反射率因子的时空变化图；

图2是图1的等值线图；

图3是集沙仪采集并经激光粒度仪分析后的2014年5月22-23日沙尘暴过程沙尘粒子概率分布，a图表示47m高度的沙尘粒子概率分布；b图表示63m高度的沙尘粒子概率分布；c图表示80m高度的沙尘粒子概率分布；

图4是不同时期沙尘谱的垂直分布，a图表示5月22日18:30时的扬沙期沙尘谱的垂直分布图；b图表示5月22日20:30时沙尘暴沙尘谱的垂直分布图；c图表示5月22日05:40浮尘期沙尘谱的垂直分布图；

图5是根据式(3)计算得到沙尘质量浓度M的时空变化；

图6是图5的等值线图；

图7是2014.05.22 14:00～05.23 10:00塔克拉玛干沙漠塔中地区采用本申请的Z-M关系式计算的700m高度质量浓度与采用沙尘谱反演的700m高度质量浓度对比示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

本申请所述基于风廓线雷达的沙尘暴沙尘质量浓度定量反演估算方法，所述方法包括：

第一步，计算沙尘暴过程雷达反射率因子Z

当风廓线雷达发射的电磁波遇到沙尘粒子群时，被沙尘粒子群散射形成回波，其中，后向散射回波可以被雷达接收系统接收。根据雷达气象方程得到沙尘粒子群的散射回波功率可以表示为公式(1)，按照公式(1)计算沙尘暴过程雷达反射率因子Z；

$P_r=\frac{{\pi }^3}{1024*ln2}*\frac{P_t{\mathrm{h\tau cG}}^2\theta \psi }{2{\lambda }^2{R}^{2}L}*|\frac{m^2-1}{m^2+1}{|}^2*Z---\left(1\right);$

在式(1)中，P_r是接收回波功率，P_t是发射功率，c为电磁波传输速度，τ是脉冲宽度，h是雷达的有效探测高度，G是天线增益，L是双程损耗，θ是水平波速宽度，Ψ是垂直波速宽度，λ是发射的电磁波波长，R是雷达到探测目标物的距离，m是沙尘的复折射指数；

第二步，计算沙尘粒子概率密度函数

沙尘谱表示在单位体积中沙尘粒子个数随粒子直径的分布；故为了计算第三步中的沙尘谱首先要确定沙尘粒子的概率密度函数，而，在1997年，董庆生对中国典型沙区的研究中发现实际沙尘粒子的概率密度函数符合对数正态分布(0-2)：

式(0-2)中p(D)表示概率密度函数、常数a是调整因子，σ表示标准方差和E表示期望值。对式(0-2)离散变形得到公式(2)。

$p\left(D_j\right)=\underset{j=1}{\overset{N}{\Sigma }}\frac{a}{\sqrt{2\pi }{\mathrm{\sigma D}}_j}\exp (-\frac{{({\mathrm{lnD}}_j-E)}^2}{2{\sigma }^2})\Delta D---\left(2\right);$

在式(2)中，p(D_j)表示直径为D_j的粒子的概率值，N表示离散化的总个数，a是调整因子，σ表示标准方差，D_j表示第j个沙尘粒子的直径，j表示积累变量，E表示期望值，ΔD表示离散化粒子间距；

基于沙尘暴过程近地层47m、63m和80m高度实测沙尘粒子概率分布资料，采用对数正态分布的概率函数(公式2)对47m、63m、80m高度的沙尘样品进行最小二乘法拟合，可得到47m、63m、80m高度的期望值E和标准方差σ。

在2010年，Zhibao Dong、Duoqing Mane对沙尘粒子分布特征研究中指出：沙尘粒子概率密度函数的标准方差和期望值随高度的变化符合指数分布。根据此，利用47m、63m和80m高度的期望值E和标准方差σ可推演得到沙尘暴过程高空不同高度的期望值E和标准方差σ；E＝0.434e^(-1.827R)；σ＝0.081e^(0.847R)；R是雷达探测目标物距离。

在此基础上，可利用公式(2)估算高空不同高度沙尘粒子的概率p(D_j)。

第三步，计算沙尘数浓度N₀和沙尘谱N(D_j)

雷达反射率因子Z和沙尘谱N(D_j)之间的关系为公式(3)，根据公式(3)计算沙尘数浓度N₀和沙尘谱N(D_j)；

$Z=\underset{j=0}{\overset{Max}{\Sigma }}N\left(D_j\right)D_j^6\Delta D=\underset{D_j=0}{\overset{D_{max}}{\Sigma }}{N}_{0}p\left(D_j\right){D_j}^6\Delta D---\left(3\right);$

在式(3)中，Z表示雷达反射率因子，N(D_j)表示直径为D_j的沙尘粒子的沙尘粒子谱，j表示积累变量，D_j表示第j个沙尘粒子的直径，p(D_j)表示直径为D_j的沙尘粒子的概率值，Max表示变量j的最大值；

第四步，计算沙尘质量浓度M

质量浓度表示在单位体积内含有的沙尘粒子的总质量，称为沙尘的质量浓度，沙尘的质量浓度M用沙尘谱N(D_i)表示为公式(4)，利用公式(4)计算沙尘质量浓度M；

$M=\frac{1}{6}\pi \rho \underset{j=1}{\overset{N}{\Sigma }}N\left(D_j\right){D_j}^3\Delta D---\left(4\right);$

在式(4)中，ρ表示沙尘的密度，N表示粒子总个数，N(D_j)表示直径为D_j的沙尘粒子的沙尘粒子谱，ΔD表示表示离散化粒子间距。

第五步，由于反射率因子Z和沙尘质量浓度M都与沙尘谱N(D_j)存在关系，采用最小二乘法方法建立雷达反射率因子Z与沙尘质量浓度M的关系，得到公式(5)；

Z＝AM^b>

式(5)中，A和b都是常数系数，完成基于风廓线雷达的沙尘暴沙尘质量浓度定量反演估算；利用公式(5)得到的雷达反射率因子Z与沙尘质量浓度M的关系计算得到沙尘暴过程高空不同高度的质量浓度；

其中，式(5)中在浮尘时期，A＝20713.5，b＝0.995；在扬沙时期，A＝22988.3，b＝1.006；在沙尘暴时期，A＝24584.2，b＝1.013。

使用并验证本发明所述方法可信性的实施例：

探测设备采用中国航天科工集团第二十三研究所研制的CFL-03型移动边界层风廓线雷达，所述探测设备的参数见表1。试验地点位于北半球中纬度欧亚大陆的塔克拉玛干沙漠腹地，坐落于中国新疆塔里木盆地中央，地理位置为39°00′N,83°40′E，海拔高度为1099.3m的塔中地区。

表2 CFL-03型移动边界层风廓线雷达参数

2015年5月22日-5月23日，受西西伯利亚冷空气影响，塔克拉玛干沙漠出现了一次明显的沙尘暴天气过程。根据地面观测资料表明：在5月22日14:00～18:50、21:12-21:46为扬沙，18:50～21:12、21:46～23:59为沙尘暴；5月23日00:00～03:40为沙尘暴，03:40～04:20为扬沙，04:20～07:00为浮尘。

利用风廓线雷达探测的垂直波束的速度功率谱数据计算得到5月22日14:00～5月23日10:00反射率因子的时空变化图1。在图1中可以得到：

①、在800m高度以下时：

扬沙阶段(14:00～18:50、21:12～21:46和03:40～04:20)反射率因子一般在8dBZ～10dBZ范围变化；

沙尘暴阶段(18:50～21:12、21:46～03:40和00:00～03:40)反射率因子一般大于10dBZ；

浮尘阶段(04:20～07:00)反射率因子在5dBZ～8dBZ范围变化；

②、在800m～1000m高度时：反射率因子明显下降，扬沙阶段降到了3dBZ左右，沙尘暴阶段降到了6dBZ左右，浮尘阶段降到了2dBZ以下；

③、在1000m～1200m高度时：在扬沙阶段、沙尘暴阶段和浮尘阶段的反射率因子均在4dBZ以下。

利用80米气象观测塔，在47米、63米和80米高度安装集沙仪，收集了2015年5月22日-5月23日沙尘暴过程三个高度层的沙尘样品，经激光粒度仪对沙尘样品进行粒子谱分析后得到沙尘粒子的概率分布特点为图3，在图3中可以看到沙尘粒子概率分布基本符合对数正态分布。

基于此次沙尘暴过程近地层47m、63m和80m高度实测沙尘粒子概率分布资料，采用对数正态分布的概率函数(公式2)对47m、63m、80m高度的沙尘样品进行最小二乘法拟合，可得到47m、63m、80m高度的期望值E和标准方差σ。

表1 表示47m、63m、80m高度的期望值E和标准方差σ

在2010年，Zhibao Dong、Duoqing Mane对沙尘粒子分布特征研究中指出：沙尘粒子概率密度函数的标准方差和期望值随高度的变化符合指数分布。根据此，利用47m、63m和80m高度的期望值E和标准方差σ可推演得到沙尘暴过程高空不同高度的期望值E和标准方差σ；E＝0.434e^(-1.827R)；σ＝0.081e^(0.847R)；R是雷达探测目标物距离。

在此基础上，可利用公式(2)估算出此次沙尘暴过程高空不同高度沙尘粒子的概率p(D_j)。

选择扬沙阶段5月22日18:30、沙尘暴阶段5月22日20:30和浮尘阶段5月23日05:40三个时刻，对不同高度的速度功率谱数据反演沙尘谱得到图4。由图4可以得到：在800m高度以下时，随高度的增大沙尘数浓度变化不大；在800m～1000m高度时，随高度的增大沙尘数浓度迅速减小了10^1.5m^-3mm^-1；当大于1000m高度时，随高度的增大沙尘数浓度缓慢减少。从图4c浮尘阶段到图4a扬沙阶段再到图4b沙尘暴阶段可以得到：相同粒子直径的沙尘数浓度依次增大了10^0.5m^-3mm^-1；随粒子直径的增大，沙尘数浓度减少的速度越来越小。

本申请在根据公式(2)和公式(3)反演得到的沙尘谱N(D_j)的基础上，根据公式(4)计算得到沙尘质量浓度M的时空变化为图5。结合反射率因子图1，由图5中可以得到：在整体上，沙尘质量浓度和反射率因子的变化趋势是一致，具体为：

①、在800m高度以下时：扬沙阶段的沙尘质量浓度在400μg/m³～700μg/m³范围变化，沙尘暴阶段的沙尘质量浓度一般在700μg/m³～1800μg/m³范围变化，最大时可以达到2000μg/m³，浮尘阶段的沙尘质量浓度在100μg/m³～400μg/m³范围变化；

②、在800m～1000m高度时：扬沙阶段沙尘质量浓度降到了100μg/m³左右，沙尘暴阶段沙尘质量浓度降到300μg/m³左右，浮尘阶段沙尘质量浓度在50μg/m³以下；

③、当大于1000m高度时，扬沙阶段沙尘质量浓度在50μg/m³以下，沙尘暴阶段沙尘质量浓度在150μg/m³以下，浮尘阶段的沙尘质量浓度基本为零。

如本领域人员公知：董庆生等人对腾格里沙漠的沙尘暴探测时得到近地面1.4m～1.6m高度的沙尘质量浓度在¹⁰⁴μg/m³～10⁶μg/m³范围变化；游来光和马培民等人利用飞机对沙尘暴探测时得到在1200m～1500m高度时沙尘质量浓度为在10²μg/m³～10³μg/m³范围变化。由图5中可以得到，本发明反演的600m～1200m高度的沙尘质量浓度在10²μg/m³～10⁴μg/m³范围变化。对比前人的研究，本发明反演的沙尘质量浓度的量级在合理的范围内，反演的结果具有可信性。

基于发生在塔克拉玛干沙漠2010年4月11日、2014年4月23日、2014年5月22日-5月23日三次沙尘暴过程风廓线雷达数据。利用反射率因子和反演对应时刻的沙尘质量浓度，通过非线性最小二乘法拟合得到浮尘阶段、扬沙阶段和沙尘暴阶段时的参数A和b为表3。

表3 不同天气Z-M系数取值

在表2中可以得到：A的值都大于20000，b的值在1的左右；从浮尘到扬沙再到沙尘暴，参数A和b的值都在增大。选择2014年5月22日-5月23日700m高度按照时间顺序利用反射率因子Z根据Z-M关系式计算的沙尘质量浓度和由沙尘谱计算的沙尘质量浓度的比较图7。

从图7中可以得到：利用本发明所述的Z-M关系计算得到的沙尘质量浓度与由沙尘谱计算得到的沙尘质量浓度随时间的变化趋势是一致的，并且两者的数值相差不大。所以利用表2的参数建立的浮尘、扬沙和沙尘暴的Z-M关系具有一定可信性。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：本发明利用风廓线雷达探测的功率谱数据计算得到了沙尘数浓度、沙尘谱和沙尘质量浓度，有效的解决了卫星遥感和微脉冲激光雷达探测沙尘暴天气的不足，实现了实时定量监测和分析沙尘暴天气的精细垂直结构和时空演变特点。本发明建立的反射率因子和沙尘质量浓度的关系式丰富了雷达气象学理论，可以进行扩展，具有很大的应用价值。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。