[go: up one dir, main page]

CN109459753B - 天气雷达数据坐标转换快速插值方法 - Google Patents

天气雷达数据坐标转换快速插值方法 Download PDF

Info

Publication number
CN109459753B
CN109459753B CN201710959444.1A CN201710959444A CN109459753B CN 109459753 B CN109459753 B CN 109459753B CN 201710959444 A CN201710959444 A CN 201710959444A CN 109459753 B CN109459753 B CN 109459753B
Authority
CN
China
Prior art keywords
interpolation
dimensional
value
sinc
dimensional matrix
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.)
Active
Application number
CN201710959444.1A
Other languages
English (en)
Other versions
CN109459753A (zh
Inventor
黄平平
李铁建
魏加华
谭维贤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tsinghua University
Inner Mongolia University of Technology
Original Assignee
Tsinghua University
Inner Mongolia University of Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Tsinghua University, Inner Mongolia University of Technology filed Critical Tsinghua University
Priority to CN201710959444.1A priority Critical patent/CN109459753B/zh
Publication of CN109459753A publication Critical patent/CN109459753A/zh
Application granted granted Critical
Publication of CN109459753B publication Critical patent/CN109459753B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • G01S13/958—Theoretical aspects
    • 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
    • G01S7/418—Theoretical aspects
    • 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

本发明提供了一种天气雷达数据坐标转换快速插值方法，根据预设的插值核点数和量化位移预先计算sinc插值核表格；计算待转换到的地理坐标系的均匀三维网格中的每个网络点的网格点坐标值所对应的以当前雷达自身为中心的极坐标系下的坐标，计算其在均匀极坐标离散网格中所处的位置；根据位置值的整数部分，从给定的离散化保存的天气雷达体积扫描数据中抽取一个三维矩阵数据块；根据位置值的小数部分从sinc插值核表格中查询并得到表格中的三行元素以分别组成列向量
Figure DDA0001434865570000011
和
Figure DDA0001434865570000012
利用三维矩阵数据块和
Figure DDA0001434865570000013
得到二维矩阵数据块；利用二维矩阵数据块和
Figure DDA0001434865570000014
得到一维列向量；利用一维列向量与
Figure DDA0001434865570000015
得到当前网格值的插值结果。本发明提升了插值速度。

Description

天气雷达数据坐标转换快速插值方法
技术领域
本发明涉及天气雷达信号与数据处理领域，尤其涉及一种用于天气雷达数据坐标转换的快速插值方法。
背景技术
天气雷达得到的数据产品通常处于以雷达自身为中心的极坐标，包含气象目标的距离，方位和俯仰坐标。对于多部天气雷达组网观测的应用，需要构建统一的坐标系，例如由经度、维度和高度组成的地理坐标系作为基准，各雷达得到的数据产品均通过插值方法由自身极坐标系变换至新的坐标系下，从而方便天气雷达数据的组网拼图。由于大部分气象目标，例如云在空间上具有连续性，并且有很多细尺度结构，因而希望插值后的散射率场在空间上要保持连续性，同时在内插过程中应最大限度地保留雷达资料中存在的原始回波结构特征。常见的插值方法包括：最近邻居法、线性内插法和sinc核插值法。其中，最近邻居法虽然速度快，但是精度过低，易导致插值后的雷达散射率在空间上不连续；线性内插法对原始插值前数据引入较大的平滑，因而易损失原有散射率细尺度结构特征。sinc核插值法是根据奈奎斯特采样定理进行插值的方法，精度最高，但是通常情况下计算复杂度较高。
发明内容
鉴于上述技术问题，本发明提供了一种用于天气雷达数据坐标转换的快速插值方法。
本发明中的天气雷达数据坐标转换快速插值方法，包括：步骤1，根据预设的插值核点数和量化位移预先计算sinc插值核表格；步骤2，通过映射关系计算出待转换到的地理坐标系的均匀三维网格中的每个网络点的网格点坐标值所对应的以当前雷达自身为中心的极坐标系下的坐标；步骤3，对于步骤2计算得到的某一个具体的距离、俯仰和方位坐标，计算其在天气雷达体积扫描数据的均匀极坐标离散网格中所处的位置；步骤4，根据步骤3得到的所述位置的位置值的整数部分，从给定的离散化保存的天气雷达体积扫描数据中抽取一个三维矩阵数据块；步骤5，分别根据步骤3得到的所述位置的位置值的小数部分，从sinc插值核表格中查询并得到表格中的三行元素以分别组成列向量
Figure GDA0003743895550000021
和
Figure GDA0003743895550000022
步骤6，利用步骤4得到的三维矩阵数据块和步骤5得到的列向量
Figure GDA0003743895550000023
进行加权运算，得到二维矩阵数据块；步骤7，利用步骤6得到的二维矩阵数据块和步骤5得到的列向量
Figure GDA0003743895550000024
进行加权运算，得到一维列向量；步骤8，利用步骤7得到的一维列向量，与步骤5得到的列向量
Figure GDA0003743895550000025
进行加权运算，得到当前网格值的插值结果。
优选地，所述步骤1中，所述的插值核点数P的典型取值是6～16间的任一偶数，量化位移1/L中的L的典型取值为大于等于10的偶数；
所述的sinc插值核表格是L+1行，P列的数值表格；
所述sinc插值核表格的第i行，第j列的元素wi,j的值通过以下公式计算得到
Figure GDA0003743895550000026
Figure GDA0003743895550000027
其中，i＝1,2,...L+1，j＝1,2,...P，sin c(x)＝sin(πx)/(πx)表示sinc函数。
优选地，所述步骤2中，由地理坐标系的网格点坐标值(xlat,m′,ylon,n′,hk′)到以当前雷达自身为中心的极坐标系坐标的映射关系的表达式为：
Figure GDA0003743895550000031
其中，(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中，由某一个具体的距离、俯仰和方位坐标
Figure GDA0003743895550000032
计算其在天气雷达体积扫描数据的均匀极坐标离散网格
Figure GDA0003743895550000033
中所处的位置(x1,x2,x3)的方法为
Figure GDA0003743895550000034
Figure GDA0003743895550000035
Figure GDA0003743895550000041
优选地，所述步骤4包括：
步骤41，分别对x1,x2,x3取整，得到
Figure GDA0003743895550000042
Figure GDA0003743895550000043
其中
Figure GDA0003743895550000044
表示取整算子；
步骤S42，在给定的离散化保存的天气雷达体积扫描数据
Figure GDA0003743895550000045
中，其中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的小数部分
Figure GDA0003743895550000046
将其除以量化位移1/L的值并四舍五入为整数
Figure GDA0003743895550000047
其中
Figure GDA0003743895550000048
表示取整算子，查询sinc插值表格并选定插值表中第L+1-m1行的元素作为加权值组成行向量
Figure GDA0003743895550000049
round()表示四舍五入取整算子；
计算x2的小数部分
Figure GDA0003743895550000051
将其除以量化位移1/L的值并四舍五入为整数
Figure GDA0003743895550000052
查询sinc插值表格并选定插值表中第L+1-m2行的元素作为加权值组成行向量
Figure GDA0003743895550000053
计算x3的小数部分
Figure GDA0003743895550000054
将其除以量化位移1/L的值并四舍五入为整数
Figure GDA0003743895550000055
查询sinc插值表格并选定插值表中第L+1-m3行的元素作为加权值组成行向量
Figure GDA00037438955500000510
优选地，所述步骤6包括：
在三维矩阵数据块sP×P×P中，针对某个固定的索引对(j,l)，将三维矩阵数据块sP×P×P的P个数据元素s(1,j,l)，s(2,j,l)，...,s(P,j,l)组成列向量
Figure GDA0003743895550000056
然后求取该向量与列向量
Figure GDA0003743895550000057
的内积
Figure GDA0003743895550000058
作为一个P×P二维矩阵数据块s′P×P×P的第j行，第l列的数据元素，其中，上标T表示矩阵或向量转置。对所有的索引对(j,l)，j,l＝1,2,...P进行遍历，直到计算得到P×P二维矩阵数据块s′P×P×P。
优选地，所述步骤7中的一维列向量的表达式为：
Figure GDA0003743895550000059
其中，
Figure GDA0003743895550000061
为一维列向量，s′P×P×P为二维矩阵数据块，
Figure GDA0003743895550000062
为列向量，T表示矩阵或向量转置。
优选地，所述步骤8中，当前网格值的插值结果为：
Figure GDA0003743895550000063
其中，T表示矩阵或向量转置。
从上述技术方案可以看出，本发明具有以下有益效果：
(1)本发明技术方案中天气雷达数据坐标变换使用sinc核插值，插值精度较高，有助于保留气象目标散射率的细尺度结构特征；
(2)本发明技术方案中插值的实现是通过查询预先计算好的sinc插值核表格获得加权值，避免了实时计算复杂度较高的sinc函数的过程，提升了插值速度。
附图说明
图1为本发明实施例天气雷达进行体积扫描示意图；
图2为本发明实施例中使用的sinc核插值原理示意图。
图3为本发明实施例一种用于天气雷达数据坐标转换的快速插值方法的流程图。
具体实施方式
为使本领域技术人员更好的理解本发明的技术方案，下面结合附图和具体实施方式对本发明作详细说明。
本发明技术方案中天气雷达数据坐标变换使用sinc核插值，插值精度较高，有助于保留气象目标散射率的细尺度结构特征；插值的实现是通过查询预先计算好的sinc插值核表格获得加权值，避免了实时计算复杂度较高的sinc函数的过程，提升了插值速度。
本发明的实施方式是针对天气雷达数据由自身极坐标数据变换到其它坐标系如由经纬高表示的地理坐标系进行的。在详细介绍本发明的实施方式的细节之前，首先简要介绍利用数据插值实现坐标变换的原理。
请参考图1所示，通常，天气雷达进行体积扫描，可获取气象目标在极坐标系下的距离r、方位角
Figure GDA0003743895550000071
俯仰角θ坐标的雷达散射率因子的基本数据产品，设其表示为
Figure GDA0003743895550000072
设数据实施坐标变换到地理坐标系后的数据应当为Z′(xlat,ylon,h)，其中xlat表示经度坐标，ylon表示维度坐标，h表示高度坐标。对于某个具体的坐标值(xlat,ylon,h)处的气象目标，坐标变换后的数据Z′(xlat,ylon,h)是由
Figure GDA0003743895550000073
通过坐标映射得到的：
Figure GDA0003743895550000074
其中，r(xlat,ylon,h),θ(xlat,ylon,h),
Figure GDA0003743895550000075
可通过大圆几何学理论求得，设天气雷达自身位置的经纬高坐标为(xr,yr,hr)，则
Figure GDA0003743895550000076
其中，R表示地球半径，s的表达式为
s＝R×arccos[sin(xlat)sin(xr)+cos(xlat)cos(xr)cos(ylat-yr)]
通常情况下，天气雷达体扫描获取的数据
Figure GDA0003743895550000081
是在极坐标系下均匀网格分布的，离散化保存的数据
Figure GDA0003743895550000082
的，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),
Figure GDA0003743895550000083
不一定正好处于极坐标系上的整数网格点上，需要利用其周围网格点上的数据通过插值方法得到变换后的数据Z′(xlat,ylon,h)。
常见的插值方法包括：最近邻居法、线性内插法和sinc核插值法。其中，最近邻居法虽然速度快，但是精度过低，易导致插值后的雷达散射率在空间上不连续；线性内插法对原始插值前数据引入较大的平滑，因而易损失原有散射率细尺度结构特征。sinc核插值法是根据奈奎斯特采样定理进行插值的方法，精度最高，但是通常情况下计算复杂度较高。
本发明的实施方式中插值是基于sinc核插值进行的。在详细介绍本发明的实施方式的细节之前，还需先简单描述sinc核插值的一些概念和原理。
通常，在满足信号带限和采样率大于奈奎斯特采样频率的前提下，离散序列s[i]，i＝1,2,...可以通过sinc核插值得到任意非整数采样点x上的高精度的插值结果s′(x)，表示如下：
Figure GDA0003743895550000091
其中，sin c(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的值通过以下公式计算得到
Figure GDA0003743895550000101
Figure GDA0003743895550000102
其中，i＝1,2,...L+1，j＝1,2,...P，sin c(x)＝sin(πx)/(πx)表示sinc函数。
表1为在P＝8，L＝12时的13行8列的插值核表格
0 0 0 0 1 0 0 0
-0.021 0.028 -0.043 0.090 0.990 -0.076 0.040 -0.027
-0.042 0.057 -0.087 0.192 0.961 -0.137 0.074 -0.051
-0.061 0.083 -0.130 0.304 0.912 -0.182 0.101 -0.070
-0.077 0.105 -0.169 0.422 0.843 -0.211 0.120 -0.084
-0.088 0.122 -0.199 0.540 0.756 -0.222 0.130 -0.092
-0.093 0.131 -0.218 0.653 0.653 -0.218 0.131 -0.093
-0.092 0.130 -0.222 0.756 0.540 -0.199 0.122 -0.088
-0.084 0.120 -0.211 0.843 0.422 -0.169 0.105 -0.077
-0.0717-1059 0.101 -0.182 0.912 0.304 -0.130 0.083 -0.061
-0.051 0.074 -0.137 0.961 0.192 -0.087 0.057 -0.042
-0.027 0.040 -0.076 0.990 0.090 -0.043 0.028 -0.021
0 0 0 1 0 0 0 0
步骤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计算出其对应的以当前雷达自身为中心的极坐标系下的坐标
Figure GDA0003743895550000111
然后都重复执行下述步骤S3～S8，其中r′,θ′,
Figure GDA0003743895550000112
分别表示映射的极坐标系下的距离、俯仰角和方位角坐标；
所述的由地理坐标系的网格点坐标值(xlat,m′,ylon,n′,hk′)到以当前雷达自身为中心的极坐标系坐标的映射关系f的表达式为：
Figure GDA0003743895550000113
其中，(xr,yr,hr)为天气雷达自身位置的经纬高坐标，R表示地球半径，s的表达式为
s＝R×arccos[sin(xlat,m′)sin(xr)+cos(xlat,m′)cos(xr)cos(ylon,n′-yr)]
步骤S3，对于步骤S2计算得到的某一个具体的距离、俯仰和方位坐标
Figure GDA0003743895550000121
计算其在天气雷达体积扫描数据的均匀极坐标离散网格
Figure GDA0003743895550000122
中所处的位置(x1,x2,x3)，其表示距离坐标r′位于rm序列中的第x1个样点，m＝1,2,...M；俯仰坐标θ′位于θn序列中的第x2个样点，n＝1,2,...N；方位坐标
Figure GDA0003743895550000123
位于
Figure GDA0003743895550000124
序列中的第x3个样点，k＝1,2,...K；x1,x2,x3为不小于1的整数或非整数；
所述的由某一个具体的距离、俯仰和方位坐标
Figure GDA0003743895550000125
计算其在天气雷达体积扫描数据的均匀极坐标离散网格
Figure GDA0003743895550000126
中所处的位置(x1,x2,x3)的方法为：
Figure GDA0003743895550000127
Figure GDA0003743895550000128
Figure GDA0003743895550000129
步骤S4，根据步骤S3得到的x1,x2,x3位置值的整数部分，从给定的离散化保存的天气雷达体积扫描数据
Figure GDA0003743895550000131
Figure GDA0003743895550000132
中，抽取出一个P×P×P维的三维矩阵数据块s(i,j,l)，i,j,l＝1,2,...P,其中整数P为所述步骤s1所示的预设的插值核点数P；
所述的根据x1,x2,x3位置值的整数部分，从给定的离散化保存的天气雷达体积扫描数据
Figure GDA0003743895550000133
Figure GDA0003743895550000134
中，抽取出一个P×P×P维的三维矩阵数据块s(i,j,l)，i,j,l＝1,2,...P的步骤又可以包含：
子步骤S41:分别对x1,x2,x3取整，得到
Figure GDA0003743895550000135
Figure GDA0003743895550000136
其中
Figure GDA0003743895550000137
表示取整算子；
子步骤S42:对于给定的离散化保存的天气雷达体积扫描数据
Figure GDA0003743895550000138
中，按第一个维度的索引从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插值核表格中查询并得到表格中的三行元素，分别组成列向量
Figure GDA0003743895550000139
和
Figure GDA00037438955500001310
所述的根据x1,x2,x3位置值的小数部分，从sinc插值核表格中查询并得到表格中的三行元素，分别组成列向量
Figure GDA0003743895550000141
和
Figure GDA0003743895550000142
的方法是：
计算x1的小数部分
Figure GDA0003743895550000143
将其除以量化位移1/L的值并四舍五入为整数
Figure GDA0003743895550000144
其中
Figure GDA0003743895550000145
表示取整算子，查询sinc插值表格并选定插值表中第L+1-m1行的元素作为加权值组成行向量
Figure GDA0003743895550000146
round()表示四舍五入取整算子；
计算x2的小数部分
Figure GDA0003743895550000147
将其除以量化位移1/L的值并四舍五入为整数
Figure GDA0003743895550000148
查询sinc插值表格并选定插值表中第L+1-m2行的元素作为加权值组成行向量
Figure GDA0003743895550000149
计算x3的小数部分
Figure GDA00037438955500001410
将其除以量化位移1/L的值并四舍五入为整数
Figure GDA00037438955500001411
查询sinc插值表格并选定插值表中第L+1-m3行的元素作为加权值组成行向量
Figure GDA00037438955500001412
步骤S6，利用步骤S4得到的三维矩阵数据块sP×P×P和步骤S5得到的列向量
Figure GDA00037438955500001413
进行加权运算，得到P×P二维矩阵数据块s′P×P×P；
所述的利用三维矩阵数据块sP×P×P和列向量
Figure GDA00037438955500001414
进行加权运算的过程为：
在三维矩阵数据块sP×P×P中，针对某个固定的索引对(j,l)，将三维矩阵数据块sP×P×P的P个数据元素s(1,j,l)，s(2,j,l)，...,s(P,j,l)组成列向量
Figure GDA0003743895550000151
然后求取该向量与列向量
Figure GDA0003743895550000152
的内积
Figure GDA0003743895550000153
作为一个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得到的列向量
Figure GDA0003743895550000154
进行加权运算，得到包含P个元素的一维列向量
Figure GDA0003743895550000155
所述的利用二维矩阵数据块s′P×P×P和列向量
Figure GDA0003743895550000156
进行加权运算，利用的是向量和矩阵的乘法，其表达式为：
Figure GDA0003743895550000157
得到一个包含P个元素的一维列向量
Figure GDA0003743895550000158
其中，上标T表示矩阵或向量转置；
步骤S8，利用步骤S7得到的包含P个元素的一维向量
Figure GDA0003743895550000159
与步骤S5得到的列向量
Figure GDA00037438955500001510
进行加权运算，得到当前网格值的插值结果，表达式为
Figure GDA00037438955500001511
其中，上标T表示矩阵或向量转置。
以上实施例仅为本发明的示例性实施例，不用于限制本发明，本发明的保护范围由权利要求书限定。本领域技术人员可以在本发明的实质和保护范围内，对本发明做出各种修改或等同替换，这种修改或等同替换也应视为落在本发明的保护范围内。

Claims (9)

1.一种天气雷达数据坐标转换快速插值方法，其特征在于，包括：
步骤1，根据预设的插值核点数和量化位移预先计算sinc插值核表格；
步骤2，通过映射关系计算出待转换到的地理坐标系的均匀三维网格中的每个网络点的网格点坐标值所对应的以当前雷达自身为中心的极坐标系下的坐标；
步骤3，对于步骤2计算得到的某一个具体的距离、俯仰和方位坐标，计算其在天气雷达体积扫描数据的均匀极坐标离散网格中所处的位置；
步骤4，根据步骤3得到的所述位置的位置值的整数部分，从给定的离散化保存的天气雷达体积扫描数据中抽取一个三维矩阵数据块；
步骤5，分别根据步骤3得到的所述位置的位置值的小数部分，从sinc插值核表格中查询并得到表格中的三行元素以分别组成列向量
Figure FDA0003728854700000011
Figure FDA0003728854700000012
和
Figure FDA0003728854700000013
步骤6，利用步骤4得到的三维矩阵数据块和步骤5得到的列向量
Figure FDA0003728854700000014
进行加权运算，得到二维矩阵数据块；
步骤7，利用步骤6得到的二维矩阵数据块和步骤5得到的列向量
Figure FDA0003728854700000015
进行加权运算，得到一维列向量；
步骤8，利用步骤7得到的一维列向量，与步骤5得到的列向量
Figure FDA0003728854700000016
进行加权运算，得到当前网格值的插值结果。
2.根据权利要求1所述的方法，其特征在于，所述步骤1中，所述的插值核点数P的典型取值是6～16间的任一偶数，量化位移1/L中的L的典型取值为大于等于10的偶数；
所述的sinc插值核表格是L+1行，P列的数值表格；
所述sinc插值核表格的第i行，第j列的元素wi,j的值通过以下公式计算得到
Figure FDA0003728854700000017
Figure FDA0003728854700000021
其中，i＝1,+，j＝1,2,...P，sinc(x)＝sin(πx)/(πx)表示sinc函数。
3.根据权利要求2所述的方法，其特征在于，所述步骤2中，由地理坐标系的网格点坐标值(xlat,m′,ylon,n′,hk′)到以当前雷达自身为中心的极坐标系坐标的映射关系的表达式为：
Figure FDA0003728854700000022
其中，(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′个高度坐标。
4.根据权利要求3所述的方法，其特征在于，所述步骤3中，由某一个具体的距离、俯仰和方位坐标
Figure FDA0003728854700000031
计算其在天气雷达体积扫描数据的均匀极坐标离散网格
Figure FDA0003728854700000032
中所处的位置(x1,x2,x3)的方法为
Figure FDA0003728854700000033
Figure FDA0003728854700000034
Figure FDA0003728854700000035
其中，
M，N和K分别表示插值前极坐标系上天气雷达数据在距离、俯仰和方位方向的数据维度。
5.根据权利要求4所述的方法，其特征在于，所述步骤4包括：
步骤41，分别对x1,x2,x3取整，得到
Figure FDA0003728854700000036
Figure FDA0003728854700000037
其中
Figure FDA0003728854700000038
表示取整算子；
步骤S42，在给定的离散化保存的天气雷达体积扫描数据
Figure FDA0003728854700000039
中，其中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。
6.根据权利要求5所述的方法，其特征在于，所述步骤5包括：
计算x1的小数部分
Figure FDA0003728854700000041
将其除以量化位移1/L的值并四舍五入为整数
Figure FDA0003728854700000042
其中
Figure FDA0003728854700000043
表示取整算子，查询sinc插值表格并选定插值表中第L+1-m1行的元素作为加权值组成行向量
Figure FDA0003728854700000044
round()表示四舍五入取整算子；
计算x2的小数部分
Figure FDA0003728854700000045
将其除以量化位移1/L的值并四舍五入为整数
Figure FDA0003728854700000046
查询sinc插值表格并选定插值表中第L+1-m2行的元素作为加权值组成行向量
Figure FDA0003728854700000047
计算x3的小数部分
Figure FDA0003728854700000048
将其除以量化位移1/L的值并四舍五入为整数
Figure FDA0003728854700000049
查询sinc插值表格并选定插值表中第L+1-m3行的元素作为加权值组成行向量
Figure FDA00037288547000000410
7.根据权利要求6所述的方法，其特征在于，所述步骤6包括：
在三维矩阵数据块sP×P×P中，针对某个固定的索引对(j,l)，将三维矩阵数据块sP×P×P的P个数据元素s(1,j,l)，s(2,j,l)，...,s(P,j,l)组成列向量
Figure FDA00037288547000000411
然后求取该向量与列向量
Figure FDA0003728854700000051
的内积
Figure FDA0003728854700000052
作为一个P×P二维矩阵数据块s′P×P×P的第j行，第l列的数据元素，其中，上标T表示矩阵或向量转置；对所有的索引对(j,l)，j,l＝1,2,...P进行遍历，直到计算得到P×P二维矩阵数据块s′P×P×P。
8.根据权利要求7所述的方法，其特征在于，所述步骤7中的一维列向量的表达式为：
Figure FDA0003728854700000053
其中，
Figure FDA0003728854700000054
为一维列向量，s′P×P×P为二维矩阵数据块，
Figure FDA0003728854700000055
为列向量，T表示矩阵或向量转置。
9.根据权利要求8所述的方法，其特征在于，所述步骤8中，当前网格值的插值结果为：
Figure FDA0003728854700000056
其中，T表示矩阵或向量转置。
CN201710959444.1A 2017-10-16 2017-10-16 天气雷达数据坐标转换快速插值方法 Active CN109459753B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710959444.1A CN109459753B (zh) 2017-10-16 2017-10-16 天气雷达数据坐标转换快速插值方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710959444.1A CN109459753B (zh) 2017-10-16 2017-10-16 天气雷达数据坐标转换快速插值方法

Publications (2)

Publication Number Publication Date
CN109459753A CN109459753A (zh) 2019-03-12
CN109459753B true CN109459753B (zh) 2022-10-11

Family

ID=65606156

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710959444.1A Active CN109459753B (zh) 2017-10-16 2017-10-16 天气雷达数据坐标转换快速插值方法

Country Status (1)

Country Link
CN (1) CN109459753B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110261857B (zh) * 2019-07-17 2022-04-15 南京信息工程大学 一种天气雷达空间插值方法
CN110412551B (zh) * 2019-07-20 2021-02-26 中国船舶重工集团公司第七二四研究所 一种超视距探测目标信息跨平台交接坐标转换方法
CN110940978A (zh) * 2019-12-09 2020-03-31 上海眼控科技股份有限公司 一种雷达ppi图像显示方法和装置、电子设备、存储介质

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6384766B1 (en) * 1997-06-18 2002-05-07 Totalförsvarets Forskningsinstitut Method to generate a three-dimensional image of a ground area using a SAR radar
JP2010278873A (ja) * 2009-05-29 2010-12-09 Victor Co Of Japan Ltd 色変換装置
CN102117227A (zh) * 2011-03-09 2011-07-06 南京恩瑞特实业有限公司 天气雷达数据的多核并行计算方法
CN102393520A (zh) * 2011-09-26 2012-03-28 哈尔滨工程大学 基于目标回波多普勒特性的声纳运动目标成像方法
CN103197299A (zh) * 2013-03-25 2013-07-10 南京信息工程大学 天气雷达径向风信息提取及量化分析系统
CN103530627A (zh) * 2013-10-23 2014-01-22 东南大学 基于二维散射中心集网格模型的isar图像恢复方法
CN103630901A (zh) * 2013-03-29 2014-03-12 中国科学院电子学研究所 机载下视阵列3-d sar成像的方法
CN105204005A (zh) * 2015-10-19 2015-12-30 中国电子科技集团公司第二十八研究所 一种基于地理坐标系的vts系统雷达回波视频显示方法
CN105701859A (zh) * 2016-02-22 2016-06-22 武汉华信联创技术工程有限公司 一种雷达单站极坐标数据的三维格点化处理方法和系统
CN107180014A (zh) * 2017-04-28 2017-09-19 华讯方舟科技有限公司 一种快速sinc插值方法及系统

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6384766B1 (en) * 1997-06-18 2002-05-07 Totalförsvarets Forskningsinstitut Method to generate a three-dimensional image of a ground area using a SAR radar
JP2010278873A (ja) * 2009-05-29 2010-12-09 Victor Co Of Japan Ltd 色変換装置
CN102117227A (zh) * 2011-03-09 2011-07-06 南京恩瑞特实业有限公司 天气雷达数据的多核并行计算方法
CN102393520A (zh) * 2011-09-26 2012-03-28 哈尔滨工程大学 基于目标回波多普勒特性的声纳运动目标成像方法
CN103197299A (zh) * 2013-03-25 2013-07-10 南京信息工程大学 天气雷达径向风信息提取及量化分析系统
CN103630901A (zh) * 2013-03-29 2014-03-12 中国科学院电子学研究所 机载下视阵列3-d sar成像的方法
CN103530627A (zh) * 2013-10-23 2014-01-22 东南大学 基于二维散射中心集网格模型的isar图像恢复方法
CN105204005A (zh) * 2015-10-19 2015-12-30 中国电子科技集团公司第二十八研究所 一种基于地理坐标系的vts系统雷达回波视频显示方法
CN105701859A (zh) * 2016-02-22 2016-06-22 武汉华信联创技术工程有限公司 一种雷达单站极坐标数据的三维格点化处理方法和系统
CN107180014A (zh) * 2017-04-28 2017-09-19 华讯方舟科技有限公司 一种快速sinc插值方法及系统

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
"Efficient Pseudopolar Format Algorithm for Down-Looking Linear-Array SAR 3-D Imaging";Han Kuoye 等;《IEEE Geoscience and Remote Sensing Letters》;20151231;第572-576页 *
"SAR极坐标格式处理波前弯曲补偿方法研究";芜晓丹;《硕士电子期刊》;20130415;全文 *
"The Polar Format Imaging Algorithm for Forward-looking Bistatic SAR";Jin-ping SUN 等;《7th European Conference on Synthetic Aperture Radar》;20081231;第1-4页 *
"高分辨机载SAR两维自聚焦处理及FPGA实现";郭江哲;《万方数据》;20170405;全文 *

Also Published As

Publication number Publication date
CN109459753A (zh) 2019-03-12

Similar Documents

Publication Publication Date Title
CN109459753B (zh) 天气雷达数据坐标转换快速插值方法
CN106291473B (zh) 嵌套式天线阵列设置方法
CN106556874A (zh) 一种近距离微波成像方法及系统
JP2023513650A (ja) 相互相関テンソルに基づく三次元の互いに素のキュービックアレイの到来方向推定方法
CN105654483A (zh) 三维点云全自动配准方法
CN108802669A (zh) 二维波达方向估计方法、二维波达方向估计装置及终端
CN113936046B (zh) 一种物体定位方法、装置、电子设备及计算机可读介质
CN104615880A (zh) 一种三维激光雷达点云匹配的快速icp方法
CN112424628A (zh) 定位设备
CN112363122A (zh) 高频地波雷达电离层噪声中弱谐波信号的提取方法及应用
CN106157258B (zh) 一种星载sar图像几何校正方法
CN104778260B (zh) 一种动态雷达环境知识库建模方法
CN112782647B (zh) 信息联合的二次等式约束最小二乘辐射源定位方法
CN106908760B (zh) 基于阵列自相关矩阵的单站无源定位方法
CN105549010A (zh) 频域合成孔径雷达成像方法
CN109507634B (zh) 一种任意传感器阵列下的基于传播算子的盲远场信号波达方向估计方法
CN113447887B (zh) 全空间定位方法、装置、设备与计算机可读存储介质
CN115529073B (zh) 基于无人机节点通信的分布式协同干扰定位方法及装置
CN117232523A (zh) 果园机器人导航线规划方法、装置、设备及介质
CN116229005A (zh) 三维巷道模型的测地线确定方法和装置
CN116486259A (zh) 遥感图像中的点目标的提取方法和装置
CN103955602A (zh) 一种综合孔径微波辐射计阵列因子成型方法
CN106154245A (zh) 基于等效阵列方向图的集中式mimo雷达阵列设计方法
CN106650645B (zh) 遥感图像的目标特征提取方法和装置
CN113835064A (zh) 一种协同校正源观测信息的加权多维标度tdoa定位方法

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
GR01 Patent grant
GR01 Patent grant