天气雷达数据坐标转换快速插值方法与流程

本发明涉及天气雷达信号与数据处理领域，尤其涉及一种用于天气雷达数据坐标转换的快速插值方法。
背景技术：
：天气雷达得到的数据产品通常处于以雷达自身为中心的极坐标，包含气象目标的距离，方位和俯仰坐标。对于多部天气雷达组网观测的应用，需要构建统一的坐标系，例如由经度、维度和高度组成的地理坐标系作为基准，各雷达得到的数据产品均通过插值方法由自身极坐标系变换至新的坐标系下，从而方便天气雷达数据的组网拼图。由于大部分气象目标，例如云在空间上具有连续性，并且有很多细尺度结构，因而希望插值后的散射率场在空间上要保持连续性，同时在内插过程中应最大限度地保留雷达资料中存在的原始回波结构特征。常见的插值方法包括：最近邻居法、线性内插法和sinc核插值法。其中，最近邻居法虽然速度快，但是精度过低，易导致插值后的雷达散射率在空间上不连续；线性内插法对原始插值前数据引入较大的平滑，因而易损失原有散射率细尺度结构特征。sinc核插值法是根据奈奎斯特采样定理进行插值的方法，精度最高，但是通常情况下计算复杂度较高。技术实现要素：鉴于上述技术问题，本发明提供了一种用于天气雷达数据坐标转换的快速插值方法。本发明中的天气雷达数据坐标转换快速插值方法，包括：步骤1，根据预设的插值核点数和量化位移预先计算sinc插值核表格；步骤2，通过映射关系计算出待转换到的地理坐标系的均匀三维网格中的每个网络点的网格点坐标值所对应的以当前雷达自身为中心的极坐标系下的坐标；步骤3，对于步骤2计算得到的某一个具体的距离、俯仰和方位坐标，计算其在天气雷达体积扫描数据的均匀极坐标离散网格中所处的位置；步骤4，根据步骤3得到的所述位置的位置值的整数部分，从给定的离散化保存的天气雷达体积扫描数据中抽取一个三维矩阵数据块；步骤5，分别根据步骤3得到的所述位置的位置值的小数部分，从sinc插值核表格中查询并得到表格中的三行元素以分别组成列向量和步骤6，利用步骤4得到的三维矩阵数据块和步骤5得到的列向量进行加权运算，得到二维矩阵数据块；步骤7，利用步骤6得到的二维矩阵数据块和步骤5得到的列向量进行加权运算，得到一维列向量；步骤8，利用步骤7得到的一维列向量，与步骤5得到的列向量进行加权运算，得到当前网格值的插值结果。优选地，所述步骤1中，所述的插值核点数p的典型取值是6～16间的任一偶数，量化位移1/l中的l的典型取值为大于等于10的偶数；所述的sinc插值核表格是l+1行，p列的数值表格；所述sinc插值核表格的第i行，第j列的元素wi,j的值通过以下公式计算得到其中，i＝1,2,...l+1，j＝1,2,...p，sinc(x)＝sin(πx)/(πx)表示sinc函数。优选地，所述步骤2中，由地理坐标系的网格点坐标值(xlat,m′,ylon,n′,hk′)到以当前雷达自身为中心的极坐标系坐标的映射关系的表达式为：其中，(xr,yr,hr)为天气雷达自身位置的经纬高坐标，r表示地球半径，s＝r×arccos[sin(xlat,m′)sin(xr)+cos(xlat,m′)cos(xr)cos(ylon,n′-yr)]xlat,m′,ylon,n′,hk′分别表示网格所代表的第m′个纬度、第n′个经度和第k′个高度坐标。优选地，所述步骤3中，由某一个具体的距离、俯仰和方位坐标计算其在天气雷达体积扫描数据的均匀极坐标离散网格中所处的位置(x1,x2,x3)的方法为优选地，所述步骤4包括：步骤41，分别对x1,x2,x3取整，得到其中表示取整算子；步骤s42，在给定的离散化保存的天气雷达体积扫描数据中，其中m＝1,2,...m，n＝1,2,...n，k＝1,2,...k，按第一个维度的索引从n1-(p/2-1)、n1-p/2、n1-p/2+1、…、n1+p/2，第二个维度的索引从n2-(p/2-1)、n2-p/2、n2-p/2+1、…、n2+p/2，第三个维度的索引从n3-(p/2-1)、n3-p/2、n3-p/2+1、…、n3+p/2，抽取出一个p×p×p维的三维矩阵数据块s(i,j,l)，i,j,l＝1,2,...p。优选地，所述步骤5包括：计算x1的小数部分将其除以量化位移1/l的值并四舍五入为整数其中表示取整算子，查询sinc插值表格并选定插值表中第l+1-m1行的元素作为加权值组成行向量round()表示四舍五入取整算子；计算x2的小数部分将其除以量化位移1/l的值并四舍五入为整数查询sinc插值表格并选定插值表中第l+1-m2行的元素作为加权值组成行向量计算x3的小数部分将其除以量化位移1/l的值并四舍五入为整数查询sinc插值表格并选定插值表中第l+1-m3行的元素作为加权值组成行向量优选地，所述步骤6包括：在三维矩阵数据块sp×p×p中，针对某个固定的索引对(j,l)，将三维矩阵数据块sp×p×p的p个数据元素s(1,j,l)，s(2,j,l)，...,s(p,j,l)组成列向量然后求取该向量与列向量的内积作为一个p×p二维矩阵数据块s′p×p×p的第j行，第l列的数据元素，其中，上标t表示矩阵或向量转置。对所有的索引对(j,l)，j,l＝1,2,...p进行遍历，直到计算得到p×p二维矩阵数据块s′p×p×p。优选地，所述步骤7中的一维列向量的表达式为：其中，为一维列向量，s′p×p×p为二维矩阵数据块，为列向量，t表示矩阵或向量转置。优选地，所述步骤8中，当前网格值的插值结果为：其中，t表示矩阵或向量转置。从上述技术方案可以看出，本发明具有以下有益效果：(1)本发明技术方案中天气雷达数据坐标变换使用sinc核插值，插值精度较高，有助于保留气象目标散射率的细尺度结构特征；(2)本发明技术方案中插值的实现是通过查询预先计算好的sinc插值核表格获得加权值，避免了实时计算复杂度较高的sinc函数的过程，提升了插值速度。附图说明图1为本发明实施例天气雷达进行体积扫描示意图；图2为本发明实施例中使用的sinc核插值原理示意图。图3为本发明实施例一种用于天气雷达数据坐标转换的快速插值方法的流程图。具体实施方式为使本领域技术人员更好的理解本发明的技术方案，下面结合附图和具体实施方式对本发明作详细说明。本发明技术方案中天气雷达数据坐标变换使用sinc核插值，插值精度较高，有助于保留气象目标散射率的细尺度结构特征；插值的实现是通过查询预先计算好的sinc插值核表格获得加权值，避免了实时计算复杂度较高的sinc函数的过程，提升了插值速度。本发明的实施方式是针对天气雷达数据由自身极坐标数据变换到其它坐标系如由经纬高表示的地理坐标系进行的。在详细介绍本发明的实施方式的细节之前，首先简要介绍利用数据插值实现坐标变换的原理。请参考图1所示，通常，天气雷达进行体积扫描，可获取气象目标在极坐标系下的距离r、方位角俯仰角θ坐标的雷达散射率因子的基本数据产品，设其表示为设数据实施坐标变换到地理坐标系后的数据应当为z′(xlat,ylon,h)，其中xlat表示经度坐标，ylon表示维度坐标，h表示高度坐标。对于某个具体的坐标值(xlat,ylon,h)处的气象目标，坐标变换后的数据z′(xlat,ylon,h)是由通过坐标映射得到的：其中，r(xlat,ylon,h),θ(xlat,ylon,h),可通过大圆几何学理论求得，设天气雷达自身位置的经纬高坐标为(xr,yr,hr)，则其中，r表示地球半径，s的表达式为s＝r×arccos[sin(xlat)sin(xr)+cos(xlat)cos(xr)cos(ylat-yr)]通常情况下，天气雷达体扫描获取的数据是在极坐标系下均匀网格分布的，离散化保存的数据的，m＝1,2,...m，n＝1,2,...n，k＝1,2,...k。m，n和k分别表示插值前极坐标系上天气雷达数据在距离、俯仰和方位方向的数据维度。某个具体的地理坐标系坐标值(xlat,ylon,h)所对应映射的极坐标r(xlat,ylon,h),θ(xlat,ylon,h),不一定正好处于极坐标系上的整数网格点上，需要利用其周围网格点上的数据通过插值方法得到变换后的数据z′(xlat,ylon,h)。常见的插值方法包括：最近邻居法、线性内插法和sinc核插值法。其中，最近邻居法虽然速度快，但是精度过低，易导致插值后的雷达散射率在空间上不连续；线性内插法对原始插值前数据引入较大的平滑，因而易损失原有散射率细尺度结构特征。sinc核插值法是根据奈奎斯特采样定理进行插值的方法，精度最高，但是通常情况下计算复杂度较高。本发明的实施方式中插值是基于sinc核插值进行的。在详细介绍本发明的实施方式的细节之前，还需先简单描述sinc核插值的一些概念和原理。通常，在满足信号带限和采样率大于奈奎斯特采样频率的前提下，离散序列s[i]，i＝1,2,...可以通过sinc核插值得到任意非整数采样点x上的高精度的插值结果s′(x)，表示如下：其中，sinc(x)＝sin(πx)/(πx)称为sinc插值核，上述插值公式表示，任意非整数采样点x上信号值可以通过对离散序列s[i]的加权和得到，而各离散序列样点的权值是sinc插值核。图2为本发明实施例中使用的sinc核插值原理示意图，一般情况下，无需对所有样点均计算sinc核权值，而是使用非整数采样点x周围的8～16点的离散序列s[i]的样本值参与插值运算，例如图2中，需要插值计算非整数样点x处的值，则使用其左右各4个整数样点值(x处左右各4个“o”符号表示的样点值)和对应的sinc函数权值(用“□”表示的sinc函数样点值)进行加权求和。进行sinc核插值公式计算时，由于根据输入的x的值计算插值核sinc(x)的值的过程涉及到三角函数运算，计算复杂度较高，为克服这个缺点，本发明专利将预先计算好的插值核以表格方式存储以避免实时计算，从而提升插值速度。图1是根据本发明实施例的一种用于天气雷达数据坐标转换的快速插值方法的流程图。请参考图1，在本发明的一个实施例中，提供了一种用于天气雷达数据坐标转换的快速插值方法，该方法可以包括：步骤s1，根据预设的插值核点数p和量化位移1/l预先计算sinc插值核表格，如表1所示；所述的插值核点数p的典型取值是6～16间的任一偶数，量化位移1/l中的l的典型取值为大于等于10的偶数；所述的sinc插值核表格是l+1行，p列的数值表格；所述sinc插值核表格的第i行，第j列的元素wi,j的值通过以下公式计算得到其中，i＝1,2,...l+1，j＝1,2,...p，sinc(x)＝sin(πx)/(πx)表示sinc函数。表1为在p＝8，l＝12时的13行8列的插值核表格00001000-0.0210.028-0.0430.0900.990-0.0760.040-0.027-0.0420.057-0.0870.1920.961-0.1370.074-0.051-0.0610.083-0.1300.3040.912-0.1820.101-0.070-0.0770.105-0.1690.4220.843-0.2110.120-0.084-0.0880.122-0.1990.5400.756-0.2220.130-0.092-0.0930.131-0.2180.6530.653-0.2180.131-0.093-0.0920.130-0.2220.7560.540-0.1990.122-0.088-0.0840.120-0.2110.8430.422-0.1690.105-0.077-0.0717-10590.101-0.1820.9120.304-0.1300.083-0.061-0.0510.074-0.1370.9610.192-0.0870.057-0.042-0.0270.040-0.0760.9900.090-0.0430.028-0.02100010000步骤s2，设要转换到的地理坐标系的均匀三维网格为ω＝{(xlat,m′,ylon,n′,hk′)|m′＝1,2,...m′，n′＝1,2,...n′，k′＝1,2,...k′}，其中xlat,m′,ylon,n′,hk′分别表示网格所代表的第m′个纬度、第n′个经度和第k′个高度坐标，m，n和k分别表示网格的尺寸维度，对于所有网格点中每一个具体的网格点坐标值(xlat,m′,ylon,n′,hk′)，首先通过映射关系f计算出其对应的以当前雷达自身为中心的极坐标系下的坐标然后都重复执行下述步骤s3～s8，其中r′,θ′,分别表示映射的极坐标系下的距离、俯仰角和方位角坐标；所述的由地理坐标系的网格点坐标值(xlat,m′,ylon,n′,hk′)到以当前雷达自身为中心的极坐标系坐标的映射关系f的表达式为：其中，(xr,yr,hr)为天气雷达自身位置的经纬高坐标，r表示地球半径，s的表达式为s＝r×arccos[sin(xlat,m′)sin(xr)+cos(xlat,m′)cos(xr)cos(ylon,n′-yr)]步骤s3，对于步骤s2计算得到的某一个具体的距离、俯仰和方位坐标计算其在天气雷达体积扫描数据的均匀极坐标离散网格中所处的位置(x1,x2,x3)，其表示距离坐标r′位于rm序列中的第x1个样点，m＝1,2,...m；俯仰坐标θ′位于θn序列中的第x2个样点，n＝1,2,...n；方位坐标位于序列中的第x3个样点，k＝1,2,...k；x1,x2,x3为不小于1的整数或非整数；所述的由某一个具体的距离、俯仰和方位坐标计算其在天气雷达体积扫描数据的均匀极坐标离散网格中所处的位置(x1,x2,x3)的方法为：步骤s4，根据步骤s3得到的x1,x2,x3位置值的整数部分，从给定的离散化保存的天气雷达体积扫描数据m＝1,2,...m，n＝1,2,...n，k＝1,2,...k中，抽取出一个p×p×p维的三维矩阵数据块s(i,j,l)，i,j,l＝1,2,...p,其中整数p为所述步骤s1所示的预设的插值核点数p；所述的根据x1,x2,x3位置值的整数部分，从给定的离散化保存的天气雷达体积扫描数据m＝1,2,...m，n＝1,2,...n，k＝1,2,...k中，抽取出一个p×p×p维的三维矩阵数据块s(i,j,l)，i,j,l＝1,2,...p的步骤又可以包含：子步骤s41:分别对x1,x2,x3取整，得到其中表示取整算子；子步骤s42:对于给定的离散化保存的天气雷达体积扫描数据m＝1,2,...m，n＝1,2,...n，k＝1,2,...k中，按第一个维度的索引从n1-(p/2-1)、n1-p/2、n1-p/2+1、…、n1+p/2，第二个维度的索引从n2-(p/2-1)、n2-p/2、n2-p/2+1、…、n2+p/2，第三个维度的索引从n3-(p/2-1)、n3-p/2、n3-p/2+1、…、n3+p/2，抽取出一个p×p×p维的三维矩阵数据块s(i,j,l)，i,j,l＝1,2,...p；步骤s5：分别根据步骤s3得到的x1,x2,x3位置值的小数部分，从sinc插值核表格中查询并得到表格中的三行元素，分别组成列向量和所述的根据x1,x2,x3位置值的小数部分，从sinc插值核表格中查询并得到表格中的三行元素，分别组成列向量和的方法是：计算x1的小数部分将其除以量化位移1/l的值并四舍五入为整数其中表示取整算子，查询sinc插值表格并选定插值表中第l+1-m1行的元素作为加权值组成行向量round()表示四舍五入取整算子；计算x2的小数部分将其除以量化位移1/l的值并四舍五入为整数查询sinc插值表格并选定插值表中第l+1-m2行的元素作为加权值组成行向量计算x3的小数部分将其除以量化位移1/l的值并四舍五入为整数查询sinc插值表格并选定插值表中第l+1-m3行的元素作为加权值组成行向量步骤s6，利用步骤s4得到的三维矩阵数据块sp×p×p和步骤s5得到的列向量进行加权运算，得到p×p二维矩阵数据块s′p×p×p；所述的利用三维矩阵数据块sp×p×p和列向量进行加权运算的过程为：在三维矩阵数据块sp×p×p中，针对某个固定的索引对(j,l)，将三维矩阵数据块sp×p×p的p个数据元素s(1,j,l)，s(2,j,l)，...,s(p,j,l)组成列向量然后求取该向量与列向量的内积作为一个p×p二维矩阵数据块s′p×p×p的第j行，第l列的数据元素，其中，上标t表示矩阵或向量转置。对所有的索引对(j,l)，j,l＝1,2,...p进行遍历，直到计算得到p×p二维矩阵数据块s′p×p×p；步骤s7，利用步骤s6得到的二维矩阵数据块s′p×p×p和步骤s5得到的列向量进行加权运算，得到包含p个元素的一维列向量所述的利用二维矩阵数据块s′p×p×p和列向量进行加权运算，利用的是向量和矩阵的乘法，其表达式为：得到一个包含p个元素的一维列向量其中，上标t表示矩阵或向量转置；步骤s8，利用步骤s7得到的包含p个元素的一维向量与步骤s5得到的列向量进行加权运算，得到当前网格值的插值结果，表达式为其中，上标t表示矩阵或向量转置。以上实施例仅为本发明的示例性实施例，不用于限制本发明，本发明的保护范围由权利要求书限定。本领域技术人员可以在本发明的实质和保护范围内，对本发明做出各种修改或等同替换，这种修改或等同替换也应视为落在本发明的保护范围内。当前第1页12