CN116879899A - 一种基于空中降水粒子谱反演的方法 - Google Patents
一种基于空中降水粒子谱反演的方法 Download PDFInfo
- Publication number
- CN116879899A CN116879899A CN202310816919.7A CN202310816919A CN116879899A CN 116879899 A CN116879899 A CN 116879899A CN 202310816919 A CN202310816919 A CN 202310816919A CN 116879899 A CN116879899 A CN 116879899A
- Authority
- CN
- China
- Prior art keywords
- raindrop
- spectrum
- power spectrum
- ground
- precipitation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000001228 spectrum Methods 0.000 title claims abstract description 200
- 238000001556 precipitation Methods 0.000 title claims abstract description 84
- 238000000034 method Methods 0.000 title claims abstract description 63
- 239000002245 particle Substances 0.000 title claims abstract description 63
- 230000008859 change Effects 0.000 claims abstract description 12
- 230000033001 locomotion Effects 0.000 claims description 36
- 238000002310 reflectometry Methods 0.000 claims description 27
- 238000005070 sampling Methods 0.000 claims description 16
- 238000001514 detection method Methods 0.000 claims description 15
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 13
- 230000002159 abnormal effect Effects 0.000 claims description 12
- 238000004422 calculation algorithm Methods 0.000 claims description 10
- 238000012937 correction Methods 0.000 claims description 9
- 239000007788 liquid Substances 0.000 claims description 9
- 238000012545 processing Methods 0.000 claims description 8
- 238000003908 quality control method Methods 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 7
- 230000011218 segmentation Effects 0.000 claims description 7
- 238000010606 normalization Methods 0.000 claims description 6
- 230000003595 spectral effect Effects 0.000 claims description 6
- 238000004458 analytical method Methods 0.000 claims description 5
- 238000001914 filtration Methods 0.000 claims description 4
- 238000006467 substitution reaction Methods 0.000 claims description 4
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 238000005315 distribution function Methods 0.000 claims description 3
- 239000012798 spherical particle Substances 0.000 claims description 3
- 239000000463 material Substances 0.000 claims description 2
- 230000008569 process Effects 0.000 abstract description 15
- 238000011161 development Methods 0.000 abstract description 5
- 230000015572 biosynthetic process Effects 0.000 abstract description 4
- 230000007246 mechanism Effects 0.000 abstract description 4
- 230000000694 effects Effects 0.000 description 8
- 230000018109 developmental process Effects 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 230000000630 rising effect Effects 0.000 description 2
- 230000035945 sensitivity Effects 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 241000533950 Leucojum Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008033 biological extinction Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000003475 lamination Methods 0.000 description 1
- 238000002156 mixing Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 239000005436 troposphere Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/95—Radar or analogous systems specially adapted for specific applications for meteorological use
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种基于空中降水粒子谱反演的方法，属于反演技术领域，解决现有技术无法精确得到地面和空中微物理参数反演参数，即无法了解大气云降水的形成、发展和变化机理，从而无法准确地预测和解释大气过程和气象灾害。包括基于C波段调频连续波雷达获取观测功率谱，并采用去噪方法去除噪声电平，其中，观测功率谱即指多普勒功率谱；基于风廓线雷达获取水平风构建大气功率谱；基于地面雨滴谱仪获取地面雨滴数据和标准化gamma分布构建雨滴谱；基于上述结果进行地面微物理参数反演和空中微物理参数反演。本发明用于反演获取地面和空中微物理参数。
Description
技术领域
一种基于空中降水粒子谱反演的方法，用于反演获取地面和空中微物理参数，属于反演技术领域。
背景技术
空中降水粒子谱通常指的是大气中的雨滴、雪花、冰粒等降水粒子的大小分布情况。空中降水粒子谱的变化是对云中复杂的降水粒子相互作用的最直接描述。同时降水粒子作为天气雷达的最直接观测对象，由于雷达探测原理的限制，不同探测量对于不同尺寸雨滴的敏感度不同。降水粒子谱的分布差异也将直接影响到雷达降水估计的精度。但现有技术存在如下技术问题：
无法精确得到地面和空中微物理参数反演参数，即无法了解大气云降水的形成、发展和变化机理，从而无法准确地预测和解释大气过程和气象灾害。
发明内容
针对上述研究的问题，本发明的目的在于提供一种基于空中降水粒子谱反演的方法，解决现有技术无法精确得到地面和空中微物理参数反演参数，即无法了解大气云降水的形成、发展和变化机理，从而无法准确地预测和解释大气过程和气象灾害。
为了达到上述目的，本发明采用如下技术方案：
一种基于空中降水粒子谱反演的方法，包括如下步骤：
步骤1、基于C波段调频连续波雷达获取观测功率谱，并采用去噪方法去除噪声电平，其中，观测功率谱即指多普勒功率谱；
步骤2、基于风廓线雷达获取水平风构建大气功率谱；
步骤3、基于地面雨滴谱仪获取地面雨滴数据和标准化gamma分布构建雨滴谱；
步骤4、基于步骤1-3得到的结果进行地面微物理参数反演和空中微物理参数反演。
进一步，所述步骤1的具体步骤为：
步骤1.1、基于C波段调频连续波雷达获取观测功率谱，即是卷积运算后的降水粒子和大气湍流运动的多普勒功率谱，多普勒功率谱的公式为：
Sobs(v，r)＝Sdsd(v-w，r)*Sair(v，r)+Snoise(v，r)
其中，Sobs(v，r)、Sair(v，r)、Sdsd(v-w，r)和Snoise(v，r)分别是实测的多普勒功率谱、大气湍流运动多普勒功率谱、降水粒子多普勒功率谱以及实测的多普勒功率谱的噪声电平，v、r、w分别表示径向速度、库数和大气垂直运动速度；
步骤1.2、采用去噪方法去除观测功率谱中的噪声电平，去噪方法包括分段法、小波去噪法和局部时频分析法；
分段法是设C波段调频连续波雷达噪声服从中心χ分布，自由度为2N/k，将观测功率谱划分成k段，统计每段的平均功率谱，最小值作为噪声电平，并去除，即得到去噪后的观测功率谱；
小波去噪法是通过小波分解和滤波去除观测功率谱中包含大于给定噪声数量的小波系数，再将得到的系数重新组合成时间域信号，即得到去噪后的观测功率谱；
局部时频分析法是对观测功率谱中LT区域的特征直方图和累计分布函数进行分析，对区分噪声和信号明显的时频小区域进行噪声过滤，即得到去噪后的观测功率谱，其中，LT是指局部时频。
进一步，所述步骤2的具体步骤为：
步骤2.1、风廓线雷达采用5波束探测模式得到垂直波束，倾斜角度均为15°的正北、正东、正南、正西方向的波束，其中，垂直波束包括得到垂直于地面的风速和风向分量，正北、正东、正南、正西方向的波束分别得到风速和风向在对应方向上的分量；
步骤2.2、基于获取的样本数据的平均值，基于步骤2.1得到水平风，比较相邻时间或相邻高度上的风速和风向的分量，将差异超过预定阈值的数据视为异常值，然后用平均值进行替换，即是通过判断垂直波束和4个倾斜方向波束的垂直速度差异和两组对应水平风差异是否超过给定阈值，如果差异超过预先设定的阈值，则限制或排除异常数据，否则将其保留；
步骤2.3、基于步骤2.2得到质控后的水平风构建大气功率谱。
进一步，所述步骤3的具体步骤为：
步骤3.1、地面雨滴谱仪获取雨滴数据，包括地面雨滴大小、下落速度和单位采样时间内经过采样面积的雨滴个数；
步骤3.2、基于获取的地面雨滴大小对其进行修正，修正公式为：
式中：D为订正后的地面雨滴的直径，Dpar为实测的地面雨滴的直径；
步骤3.3、去除修正直径为6-8mm的雨滴后，基于标准化gamma分布构建雨滴谱，公式为：
N(D)＝N0Dμexp(-ΛD)
其中，N0表示截距参数，μ表示形状因子，人表示斜率参数。
进一步，所述步骤4的具体步骤为：
步骤4.1、基于雨滴数据和雨滴谱进行地面微物理参数反演，得到地面微物理参数，包括地面雨滴数量浓度、地面雨强、地面液态含水量和地面雷达反射率因子；
步骤4.2、基于多普勒功率谱、大气功率谱和雨滴谱进行空中物理参数反演，得到空中物理参数，包括空中雨滴数量浓度、空中雨强、空中液态含水量、空中雷达反射率因子、垂直大气运动和空中雨滴谱。
进一步，所述步骤4.1的具体步骤为：
步骤4.11、基于地面雨滴下落速度和单位采样时间内经过采样面积的雨滴个数得到雨滴数浓度，公式为：
式中：nij为修正后的雨滴直径位于第i个尺度通道且下落速度位于第j个速度通道的雨滴数，A为仪器采样面积，Δt为采样时间间隔，Di和ΔDi为修正后的地面雨滴在第i个尺度通道的中心直径和该区间的尺度间隔，Vj为第j个速度通道的中心速度；
步骤4.12、基于雨滴数浓度计算得到地面雨强R1、地面液态含水量LWC1和地面雷达反射率因子Z1，计算公式为：
式中：ρ为水的密度；
步骤4.13、基于雨滴谱采用阶矩法对标准化gamma分布参数进行拟合，拟合时引入雨滴谱的阶距Mn，其中，标准化gamma分布的参数N0、μ、Λ用雨滴谱的2、4、6阶距估算得到，公式为：
其中，n表示第n阶矩；
步骤4.14、基于阶距Mn得到地面雨滴的质量加权平均直径D1m，公式为：
步骤4.15、基于质量加权平均直径D1m，得到归一化的数量浓度N1w，计算公式为：
其中，w表示大气垂直运动速度。
进一步，所述步骤4.2的具体步骤为：
步骤4.21、基于大气功率谱，采用傅里叶反卷积算法消除去噪后的多普勒功率谱中的大气谱展宽；
步骤4.22、基于雨滴谱构建初始高度降水粒子功率谱，即在静稳大气条件下，降水粒子多普勒功率谱Sdsd(v，r)用雨滴谱分布Ndsd(D)和后向散射截面表示为：
其中，Sdsd(v，r)表示降水粒子多普勒功率谱，即初始高度降水粒子功率谱，λ为雷达波长，为5.42cm，|Kw|2为介电常数，表示下落末速度与雨滴直径关系，D为修正后的雨滴；
当波长为5.42cm时，此时降水粒子发生瑞利散射，若将其视为球形粒子，此时，雷达后向散射截面与粒子直径的6次方成正比，后向散射截面可以表示为：
其中，「表示欧拉函数，μ表示形状因子，Λ表示斜率参数，其中的积分从0延伸到无穷大，反演过程中雨滴谱受到最小和最大雨滴直径的限制，Dm表示地面雨滴的质量加权平均直径；
空中雨强R和空中液态水含水量LWC均是空中雨滴大小分布的三阶矩：
其中，ΔD表示空中雨滴直径的变化率，V(D)表示空中雨滴下降速度与直径之间关系，p0代表地面的大气密度，ph表示地面的大气密度，h为观测高度；
步骤4.23、基于步骤4.21得到的结果获取初始高度大气运动，并剔除初始高度降水粒子功率谱中相邻的分布差异；
步骤4.24、空中雷达反射率因子Z是雨滴谱的零阶矩，公式为：
代入：
N(D；Nw，Dm，μ)＝Nwf(D；Dm，μ)
得到表示纯降水粒子在瑞利散射时的反射率，即空中雷达反射率因子Z为：
其中，Δv分别表示空中雨滴下落速度的变化率；
步骤4.25、基于步骤4.23得到的结果，对每一距离库的雨滴谱逐一进行归一化处理，并通过最小化代价函数得到垂直大气运动；
归一化公式为：
最小化代价函数为：
步骤4.26、基于步骤4.23和步骤4.25得到的结果，反卷积得到空中雨滴谱。
本发明同现有技术相比，其有益效果表现在：
本发明基于订正了雨滴后得到的雨滴谱数据，使用标准化gamma雨滴谱分布，实时建立雨滴谱分布模型为空中降水微物理参数反演提供初始值，基于此前提下，利用排除谱展宽后的垂直指向雷达数据，通过对反射率因子分解提取因空中雨滴谱演变带来的下落速度差异在雷达功率谱上造成的相应差异，进而基于归一化功率谱得到大气垂直速度(即指垂直大气速度)，基于剔除大气垂直速度的降水功率谱，利用粒子直径与速度的关系，得到空中雨滴谱分布，本发明能精确得到地面和空中微物理参数反演参数，即能利用这些参数有效了解大气云降水的形成、发展和变化机理，从而准确地预测和解释大气过程和气象灾害。
附图说明
图1为本发明中空中微物理参数反演流程；
图2为本发明中层云个例反演示例示意图；
图3为本发明中积层混合云个例反演示例；
图4为本发明中对流云个例反演示例；
图5为本发明中的对比验证散点图。
具体实施方式
下面将结合附图及具体实施方式对本发明作进一步的描述。
观测设备：
C波段调频连续波雷达(C-band FMCW：Frequency Modulation Continuous Waveradar)有别于其他气象雷达，采用连续波体制调制发射信号频率。通过频率调制发射电磁波，将发射频率从低到高逐渐增加，这个过程是通过向频率源注入线性调频信号来实现的。当电磁波信号遇到目标并反射回来时，FMCW雷达接收到的信号频率和发射频率之差就是目标相对速度所对应的多普勒频移。通过测量这个多普勒频移，可以确定目标的速度。同时，FMCW雷达还可以利用接收到的多普勒频移的变化来确定目标的位置。具体来说，当雷达向目标发射信号时，会同时记录下发射信号的时间。当接收到目标反射回来的信号时，也会记录下接收到信号的时间。通过计算发射信号和接收信号之间的时间差，可以确定信号传播的时间，从而确定目标与雷达之间的距离。综上所述，FMCW雷达可以实现对目标的高精度探测，它的优点是测距和测速的性能都非常出色，同时可以在长距离和弱信号环境下实现探测，并且可以实时采集到目标的位置、速度、方位等信息，对于对流层、强对流等气象现象的探测具有较高的分辨率。同时C波段的频率范围较小，波长较短，在大气层中能够较好地穿透云层，可以实现云层内部物理过程的探测，探测效果较好，这使针对复杂的真实云中微物理过程探测成为可能。
研究使用的雷达具有较高的时空间分辨率分别为3s和30m，探测高度范围0.03～15km，基本覆盖了整个华南降水过程形成发展的全高度且探测无盲区，足以保证对降水过程的全高度连续观测。此外，雷达的最小可检测信号功率为-166dBm，相较于常规天气雷达的-115dBm更优秀，能够更好地反映出降水过程中的细节信息。
风廓线雷达(Wind Profiling Radar)是一种利用雷达技术测量气象和气候要素的仪器。它可以实时监测和测量不同高度范围内的水平和垂直风速，同时还能反演温度、湿度和涡度等大气参数。相较于传统中的探空气球能提供更加充足的各时段信息。风廓线雷达的原理是利用多普勒效应来测量大气中的风速。当雷达发射出一个频率稳定的电磁波向大气上空发送时，受到风场的影响，回波的频率就会发生多普勒频移，从而反映出大气中的风速情况。同样地，风廓线雷达可以利用多普勒效应通过反演其他大气参数。
风廓线雷达其高度分辨率在不同的探测模式下分别为60米和120米。同时，其时间分辨率为6分钟，探测范围分别覆盖了0.1至2.86千米和0.99至5.55千米的范围。除此之外，该雷达还能够提供多项基本观测数据产品，其中包括水平风、垂直风、径向风以及折射率结构常数。这些数据能够被用于实时监测一个较大范围内的天气变化，并提供具有高水平和垂直分辨率的详细大气参数信息。
地面雨滴谱仪工作原理是，利用两个激光束打在雨水粒子上，雨水粒子反射回的光被接收后，通过时间分辨技术和回波强度来分析雨水粒子的质量和速度。因此，PARSIVEL激光雨滴谱仪不仅可以实时监测雨滴的大小、速度和数量等信息，还可以得到准确的雨滴谱图。它通过同时监测雨滴的大小、速度和数量等信息，可以得到更为精细的雨滴谱分布，对于天气预报、气候研究等方面具有很重要的应用价值。采样对象为区域(18cm×3cm)内所有粒子的直径和速度信息，时间分辨率为1分钟。该仪器的信号处理器将粒子的直径和速度分为32个通道，反演出雨量、雨强、反射率强度、质量加权平均直径(质量/体积权重平均直径)Dm、归一化截距参数Nw等参数。
数据质控：
C波段调频连续波雷达的数据质控主要聚焦在多普勒功率谱数据的质控。多普勒功率谱数据是根据雷达回波信号的多普勒频移(即气象目标速度引起的频率变化)进行处理得到的。在雷达多普勒模式下，雷达将发射一段连续的脉冲，这些脉冲将对气象目标进行扫描，得到回波信号后，通过多普勒处理，提取出气象目标的速度信息和功率谱信息。因此，气象雷达多普勒功率谱数据是雷达的最基础数据，回波强度、径向速度、速度谱宽等具体参数都可以通过多普勒功率谱数据计算得到。对于多普勒功率谱的数据质控主要分为两方面，噪声电平的去除与功率谱展宽的订正。
噪声电平是多普勒功率谱中需要去除的主要部分。常用的去噪方法包括分段法、小波去噪和局部时频噪声过滤等。分段法是设C波段调频连续波雷达噪声服从中心χ分布，自由度为2N/k，将观测功率谱划分成k段，统计每段的平均功率谱，最小值作为噪声电平，并去除。小波去噪法是通过小波分解和滤波去除观测功率谱中包含大于给定噪声数量的小波系数，再将得到的系数重新组合成时间域信号，即得到去噪后的观测功率谱。局部时频噪声处理结合了时域、频域和局部区域信息，能更好地降低噪声电平。局部时频分析法是对观测功率谱中LT区域的特征直方图和累计分布函数进行分析，对区分噪声和信号明显的时频小区域进行噪声过滤，即得到去噪后的观测功率谱，其中，LT是指局部时频。在三种方法之中，分段法由于其高稳定性和准确性以及快速处理效率被采用。该方法在去除噪声的同时保留了全部信息。
为了获取更加精确的雨滴谱分布，消除大气谱展宽对多普勒功率谱形状的影响是必不可少的。大气湍流功率谱会对降水功率谱产生卷积作用，从而导致观测所得的多普勒功率谱方差增大，进而引起谱展宽。采用了直接反卷积算法来进行谱展宽订正。我们基于质控后的风廓线雷达观测水平风数据，建立了大气功率谱，并使用傅里叶反卷积算法来消除C波段调频连续波雷达观测多普勒功率谱中由大气谱展宽导致的影响。通过采用这种方法，成功地消除了大气湍流功率谱对谱展宽造成的影响，进而可以获取更加准确的雨滴谱分布。
风廓线雷达数据可以用于获取大气中的高度分辨率水平风数据，但这些数据可能会被环境噪声、干扰和雷达本身的问题影响。尤其是在降水期间，风廓线雷达的数据质量会受到许多影响，包括大气湍流、降水、湍流剪切等，因此需要进行更严格的数据质控。研究中的风廓线雷达采用常用的5波束探测模式，该模式在空间分辨率和时间分辨率上具有一定优势。具体来说，5波束探测模式包括以下5个波束：垂直波束可以得到垂直于地面的风速和风向信息，是其他波束的基础。正北、正东、正南、正西方向波束(倾斜角度均为15°)，分别得到风速和风向在对应方向上的分量。通过同时收集这5个波束的数据，并进行合理的处理和分析，可以得到更加全面、准确的风场信息。
风廓线雷达数据质控采用一致性平均方法，它的基本思想是使用样本数据的平均值作为参考，比较相邻时间或相邻高度上的风速、风向数据，将差异超过预定阈值的数据视为异常值，然后用平均值替换这些异常值。即是通过判断垂直波束和4个倾斜方向波束的垂直速度差异和两组对应水平风差异是否超过给定阈值。如果差异超过预先设定的阈值，则限制或排除异常数据，否则将其保留。重复以上步骤，直到所有数据都经过质控为止。
地面降水测量使用0TT Parsivel，可测量雨滴的大小和下落速度。然而，雨滴在下落过程中会受到形变的影响。一个孤立的、相对较大的雨滴在下落末速度落下时会呈现出明显的扁平化，并在整体上呈椭球形，而不是流线型的圆锥形状，通常描述为“泪滴”的形状。这会导致雨滴测量偏大。因此，需要进行形变修正，其公式如下：
式中：D为订正后的等效直径(mm)；Dpar为实测的直径(mm)。同时在自然界中通常认为直径6-8mm的雨滴(具体取决于湍流强度)是不稳定的，故而也将舍弃这部分数据。
地面微物理参数反演
地面微物理参数反演基于Parsivel的观测数据。观测数据为单位采样时间内经过采样面积的雨滴个数。因此，雨滴数浓度N(m-3·mm-1)计算如下：
式中：nij为雨滴直径位于第i个尺度通道且下落速度位于第j个速度通道的雨滴数；A为仪器采样面积，54em2；Δt为采样时间间隔，60s；Di和ΔDi为第i个尺度通道的中心直径(mm)和该区间的尺度间隔(mm)；Vj为第j个速度通道的中心速度(m·s-1)。
地面雨强R1(mm·h-1)、地面液态含水量LWC1(mg·m-3)、地面雷达反射率因子Z1可以通过以下公式计算：
式中：ρ为水的密度(g·cm-1)。
目前gamma分布被广泛地应用于描述地面以及空中的雨滴谱，基于雨滴谱数据，公式可以表达为：
N(D)＝N0Dμexp(-ΛD) (3-5)
N0是截距参数(m-3·mm-μ-1)；μ是形状因子；人是斜率参数(mm-1)。采用阶矩法对gamma分布参数进行拟合，拟合时引入雨滴谱阶距Mn其下标为第n阶矩，gamma分布的参数(N0，μ，A)可用观测值的2、4、6阶距估算，计算公式如下：
进而质量加权平均直径Dm(mm)可表示为：
归一化的数量浓度NW(mm-1·m-3)为
空中微物理参数反演
各种降水探测设备在实际应用中获取的返回信号实际上是由降水粒子的散射信号和大气湍流散射信号的综合影响引起的。采用C-FMCW雷达技术进行测量，获得的实测功率谱，其本质是卷积运算后的降水粒子和大气湍流运动的多普勒功率谱。以下是对应的数学公式[37]：
Sobs(v，r)＝Sdsd(v-w，r)*Sair(v，r)+Snoise(v，r) (3-9)
式中：Sobs(v，r)[(mm6m-3)(ms-1)-1]、Sair(v，r)、Sdsd(v-w，r)与Snoise(v，r)分别是实测多普勒功率谱、大气湍流运动多普勒功率谱、降水粒子多普勒功率谱以及实测多普勒功率谱噪声。变量v、r、w分别表示径向速度和库数与大气垂直运动速度。对于观测到的多普勒功率谱的直流分量和噪声电平，我们采用了预处理方法。其中，噪声电平指的是多普勒功率谱中所有雷达噪声的平均功率，而这个数值会随着每个距离库的变化而发生变化。此外，由于大气湍流和大气垂直运动的影响，多普勒功率谱具有左右移动的倾向，这可能导致我们检索到的雨滴谱被低估或高估。例如，详细的误差统计分析表明，当湍流的垂直速度偏差为1ms-1时，在雨滴谱反演中，雨率的误差可能高达50％。因此，我们无法直接获得雨滴谱的分布并进行定量分析。
本案提供了一种用于进行大气垂直运动反演的算法。在该过程中，我们需要考虑到降水的影响，因为强烈的大气垂直运动会促进降水粒子各种演变过程，从而影响降水粒子功率谱在多普勒速度轴上的分布，导致反演得到的雨滴谱可能会低估或高估。因此，我们采用了一种新的原理来提取大气垂直速度——通过比较相邻两个距离库的降水粒子功率谱的差异来分析降水粒子功率谱的变化。这种差异是由大气垂直速度和雨滴谱在下落过程中的演变所导致的。通过提取大气垂直速度，我们可以减轻这种影响并提高雨滴谱的反演精度。
在静稳大气条件下，降水粒子多普勒功率谱Sdsd(v，r)可用雨滴谱分布Ndsd(D)和后向散射截面表示为：
式中：λ为雷达波长，5.42cm；|Kw|2为介电常数，0.93；表示下落末速度与雨滴直径关系。
当波长为5.42cm时，此时降水粒子发生瑞利散射，若将其视为球形粒子。此时，雷达后向散射截面与粒子直径的6次方成正比，后向散射截面可以表示为：
雨滴谱若用三参数的gamma分布来表示：
N(D)＝N0Dμexp(-ΛD) (3-12)
其中，μ(mm-1)、Λ(mm-1)和N0(mm-1-μm-3)分别为gamma分布中的形状、斜率和截距参数。
Willis(1984)给出了DSD的另一种表达式，用于比较具有不同液态水含量值的雨滴谱，称为标准化gamma雨滴谱分布，表示为：
式中为Nw归一化数浓度参数，Dm(mm)为质量加权平均直径，μ为gamma形状因子。「表示欧拉函数。其中的积分从0延伸到无穷大。实际上，反演过程中雨滴谱受到最小(Dmin)和最大(Dmax)雨滴直径的限制。
降水参数可以用雨滴的大小分布来计算。空中雨率R(mm·h-1)和空中液态水含水量LWC(g·m-3)均是液滴大小分布的三阶矩：
在理想的瑞利后向散射情况下，不同参数对于雨滴大小分布的敏感性是各不相同的。例如，反射率更容易受到大雨滴的影响，但是相较于同样反射率条件下的小雨滴，它们对LWC的贡献较小；同时，云滴对LWC的贡献也是显著的[39]。
此时雨滴谱也可以使用归一化gamma分布建模，雨滴谱分布也可以写成，
N(D；Nw，Dm，μ)＝Nwf(D；Dm，μ) (3-18)
其中
利用公式(3-18)利用雨滴大小分布计算降水粒子功率谱：
式中，D6表示假设瑞利散射的雨滴背散射截面，ΔD和Δv分别表示雨滴直径和下落速度的变化率。
此时，空中雨滴下降速度与直径之间关系可以表达为：
式中，在观测高度为h处，地面的大气密度为ρh，而ρ0则代表地面的大气密度。两者的比值受到许多因素的影响，包括地面温度、大气温度随高度变化的速率和观测高度。在低空条件下，可以将这个比值简化成只与高度有关的函数。
此时空中雨滴的下降速度是纯雨滴谱的一阶矩，可以表示为：
同时，反射率Z是雨滴谱的零阶矩，可以表示为：
代入公式(3-18)、公式(3-19)、公式(3-20)、可转化为：
上式表示了纯降水粒子在瑞利散射时的反射率。
因为反射率因子Z(r)(mm-6m-3)和Dm(mm)使用不同的单位来测量的。取公式(3-25)两侧的对数，ZdB(r)(dBZ)可以分解为数浓度函数NwdB(dB)和粒子谱形状数即以下形式：
为反射率形状因子，等价于归一化的后向散射反射率，即Nw＝1的反射率。由于Nw在公式(3-29)中被抵消，因此形状函数和雨滴下降速度仅取决于形状因子μ和平均质量加权直径Dm。/>
也可以表达为如下形式：
其中，f为雷达工作频率。
由于相邻距离库的多普勒速度谱差异主要受下降过程中湍流和雨滴谱变化影响，可通过该原理提取大气的垂直速度。对两个多普勒功率谱进行归一化处理后，计算它们间的差异可得到谱在多普勒速度轴上的位置。该位置反映了垂直大气运动和雨滴速度差的总和。
对每一距离库的雨滴谱逐一进行归一化处理。
确认雨滴速度差后，知/>之间的差异应取决于垂直大气运动。垂直大气运动可以通过最小化代价函数(Cost Function，CF)来得到。
当CF值达到最小时可以认为垂直大气运动被完全消除。垂直大气运动的绝对值可以通过移动次数和速度分辨率的乘积来计算。通过将功率谱向右移动以求得最小CF，可以确定气流为上升气流，否则为下沉气流。
为了更直观地展示反演结果，我们从2020年6月华南的一次累计四天的持续性降水事件中典型的降水云体中选取了层云、积层混合云和对流云，并分别选取一个径向进行反演示例。在此过程中，我们采用向下为正的方向来描述大气运动，如图2所示。
层云个例选取2020年6月8日14：00：14 LCT的多普勒功率谱。如图2中的a所示，在5km海拔处有较为明显的反射率增强处标识了亮带所在高度，在亮带高度以下反射率均较低表现为明显的层云特征。2中的b为去除垂直大气运动后的降水功率谱，用黑色和红线分别标识了平均垂直大气运动和平均粒子下降速度。整体而言，平均垂直大气运动较为稳定整体呈现较弱下沉，符合层云出现在降水***消亡阶段整体呈现较弱下沉的认知，雨滴谱随高度分布的情况如图2中的c所示，整体较为稳定。
积层混合云个例选取2020年6月9日07：00：21 LCT的多普勒功率谱。如图3中的a所示，在5km海拔处有较为明显的反射率增强处标识了亮带所在高度，在亮带高度以下反射率依旧较强，既有对流云又有层状云特征。图3中的b为去除垂直大气运动后的降水功率谱，整体而言，平均垂直大气运动较为稳定整体呈现较强下沉，雨滴谱随高度分布的情况如图3中的c所示，整体较为数浓度对比层云更大。
对流个例选取2020年6月9日06：05：07 LCT的多普勒功率谱。如图4中的a所示，整体无明显亮带结构存在，对比其他个例在亮带高度下反射率无明显减弱表现为明显的对流特征。图4中的b为去除垂直大气运动后的降水功率谱。整体而言，平均垂直大气运动上层呈现上升，下层呈现下沉结合***整体的倾斜性，表现出了对流中上升下沉气流交替出现的特点，符合前人的认知。雨滴谱随高度分布的情况如图4中的c所示，整体来看大雨滴数浓度较高符合Raut的结论。
该研究采用反演算法，利用雷达频谱高度连续性，消除垂直大气运动对结果的影响。地面雨滴谱仪可测量不同大小的降水粒子在地面的分布，而空中微物理参数反演则通过雷达频谱在高度上的连续性获得不同高度处的微物理参数。为验证结论，我们对比分析了地面雨滴谱仪数据和空中150米高度数据，选取反射率因子和Nw进行对比，以验证反演方法在微物理参数反演过程中的合理性和雷达数据可用性。公式推导表明Nw将在公式(3-29)中被抵消，因此选用Nw进行对比是合理的。图3-5显示反射率因子和logNw两个参数的对比结果，反射率因子的相关性为0.95，均方根误差为2.12，logNw的相关性为0.92，均方根误差为0.073。总体而言，这些结果符合预期，表明所采用的前期的数据订正和反演算法是合理的，并验证了其可靠性。
此外，观察到logNw的反演结果在小数浓度处低估，在高数浓度处高估，这可能是因为采用的雨滴直径范围(0.25mm-6mm)不当。然而，由于在自然界中，直径大于6-8毫米的雨滴是不稳定的，具体临界值取决于气流的湍流水平[44]，因此难以给出准确的边界。
地面及空中的微物理参数反演算法，着重介绍了基于地面雨滴谱数据、风廓线雷达和C波段调频连续波雷达联合反演大气动力参数和微物理参数的算法。其中地面雨滴谱数据的置入使得算法脱离了近地面大气垂直运动速度等于零的假设，进一步提高了数据的可信度。同时风廓线雷达水平风数据仿真大气湍流功率谱，对C波段雷达的观测的降水功率谱进行反卷积剔除大气谱展宽影响。对比验证数据来源于2020年6月华南的一次累计四天的持续性降水，结果表明根据雷达资料和地面雨滴谱仪数据，对空中雨滴谱实时反演具有一定的可行性。这无疑对于雷达估测降水和模式参数化研究有一定促进作用。
以上仅是本发明众多具体应用范围中的代表性实施例，对本发明的保护范围不构成任何限制。凡采用变换或是等效替换而形成的技术方案，均落在本发明权利保护范围之内。
Claims (7)
1.一种基于空中降水粒子谱反演的方法，其特征在于，包括如下步骤：
步骤1、基于C波段调频连续波雷达获取观测功率谱，并采用去噪方法去除噪声电平，其中，观测功率谱即指多普勒功率谱；
步骤2、基于风廓线雷达获取水平风构建大气功率谱；
步骤3、基于地面雨滴谱仪获取地面雨滴数据和标准化gamma分布构建雨滴谱；
步骤4、基于步骤1-3得到的结果进行地面微物理参数反演和空中微物理参数反演。
2.根据权利要求1所述的一种基于空中降水粒子谱反演及雷达估测降水数据的方法，其特征在于，所述步骤1的具体步骤为：
步骤1.1、基于C波段调频连续波雷达获取观测功率谱，即是卷积运算后的降水粒子和大气湍流运动的多普勒功率谱，多普勒功率谱的公式为：
Sobs(v，r)＝Sdsd(v-w，r)*Sair(v，r)+Snoise(v，r)
其中，Sobs(v，r)、Sdsd(v，r)、Sdsd(v-w，r)和Snoise(v，r)分别是实测的多普勒功率谱、大气湍流运动多普勒功率谱、降水粒子多普勒功率谱以及实测的多普勒功率谱的噪声电平，v、r、w分别表示径向速度、库数和大气垂直运动速度；
步骤1.2、采用去噪方法去除观测功率谱中的噪声电平，去噪方法包括分段法、小波去噪法和局部时频分析法；
分段法是设C波段调频连续波雷达噪声服从中心χ分布，自由度为2N/k，将观测功率谱划分成k段，统计每段的平均功率谱，最小值作为噪声电平，并去除，即得到去噪后的观测功率谱；
小波去噪法是通过小波分解和滤波去除观测功率谱中包含大于给定噪声数量的小波系数，再将得到的系数重新组合成时间域信号，即得到去噪后的观测功率谱；
局部时频分析法是对观测功率谱中LT区域的特征直方图和累计分布函数进行分析，对区分噪声和信号明显的时频小区域进行噪声过滤，即得到去噪后的观测功率谱，其中，LT是指局部时频。
3.根据权利要求1所述的一种基于空中降水粒子谱反演及雷达估测降水数据的方法，其特征在于，所述步骤2的具体步骤为：
步骤2.1、风廓线雷达采用5波束探测模式得到垂直波束，倾斜角度均为15°的正北、正东、正南、正西方向的波束，其中，垂直波束包括得到垂直于地面的风速和风向分量，正北、正东、正南、正西方向的波束分别得到风速和风向在对应方向上的分量；
步骤2.2、基于获取的样本数据的平均值，基于步骤2.1得到水平风，比较相邻时间或相邻高度上的风速和风向的分量，将差异超过预定阈值的数据视为异常值，然后用平均值进行替换，即是通过判断垂直波束和4个倾斜方向波束的垂直速度差异和两组对应水平风差异是否超过给定阈值，如果差异超过预先设定的阈值，则限制或排除异常数据，否则将其保留；
步骤2.3、基于步骤2.2得到质控后的水平风构建大气功率谱。
4.根据权利要求3所述的一种基于空中降水粒子谱反演及雷达估测降水数据的方法，其特征在于，所述步骤3的具体步骤为：
步骤3.1、地面雨滴谱仪获取雨滴数据，包括地面雨滴大小、下落速度和单位采样时间内经过采样面积的雨滴个数；
步骤3.2、基于获取的地面雨滴大小对其进行修正，修正公式为：
式中：D为订正后的地面雨滴的直径，Dpar为实测的地面雨滴的直径：
步骤3.3、去除修正直径为6-8mm的雨滴后，基于标准化gamma分布构建雨滴谱，公式为：
N(D)＝N0Dμcxp(-ΛD)
其中，N0表示截距参数，μ表示形状因子，Λ表示斜率参数。
5.根据权利要求4所述的一种基于空中降水粒子谱反演及雷达估测降水数据的方法，其特征在于，所述步骤4的具体步骤为：
步骤4.1、基于雨滴数据和雨滴谱进行地面微物理参数反演，得到地面微物理参数，包括地面雨滴数量浓度、地面雨强、地面液态含水量和地面雷达反射率因子；
步骤4.2、基于多普勒功率谱、大气功率谱和雨滴谱进行空中物理参数反演，得到空中物理参数，包括空中雨滴数量浓度、空中雨强、空中液态含水量、空中雷达反射率因子、垂直大气运动和空中雨滴谱。
6.根据权利要求5所述的一种基于空中降水粒子谱反演及雷达估测降水数据的方法，其特征在于，所述步骤4.1的具体步骤为：
步骤4.11、基于地面雨滴下落速度和单位采样时间内经过采样面积的雨滴个数得到雨滴数浓度，公式为：
式中：nij为修正后的雨滴直径位于第i个尺度通道且下落速度位于第j个速度通道的雨滴数，A为仪器采样面积，Δt为采样时间间隔，Di和ΔDi为修正后的地面雨滴在第i个尺度通道的中心直径和该区间的尺度间隔，Vj为第j个速度通道的中心速度；
步骤4.12、基于雨滴数浓度计算得到地面雨强R1、地面液态含水量LWC1和地面雷达反射率因子Z1，计算公式为：
式中：ρ为水的密度；
步骤4.13、基于雨滴谱采用阶矩法对标准化gamma分布参数进行拟合，拟合时引入雨滴谱的阶距Mn，其中，标准化gamma分布的参数N0、μ、人用雨滴谱的2、4、6阶距估算得到，公式为：
其中，n表示第n阶矩；
步骤4.14、基于阶距Mn得到地面雨滴的质量加权平均直径D1m，公式为：
步骤4.15、基于质量加权平均直径Dlm，得到归一化的数量浓度N1w，计算公式为：
其中，w表示大气垂直运动速度。
7.根据权利要求6所述的一种基于空中降水粒子谱反演及雷达估测降水数据的方法，其特征在于，所述步骤4.2的具体步骤为：
步骤4.21、基于大气功率谱，采用傅里叶反卷积算法消除去噪后的多普勒功率谱中的大气谱展宽；
步骤4.22、基于雨滴谱构建初始高度降水粒子功率谱，即在静稳大气条件下，降水粒子多普勒功率谱Sdsd(v，r)用雨滴谱分布Ndsd(D)和后向散射截面表示为：
其中，Sdsd(v，r)表示降水粒子多普勒功率谱，即初始高度降水粒子功率谱，λ为雷达波长，为5.42cm，|Kw|2为介电常数，表示下落末速度与雨滴直径关系，D为修正后的雨滴；
当波长为5.42cm时，此时降水粒子发生瑞利散射，若将其视为球形粒子，此时，雷达后向散射截面与粒子直径的6次方成正比，后向散射截面可以表示为：
其中，Γ表示欧拉函数，μ表示形状因子，人表示斜率参数，其中的积分从0延伸到无穷大，反演过程中雨滴谱受到最小和最大雨滴直径的限制，Dm表示地面雨滴的质量加权平均直径；
空中雨强R和空中液态水含水量LWC均是空中雨滴大小分布的三阶矩：
其中，ΔD表示空中雨滴直径的变化率，V(D)表示空中雨滴下降速度与直径之间关系，ρ0代表地面的大气密度，ρh表示地面的大气密度，h为观测高度；
步骤4.23、基于步骤4.21得到的结果获取初始高度大气运动，并剔除初始高度降水粒子功率谱中相邻的分布差异；
步骤4.24、空中雷达反射率因子Z是雨滴谱的零阶矩，公式为：
代入：
N(D；Nw，Dm，μ)＝Nwf(D；Dm，μ)
得到表示纯降水粒子在瑞利散射时的反射率，即空中雷达反射率因子Z为：
其中，Δv分别表示空中雨滴下落速度的变化率；
步骤4.25、基于步骤4.23得到的结果，对每一距离库的雨滴谱逐一进行归一化处理，并通过最小化代价函数得到垂直大气运动；
归一化公式为：
最小化代价函数为：
步骤4.26、基于步骤4.23和步骤4.25得到的结果，反卷积得到空中雨滴谱。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310816919.7A CN116879899A (zh) | 2023-07-05 | 2023-07-05 | 一种基于空中降水粒子谱反演的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310816919.7A CN116879899A (zh) | 2023-07-05 | 2023-07-05 | 一种基于空中降水粒子谱反演的方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116879899A true CN116879899A (zh) | 2023-10-13 |
Family
ID=88270807
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310816919.7A Pending CN116879899A (zh) | 2023-07-05 | 2023-07-05 | 一种基于空中降水粒子谱反演的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116879899A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117706512A (zh) * | 2023-12-13 | 2024-03-15 | 安徽省气象台 | 一种融合温度判定与背景依赖的水凝物反演方法和*** |
-
2023
- 2023-07-05 CN CN202310816919.7A patent/CN116879899A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117706512A (zh) * | 2023-12-13 | 2024-03-15 | 安徽省气象台 | 一种融合温度判定与背景依赖的水凝物反演方法和*** |
| CN117706512B (zh) * | 2023-12-13 | 2024-05-10 | 安徽省气象台 | 一种融合温度判定与背景依赖的水凝物反演方法和*** |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Bianco et al. | Convective boundary layer depth: Improved measurement by Doppler radar wind profiler using fuzzy logic methods | |
| US8098188B2 (en) | Method of characterizing the convection intensity of a cloud, by a meteorological radar | |
| Collis et al. | Statistics of storm updraft velocities from TWP-ICE including verification with profiling measurements | |
| Protat et al. | The accuracy of radar estimates of ice terminal fall speed from vertically pointing Doppler radar measurements | |
| KR101531224B1 (ko) | 이중편파 레이더 기반의 강수 추정 시스템 및 그 방법 | |
| CN115526272A (zh) | 基于多源资料融合的低云低能见度天气识别及反演方法 | |
| Giangrande et al. | Automated retrievals of precipitation parameters using non-Rayleigh scattering at 95 GHz | |
| Gage et al. | Use of two profilers during MCTEX for unambiguous identification of Bragg scattering and Rayleigh scattering | |
| Szyrmer et al. | Snow studies. Part IV: Ensemble retrieval of snow microphysics from dual-wavelength vertically pointing radars | |
| CN116879899A (zh) | 一种基于空中降水粒子谱反演的方法 | |
| CN115825991A (zh) | 基于测风激光雷达的风切变识别方法及设备、储存介质 | |
| Frehlich | Scanning doppler lidar for input into short-term wind power forecasts | |
| Koistinen et al. | Operational measurement of precipitation in cold climates | |
| Lin et al. | Correction to raindrop size distributions measured by PARSIVEL disdrometers in strong winds | |
| Waitz et al. | PHIPS-HALO: the airborne Particle Habit Imaging and Polar Scattering probe–Part 3: Single-particle phase discrimination and particle size distribution based on the angular-scattering function | |
| Williams et al. | Cluster analysis techniques to separate air motion and hydrometeors in vertical incident profiler observations | |
| Lin et al. | Evaluation of the raindrop size distribution representation of microphysics schemes in Typhoon Lekima using disdrometer network observations | |
| Melchionna et al. | A new algorithm for the extraction of cloud parameters using multipeak analysis of cloud radar data. First application and results | |
| Ruan et al. | Analysis of the microphysical properties of a stratiform rain event using an L-Band profiler radar | |
| Bonin et al. | Detection of range-folded returns in Doppler lidar observations | |
| Pang et al. | Estimating raindrop size distributions and vertical air motions with spectral difference using vertically pointing radar | |
| Kumar et al. | The estimation of convective mass flux from radar reflectivities | |
| Deshpande et al. | UHF wind profiler observations during a tropical pre-monsoon thunderstorm—A case study | |
| Cohn et al. | Radial velocity and wind measurement with NIMA–NWCA: Comparisons with human estimation and aircraft measurements | |
| Orikasa et al. | Ice crystal concentration in midlatitude cirrus clouds: In situ measurements with the balloonborne hydrometeor videosonde (HYVIS) |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |