多普勒雷达风速资料的填补方法技术

本发明专利技术公开了一种多普勒雷达风速资料的填补方法，其方法包括以下步骤：第一步：多普勒雷达以低仰角α扫描收集覆盖雷达站探测范围的大面积降水过程或睛空回波的雷达样本数据，所述雷达样本数据为极坐标下某一距离圈上随方位角分布的径向速度Vr0(θ)，该径向速度Vr0(θ)由回波区A的径向速度V/r0(θ)和缺测的无回波区B的径向速度V//r0(θ)所组成；第二步：利用迭代法根据所述回波区A的径向速度V/r0(θ)对所述某一距离圈上的缺测无回波区B的多普勒速度资料进行填补。本发明专利技术对于不大于90°连续累积的多普勒速度资料缺测区和小于120°非连续累积缺测区都能进行有效填补，填补后的速度图比填补前更能反应出大范围的速度场特征。

【技术实现步骤摘要】

本专利技术属于气象雷达资料数据处理领域，尤其涉及一种非线性风场下多普勒雷达风速资料的填补方法。
技术介绍
多普勒天气雷达所探测到的风场信息缺测问题，一直是影响多普勒天气雷达风场信息分析应用和后续速度产品开发的一个关键，目前尚无有效技术解决这个问题。对于风场信息缺失的填补，目前最常用的方法是插值法，如台北台湾大学2002出版的汪家慧《应用单多普勒雷达观测进行台风涡旋初始化》中用高仰角的风场数据对相应低仰角的风场数据进行填补，但这只适用于高低层风场较为一致的情况；应用气象学报2002,13(5):591-599出版的梁海河,张沛源,葛润生《多普勒天气雷达风场退模糊方法的研究》中提出K-邻域频数法来消除风场中的噪声和填补缺测，仅针对缺测点；南京气象学院2003出版的沃伟峰《多普勒天气雷达二次产品开发平台的设计和实现 》则通过与缺测区域相邻的数据进行线性插值来填补缺测，但只适用于连续缺测方位角小于8°的情况；气象2010,36(5):1-12出版的邓勇、尹丽云、许迎杰等利用VAD方法和迭代填补技术进行缺测区域迭代填补试验，这些反演方法都是建立在风场在一定空间呈线性分布的假设上，总而言之，目前还没有能有效填补非线性风场下多普勒风速资料的方法。
技术实现思路
本专利技术的目的是提供一种在非线性风场中能有效地填补多普勒风速资料的方法。为达到上述目的，本专利技术提供了一种多普勒雷达风速资料的填补方法，包括步骤：第一步：多普勒雷达以低仰角α扫描收集覆盖雷达站探测范围的大面积降水过程或睛空回波的雷达样本数据，所述雷达样本数据为极坐标下某一距离圈上随方位角分布的径向速度Vr0(θ)，该径向速度Vr0(θ)由回波区A的径向速度V/r0(θ)和缺测的无回波区B的径向速度V//r0(θ)所组成；第二步：利用迭代法根据所述回波区A的径向速度V/r0(θ)对所述某一距离圈上的缺测无回波区B的多普勒速度资料进行填补；其特征在于：使用迭代法对所述某一距离圈上的缺测无回波区B的径向速度进行填补包括如下步骤：S1：根据VAD理论推导出用傅立叶系数来计算极坐标系下的非线性径向风速Vr(θ)的公式（1）； （1）S2：将公式（1）进行三阶傅里叶系数展开得到公式（2）， （2），其中，ap为傅立叶系数，fp(θ)为基函数，f0(θ)=1，f1(θ)=sin(θ)，f2(θ)= cos(θ)，f3(θ)= sin(2θ)，f4(θ)= cos(2θ)，f5(θ)= sin(3θ)，f6(θ)= cos(3θ)；S3：将回波区A的径向速度V/r0(θ)分别代入公式（2）形成一组径向速度V/r0(θ)的傅里叶展开方程组（3）， （3）M表示回波区A径向速度V/r0(θ)样本数据的总个数；S4：根据最小二乘法的性质，对回波区A原雷达探测数据V/r0(θ)与经VAD反演得到的径向速度的误差函数求极值，获得求解傅立叶系数P＊=(a0，a1，a2，a3，a4，a5，a6)最优解的迭代公式（11）： （11）其中，μk为阻尼因子，I为单位矩阵，J（PK）为傅立叶系数Ｐｋ(a0，a1，a2，a3，a4，a5，a6)的雅克比矩阵，JT（PK）为J（PK）的转置矩阵；S5：根据迭代公式（1１）逐步迭代求解即可得出傅立叶系数P＊(a0，a1，a2，a3，a4，a5，a6)的最优解；S6：将步骤S５求得的傅立叶系数P＊=(a0，a1，a2，a3，a4，a5，a6) 最优解代入公式（1）中，获得一组新的覆盖有回波区A和无回波区B的径向速度Vr1(θ)， 所述Vr1(θ)中对应于 Vr0(θ)中有回波A的部分为V/r1(θ)，对应于Vr0(θ)中无回波B的部分为V//r1(θ)，用V//r1(θ)去填补无回波区B中的V//r0(θ) ，在有回波区A中仍保留原来的雷达探测值V/r0(θ)，由此得到已消除原无回波区B资料空缺的新的径向速度Vr随方位角θ的分布，定义为Vr2(θ)；S7：计算有回波区A中原雷达探测数据V/r0(θ)与经VAD反演及迭代求解得到的数据V/r1(θ)之间的均方根误差RMS；S8：重复步骤S１至S７,当所述均方根误差RMS的值或连续多组均方根误差RMS的滑动平均值小到符合预设的约束条件时迭代填补过程结束。所述步骤S４推导迭代公式（1１）包括如下步骤：S4-1：定义实测回波与反演回波数值误差，记为 (４)其中，j=1，2，…，M；S4-2：假设误差目标函数: （5）其中P表示傅立叶系数，P=(a0，a1，a2，a3，a4，a5，a6);S4-3：为使误差目标函数最小，即要求出最优解的傅立叶系数P＊使得Φ(P)最小，根据多项式极值必要条件: （6）即 （7）其中J(P＊)为如下雅克比矩阵：（8）。S4-4：根据泰勒公式把JT(P)和dif(P)在初值P=P＊处展开，将其转化为线性方程组并以适当步长Sk（其中迭代次数k=0，1，…）逐步迭代求解，经正规化处理后，迭代方程和法方程分别如下： （9） （10）其中，μk为阻尼因子，I为单位矩阵，J（PK）为傅立叶系数Ｐｋ(a0，a1，a2，a3，a4，a5，a6)的雅克比矩阵，JT（PK）为J（PK）的转置矩阵；S4-5：把公式（10）代入（9）即可获得求解最优解的傅立叶系数P＊=(a0，a1，a2，a3，a4，a5，a6)的迭代公式（11）。所述步骤S4－4的步长Sk具体为１度。所述步骤S7中均方根误差RMS的计算公式如下：所述连续多组均方根误差RMS的滑动平均值是指连续三次或三次以上迭代的均方根误差RMS的滑动平均值。所述多普勒雷达低仰角α扫描是指雷达仰角α小于10度。在所述第一步收集雷达样本数据时,若所述径向速度Vr0(θ)存在速度模糊现象，则需进行退速度模糊处理。与现有技术相比较，本专利技术具有如下优点：本专利技术对于不大于90°连续累积的多普勒速度资料缺测区和小于120°非连续累积缺测区都能进行有效填补，填补后的多普勒径向速度资料所产生的图像与实测图像能很好的融为一体，填补效果极佳。附图说明参照附图，本专利技术将得到更好的理解。附图如下：图1是本专利技术多普勒雷达风速资料填补方法流程图；图2是本专利技术多普勒雷达风速资料填补方法第二步使用迭代法对缺测无回波区B的径向速度进行填补的流程图：图3是本专利技术多普勒雷达风速资料填补方法的流程图；图4是利用本专利技术多普勒雷达风速资料填补方法填补速度资料后的效果示意图；图5是利用本专利技术多普勒雷达风速资料填补方法对非连续累积缺口为60°时的雷达风速场资料填补前后效果对比示意图；图6是利用本专利技术多普勒雷达风速资料填补方法对非连续累积缺口为120°时的雷达风速场资料填补前后效果对比示意图；图7是利用本专利技术多普勒雷达风速资料填补方法对连续累积缺口为30°时的雷达风速场资料填补前后效果对比示意图本文档来自技高网...

【技术保护点】
一种多普勒雷达风速资料的填补方法，包括步骤：第一步：多普勒雷达以低仰角α扫描收集覆盖雷达站探测范围的大面积降水过程或睛空回波的雷达样本数据，所述雷达样本数据为极坐标下某一距离圈上随方位角分布的径向速度Vr0(θ)，该径向速度Vr0(θ)由回波区A的径向速度V/r0(θ)和缺测的无回波区B的径向速度V//r0(θ)所组成；第二步：利用迭代法根据所述回波区A的径向速度V/r0(θ)对所述某一距离圈上的缺测无回波区B的多普勒速度资料进行填补；其特征在于：使用迭代法对所述某一距离圈上的缺测无回波区B的径向速度进行填补包括如下步骤：S1：根据VAD理论推导出用傅立叶系数来计算极坐标系下的非线性径向风速Vr(θ)的公式（1）；（1）S2：将公式（1）进行三阶傅里叶系数展开得到公式（2），（2），其中，ap为傅立叶系数，fp(θ)为基函数，f0(θ)=1，f1(θ)=sin(θ)，f2(θ)= cos(θ)，f3(θ)= sin(2θ)，f4(θ)= cos(2θ)，f5(θ)= sin(3θ)，f6(θ)= cos(3θ)；S3：将回波区A的径向速度V/r0(θ)分别代入公式（2）形成一组径向速度V/r0(θ)的傅里叶展开方程组（3），（3）M表示回波区A径向速度V/r0(θ)样本数据的总个数；S4：根据最小二乘法的性质，对回波区A原雷达探测数据V/r0(θ)与经VAD反演得到的径向速度的误差函数求极值，获得求解傅立叶系数P＊=(a0，a1，a2，a3，a4，a5，a6)最优解的迭代公式（11）：（11）其中，μk为阻尼因子，I为单位矩阵，Ｊ（Ｐｋ）为傅立叶系数Ｐｋ(a0，a1，a2，a3，a4，a5，a6)的雅克比矩阵，JT（Ｐｋ）为Ｊ（Ｐｋ）的转置矩阵；S5：根据迭代公式（11）逐步迭代求解即可得出傅立叶系数P＊(a0，a1，a2，a3，a4，a5，a6)的最优解；S6：将步骤S５求得的傅立叶系数P＊=(a0，a1，a2，a3，a4，a5，a6) 最优解代入公式（1）中，获得一组新的覆盖有回波区A和无回波区B的径向速度Vr1(θ)， 所述Vr1(θ)中对应于 Vr0(θ)中有回波A的部分为V/r1(θ)，对应于Vr0(θ)中无回波B的部分为V//r1(θ)，用V//r1(θ)去填补无回波区B中的V//r0(θ) ，在有回波区A中仍保留原来的雷达探测值V/r0(θ)，由此得到已消除原无回波区B资料空缺的新的径向速度Vr随方位角θ的分布，定义为Vr2(θ)；S7：计算有回波区A中原雷达探测数据V/r0(θ)与经VAD反演及迭代求解得到的数据V/r1(θ)之间的均方根误差RMS；S8：重复步骤S１至S７,当所述均方根误差RMS的值或连续多组均方根误差RMS的滑动平均值小到符合预设的约束条件时迭代填补过程结束。...

【技术特征摘要】
1.一种多普勒雷达风速资料的填补方法，包括步骤：第一步：多普勒雷达以低仰角α扫描收集覆盖雷达站探测范围的大面积降水过程或睛空回波的雷达样本数据，所述雷达样本数据为极坐标下某一距离圈上随方位角分布的径向速度Vr0(θ)，该径向速度Vr0(θ)由回波区A的径向速度V/r0(θ)和缺测的无回波区B的径向速度V//r0(θ)所组成；第二步：利用迭代法根据所述回波区A的径向速度V/r0(θ)对所述某一距离圈上的缺测无回波区B的多普勒速度资料进行填补；其特征在于：使用迭代法对所述某一距离圈上的缺测无回波区B的径向速度进行填补包括如下步骤：S1：根据VAD理论推导出用傅立叶系数来计算极坐标系下的非线性径向风速Vr(θ)的公式（1）； （1）S2：将公式（1）进行三阶傅里叶系数展开得到公式（2）， （2），其中，ap为傅立叶系数，fp(θ)为基函数，f0(θ)=1，f1(θ)=sin(θ)，f2(θ)= cos(θ)，f3(θ)= sin(2θ)，f4(θ)= cos(2θ)，f5(θ)= sin(3θ)，f6(θ)= cos(3θ)；S3：将回波区A的径向速度V/r0(θ)分别代入公式（2）形成一组径向速度V/r0(θ)的傅里叶展开方程组（3）， （3）M表示回波区A径向速度V/r0(θ)样本数据的总个数；S4：根据最小二乘法的性质，对回波区A原雷达探测数据V/r0(θ)与经VAD反演得到的径向速度的误差函数求极值，获得求解傅立叶系数P＊=(a0，a1，a2，a3，a4，a5，a6)最优解的迭代公式（11）： （11）其中，μk为阻尼因子，I为单位矩阵，Ｊ（Ｐｋ）为傅立叶系数Ｐｋ(a0，a1，a2，a3，a4，a5，a6)的雅克比矩阵，JT（Ｐｋ）为Ｊ（Ｐｋ）的转置矩阵；S5：根据迭代公式（11）逐步迭代求解即可得出傅立叶系数P＊(a0，a1，a2，a3，a4，a5，a6)的最优解；S6：将步骤S５求得的傅立叶系数P＊=(a0，a1，a2，a3，a4，a5，a6) 最优解代入公式（1）中，获得一组新的覆盖有回波区A和无回波区B的径向速度Vr1(θ)， 所述Vr1(θ)中对应于 Vr0(θ)中有回波A的部分为V/r1(θ)，对应于Vr0(θ)中无回波B的部分为V//r1(θ)，用V//r1(θ)去填补无回波区B中的V//r0(θ) ，在有回波区A中仍保留原来的雷达探测值V/r0(θ...

【专利技术属性】
技术研发人员：汤达章，黄铃光，尹丽云，
申请(专利权)人：福建省气象台，
类型：发明
国别省市：福建;35