CN111693999B - 基于雷达测风组合策略的多传感器融合风速风向估计方法 - Google Patents
基于雷达测风组合策略的多传感器融合风速风向估计方法 Download PDFInfo
- Publication number
- CN111693999B CN111693999B CN202010481799.6A CN202010481799A CN111693999B CN 111693999 B CN111693999 B CN 111693999B CN 202010481799 A CN202010481799 A CN 202010481799A CN 111693999 B CN111693999 B CN 111693999B
- Authority
- CN
- China
- Prior art keywords
- wind
- wind speed
- sensor
- wind direction
- value
- 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
Links
- 238000005259 measurement Methods 0.000 title claims abstract description 72
- 238000000034 method Methods 0.000 title claims abstract description 54
- 230000004927 fusion Effects 0.000 title claims abstract description 45
- 239000013598 vector Substances 0.000 claims abstract description 21
- 238000012216 screening Methods 0.000 claims abstract description 10
- 238000004364 calculation method Methods 0.000 claims description 10
- 238000002474 experimental method Methods 0.000 claims description 8
- 238000000354 decomposition reaction Methods 0.000 claims description 4
- 230000000694 effects Effects 0.000 abstract description 12
- 238000012937 correction Methods 0.000 abstract description 4
- 238000007500 overflow downdraw method Methods 0.000 abstract description 3
- 238000004422 calculation algorithm Methods 0.000 description 28
- 239000012530 fluid Substances 0.000 description 16
- 238000001228 spectrum Methods 0.000 description 14
- 230000003068 static effect Effects 0.000 description 11
- 238000011161 development Methods 0.000 description 7
- 238000005516 engineering process Methods 0.000 description 5
- 238000001914 filtration Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 238000005096 rolling process Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 238000012935 Averaging Methods 0.000 description 2
- 238000009434 installation Methods 0.000 description 2
- 238000013178 mathematical model Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- XUIMIQQOPSSXEZ-UHFFFAOYSA-N Silicon Chemical group [Si] XUIMIQQOPSSXEZ-UHFFFAOYSA-N 0.000 description 1
- 238000003723 Smelting Methods 0.000 description 1
- 241000538662 Tangaroa Species 0.000 description 1
- 244000153888 Tung Species 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 239000003973 paint Substances 0.000 description 1
- 238000005192 partition Methods 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 238000012552 review Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 230000005477 standard model Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
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/956—Radar or analogous systems specially adapted for specific applications for meteorological use mounted on ship or other platform
-
- 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)
- Ocean & Marine Engineering (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明提供的是一种基于雷达测风组合策略的多传感器融合风速风向估计方法。包括获取传感器权值;根据雷达测量的自由来流风向值作为相对风向的参考约束,对风向角度范围约束下的传感器信息进行校验,筛选对应分组的传感器权值,进而估计出不同传感器权值下的风向、风速值;将得到的每组风速分量矢量相加,计算出每组风向和风速的估计值,引入代价函数,在加入雷达测量风向角度误差范围内,输入每组风向和风速的估计值,输出代价函数最小的结果作为风向和风速的最优估计值。相比于传统的最优加权融合方法,本发明的方法大大提高了校正后的风速、风向测量精度,比传统的方法融合效果更好。
Description
技术领域
本发明涉及的是一种风速风向估计方法,特别是一种多传感器融合风速风向估计方法。
背景技术
现如今,随着科技的发展,国家越来越重视海洋的开发与利用,海洋中存在着丰富的资源。海上的船只设备的发展必须紧跟时代的发展。对于一些可以提供飞机起飞或者降落的船舶,船面气流场的准确空间分布情况测量和估计是非常重要的一个环节。例如在直升机上升或降落时,需要为飞行员实时提供船面稳态自由来流风速、风向测量或估计参数,以便飞行员控制飞机安全起飞或降落[2],如果不能准确得到这些参数,则会对飞机起飞、降落过程造成隐形的威胁因素。见参考文献[1-2](洪术华,宋雍,叶景波,陈荣.海洋工程发展现状与跨越发展战略[J].船舶工程,2019,41(S2):264-268.马欣瑞.航母舰岛位置及构型对舰载机着舰安全影响研究[C].中国力学学会流体力学专业委员会.2019年全国工业流体力学会议摘要集.中国力学学会流体力学专业委员会:北京航空航天大学陆士嘉实验室,2019:21.)。
在舰面风速、风向传感器技术方面,常用的风车、风杯式传感器装置大多数使用在海上的大小船只中,它通过尾部的偏转角来测量风向,但这类传感器测量精度不高,可靠性也比较差[3],因此,它一般被用于普通船只。近年来随着科技的发展出现了固态风速风向传感器[4],其利用红外线、超声波或硅芯等测量风场信息,克服了传统风传感器的缺点,有较高的测量精度[5]。但由于受船舶甲板及船上建筑物的影响,涡流区安装的传感器所测风场与实际真实自由来流风场还是存在较大的误差。见参考文献[3-5](PAPADOPOULOS K H,STEFANTOS N C,PAULSEN US,et al.Effects of turbulence and flow inclination onthe performance of cup anemometers in the field[J].Boundary-LayerMeteorology,2001,101(1):77-107.张燕波,沈广平,董自强等.基于微控制器的风速风向传感器系统设计[J].仪器仪表学报,2009,30(10):2144-2149.王国峰,赵永生,范云生.风速风向测量误差补偿算法的研究[J].仪器仪表学报,2013,(04):786-790.)
船体上通常会安装两个以上的风速风向传感器测量风速风向。大致分布位置为左右对称或者在船首位置。传统的测风方法为取船体所有风速风向传感器测量值中风速最大的一个传感器测量值作为风速风向的测量值。这种测量方法并不能准确地得到风速风向值,只是排除了一些所受影响较大的测量值,最后得到的值仍然存在误差,所以需要找到一种可以更加准确估计舰面稳态风速风向值的方法。1982年,刘连吉等用实验的方法研究了船上测风仪测量的风场数据与浮标测量的风场数据之间的误差,并提出了降低测风仪测量误差的改进意见[6]。2004年,Popinet[7]研究了R/V Tangaroa对平均气流和湍流尾流特性的影响。发现风速只弱依赖于船舶速度、船舶运动和海洋状态,强依赖于相对风向。2012年,刘慧和胡松等人通过船测风场数据与QuikSCAT卫星遥感所测风场数据的比对分析,研究了船测风场的偏差统计特征[8]。2014年,周亦武和王国锋等人分析了船舶在运动时风速仪测量的误差,通过大量试验数据的分析,提出了一种多变量非线性拟合的补偿算法,降低了船舶在运动状态下舰船测风仪的风场测量误差[9]。2014年,郜冶、谭大力[10]等以CFD计算软件FLUNENT13.0为平台,仿真了不同工况下风速风向仪周边气流场的失调状况。结果表明,风速仪安装位置向外平移有利于减小风向角和风速值测量误差。此外,桅杆的近风侧传感器测量误差更小,且与迎风角度成反比,而远风侧测量误差则正好相反。2018年,胡桐[11]等对船舶横纵摇工况下船体周围气流场的分布情况进行数值计算,并使用实船实验数据对数值计算进行了验证。结果表明,船舶的横摇与纵摇会影响风速风向仪对自由来流的风速风向测量,且在不同横摇纵摇角度速度下,只对传感器数据进行姿态校正是不全面的,还需要对船体上层建筑对气流的影响程度进行分析与校正。2018年,李盼飞[12]利用数值模拟,仿真了舰船的气流分布情况,分析发现风速仪的误差与位置有关,进而提出了一种基于矢量差的风速仪的最优布置方案,并使用多传感器的最优加权融合算法与神经网络对舰面风速风向进行估计。见参考文献[6-12](刘连吉.对现用海洋水文气象要素调查规范的改进意见[J].海洋技术,1982,(02):70-73.Popinet S,Smith M,Stevens C.Experimental andnumerical study of the turbulence characteristics of air flow around aresearch vessel[J].Journal of Atmospheric and Oceanic Technology,2004,12(10):1575–1589.刘慧,胡松,邹晓荣.船测资料与智利外海QuikSCAT风场比较分析[J].遥感技术与应用,2012,27(05):763-769.周亦武,王国锋,赵永生.船舶摇摆状态下风速测量误差分析与补偿研究[J].仪器仪表学报,2014,35(06):1239-1245.郜冶,谭大力,李海旭,王金玲,刘长猛.舰载风速仪测量误差与安装位置的关系研究[J].哈尔滨工程大学学报,2014,35(10):1195-1200.胡桐,漆随平,王东明.船体纵横摇工况下海面风测量数据偏差分析[J].海洋技术学报,2018,37(3).李盼飞.基于多传感器融合的舰面稳态风场估计技术研究[D].哈尔滨工程大学,2018.)
现有的传感器对于风速、风向的估计精度不佳,由于船上传感器测风的误差与自由来流风向/风速均呈非线性关系,传统的最优加权融合算法所估计出的风并不会在所有风向角度下保持高精度。船上测风的另一个手段是利用航海雷达图像反演距离本船数公里外海面的风向值,其优势在于测量区域远离本船、自由流场没有被船体等结构物扰流破坏,但其缺点是精度较低。
发明内容
本发明的目的在于提供一种能够提高校正后的风速、风向测量精度的基于雷达测风组合策略的多传感器融合风速风向估计方法。
本发明的目的是这样实现的:
步骤1,获取传感器权值:
离线展开获取风速风向测量值实验,对桅杆部分建模后的结构进行网格划分以及边界条件、入口条件的设定,通过CFD计算软件计算出在不同传感器布置位置处的风向、风速值,将得到的不同传感器布置位置处的传感器风速风向值作为传感器融合的基本数据,在规定风向范围内对传感器风速值以固定步长逐一进行分组,对于多个传感器的条件下,对应求出分组后的每组传感器所占权值;
步骤2,引入雷达风向进行数据筛选:
根据雷达测量的自由来流风向值作为相对风向的参考约束,对风向角度范围约束下的传感器信息进行校验,筛选对应分组的传感器权值,进而估计出不同传感器权值下的风向、风速值;
步骤3,多传感器融合:
将得到的每组风速分量矢量相加,计算出每组风向和风速的估计值,引入代价函数,在加入雷达测量风向角度误差范围内,输入每组风向和风速的估计值,输出代价函数最小的结果作为风向和风速的最优估计值。
本发明还可以包括:
1.步骤1进一步包括:
步骤1.1,离线通过CFD计算软件开展获取不同传感器布置位置处的传感器风速风向测量值的实验,对桅杆部分建模,对桅杆模型结构进行网格划分,设定整个模型的边界条件,根据设定好的物理模型与方程计算出在不同传感器布置位置处的风向、风速值;
步骤1.2,根据步骤1.1得到不同传感器布置位置处的传感器风速风向值作为传感器融合的基本数据,将自由来流真风速vref与每个传感器布置位置风速vsensor沿X轴、Y轴进行正交分解;
步骤1.3,根据步骤1.2确定的真风速与所有传感器风速分量数据,在0°~180°的风向范围内对传感器风速值以步长δ逐一进行分组;
步骤1.4,对步骤1.3得到的所有分组,计算出每组角度下每个传感器对应的权值,计算方法如下:
设每个传感器之间是相互独立的,且其测量值是真实值的无偏估计,设各个传感器的X轴风速分量为Xxi,Y轴风速分量为Xyi,其中i=1,2,…,n,n为传感器个数;ωxi代表第i个传感器X轴风速分量的权重值,ωyi代表第i个传感器Y轴风速分量的权重值,为X轴分量最后测量的估计值,为Y轴分量最后测量的估计值;
对上式等号的左右两侧分别取期望,表示如下:
由于估计值与每个传感器X轴风速分量测量值Xxi、Y轴风速分量测量值Xyi均视为X轴风速分量、Y轴分速分量真实值Xx、Xy的无偏估计,上式等号左右两侧的期望值相同,所以推算出传感器权重系数的和为1,
根据上述公式,最后估计X轴风速分量、Y轴风速分量的方差值与每个传感器的权重系数、测量方差有关,若要找到最优的估计值,只需要分别找到方差值最小值,又由于权重系数满足权重系数的和为1,通过拉格朗日乘数法,找到最优的方差值、公式如下:
对上式的变量ωxi、ωyi与λ分别进行求导,使导数为0,求出估计值方差的极小值,
权重系数的求解公式为:
则传感器的估计值为:
上式中的每个传感器的测量方差通过测量的数据与参考的数据计算得到,从而求出权重系数。
2.步骤2进一步包括:
步骤2.1,引入雷达的风向测量数据θradar对传感器测量数据进行风向角度范围约束;
步骤2.2,雷达测量的风向数据存在误差,在加入雷达测量出的风向误差约束角度范围θradar±θ0内,在完成步骤1的分组后查找相对应角度范围下每个传感器的权值分布并存储;
步骤2.3,根据2.2查找得到的每个传感器对应的权值,估计每组角度范围下风速风向估计值。
3.步骤3进一步包括:
步骤3.2,根据步骤3.1得到每组风速风向估计值,引入了一个代价函数,在加入雷达测量出的风向误差约束角度范围θradar±θ0内,取其中代价函数最小的那一组风速作为输出,代价函数定义为:
本发明首先,离线展开获取风速风向测量值实验,对桅杆部分建模后的结构进行网格划分以及边界条件、入口条件的设定,根据设定好的物理模型与方程通过CFD计算软件计算出在不同传感器布置位置处的风向、风速值,将得到的不同传感器布置位置处的传感器风速风向值作为传感器融合的基本数据,在规定风向范围内对传感器风速值以固定步长逐一进行分组,对于多个传感器的条件下,对应求出分组后的每组传感器的权值;然后,根据雷达测量的自由来流风向值作为相对风向的参考约束,对该风向角度范围约束下的传感器信息进行校验,筛选对应分组的传感器权值,进而估计出不同传感器权值下的风向、风速值;最后,将得到的每组风速分量矢量相加,计算出每组风向和风速的估计值。引入代价函数,在加入雷达测量风向角度误差范围内,输入每组风向和风速的估计值,输出代价函数最小的结果作为风向和风速的最优估计值。将传统加权融合算法与本发明进行了效果对比,并对比加入传感器噪声后的效果,结果证明,相比于传统的最优加权融合方法,该方法大大提高了校正后的风速、风向测量精度,比传统的方法融合效果更好。
本发明提供了一种基于雷达测风组合策略的多传感器融合风速风向估计方法,相较传统的最优加权融合方法,本发明提高了校正后的风速、风向测量精度。本发明是一种基于航海雷达图像测量的风向参数组合策略的多传感器融合风速风向估计方法,该方法将通过航海雷达图像反演所测得的风向参数作为先验条件,在该风向参数组合策略下,利用舰面风速风向仪的实测值进行舰面稳态自由来流的风速风向参数估计。
本发明考虑除了使用传感器测量的风速、风向数据外,还加入了雷达图像反演的风向数据作为参考项,从而实现了基于雷达测风组合策略的多传感器最优加权算法,估计真实自由来流的风参数。
附图说明
图1物理模型结构图;
图2风向角度与传感器位置示意图;
图3a-图3b网格划分图与计算模型的局部网格划分图;
图4本发明的流程图;
图5a-图5b验证实例中值滤波前后雷达图;
图6验证实例提取出的静态特征雷达图像;
图7极坐标图像区域划分和区域选取;
图8验证实例中极坐标插值笛卡尔坐标示意图;
图9笛卡尔坐标下的风场信息图像;
图10验证实例中得到风条纹功率谱图;
图11a-图11b传统最优加权算法效果图;
图12a-图12b基于雷达测风组合的最优加权算法效果图;
图13a-图13b加入传感器噪声情况下的基于雷达测风组合的最优加权算法效果图。
具体实施方式
下面举例对本发明作进一步描述。
具体实施方式一:
步骤1,获取传感器权值:
离线展开获取风速风向测量值实验,对舰船的桅杆部分建模后的结构进行网格划分以及边界条件、入口条件的设定,根据设定好的物理模型与方程通过CFD计算软件计算出在不同传感器布置位置处的风向、风速值,将得到的不同传感器布置位置处的传感器风速风向值作为传感器融合的基本数据,在规定风向范围内对传感器风速值以固定步长逐一进行分组,对于多个传感器的条件下,对应求出分组后的每组传感器所占权重。
步骤2,引入雷达风向进行数据筛选:
根据雷达测量的自由来流风向值作为相对风向的参考约束,对该风向角度范围约束下的传感器信息进行校验,筛选对应分组的传感器权值,进而估计出不同传感器权值下的风向、风速值。
步骤3,多传感器融合:
将得到的每组风速分量矢量相加,计算出每组风向和风速的估计值。引入代价函数,在加入雷达测量风向角度误差范围内,输入每组风向和风速的估计值,输出代价函数最小的结果作为风向和风速的最优估计值。
具体实施方式二:
在具体实施方式一的基础上,步骤1进一步包括以下步骤:
步骤1.1,离线通过CFD计算软件开展获取不同传感器布置位置处的传感器风速风向测量值的实验,对舰船的桅杆部分建模,对桅杆模型结构进行网格划分,设定整个模型的边界条件,根据设定好的物理模型与方程计算出在不同传感器布置位置处的风向、风速值;
步骤1.2,根据步骤1.1得到不同传感器布置位置处的传感器风速风向值作为传感器融合的基本数据,将自由来流真风速vref与每个传感器布置位置风速vsensor沿X轴、Y轴进行正交分解;
步骤1.3,根据步骤1.2确定的真风速与所有传感器风速分量数据,在0°~180°的风向范围内对传感器风速值以步长δ逐一进行分组;
步骤1.4,对步骤1.3得到的所有分组,计算出每组角度下每个传感器对应的权值,计算方法如下:
设每个传感器之间是相互独立的,且其测量值是真实值的无偏估计,设各个传感器的X轴风速分量为Xxi,Y轴风速分量为Xyi,其中i=1,2,…,n,n为传感器个数。ωxi代表第i个传感器X轴风速分量的权重值,ωyi代表第i个传感器Y轴风速分量的权重值,为X轴分量最后测量的估计值,为Y轴分量最后测量的估计值。
对上式等号的左右两侧分别取期望,如下式所示:
由于估计值与每个传感器X轴风速分量测量值Xxi、Y轴风速分量测量值Xyi均视为X轴风速分量、Y轴分速分量真实值Xx、Xy的无偏估计,上式等号左右两侧的期望值相同,所以可以推算出传感器权重系数的和为1。
根据上述公式,可以看出,最后估计X轴风速分量、Y轴风速分量的方差值与每个传感器的权重系数、测量方差有关,若要找到最优的估计值,只需要分别找到方差值最小值,又由于权重系数满足权重系数的和为1,通过拉格朗日乘数法,找到最优的方差值,公式如下:
对上式的变量ωxi、ωyi与λ分别进行求导,使导数为0,求出估计值方差的极小值。
权重系数的求解公式为:
则传感器的估计值为:
上式中的每个传感器的测量方差可以通过测量的数据与参考的数据计算得到,从而求出权重系数。
具体实施方式三:
在具体实施方式二的基础上,步骤2进一步包括以下步骤:
步骤2.1,引入雷达的风向测量数据θradar对传感器测量数据进行风向角度范围约束;
步骤2.2,雷达测量的风向数据存在误差,在加入雷达测量出的风向误差约束角度范围θradar±θ0内,在完成步骤1的分组后查找相对应角度范围下每个传感器的权值分布并存储;
步骤2.3,根据2.2查找得到的每个传感器对应的权值,估计每组角度范围下风速风向估计值。
具体实施方式四:
在具体实施方式三的基础上,步骤3进一步包括以下步骤:
步骤3.2,根据步骤3.1得到每组风速风向估计值,引入了一个代价函数,在加入雷达测量出的风向误差约束角度范围θradar±θ0内,取其中代价函数最小的那一组风速作为输出,代价函数定义为:
下面将结合附图,通过验证实例,对本发明的一种基于雷达测风组合的多传感器融合风速风向估计方法作进一步的详细说明,并验证其技术效果。本发明流程图见附图4,具体可以分为以下几步,第一步为获取多传感器权值,第二步为引入雷达风向进行数据筛选,第三步为多传感器融合。
本发明实施所用的风传感器型号是Model-05103,风向偏差为3°,风速偏差为0.3m/s,主要的技术参数如表一所示:
表一Model-05103型风传感器参数
结合附图1~13,本发明具体实施步骤为:
第一步为获取传感器权值:
步骤1.1,离线通过CFD计算软件开展获取不同传感器布置位置处的传感器风速风向测量值的实验,对舰船的桅杆部分使用ProE软件进行建模后得到简化的物理模型,见附图1。本例舰船安装测风传感器的位置见附图2,数量为8个。对桅杆模型结构采用ICEM进行网格划分,见附图3a图3b,其中图3a为整体的计算区域网格划分图,图3b为计算模型的局部网格划分图。设定整个模型计算的入口条件为给定的自由来流风速风向参数,计算区域的出口条件设置为压力出口,气流视为粘性绕流,计算区域为无滑移壁面,其他区域为自由滑移壁面。根据设定好的物理模型与方程计算出空间各个划分网格位置的风向、风速值,其中方程包括:
1)质量守恒方程
质量守恒定律作为自然界中三大定律之一,所有流动问题都必须遵循该定律,其质量守恒方程如下所示:
式中:ρ——密度
t——时间
u——速度矢量在X轴方向的分量
v——速度矢量在Y轴方向的分量
w——速度矢量在Z轴方向的分量
本文研究的是稳态风场在舰面绕流的问题,因此本文中研究的流体为稳定、不可压流体,故密度ρ为常数。因此上式可改写为:
2)动量守恒方程
所有流动问题都必须遵循牛顿第二定律,其动量守恒方程如下所示:
式中:ρ——流体密度
t——时间
u——速度矢量
u——速度矢量在X轴方向的分量
v——速度矢量在Y轴方向的分量
w——速度矢量在Z轴方向的分量
p——流体微元体上的压力
τxy等——粘性应力τ的分量
Fx等——微元体上的体力
3)能量守恒方程
所有含有热交换的流动问题都必须遵循该定律,其方程如下所示:
式中:ρ——流体密度
T——温度
u——速度矢量在X轴方向的分量
v——速度矢量在Y轴方向的分量
w——速度矢量在Z轴方向的分量
k——流体的传热系数
cp——比热容
ST——流体内热源及其机械能转换为热能的能量之和
4)湍流模型控制方程
湍流是一种非常复杂的、由粘性力引起的三维非稳态、带旋转的不规则流动。其中雷诺数是判断流体为层流或湍流的一个重要依据,其公式如下所示:
式中:Re——雷诺数
ρ——流体密度
v——流体特征速度
L——流体特征长度
u——动力粘性系数
5)RNG k-ε模型
与标准k-ε模型相比,RNG k-ε模型有很多改进:第一,RNG k-ε模型考虑了湍流涡旋,且将标准模型中用户给定的Prandtl常数更改为了一个解析公式。第二,RNG k-ε模型中增添了低雷诺数流动粘性解析公式。
本例采用RNG k-ε模型来对计算模型周围的气流场进行数值模拟,其运输方程见下式:
式中:ρ——流体密度
k——湍动能
ε——耗散率
ui——速度矢量在xi方向的分量
μ——流体粘性系数
μt——湍动粘度
Gk——湍动能k的产生项
σk——Prandtl数,取值1.0
σε——Prandkl数,取值1.3
C1ε——经验常数,取值1.44
C2ε——经验常数,取值1.92
步骤1.2,根据步骤1.1得到8个不同传感器布置位置处的传感器风速风向值作为传感器融合的基本数据,将自由来流真风速vref与8个传感器布置位置风速vsensor沿X轴、Y轴进行正交分解。
步骤1.3,根据步骤1.2确定的真风速与所有传感器风速分量数据,对0°~360°的风向范围以每5°为步长逐一进行分组,本例的分组数量为72。
步骤1.4,对步骤1.3得到的所有72个分组,计算出每组下每个传感器对应的权值,计算方法如下:
设每个传感器之间是相互独立的,且其测量值是真实值的无偏估计,设各个传感器的X轴风速分量为Xxi,Y轴风速分量为Xyi,其中i=1,2,…,n,n为传感器个数、本例为8。ωxi代表第i个传感器X轴风速分量的权重值,ωyi代表第i个传感器Y轴风速分量的权重值,为X轴分量最后测量的估计值,为Y轴分量最后测量的估计值。
对上式等号的左右两侧分别取期望,如下式所示:
由于估计值与每个传感器X轴风速分量测量值Xxi、Y轴风速分量测量值Xyi均视为X轴风速分量、Y轴分速分量真实值Xx、Xy的无偏估计,上式等号左右两侧的期望值相同,所以可以推算出传感器权重系数的和为1。
根据上述公式,可以看出,最后估计X轴风速分量、Y轴风速分量的方差值与每个传感器的权重系数、测量方差有关,若要找到最优的估计值,只需要分别找到方差值最小值,又由于权重系数满足权重系数的和为1,通过拉格朗日乘数法,找到最优的方差值,公式如下:
对上式的变量ωxi、ωyi与λ分别进行求导,使导数为0,求出估计值方差的极小值。
权重系数的求解公式为:
则传感器的估计值为:
上式中的每个传感器的测量方差可以通过测量的数据与参考的数据计算得到,从而求出权重系数。
第二步为引入雷达风向进行数据筛选。包括以下步骤:
步骤2.1,引入雷达的风向测量数据θradar对传感器测量数据进行风向角度范围约束,本例应用导航雷达图像反演海面风向θradar,步骤如下:
第一步,对实测雷达图像序列应用3×3模板的2-D非线性平滑中值滤波,抑制同频对海面风场研究的影响。
其中g(s,t)为雷达图像像元点在极坐标位置(s,t)的图像回波强度值;f′(r,θ)为滤波后图像在极坐标位置(r,θ)的灰度值;N(r,θ)为中心点在(r,θ)处的像元点,(s,t)取以(r,θ)为中心的8个点。
将3×3模板中值滤波器的模板中心与极坐标图像的某个像元位置重合,将其与周围8个相邻像元点的回波强度值进行排列,取中间的回波强度值赋给中心位置的像元,模板遍历全幅雷达图像获得中值滤波后的图像序列,附图5a图5b为本例中值滤波前后雷达图像。
第二步,极坐标图像归一化,使图像序列中的每幅图像的线数固定,从而达到像元点数和像元位置固定。
本例应用的X波段导航雷达图像是由大约3600条射线状线(3600个角度值)和600个同心圆组成,间隔角度大约是0.1°,由于雷达工作时脉冲数量不定和天线旋转时受到各种外界环境的干扰,所以线数是不固定的。为了使图像的线数固定即间隔角度固定使用极坐标图像归一化算法,步骤如下:
①建立一个1800条线(1800个角度值)和600个同心圆组成的极坐标新图像;
②将原图像角度值与新角度值相等,或大于原图像角度值的第一条线上的雷达回波强度赋到新图像的线上;
③不断重复②直到新图像上的N条线上都具有原雷达回波强度值,从而获得归一化的新极坐标图像。
第三步,应用雷达图像序列建立全局低通滤波器,滤除海浪等高频信号,仅保留静态或准静态频率特征信号,主要为风场导致的风条纹信号。具体步骤如下:
①根据海面风条纹属于静态频率特性,在单幅雷达图像上无法获得风条纹,对归一化的雷达图像序列在同一位置的像元点上按时间(80s)进行积分平均,构建一个图像空间全局低通滤波器:
式中f(θ,r)为风条纹图像,f′(θ,r,t)雷达图像序列时间为t的单幅雷达图像,Nt为时间序列;
②构建1800条线和600个同心圆的二维极坐标图像;
③将获得的1800×600个低通滤波后的新像元灰度值按位置赋予新构建的二维极坐标图像上,得到二维极坐标海面静态特征图像,如附图6所示。
第四步,极坐标图像区域划分,将二维极坐标图像的像元点按照128*128的数量将图像分成14个小区域,如附图7所示。计算每个小区域的海面风向,最终海面平均风向为每个小区域风向的平均值。图7中选框区域的中心角度为203°,像元灰度值为f(ri,θj)(i,j=1,2,…,128)。
第五步,图7中选框区域的极坐标应用最近点插值为笛卡尔坐标,主要是将原区域的雷达回波强度值应用位置最近的方式,将值赋予到已建立的笛卡尔坐标下,具体步骤如下:
②选取区域所在的中心角度θc和半径rc,计算选取极坐标图像中中心点在全幅极坐标图像中的笛卡尔坐标位置(xc,yc),以中心点所在位置为中心计算选取区域内的像元点在全幅极坐标图像中的笛卡尔坐标的x-y轴分量(xi,yi);
③应用最近点插值,当两个坐标(xi,yi)和位置最近时将f(ri,θj)赋到上,将选取的极坐标区域图像上的灰度值都插值到新生成的笛卡尔坐标上,横纵坐标为像元点所在位置,这样就得到笛卡尔坐标下的二维静态特征图像,插值方法如附图8。图9为图7中红框区域插值结果图,横纵坐标为128*128像元点所在位置。
第六步,获取海面静态特征二维能量谱图,对二维图像进行二维离散快速傅里叶变换(2D FFT)及可得到图像的能量谱图。本例对二维笛卡尔坐标下的海面静态特征图像应用二维离散快速傅里叶变换(2D FFT)获得海面静态特征图像能量谱,数学模型为:
F(kx,ky)为海面静态特征图像的傅里叶系数,其中复数指数项可以写成:
其中,
(kx,ky)为f(x,y)在能量谱T的坐标,dx、dy为雷达图像分辨率,本例中公式dx=dy=7.5m。
第七步,对静态特征图像能量谱应用波数能量谱尺度分离,将风条纹信号的谱能量从静态特征图像能量谱中分离出来,主要是依据风条纹波长L和频域坐标尺度k的关系来提取。由于风条纹的尺度为100~500m,计算得到风条纹能量谱波数下限kd:
风条纹能量谱波数上限为kt:
根据能量谱波数上下限设计能量谱带通滤波器,将风条纹信号能量谱提取出来,数学模型为:
I(kx,ky)为风条纹能量谱,图10波数尺度分离后I(kx,ky)在(kx,ky)坐标下得到的风条纹等高线能量谱图,虚线为(kx,ky)坐标下的波长范围。
第八步,获得海面方向信息。由于傅里叶变换的周期性,得到的风条纹能量谱关于一三或二四象限镜像对称,两个能量集中区域连线方向为风条纹的垂直向,而连线的垂直方向则为风条纹平行方向,由于风条纹方向与风向平行,即可得到海面方向。
由于选取区域和船艏向的影响,导致计算的风向需要校正后才能得到相对正北方向的海表面风向,校正公式为:
Nw=|θc|+|α|+|β|
Nw为相对正北向的海面风向;θc为选取区域的中心角,选取为110°;β为笛卡尔坐标下计算得到的风条纹方向;α为船艏向,为93°。图9中得到Nw为33°或213°。
由于单点风向标测得的海面风向为36°属于第一个象限,将计算出来两个方向与第一个象限比较,保留第一个象限角度,去除第三象限的角度,从而获得的海面风向为33°。
θradar为最终计算得到的海面主风向,区域划分的越多则平均风向精度会有所增加,但要保证在进行二维离散快速傅里叶变换时像元的数量,不然在谱图像上无法获得风条纹尺度的波数k,这样就无法获得准确的海面风向信息。
步骤2.2,本例考虑雷达测量风向存在±20°的误差,所以在雷达测量出的风向约束角度θradar±20°范围内,对风向数据每5°作为一个范围,将该风向下所有风速的数据做多传感器的最优加权融合,查找到传感器对应的权值分布如下式所示:
式中,i代表传感器的序号,n为传感器的个数、本例为8,f(θradar)j代表查找得到的j组传感器对应权值,j为根据雷达测量风向分组序号,j=1,2,…,9,本例为9。wxij为第j组第i个传感器的在雷达测量角度θradar下的X轴分量的最优加权权值,wyij为第j组第i个传感器在雷达测量角度θradar下的Y轴分量的最优加权权值。
步骤2.3,根据2.2查找得到的9组每个传感器对应的不同权值,估计每组角度范围下风速风向估计值。
第三步为多传感器融合。包括以下步骤:
步骤3.2,根据步骤3.1得到9组风速风向估计值,引入了一个代价函数,在θradar±20°范围内,取其中代价函数最小的那一组风速作为输出,代价函数定义为:
通过仿真结果表明,本发明相较于传统加权融合算法的结果更好,因此利用传统加权融合与本发明进行对比、加入传感器噪声条件下传统加权融合算法与加入传感器噪声条件下本发明进行了对比,对于估计的风向使用风向绝对误差MAE作为标志,对风速的估计则选择风速相对误差MRE作为标志,以验证算法性能。
附图11a-图11b为传统最优加权融合算法效果图,图中共描绘了五种工况下的风速相对误差和风向绝对误差图。其中工况A的标准风速为3m/s,工况B的标准风速为6m/s,工况C的标准风速为9m/s,工况D的标准风速为12m/s,工况E的标准风速为15m/s。横坐标代表真实风向值。图11a为风速相对误差,图中可以看到风向角度范围在0°到180°的范围,整个误差情况总体较高,大部分的风速相对误差值分布在5%到15%之间。图11b为风向绝对误差,同样可以观察到,整个图中的大部分风向角度下的风向绝对误差为5°以上,只有中间一小部分的风向绝对误差较低。
附图12a-图12b为基于雷达测风组合的最优加权融合算法效果图,图12a为风速相对误差图,可以观察到,此时的风速相对误差值大部分分布在5%以下,只有少数部分的点超过5%。图12b为风向绝对误差,图中的风向误差分布大部分在0°到3°之间,只有中间一部分误差较大,但大部分在5°以下,误差最大值只有7°。
表二为传统加权融合算法和基于雷达测风组合的多传感器融合风速风向估计算法误差校正误差对比,从表中可以看出,加入了雷达测风组合的多传感器融合风速风向估计算法的风向误差均值比传统的加权融合算法误差均值降低了3.3°,而风速相对误差则整体提升了4.01%。
表二传统加权融合算法和基于雷达测风组合的多传感器融合风速风向估计算法误差对比
附图13a-图13b为加入传感器噪声情况下的基于雷达测风组合的多传感器融合风速风向估计算法效果图,图图13a为风速相对误差图,可以看到,风速相对误差的分布范围集中在0%到10%之间,并且绝大部分在5%之下。图13b为风向绝对误差,风向相对误差的分布范围则整体在0°到8°之间,并且绝大部分的风向相对误差在5°之下。
表三为传感器噪声下的传统算法与基于雷达测风组合的多传感器融合风速风向估计算法误差对比,从表中可以看出,改进的基于雷达测风组合的多传感器融合风速风向估计算法的风向绝对误差为2.51°,比传统的最优加权融合算法精确了2.45°。同样地,风速相对误差值比传统的方法精确度提升了0.89%。
表三传感器噪声下的传统算法和基于雷达测风组合的多传感器融合风速风向估计算法误差对比
本发明所提出的基于测风组合的多传感器融合风速风向估计方法,相比于传统的最优加权融合方法,该方法大大提高了校正后的风速、风向测量精度,比传统的方法融合效果更好。
Claims (3)
1.一种基于雷达测风组合策略的多传感器融合风速风向估计方法,其特征是:
步骤1,获取传感器权值:
离线展开获取风速风向测量值实验,对桅杆部分建模后的结构进行网格划分以及边界条件、入口条件的设定,通过CFD计算软件计算出在不同传感器布置位置处的风向、风速值,将得到的不同传感器布置位置处的传感器风速风向值作为传感器融合的基本数据,在规定风向范围内对传感器风速值以固定步长逐一进行分组,对于多个传感器的条件下,对应求出分组后的每组传感器所占权值;
步骤1.1,离线通过CFD计算软件开展获取不同传感器布置位置处的传感器风速风向测量值的实验,对桅杆部分建模,对桅杆模型结构进行网格划分,设定整个模型的边界条件,根据设定好的物理模型与方程计算出在不同传感器布置位置处的风向、风速值;
步骤1.2,根据步骤1.1得到不同传感器布置位置处的传感器风速风向值作为传感器融合的基本数据,将自由来流真风速vref与每个传感器布置位置风速vsensor沿X轴、Y轴进行正交分解;
步骤1.3,根据步骤1.2确定的真风速与所有传感器风速分量数据,在0°~180°的风向范围内对传感器风速值以步长δ逐一进行分组;
步骤1.4,对步骤1.3得到的所有分组,计算出每组角度下每个传感器对应的权值,计算方法如下:
设每个传感器之间是相互独立的,且其测量值是真实值的无偏估计,设各个传感器的X轴风速分量为Xxi,Y轴风速分量为Xyi,其中i=1,2,…,n,n为传感器个数;ωxi代表第i个传感器X轴风速分量的权重值,ωyi代表第i个传感器Y轴风速分量的权重值,为X轴分量最后测量的估计值,为Y轴分量最后测量的估计值;
对上式等号的左右两侧分别取期望,表示如下:
由于估计值与每个传感器X轴风速分量测量值Xxi、Y轴风速分量测量值Xyi均视为X轴风速分量、Y轴分速分量真实值Xx、Xy的无偏估计,上式等号左右两侧的期望值相同,所以推算出传感器权重系数的和为1,
根据上述公式,最后估计X轴风速分量、Y轴风速分量的方差值与每个传感器的权重系数、测量方差有关,若要找到最优的估计值,只需要分别找到方差值最小值,又由于权重系数满足权重系数的和为1,通过拉格朗日乘数法,找到最优的方差值、公式如下:
对上式的变量ωxi、ωyi与λ分别进行求导,使导数为0,求出估计值方差的极小值,
权重系数的求解公式为:
则传感器的估计值为:
上式中的每个传感器的测量方差通过测量的数据与参考的数据计算得到,从而求出权重系数;
步骤2,引入雷达风向进行数据筛选:
根据雷达测量的自由来流风向值作为相对风向的参考约束,对风向角度范围约束下的传感器信息进行校验,筛选对应分组的传感器权值,进而估计出不同传感器权值下的风向、风速值;
步骤3,多传感器融合:
将得到的每组风速分量矢量相加,计算出每组风向和风速的估计值,引入代价函数,在加入雷达测量风向角度误差范围内,输入每组风向和风速的估计值,输出代价函数最小的结果作为风向和风速的最优估计值。
2.根据权利要求1所述的基于雷达测风组合策略的多传感器融合风速风向估计方法,其特征是步骤2进一步包括:
步骤2.1,引入雷达的风向测量数据θradar对传感器测量数据进行风向角度范围约束;
步骤2.2,雷达测量的风向数据存在误差,在加入雷达测量出的风向误差约束角度范围θradar±θ0内,在完成步骤1的分组后查找相对应角度范围下每个传感器的权值分布并存储;
步骤2.3,根据2.2查找得到的每个传感器对应的权值,估计每组角度范围下风速风向估计值。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010481799.6A CN111693999B (zh) | 2020-05-27 | 2020-05-27 | 基于雷达测风组合策略的多传感器融合风速风向估计方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010481799.6A CN111693999B (zh) | 2020-05-27 | 2020-05-27 | 基于雷达测风组合策略的多传感器融合风速风向估计方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111693999A CN111693999A (zh) | 2020-09-22 |
CN111693999B true CN111693999B (zh) | 2023-05-05 |
Family
ID=72478747
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010481799.6A Active CN111693999B (zh) | 2020-05-27 | 2020-05-27 | 基于雷达测风组合策略的多传感器融合风速风向估计方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111693999B (zh) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112258850A (zh) * | 2020-10-16 | 2021-01-22 | 南京莱斯网信技术研究院有限公司 | 一种车路协同系统边缘侧多传感器数据融合系统 |
CN113253232B (zh) * | 2021-05-31 | 2021-09-17 | 中国人民解放军国防科技大学 | 结合机器学习与三维变分同化的二维风场反演方法和装置 |
CN114494894B (zh) * | 2022-04-18 | 2022-07-22 | 中国海洋大学 | 海洋黑涡自动识别与关键参数反演方法、装置和电子设备 |
CN116953290B (zh) * | 2023-09-18 | 2024-01-12 | 浙江力夫传感技术有限公司 | 风速变送器探杆角度调整方法及风速变送器 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105975763A (zh) * | 2016-04-29 | 2016-09-28 | 国家卫星海洋应用中心 | 一种多源海面风场的融合方法和装置 |
CN106199605A (zh) * | 2016-07-06 | 2016-12-07 | 西南技术物理研究所 | 风场误差修正方法 |
CN109100717A (zh) * | 2018-06-11 | 2018-12-28 | 广州地理研究所 | 一种多源微波遥感海面风场数据融合方法及其装置 |
CN109541963A (zh) * | 2018-11-12 | 2019-03-29 | 北京应用气象研究所 | 一种基于侧滑角信息的无人机测风建模技术 |
CN110286390A (zh) * | 2019-06-11 | 2019-09-27 | 中国科学院合肥物质科学研究院 | 一种指定路径风速测量方法、装置及测风雷达标定方法 |
EP3588128A1 (en) * | 2018-06-26 | 2020-01-01 | Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. | Method for detection and height and azimuth estimation of objects in a scene by radar processing using sparse reconstruction with coherent and incoherent arrays |
CN111025295A (zh) * | 2019-11-22 | 2020-04-17 | 青岛海狮网络科技有限公司 | 一种基于岸基雷达的多船协同感知数据融合系统及方法 |
-
2020
- 2020-05-27 CN CN202010481799.6A patent/CN111693999B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105975763A (zh) * | 2016-04-29 | 2016-09-28 | 国家卫星海洋应用中心 | 一种多源海面风场的融合方法和装置 |
CN106199605A (zh) * | 2016-07-06 | 2016-12-07 | 西南技术物理研究所 | 风场误差修正方法 |
CN109100717A (zh) * | 2018-06-11 | 2018-12-28 | 广州地理研究所 | 一种多源微波遥感海面风场数据融合方法及其装置 |
EP3588128A1 (en) * | 2018-06-26 | 2020-01-01 | Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. | Method for detection and height and azimuth estimation of objects in a scene by radar processing using sparse reconstruction with coherent and incoherent arrays |
CN109541963A (zh) * | 2018-11-12 | 2019-03-29 | 北京应用气象研究所 | 一种基于侧滑角信息的无人机测风建模技术 |
CN110286390A (zh) * | 2019-06-11 | 2019-09-27 | 中国科学院合肥物质科学研究院 | 一种指定路径风速测量方法、装置及测风雷达标定方法 |
CN111025295A (zh) * | 2019-11-22 | 2020-04-17 | 青岛海狮网络科技有限公司 | 一种基于岸基雷达的多船协同感知数据融合系统及方法 |
Also Published As
Publication number | Publication date |
---|---|
CN111693999A (zh) | 2020-09-22 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111693999B (zh) | 基于雷达测风组合策略的多传感器融合风速风向估计方法 | |
WO2021218424A1 (zh) | 基于rbf神经网络的航海雷达图像反演海面风速方法 | |
Diebold et al. | Flow over hills: a large-eddy simulation of the Bolund case | |
CN109345875B (zh) | 一种提高船舶自动识别系统测量精度的估计方法 | |
CN112711899B (zh) | 一种蒸发波导高度的融合预测方法 | |
Dobson et al. | Air-sea interaction: instruments and methods | |
CN110184885B (zh) | 一种基于智能手机测试路面平整度的方法 | |
CN107632964A (zh) | 一种平面地磁异常场向下延拓递归余弦变换法 | |
Huang et al. | Wave height estimation from X-band nautical radar images using temporal convolutional network | |
Eliassen et al. | Coherence of turbulent wind under neutral wind conditions at FINO1 | |
EP4127457A1 (en) | System and method for wind flow turbulence measurement by lidar in a complex terrain | |
CN109739088A (zh) | 一种无人船有限时间收敛状态观测器及其设计方法 | |
Fuentes-Pérez et al. | Underwater vehicle speedometry using differential pressure sensors: Preliminary results | |
Thornhill et al. | Ship anemometer bias management | |
Mateer et al. | Skin-friction measurements and calculations on a lifting airfoil | |
Hansen et al. | Vortex-containing wakes of surface obstacles | |
Yao et al. | Optimal design of hemispherical 7-hole probe tip with perpendicular holes | |
Redondo et al. | Eddy measurements, coastal turbulence and statistics in the gulf of Lions | |
CN113971350B (zh) | 风速场拟合补缺方法及装置、介质 | |
CN115683544A (zh) | 无人机旋翼扰动校正方法及装置 | |
Bykov et al. | Objective analysis of the structure of three-dimensional atmospheric fronts | |
CN117454807B (zh) | 基于光学设备防护罩的多尺度cfd数值模拟方法 | |
CN117454614A (zh) | 一种基于模糊数学的加权融合风速及风向估计方法 | |
CN108562896A (zh) | 一种基于三维正压浅海大陆架模型的深层海流反演方法 | |
Gehlert et al. | Unsteady vorticity shedding from a circular cylinder: surging, spinning and gust encounters |
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 |