CN111868533A - 风流感测系统和用于确定风流的速度场的方法 - Google Patents

风流感测系统和用于确定风流的速度场的方法 Download PDF

Info

Publication number
CN111868533A
CN111868533A CN201880091302.0A CN201880091302A CN111868533A CN 111868533 A CN111868533 A CN 111868533A CN 201880091302 A CN201880091302 A CN 201880091302A CN 111868533 A CN111868533 A CN 111868533A
Authority
CN
China
Prior art keywords
velocity
altitude
velocity field
cfd
altitudes
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.)
Granted
Application number
CN201880091302.0A
Other languages
English (en)
Other versions
CN111868533B (zh
Inventor
S·纳比
P·格罗弗
M·德纳斯
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.)
Mitsubishi Electric Corp
Original Assignee
Mitsubishi Electric Corp
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 Mitsubishi Electric Corp filed Critical Mitsubishi Electric Corp
Publication of CN111868533A publication Critical patent/CN111868533A/zh
Application granted granted Critical
Publication of CN111868533B publication Critical patent/CN111868533B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01P—MEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
    • G01P5/00—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft
    • G01P5/001—Full-field flow measurement, e.g. determining flow velocity and direction in a whole region at the same time, flow visualisation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01P—MEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
    • G01P5/00—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft
    • G01P5/26—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft by measuring the direct influence of the streaming fluid on the properties of a detecting optical wave
    • 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
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88—Lidar systems specially adapted for specific applications
    • G01S17/95—Lidar 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/48—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00
    • G01S7/4808—Evaluating distance, position or velocity data
    • 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
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/02—Systems using the reflection of electromagnetic waves other than radio waves
    • G01S17/50—Systems of measurement based on relative movement of target
    • G01S17/58—Velocity or trajectory determination systems; Sense-of-movement determination systems
    • 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)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Electromagnetism (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Environmental & Geological Engineering (AREA)
  • Multimedia (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Ecology (AREA)
  • Environmental Sciences (AREA)
  • Biodiversity & Conservation Biology (AREA)
  • Atmospheric Sciences (AREA)
  • Indicating Or Recording The Presence, Absence, Or Direction Of Movement (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

风流感测系统通过利用减小误差的加权组合的成本函数的操作参数来模拟风流的计算流体动力学(CFD)，确定在每个海拔高度处的速度场的第一近似值；从速度场的第一近似值确定在每个海拔高度处的垂直速度的水平导数；以及使用每个海拔高度的速度场、径向速度的测量结果在三维轴上的投影、以及相应速度场的垂直速度的水平导数之间的几何关系，确定速度场的第二近似值。在CFD的成本函数中，每个误差对应于海拔高度之一，并且包括在相应海拔高度处在视线点处的测量速度与通过CFD针对相应海拔高度所模拟的在视线点处的模拟速度之间的差。加权组合中的至少一些误差具有不同权重。

Description

风流感测系统和用于确定风流的速度场的方法
技术领域
本发明一般涉及遥感，并且更具体地，涉及感测通过复杂地形的风的流动。
背景技术
测量大气或风的运动对于许多应用，尤其是在气象学中以及用于诸如机场和风电场之类的场所的监视和表征来说是重要的。测量气团在广泛海拔高度范围内或在与大体积相对应的区域中的位移通常很有用。这种大体积位移的感测用诸如杯型风速仪之类的风速仪来执行可能是不切实际的，并且需要能够进行远程测量的遥感仪器。这些仪器尤其包括雷达、激光雷达(LiDAR)和声雷达(SODAR)。雷达和激光雷达系统分别使用高频和光学频率范围内的电磁波。声雷达系统使用声波。
例如，用这种仪器测量气团的运动通常如下进行。该设备沿待测区域中的传输轴连续地或作为脉冲来传输一个或更多个波束(声波和/或电磁波)。沿着不同传输轴的传输可以是同时的或顺序的。波束在大气中特别是由于遇到异质性(气溶胶、颗粒、电磁波折射率的变化或声波声阻抗的变化)而受到散射效应。当它们散布在气团或运动颗粒中时，这些波束也会由于Doppler(多普勒)效应而发生频移。
背向散射波束由根据测量轴定向的一个或更多个接收器检测。这些接收器在沿其测量轴的方向上检测被大气散射的波。然后可以例如通过测量飞行时间的方法或通过干涉法测量相移的方法，来计算沿着检测器的测量轴发生散射的距离。气团或颗粒沿着测量轴的径向速度也可以通过利用Doppler效应测量波的频移来获得。该测得的径向速度对应于散射位置的速度矢量在检测器的测量轴上的投影。
遥感仪器，尤其是包括适于测量大气下层中的风的特征的激光雷达系统，通常是单静态类型的。这表示使用相同的光学器件或相同的天线(声学或电磁)来发送和接收信号。探测到的体积通常沿着其顶点位于仪器的光学元件或天线的高度的圆锥分布。沿着圆锥的仪器的每个脉冲束对颗粒沿着与传输轴重合的测量轴运动的径向速度进行测量。这样就获得了风的径向速度的测量值，其代表了风矢量在波束传播轴上的投影。
然后，基于风的径向速度的测量结果，计算所有感兴趣体积中的风矢量，即速度场。在现有装置中，通常使用纯几何模型实施此计算。这些模型的缺点在于它们基于有时是不太现实的假设，特别是在整个样本测量期间风的空间和时间均匀性。根据该假设，在给定的海拔高度下，仪器探测到的大气的每个点处，风向矢量是相同的。
使用风的空间和时间均匀性，一些方法可以通过求解描述了风矢量与其沿着由径向速度的测量结果组成的测量结果的轴的投影之间的几何关系的、具有三个未知数的至少三个方程式的系统，从在至少三个不同方向中在一个且相同的海拔高度处测量到的径向速度的至少三个测量值，来计算在给定海拔高度下的风矢量的分量。使用几何计算的方法的示例是“Velocity Azimuth Display(速度方位角显示)”(VAD)方法，例如参见美国专利No.4,735,503。但是，此方法基于在给定海拔高度下风的空间均匀性的相同假设。
当在基本上平坦的地形(具有很少起伏或没有起伏的地形，或离岸的地形)上实施测量时，使用重构风矢量的几何技术的、用于风测量的大气遥感仪器允许精确测量风的平均速度。例如，凭借激光雷达系统，在10分钟内平均的测量值所得到的相对误差相对于由校准的杯型风速仪组成的参考而言在2％以下。
另一方面，当在诸如起伏或山区地形、森林覆盖的地形等的更复杂地形上实施测量时，确定水平和垂直速度的准确性会显著恶化。针对用这些相同的激光雷达系统实施的测量结果并且相对于校准的杯型风速仪，在复杂位置上观察到在10分钟计算出的平均值的相对误差为5％至10％的量级。
因此，当前实现几何模型的遥测装置不能在复杂的地形上足够精确地测量风的水平和/或垂直速度和方向。实际上，在复杂的地形上，在仪器探测到的大气的体积中，在给定的海拔高度下，风不再被认为是同质的。然而，在这些条件下，尤其是在风电场发展的背景下，精确测量风的特性是有用的。
一些方法(例如，US 2013/0282285)提出了使用几种不同的保真度数值模型以及优化方法，以匹配由激光雷达测量到的流速。但是，确定证明这些模型正确性的条件仍然是个难题。
因此，仍然需要一种适合于感测通过复杂地形的风流的系统和方法。
发明内容
一些实施方式基于以下认识：诸如激光雷达之类的遥感系统在视线(LOS)方向上测量盛行风的径向速度并估计3D风速场面临许多挑战。例如，使用圆锥图案的基于单个地面的激光雷达扫描的遥感系统会遭受所谓的“独眼困境”。这种现象是指以下事实：使用单个LOS测量结果无法执行任意3D速度场的精确重构。该问题可以至少部分通过以下假设来解决：在测量的圆锥内的每个水平面上，速度是均匀的和/或恒定的。换句话说，速度仅随着在垂直方向上移动而改变。这种均匀速度假设允许使用速度场与测量到的径向速度在速度场上的投影之间的几何关系来计算在每个水平面处的“平均”水平速度和垂直速度。
在远离风力发电机、建筑物以及茂密森林中超过一定高度的平坦地形上，均匀速度假设通常非常准确。但是，对于在诸如山丘或者大型建筑物或其它城市结构附近之类的复杂地形上方的风流，能够违反这种均匀速度假设。
一些实施方式基于以下认识：能够使用垂直速度的水平导数来校正在感测通过复杂地形的风流中的均匀速度假设。这种垂直速度的水平导数能够示出在给定海拔高度处垂直速度的方向和/或大小的变化。甚至示出了垂直速度的线性变化的一阶导数能够用于提高风流感测的准确度，这是因为均匀速度假设导致感测到的速度场的误差中的主要项。
一些实施方式基于以下认识：通过仅强制在流经复杂地形的风的重构速度场上的非压缩性，可能无法去除在均匀假设中出现的垂直速度的水平导数的偏差。具体地，试图将诸如空气的非压缩性之类的风流的物理特性强加于风流重构上。然而，一些实施方式基于以下认识：能够通过将非压缩性条件强加于感兴趣的域内的流体体积上，防止由均匀速度假设引起的误差的数学结构，以减少主要项。
计算流体动力学(CFD)是使用数值分析和数据结构来解决和分析涉及流体流的问题的流体力学的一个分支。使用计算机来执行模拟液体和气体与由边界条件限定的表面之间的相互作用所需的计算。一些实施方式基于一般的理解，即CFD能够用于从由激光雷达感测到的圆锥上的风的测量结果来估计风的速度场。然而，关于复杂地形上方的风流的操作参数(诸如边界条件)通常是未知的，并且这些操作参数的近似值会不期望地降低风流感测的准确度。
一些实施方式基于以下认识：尽管CFD近似值对于确定速度场可能不够准确，但是CFD近似值对于在给定海拔高度的垂直速度重构的水平导数的平均值来说能够足够精确，这进而能够用于校正由于均匀速度假设导致的偏差。为此，一些实施方式使用CFD近似值来确定垂直速度的水平导数，并且结合在期望海拔高度上的风流的径向速度测量结果使用垂直速度的水平导数来确定期望海拔高度的速度场。以这种方式，能够实现使用径向速度测量结果的速度场感测的目标准确度。
为此，一些实施方式使用风流的计算流体动力学(CFD)的模拟来产生速度场的第一近似值。然而，代替使用速度场的该第一近似值作为最终结果，实施方式从速度场的第一近似值确定在每个海拔高度处的垂直速度的水平导数，并且使用每个海拔高度的速度场、径向速度的测量结果在三维轴上的投影、以及相应速度场的垂直速度的水平导数之间的几何关系，确定速度场的第二近似值。在许多实验和测试中，各种实施方式表明，通过去除均匀速度假设的偏差而确定的速度场的第二近似值比使用风流的CFD所确定的速度场的第一近似值更准确。
一些实施方式基于以下认识：当将CFD用于垂直速度重构的水平导数的目的时，需要近似操作参数的特定集合。可以基于这样的参数对感兴趣的最终量(即，垂直速度的水平导数)的影响来选择这些操作参数。以这种方式，能够提高感兴趣的量(即，垂直速度的水平导数)的准确度。
例如，一些实施方式基于垂直速度的水平导数对这些操作参数的值的变化的灵敏度来选择操作参数。在一个实施方式中，在CFD模拟期间所近似的操作参数的基于目的的集合中选择灵敏度高于阈值的操作参数。
例如，一些实施方式选择诸如地形粗糙度、流入平均速度、流入湍流强度和大气稳定性条件之类的操作参数，以用于确定垂直速度的水平导数。
以这种方式，一些实施方式使CFD的未知操作参数适于CFD近似的目的。操作参数的这种适配减轻了计算负担，而没有降低感兴趣的量的CFD近似的准确度。
一些实施方式基于以下认识：可以使用基于直接伴随循环(DAL)的CFD框架来估计基于垂直速度的水平导数对操作参数的变化的灵敏度而选择的这些操作参数。该框架通过使估计在正向CFD模拟与可用测量结果之间的视线数据及其梯度中的误差的成本函数最小化、然后以迭代方式求解灵敏度(或伴随CFD)方程式，来同时校正用于共同目的的各种未知操作参数。用于共同目的的参数的灵敏度指示了基于DAL的CFD框架的收敛方向。一旦收敛，基于DAL的CFD框架能够产生操作参数的当前值以及针对这些操作参数而确定的速度场二者。
在DAL收敛之后，一些实施方式提取感兴趣的量(即，垂直速度梯度)，以使用在测量的圆锥上的激光雷达视线(LOS)校正在复杂地形上方的风速重建中的偏差误差。另外，与强制非压缩性来校正均匀速度假设相对比，一个实施方式强制空气的非压缩性来校正由LOS测量结果的稀疏性引起的某些类型的误差。该实施方式通过强制与在每个海拔高度的径向速度的测量结果一致的风流的非压缩性和正则化，产生在每个海拔高度的速度场的非恒定值的密集网格。
一些实施方式基于以下认识：当将CFD用于垂直速度重构的水平导数的目的时，需要将特定成本函数公式化以选择CFD的操作参数。能够基于改进最终感兴趣的量(即，垂直速度的水平导数)的估计值的目的来选择该成本函数。值得注意地是，该成本函数能够不同于当CFD的目标不同时所选择的成本函数。
具体地，一些实施方式基于以下认识：垂直速度的水平导数依据海拔高度对速度场的影响不同。为此，当CFD的目的是确定垂直速度的水平导数时，成本函数需要不同地考虑不同的海拔高度，例如，利用不同的权重。例如，在一些实施方式中，成本函数包括代表针对不同海拔高度的CFD的准确度的误差的加权组合。例如，误差包括对应于第一海拔高度的第一误差和对应于第二海拔高度的第二误差，并且在误差的加权组合中的第一误差的权重不同于在误差的加权组合中的第二误差的权重。
例如，一些实施方式选择成本函数，该成本函数是测量位置的海拔高度的递增函数。在一个实施方式中，由于均匀假设而导致的一阶误差项的分析形式被用于推导成本函数的形式。以这种方式，操作参数适于提高垂直速度的水平导数的确定准确度。
一些实施方式在成本函数的灵敏度的最大下降方向上迭代地确定操作参数。为此，多次执行(即，每次迭代执行一次)CFD模拟。为此，期望减少CFD模拟的计算成本。
一些实施方式基于以下认识：不是对每个新的测量数据集(例如，对于每个新的风向和/或新的地形)使用高保真CFD解，而是通过场反转和机器学习(FIML)方法能够使用高效和准确的低保真CFD模拟。低保真CFD模拟通过内部参数相关模型来近似流中的高频项。为此，实施方式使用包括垂直速度的水平导数的特征向量来应用机器学习，以学习内部参数对流特征的依赖性。该实施方式能够通过使用预处理步骤来减轻CFD近似的计算负担。
因此，一个实施方式公开了一种风流感测系统，该风流感测系统用于从在每个海拔高度处的径向速度的测量结果的集合来确定在一组不同海拔高度处的风流的速度场。该系统包括：输入接口，其接受针对每个海拔高度的在视线点处的径向速度的测量结果的集合；处理器，其被配置为：通过模拟减小误差的加权组合的成本函数的风流操作参数的计算流体动力学(CFD)，确定在每个海拔高度处的速度场的第一近似值，其中，每个误差对应于海拔高度之一，并且包括在相应海拔高度处在视线点处的测量速度与通过CFD针对相应海拔高度用操作参数的当前值所模拟的在视线点处的模拟速度之间的差，其中，误差包括对应于第一海拔高度的第一误差和对应于第二海拔高度的第二误差，其中，在误差的加权组合中的第一误差的权重不同于在误差的加权组合中的第二误差的权重；从每个海拔高度处的速度场的第一近似值，确定在每个海拔高度处的垂直速度的水平导数；以及使用每个海拔高度的速度场、径向速度的测量结果在三维轴上的投影、以及相应速度场的垂直速度的水平导数之间的几何关系，确定每个海拔高度处的速度场的第二近似值；以及输出接口，其呈现风流的速度场的第二近似值。
另一实施方式公开了一种风流感测方法，该风流感测方法用于从在每个海拔高度处的径向速度的测量结果的集合来确定在一组不同海拔高度处的风流的速度场，其中该方法使用与实现该方法的存储指令联接的处理器，其中指令在由处理器执行时实施该方法的步骤，步骤包括：接受针对每个海拔高度的在视线点处的径向速度的测量结果的集合；通过用减小误差的加权组合的成本函数的操作参数来模拟风流的计算流体动力学(CFD)，确定在每个海拔高度处的速度场的第一近似值，每个误差对应于海拔高度之一，并且包括在相应海拔高度处在视线点处的测量速度与通过CFD针对相应海拔高度用操作参数的当前值所模拟的在视线点处的模拟速度之间的差，其中误差包括对应于第一海拔高度的第一误差和对应于第二海拔高度的第二误差，其中在误差的加权组合中的第一误差的权重不同于在误差的加权组合中的第二误差的权重；从速度场的第一近似值，确定在每个海拔高度处的垂直速度的水平导数；以及使用每个海拔高度的速度场、径向速度的测量结果在三维轴上的投影、以及相应速度场的垂直速度的水平导数之间的几何关系，确定速度场的第二近似值；以及呈现风流的速度场的第二近似值。
又一个实施方式公开了一种非暂时性计算机可读存储介质，该非暂时性计算机可读存储介质上实现有用于执行以下方法的由处理器可执行的程序，该方法用于从在每个海拔高度处的径向速度的测量结果的集合来确定在一组不同海拔高度处的风流的速度场。该方法包括：接受针对每个海拔高度的在视线点处的径向速度的测量结果的集合；通过用减小误差的加权组合的成本函数的操作参数来模拟风流的计算流体动力学(CFD)，确定在每个海拔高度处的速度场的第一近似值，每个误差对应于海拔高度之一，并且包括在相应海拔高度处在视线点处的测量速度与通过CFD针对相应海拔高度用操作参数的当前值所模拟的在视线点处的模拟速度之间的差，其中误差包括对应于第一海拔高度的第一误差和对应于第二海拔高度的第二误差，其中在误差的加权组合中的第一误差的权重不同于在误差的加权组合中的第二误差的权重；从速度场的第一近似值，确定在每个海拔高度处的垂直速度的水平导数；以及使用每个海拔高度的速度场、径向速度的测量结果在三维轴上的投影、以及相应速度场的垂直速度的水平导数之间的几何关系，确定速度场的第二近似值，以及呈现风流的速度场的第二近似值。
附图说明
[图1]
图1示出了根据一些实施方式的风流感测系统的不同组件的硬件图。
[图2A]
图2A示出了根据一些实施方式的被配置为测量风流的径向速度的示例性遥感仪器的示意图。
[图2B]
图2B示出了由一些实施方式测量到的径向速度的几何形状的示意图。
[图3]
图3示出了一些实施方式所使用的在复杂地形上方的风的遥感的示意图。
[图4]
图4示出了一些实施方式估计风流的速度场所使用的风流的示例性参数的示意图。
[图5A]
图5A示出了由一些实施方式使用的计算流体动力学(CFD)框架的框图。
[图5B]
图5B示出了根据一个实施方式的用于确定无偏差速度场的方法的框图。
[图6A]
图6A示出了一些实施方式获得垂直速度的水平梯度而使用的基于CFD模拟的框架的框图。
[图6B]
图6B示出了由一些实施方式所确定的地形的网格的示例。
[图7]
图7示出了根据一些实施方式的用于选择操作参数的方法的框图。
[图8]
图8示出了根据一些实施方式的用于确定操作参数的当前值的方法的流程图。
[图9]
图9示出了根据一些实施方式的对成本函数中的不同项指配不同权重的示意图。
[图10]
图10示出了根据一些实施方式的用于迭代地确定CFD的操作参数的方法的框图。
[图11]
图11示出了根据一些实施方式的与风流感测有关的单个平面上的各种数据点的示例。
[图12]
图12示出了根据一个实施方式的用于确定风流的水平速度的方法的流程图。
[图13]
图13示出了根据一个实施方式的用于CFD模拟的方法的框图。
具体实施方式
图1示出了根据一些实施方式的风流感测系统100的不同组件的硬件图。风流感测系统100包括被配置为执行所存储的指令的处理器120，以及存储处理器可执行的指令的存储器140。处理器120可以是单核处理器、多核处理器、计算集群或任何数量的其它配置。存储器140可以包括随机存取存储器(RAM)、只读存储器(ROM)、闪存或任何其它合适的存储器系统。处理器120通过总线106连接到一个或更多个输入输出接口和/或装置。
存储器140中存储的这些指令实现了一种用于从在每个海拔高度处的径向速度的测量结果的集合来确定在一组不同海拔高度处的风流的速度场的方法。为此，风流感测系统100还可以包括储存装置130，该储存装置130适于存储不同模块，该不同模块存储处理器120的可执行指令。该储存装置存储：CFD模拟模块131，其被配置为通过用操作参数的当前值来模拟风流的计算流体动力学(CFD)来求在每个海拔高度处的速度场；CFD操作参数模块132，其被配置为确定降低成本函数的操作参数的值；以及水平导数模块133，其被配置为从各个海拔高度处的速度场确定在各个海拔高度处的垂直速度的水平导数；以及速度场模块，其被配置为使用垂直速度的水平导数和径向速度的测量结果来确定包括水平速度的速度场。储存装置130可以使用硬盘驱动器、光盘驱动器、拇指驱动器、驱动器阵列或其任何组合来实现。
风流感测系统100包括接受针对每个海拔高度的视线点(line-of-site points)处的径向速度的测量结果195的集合的输入接口。例如，在一些实现中，输入接口包括风流感测系统100内的将处理器120连接到键盘111和定点装置112的人机接口110，其中，定点装置112可以包括鼠标、轨迹球、触摸板、操纵杆、定点杆、手写笔或触摸屏等。附加地或另选地，输入接口可以包括适于通过总线106将风流感测系统100连接到网络190的网络接口控制器150。通过网络190，测量结果195可以被下载并存储在储存系统130内，用于存储和/或进一步处理。在一些实现中，网络190通过无线或有线连接将风流感测系统100与被配置为测量风流的径向速度的遥感仪器连接。
风流感测系统100包括呈现风流的速度场的输出接口。例如，风流感测系统100可以通过总线106链接到显示接口160，该显示接口160适于将风流感测系统100连接到显示装置165，其中显示装置165可以包括计算机监视器、相机、电视、投影仪或移动装置等。
例如，风流感测系统100可以连接至系统接口170，该系统接口170适合于将风流感测系统连接至基于重构的速度场而控制的不同系统175。附加地或另选地，风流感测系统100可以通过总线106连接到应用接口180，该应用接口180适于将风流感测系统100连接到可以基于速度场重构的结果而操作的应用装置185。
图2A示出了根据一些实施方式的被配置为测量风流的径向速度的示例性遥感仪器的示意图。图2A的示例示出了激光雷达250的示意图，该激光雷达250测量在不同海拔高度处的风流的径向速度195。不同的实施方式使用不同的遥感仪器进行远程测量。这些仪器的示例包括雷达、激光雷达和声雷达。为了清楚起见，本公开使用诸如图1中所示的激光雷达250之类的激光雷达作为示例性遥感仪器。
对象相对于给定点的径向速度是对象与该点之间的距离的变化率。也就是说，径向速度是对象速度的指向将对象和该点连接的半径方向的分量。在大气测量的情况下，该点是诸如雷达、激光雷达和声雷达之类的遥感仪器在地球上的位置，因此径向速度表示对象离开或接近接收仪器的速度。此测得的径向速度也称为视线(LOS：line-of-sight)速度。
遥感仪器通过描述气流的速度场来确定感兴趣体积中的诸如空气之类的流体的流动。例如，激光雷达250包括：对返回信号进行谱分析的激光器251或声学发送器和接收器；用于执行进一步的计算的计算机257；以及用于将发送器和/或接收器对准空间中的距所述发送器和接收器相当大距离的目标的导航仪。接收器沿着测量轴检测由于在遥感系统和所述目标之间存在污染物而散射的返回信号253。波沿着由可能瞄准方向形成的锥面200传输。由于特定的空气污染物，通过Doppler效应从频移中推论出在所述目标处感兴趣体积255处的颗粒的径向速度。
图2B示出了由一些实施方式沿着锥面200以及沿圆锥的中心线208在特定海拔高度处测量到的径向速度的几何形状的示意图。激光雷达测量结果提供了风的径向(视线)速度分量，由于所谓的“独眼巨人(cyclops)”困境，使得其难以精确地确定风的大小和方向。沿一个波束的径向速度201示出了使用位于笛卡尔坐标系204、205和206中的位置203处的激光雷达的速度矢量220的投影。
这里，θ204是从北方202顺时针测量到的水平风向，ψ207是波束的仰角，(u,v,w)是空间中每个点的风速V的x分量206、y分量204和z分量205。
在每个海拔高度处的水平速度vh定义为
Figure BDA0002682572480000101
在每个海拔高度上，将径向速度(也称为LOS速度)定义为
Figure BDA0002682572480000102
图3示出了由一些实施方式使用的在复杂地形300上的风的遥感的示意图。布置在点308(例如，在山顶)处的激光雷达250执行圆锥上的一系列视线测量，包括沿锥面的测量309、310、312、313和沿中心线311的测量。针对被示为不同平面301的不同海拔高度取得测量结果。以这种方式，对于每个海拔高度301，圆锥上的测量结果为圆上的测量结果，该圆上的测量结果包括在圆的圆周上的不同视线点处测量到的径向速度在不同角度方向上的多个测量结果以及在圆心处测量到的径向速度在垂直方向上的一个测量结果。
一个实施方式旨在确定每个海拔高度的风流的水平速度vh。给定这些测量结果，可以使用几何关系并假设风速在每个平面上是均匀的，从径向速度vR的测量结果确定水平速度vh的估计值。在此，VL＝(uL,vl,wL)是基于均匀假设的风流的估计速度。
例如，以下公式产生对应于指向北、东、南、西和顶部的波束的径向速度V1、V2、V3、V4、V5的估计速度：
Figure BDA0002682572480000111
wL＝V5
一些实施方式基于以下认识：对于诸如地形300之类的复杂地形，均匀速度假设导致水平速度的激光雷达估计中的偏差。主要误差是由于垂直速度w在垂直方向(例如，沿着山丘)上的变化而引起的。为此，一些实施方式基于以下认识：能够使用垂直速度的水平导数来校正在感测通过复杂地形的风流中的均匀速度假设。
图4示出了一些实施方式估计风流的速度场所使用的示例性参数的示意图。一些实施方式基于以下认识：在感测通过复杂地形的风流中的均匀速度假设是不正确的，但是能够使用垂直速度的水平导数来校正。图4示出了当装置402放置在山顶附近(例如，位置308)时风流的二维图示。垂直速度的水平导数能够示出在给定海拔高度处垂直速度401在方向和/或大小上的变化。在此示例中，垂直速度的水平导数表示垂直速度沿山丘的一个山坡直到点410处的山顶的增加、以及垂直速度从点410起在山丘的另一山坡上的减小。甚至示出了垂直速度线性变化的一阶导数能够用来提高风流感测的准确度，这是因为均匀速度假设导致感测速度场的误差中的主要项。
对于装置上方海拔高度z处的任何点，可以将误差或偏差写入一阶，如下所示：
Figure BDA0002682572480000112
因此，由均匀假设引起的偏差与i)装置上方的海拔高度z，ii)垂直速度的水平梯度dw/dx和dw/dy成正比。这种误差不是仰角ψ的函数，减小该角度将不会减小水平速度的偏差。
一些实施方式基于以下认识：仅基于可用的径向速度测量结果，可能无法获得dw/dx和dw/dy的估计。由于扫描波束的对称性，所得到的方程式组是不确定的。因此，不能从这样的典范扫描中获得垂直速度的任何水平梯度信息。
流的非压缩性是指在流体颗粒中材料密度恒定的流，流体颗粒是随流速而移动的无穷小体积。这种物理原理基于质量守恒。
一些实施方式基于以下认识：由均匀速度假设引起的主要误差(leading ordererrors)是非压缩的。换句话说，这能够表明由海拔高度和垂直速度的水平梯度的乘积组成的偏差项保存了质量。这意味着对于在复杂地形上方的流，在感兴趣域内的流的体积上强加非压缩条件将不会校正由均匀流假设导致的主要误差。
计算流体动力学(CFD)是使用数值分析和数据结构来解决和分析涉及流体流的问题的流体力学的一个分支。使用计算机来执行模拟液体和气体与由边界条件限定的表面之间的相互作用所需的计算。一些实施方式基于一般的理解，即CFD能够用于从由激光雷达感测到的圆锥上的风的测量结果来估计风的速度场。然而，关于复杂地形上方的风流的操作参数(诸如边界条件)通常是未知的，并且这些操作参数的近似值能够不期望地降低风流感测的准确度。
一些实施方式基于以下认识：尽管CFD近似值对于确定速度场可能不够准确，但是CFD近似值对于在给定海拔高度的垂直速度重构的水平导数的平均值能够足够精确，这进而能够用于校正由于均匀速度假设导致的偏差。为此，一些实施方式使用CFD近似值来确定垂直速度的水平导数，并且结合在期望海拔高度上的风流的径向速度测量结果使用垂直速度的水平导数来确定期望海拔高度的速度场。以这种方式，能够实现使用径向速度测量结果感测速度场的目标准确度。
图5A示出了一些实施方式出于获得在每个感兴趣海拔高度处垂直速度的水平梯度的精确测量值的目的而求解盛行的风流所使用的计算流体动力学(CFD)框架的框图。使用视线测量结果501，通过CFD模拟获得502速度场的第一近似值。例如，CFD模拟可以由图1所示的模块131执行。在许多情况下，CFD模拟需要诸如边界条件之类的操作参数的知识。这些操作参数通常是未知的。为此，一些实施方式确定减小估计的径向速度和测量的径向速度之间的差的操作参数。例如，这种估计可以由图1所示的模块132执行。
尽管CFD提供的第一近似值中的速度场对于所需目的而言是不准确的，但能够以所需的准确度提取垂直速度的水平梯度的估计值503。这种提取可以由模块133执行。CFD模拟得到在网格的离散点处的速度场。使用该速度场，使用有限差分法计算每个离散点的x导数和y导数。然后，通过在各个平面上对x和y方向上的导数取平均值，提取在每个平面处垂直速度的x和y水平导数
Figure BDA0002682572480000131
的单一值。
然后，使用垂直速度的这种水平梯度以及视线速度和风速之间的几何关系，来使用式(3a)和式(3b)基于均匀假设来校正504偏差的水平速度分量uL和vL。这种估计可以由模块134执行。
图5B示出了根据一个实施方式的用于确定无偏差速度场的方法的框图。为了确定速度场的第二近似值，实施方式针对每个海拔高度确定510速度场的均匀速度假设下的有偏差速度场，并且针对每个海拔高度使用相应海拔高度的垂直速度的水平导数
Figure BDA0002682572480000132
去除520有偏差速度场的均匀速度假设的偏差。
例如，式3a/式3b用于通过减去偏差项
Figure BDA0002682572480000133
和
Figure BDA0002682572480000134
从偏差速度场uL,vL获得无偏差速度场(u,v)。
图6A示出了一些实施方式获得垂直速度的水平梯度650而使用的基于CFD模拟的框架的框图。实施方式执行预处理步骤以定义CFD模拟的几何形状和物理边界。例如，一些实现使用计算机辅助设计(CAD)来定义模拟范围。流体所占的体积被划分为离散的单元(网格)。装置中存在提取地形的地理位置610的GPS。然后将该位置与存储在装置存储器中的可用数据集进行比较，以生成地形数据。可以使用诸如Google或NASA数据库之类的各种资源来收集地形数据。在读取地形信息之后，选择最佳半径以构造网格620。
图6B示出了由一些实施方式确定的网格620的示例。在各种实现中，网格620可以是均匀的或不均匀的、结构化的或非结构化的、由六面体、四面体、棱柱、金字塔或多面体元素的组合组成。域的最佳半径被选择为使得基于风向将重要地形结构捕获到网格中。基于所选地形上的选定半径，来生成网格。网格的分辨率由用户设置。在一些实现中，网格的默认分辨率值为25m。
在预处理期间，还指定操作参数630的值。在一些实施方式中，操作参数规定了在流体域的所有边界面处的流体行为和性质。场(速度、压力)的边界条件规定了函数本身的值，或者函数的法线导数的值，或者为法线导数和变量本身赋予值的曲线或曲面的形式、或者在给定区域上函数的值与函数的导数之间的关系。在由地形所定义固体表面的涉及流体速度的边界条件可以设置为零。流入速度基于在平坦地形上方风的方向以及相对于高度而具有对数轮廓的速度决定。
一些实施方式通过用操作参数的当前值求解640定义了风流的Navier-Stokes(纳维叶－斯托克斯)方程式的变型之一，来执行CFD的模拟。例如，CFD连同质量和能量守恒一起求解Navier-Stokes方程式。该组方程式被证明代表任何诸如空气之类的Newtonian(牛顿)流体的机械性能，并且实现模拟大气流动。Navier-Stokes方程式的离散化是以能够应用于计算流体动力学的方式对方程式进行重新公式化。数值方法可以是有限体积法、有限元素法或有限差分法、以及谱方法或谱元素方法。
控制方程式Navier-Stokes如下：
Figure BDA0002682572480000141
Figure BDA0002682572480000142
是散度算子。
Figure BDA0002682572480000143
是梯度算子，并且
Figure BDA0002682572480000144
是拉普拉斯算子。
一些实施方式将式4a和4b表示为N(p,V)＝0，流入速度和方向由Vin,Θin表示；p：空气压力[pa]或[atm]；ρ：空气密度[kg/m3]；v：运动粘滞性[m2/s]。
在CFD模拟之后，实施方式提取650垂直速度的水平梯度。
图7示出了根据一些实施方式的用于选择操作参数的方法的框图。例如，一些实施方式基于垂直速度的水平导数对这些操作参数的值的变化的灵敏度720来选择725操作参数710。在一个实施方式中，在CFD模拟期间所近似的操作参数的基于目的的集合中选择灵敏度高于阈值715的操作参数。以这种方式，一些实施方式使CFD的未知操作参数适于CFD近似的目的。操作参数的这种适配减轻了计算负担，而没有降低感兴趣的量的CFD近似的准确度。例如，一些实施方式选择操作参数，诸如地形粗糙度、流入平均速度、流入湍流强度和大气稳定性条件。
在一些实施方式中，操作参数包括流入边界条件(速度、方向)、表面粗糙度、大气稳定性。在一个实施方式中，将操作参数选择为流入边界条件(速度、方向)、表面粗糙度、流入湍流动能和耗散。这些参数的值不能直接从LOS激光雷达测量结果中获得。
Figure BDA0002682572480000151
Cμ是k-∈湍流模型中的常数
κ是von Karman(冯·卡门)常数
V*是摩擦速度[m/s]
Vref是在参考位置处选择的参考速度。参考位置可以是任意的。[m/s]
zref是参考海拔高度。[m]
z0是表面粗糙度。[m]
湍流动能k是湍流波动的每单位质量的动能。湍流耗散∈是湍流动能转换成内部热能的速率。
一些实施方式基于以下认识：在许多情况下，用于模拟CFD的操作参数是未知的。例如，对于上述情况，在流入时的式5a-5d中，Vref,zref,z0是未知操作参数，并且遥感测量结果不直接提供此类值。
图8示出了根据一些实施方式的用于确定操作参数的当前值的方法的流程图。具体地，一些实施方式确定820使在一组视线点处的径向速度的测量结果810与由CFD使用操作参数的当前值执行的在同一组视线点处的径向速度的估计结果830之间的误差最小化的操作参数。
一些实施方式基于以下认识：当CFD用于提取垂直速度的水平导数时，需要使特定成本函数840最小化以获得操作参数的估计值。具体地，一些实施方式基于以下认识：垂直速度的水平导数依据海拔高度对速度场具有不同的影响。为此，成本函数840包括误差的加权组合。每个误差对应于海拔高度之一，并且包括在相应海拔高度处在视线点处测量到的速度与由CFD针对相应海拔高度使用操作参数的当前值所模拟的在视线点处的模拟速度之间的差。另外，至少一些误差的权重是不同的。例如，误差包括对应于第一海拔高度的第一误差和对应于第二海拔高度的第二误差，其中，在误差的加权组合中的第一误差的权重与误差的加权组合中的第二误差的权重不同。
图9示出了根据一些实施方式的为成本函数中的不同项指配不同权重的过程。在一些示例中，成本函数840返回表示CFD模拟920与在各个海拔高度沿不同波束的视线的激光雷达数据910的匹配程度的数字。为此，当CFD的目的是确定垂直速度的水平导数时，成本函数需要不同地考虑不同的海拔高度，例如利用不同的权重930。例如，在一些实施方式中，成本函数包括代表不同海拔高度的CFD的准确度的误差的加权组合。
在一实施方式中，成本函数为
Figure BDA0002682572480000161
i是各个测量点，vR,i是在点i的位置处的视线速度，vR,CFD是在点i的位置处从CFD模拟计算出的径向速度，wi是加权因子。每个项中的误差与测量结果和CFD之间的径向速度之差成比例。由于在较高海拔高度处估计垂直速度梯度对由均匀假设导致的偏差有更大的贡献(参见式3)，因为为了向在较高海拔高度处估计垂直速度梯度赋予更大的权重，一些实现将加权因子wi设置为与海拔高度(即，发送器位置之上的高度)成比例。
例如，vR,CFD是从风流的CFD模拟获得的径向速度的集合，以产生降低式(6)中给出的误差的加权组合的成本函数的速度场的第一近似值。
该vR的集合表示由风流的遥感装置给出的径向(或视线)速度的测量结果。这样的值具有非常小的误差，并且用作在波束方向上的风的真实值。式(6)中用i表示的每个项对应于由于海拔高度之一引起的误差，并且包括在相应海拔高度处的视线点处测量到的速度vR与由CFD针对相应海拔高度所模拟的在视线点处的模拟速度之间的差。误差的加权组合中每个误差的权重是相应海拔高度值的递增函数。
一些实施方式基于以下认识：可以使用基于直接伴随循环(DAL)的CFD框架来估计操作参数的未知值。该框架通过使估计在正向CFD模拟与可用激光雷达测量结果之间的视线数据及其梯度中的误差的成本函数最小化，然后以迭代方式求解灵敏度(或伴随CFD)方程式，来同时校正用于共同目的的各种未知参数。用于共同目的的参数的灵敏度指示了基于DAL的CFD框架的收敛方向。同时校正减少了更新多个操作参数的计算时间。
一些实施方式用(ξ1,ξ2,…ξn)表示需要估计的操作参数的集合。然后，成本函数J对任何操作参数ξi的灵敏度可以表示为
Figure BDA0002682572480000171
图10示出了一个实施方式以迭代方式确定操作参数和CFD模拟结果所使用的DAL的实现的框图。该实施方式通过评估CFD模拟来估计操作参数的最可能值。DAL是一种优化方法，其以迭代1080方式求解CFD方程式1020和伴随(或灵敏度)方程式1030，以获得在操作参数的当前估计值下针对未知操作参数的成本函数的灵敏度1040。用操作参数的猜测来初始化1010DAL方法。例如，使用Bernoulli(伯努利)方程式来估计流入速度，并且使用均匀假设来估计角度。在每次迭代之后，使用共轭梯度下降来更新操作参数的当前值的估计1050，在成本函数的灵敏度的最大下降方向上来更新该共轭梯度下降。为此，多次执行(即，每个迭代执行一次)CFD模拟，并且如果自先前迭代的估计变化低于阈值，则认为DAL方法收敛1060。
DAL方法是通过公式化Lagrangian(拉格朗日)获得的
Figure BDA0002682572480000172
由于方程式(Navier-Stokes方程式)中的N(p,V)＝0，因此当p和V的值准确时，方程式L和J相等。考虑到ξi的变化，L的变化可以表示为
Figure BDA0002682572480000173
为了确定该项
Figure BDA0002682572480000174
选择伴随变量以满足
Figure BDA0002682572480000175
因此，DAL方法涉及新变量(Va,pa)，它们分别表示使
Figure BDA0002682572480000176
可计算的伴随速度和压力。
在一个实施方式中，未知参数被选择为Vin,Θin，即，流入速度和流入角度。因此，将寻找使J最小化的Vin,Θin的问题变换为寻找使增强目标函数L最小化的Vin,Θin的问题。例如，为确定δJ/δVin和δJ/δΘin，能够通过设置ξi＝Vin或ξi＝Θin，来使用DAL方法。
pa：伴随压力
Va：伴随速度
步骤1030中的伴随方程式为
Figure BDA0002682572480000181
算子
Figure BDA0002682572480000182
对应于速度矢量梯度的转置。
伴随变量能够用于确定成本函数对任何操作参数的灵敏度1040
Figure BDA0002682572480000183
例如，针对流入速度Vin，成本函数的灵敏度可以写为式11
Figure BDA0002682572480000184
Ain：计算域Ω的流入面积[m2]
n：Ain的单位法线向量[m2]
通过使用梯度下降算法，操作参数的估计值能够更新1050为：
Figure BDA0002682572480000185
λ是代表步长的正常数，能够使用多种标准算法来选择该正常数。
对于每次迭代，需要进行两次模拟以使用式7确定每个操作参数的灵敏度。当操作变量的数量大时，这种计算成本相当大的。使用DAL方法，无论未知参数的数量如何，每次迭代仅求解式(4)和(10)一次，并且因此降低了计算成本，并使优化问题可求解。这是伴随方法相对于通过直接测量成本函数的扰动来确定成本函数的灵敏度的方法的优点。
在DAL收敛以产生外部操作参数的当前值1070之后，一些实施方式提取感兴趣的量(即，垂直速度梯度)，以使用在测量圆锥上的激光雷达视线(LOS)来校正在复杂地形上方的风速重建中的偏差误差。
图11示出了根据一些实施方式的与风流感测有关的单个平面上的各种数据点的示例。在该示例中，圆1120上的点是测量径向速度1110的点。在一些实现中，每个海拔高度的速度场包括在圆锥(即，圆1120)内外的风的速度值。
附加地或另选地，在一些实施方式中，在每个海拔高度处的垂直速度的水平导数定义了在定义针对相应海拔高度的激光雷达的测量结果的、圆锥的圆心处的垂直速度的梯度。例如，一些实施方式将每个海拔高度的速度和/或梯度求平均，以产生圆锥的中心1150和圆1120。在那些实施方式中，通过几何关系并且使用速度的水平梯度去除偏差而获得的速度场的第二近似值提供在每个平面(或每个海拔高度)上的速度场的单一值。以这种方式，在1130和1140处的无偏差速度值被认为等于那个单一值。
为此，在一个实施方式中，速度场的第二近似值包括针对每个海拔高度的速度场的单一值。另外，该实施方式通过强迫与在每个海拔高度上的径向速度的测量结果一致的风流的非压缩性和正则化，将单一值转换成在每个海拔高度上的速度场的非恒定值的密集网格。在这样的变换之后，在圆锥内外的点(诸如，点1120、1150、1130和1140)处的水平速度能够具有不同的值。
图12示出了根据一个实施方式的用于确定风流的水平速度的方法的流程图。该实施方式从单个平面1210上的所有点处的相同无偏差值开始，强制空气1220的非压缩性来校正由LOS测量结果的稀疏性引起的误差和由于均匀假设而产生的二阶误差1220，以产生速度值的密集网格1230。在各个实施方式中，点的网格密度是用户指定的值。值得注意的是，在该实施方式中，与强制非压缩性来校正由于均匀速度假设而导致的主要误差相对比，空气的非压缩性被用于校正二阶误差。
在一些实现中，使用基于直接伴随循环的算法来确定速度场的非恒定值的密集网格。该算法开始于在每个平面中在网格上的所有离散点处插入无偏差速度值1210。DAL问题1220被公式化为强制体积的非压缩性，同时使具有以下两项的成本函数最小化：一项用于测量离散点处的最终速度场与初始速度场之间的差异，另一项是用于增加速度场的平滑度的正则项。所得到的伴随变量λ的伴随方程式为：
Figure BDA0002682572480000191
其中，Uk是DAL循环第k次迭代的速度场。在每次迭代结束时，执行以下更新：
Figure BDA0002682572480000192
当达到收敛时，算法结束。
在Navier-Stokes方程式的求解中，计算成本取决于流体的速度和粘度。对于大气流，由于风速高而空气的粘度小，因此计算成本非常大。这导致流中的不稳定惯性力明显大于稳定粘滞力的、所谓的高Reynold(雷诺)数流。为了完全解决动力学问题并避免数值不稳定性，湍流的所有空间尺度必须在计算网格中与包含大部分动能的运动相关联地、从最小的耗散尺度(Kolmogorov(柯尔莫哥洛夫)尺度)直到与域尺寸成比例的积分尺度进行分解。
大涡模拟(LES)是用于求解流体力学的控制方程式的流行类型的CFD技术。Kolmogorov的自相似理论的含义是大涡流取决于几何形状，而尺度越小则越通用。此特征使得在计算中显式地求解大涡流而通过使用子网格规模模型(SGS模型)隐式地考虑小涡流。使用LES方法的CFD模拟能够以高保真度模拟流场，但是计算成本也非常高。
一些实施方式基于以下认识：能够修改低保真模型以学习针对结果中期望准确度所需的内部模型参数，而不是对每个新的测量数据集(例如，对于每个新的风向和/或新的地形)使用高保真CFD解。
图13示出了根据一个实施方式的用于CFD模拟的方法的框图。低保真CFD模拟通过取决于一些内部参数的模型来近似流中的小尺度项。为此，实施方式使用包括垂直速度的水平导数的特征向量来应用场反转和机器学习(FIML)方法1302，以学习低保真模型1303的内部参数对高保真LES模拟1301中的流特征的依赖性。
我们使用的低保真CFD模型是Reynolds(雷诺)平均的Navier–Stokes方程式(或RANS方程式)。这些是流体流的运动的时间平均方程式。这些模型具有内部参数，以近似由于低保真度而无法分解的项。这些内部参数的正确值是特定于问题的，因此为了使RANS与LES几乎一样准确，采用了FIML框架。将RANS与FIML结合使用的显著优势是，与高保真LES模拟相比，高Reynolds数的CFD模拟的成本降低了几个数量级，同时保持了所需的准确度。一旦低保真模型的内部模型参数是离线固定的，那么如果操作参数是已知的，则能够执行基于RANS的CFD模拟。
以这种方式，在一些实施方式中，通过求解雷诺平均的Navier-Stokes(RANS)方程式来执行风流的CFD的模拟，同时使用利用包括每个海拔高度的垂直速度的水平导数的特征向量的场反转和机器学习(FIML)来确定RANS方程式的内部操作参数。
为了以流入速度和风向为焦点来寻找最可行操作参数，使用用于计算水平速度的均匀假设，激光雷达与杯型风速仪的相对误差约为8％，而使用根据本发明的方法并使用CFD和DAL，激光雷达与杯型风速仪的相对误差约为1％。此外，一些实施方式强制非压缩性假设，以重构锥形区域内外的密度场。
能够以多种方式中的任何方法来实现本发明的上述实施方式。例如，可以使用硬件、软件或其组合来实现实施方式。当以软件实现时，软件代码能够在任何合适的处理器或处理器集上执行，无论其是设置在单台计算机中还是分布在多台计算机当中。这种处理器可以实现为集成电路，在集成电路组件中具有一个或更多个处理器。但是，可以使用任何适当格式的电路来实现处理器。
另外，本公开的实施方式可以被体现为一种方法，已经提供了该方法的示例。作为该方法的一部分而执行的动作可以以任何合适的方式排序。因此，可以构造实施方式，在该实施方式中以与所例示的次序不同的次序来执行动作，这可以包括同时执行一些动作，即使这些动作在示例性实施方式中被示为顺序动作。
在权利要求中使用诸如“第一”、“第二”之类的序数术语来修饰权利要求要素本身并不意味着一个权利要求要素相对于另一权利要求要素的优先权、优先级或次序，或者执行方法动作的时间次序，而仅用作标签以将具有一定名称的一个权利要求要素与具有相同名称(除了使用序数词之外)的另一要素区分开，以区分权利要求要素。

Claims (20)

1.一种风流感测系统，该风流感测系统用于从在每个海拔高度处的径向速度的测量结果的集合来确定在一组不同海拔高度处的风流的速度场，该风流感测系统包括：
输入接口，该输入接口接受针对每个所述海拔高度的在视线点处的径向速度的所述测量结果的集合；
处理器，该处理器被配置为：
通过使用使误差的加权组合的成本函数最小化的操作参数来模拟所述风流的计算流体动力学CFD，确定在每个所述海拔高度处的所述速度场的第一近似值，其中，每个误差对应于所述海拔高度之一，并且包括在相应海拔高度处在测量点处的测量速度与通过所述CFD针对所述相应海拔高度利用所述操作参数的当前值所模拟的所述测量点处的模拟视线速度之间的差，其中，所述误差包括对应于第一海拔高度的第一误差和对应于第二海拔高度的第二误差，其中，在所述误差的加权组合中的所述第一误差的权重不同于在所述误差的加权组合中的所述第二误差的权重；
从每个所述海拔高度处的所述速度场的所述第一近似值，确定在每个所述海拔高度处的垂直速度的水平导数；以及
使用每个所述海拔高度的速度场、所述径向速度的测量结果在三维轴上的投影、以及相应速度场的所述垂直速度的水平导数之间的几何关系，确定每个所述海拔高度处的所述速度场的第二近似值；以及
输出接口，该输出接口呈现所述风流的速度场的所述第二近似值。
2.根据权利要求1所述的感测系统，其中，在所述误差的加权组合中的每个误差的权重是所述相应海拔高度的值的递增函数。
3.根据权利要求1所述的感测系统，其中，所述处理器通过执行一组迭代直到满足终止条件来确定所述速度场的所述第一近似值，其中，对于每次迭代，所述处理器被配置为：
利用所述操作参数的当前值执行所述风流的CFD的模拟，以减小所述成本函数的值；
确定使用利用所述操作参数的当前值进行的所述CFD所确定的所述成本函数的值对所述操作参数的值的变化的灵敏度；以及
根据所确定的灵敏度更新所述操作参数。
4.根据权利要求3所述的感测系统，
其中，所述处理器通过利用所述操作参数的当前值来求解定义所述风流的Navier-Stokes方程式，来执行所述CFD的模拟；
其中，所述处理器基于使用直接伴随循环的所述CFD的伴随模型来确定所述成本函数的值的灵敏度；
其中，所述处理器在所述灵敏度的最大下降方向上更新所述操作参数的当前值。
5.根据权利要求3所述的感测系统，其中，所述径向速度的所述测量结果是由基于地面的激光雷达在圆锥上测量的，其中，所述操作参数包括地形粗糙度、流入速度、流入湍流强度和大气稳定性条件，其中，所述CFD使用在所述激光雷达的基于地面的位置处的地形的形状。
6.根据权利要求3所述的感测系统，其中，所述操作参数是基于所述操作参数对所述垂直速度的水平导数的值的变化的灵敏度而选择的。
7.根据权利要求1所述的感测系统，其中，为了确定所述速度场的所述第二近似值，所述处理器被配置为：
针对每个所述海拔高度，确定在所述速度场的均匀速度假设下的有偏差速度场；以及
使用每个所述海拔高度的所述垂直速度的水平导数，去除相应海拔高度的所述有偏差速度场的所述均匀速度假设的偏差。
8.根据权利要求1所述的感测系统，其中，通过求解雷诺平均的Navier-Stokes RANS方程式来执行所述风流的CFD的模拟。
9.根据权利要求8所述的感测系统，其中，使用利用包括每个所述海拔高度的所述垂直速度的水平导数的特征向量的场反转和机器学习FIML，来确定所述RANS方程式的内部操作参数。
10.根据权利要求1所述的感测系统，其中，所述径向速度的所述测量结果是由基于地面的激光雷达在圆锥上测量的，使得对于每个海拔高度，所述圆锥上的测量结果为圆上的测量结果，该圆上的测量结果包括在该圆的圆周上的不同视线点处测量到的所述径向速度在不同角度方向上的多个测量结果、以及在所述圆的中心处测量到的所述径向速度在垂直方向上的一个测量结果。
11.根据权利要求10所述的感测系统，其中，在每个所述海拔高度处的所述垂直速度的水平导数定义了在定义所述激光雷达针对相应海拔高度的测量结果的、所述圆锥的所述圆的中心处的所述垂直速度的梯度。
12.根据权利要求11所述的感测系统，其中，针对每个所述海拔高度的所述速度场包括所述圆锥内外的风的速度值。
13.根据权利要求1所述的感测系统，其中，所述速度场的所述第二近似值包括针对每个所述海拔高度的所述速度场的单一值，并且其中，所述处理器通过强制与在每个所述海拔高度处的径向速度的测量结果一致的所述风流的非压缩性和正则化，将所述单一值变换为在每个所述海拔高度处所述速度场的非恒定值的密集网格。
14.根据权利要求13所述的感测系统，其中，所述速度场的非恒定值的所述密集网格是使用直接伴随循环来确定的。
15.根据权利要求13所述的感测系统，该感测系统还包括：
控制器，该控制器基于在每个所述海拔高度处所述速度场的非恒定值的所述密集网格来控制风敏系统。
16.根据权利要求15所述的感测系统，其中，所述风敏系统是风力涡轮机。
17.根据权利要求1所述的感测系统，该感测系统还包括：
控制器，所述控制器基于所述风流的速度场的所述第二近似值来控制风敏系统。
18.一种风流感测方法，该风流感测方法用于从在每个海拔高度处的径向速度的测量结果的集合来确定在一组不同海拔高度处的风流的速度场，其中，该方法使用与实现该方法的存储指令联接的处理器，其中，所述指令在由所述处理器执行时实施该方法的步骤，所述步骤包括：
接受针对每个海拔高度的在视线点处的径向速度的所述测量结果的集合；
通过使用减小误差的加权组合的成本函数的操作参数来模拟所述风流的计算流体动力学CFD，确定在每个所述海拔高度处的所述速度场的第一近似值，每个误差对应于所述海拔高度之一，并且包括在相应海拔高度处在所述视线点处的测量速度与通过所述CFD针对相应海拔高度利用所述操作参数的当前值所模拟的在所述视线点处的模拟速度之间的差，其中，所述误差包括对应于第一海拔高度的第一误差和对应于第二海拔高度的第二误差，其中，在所述误差的加权组合中的所述第一误差的权重不同于在所述误差的加权组合中的所述第二误差的权重；
从所述速度场的所述第一近似值，确定在每个所述海拔高度处的垂直速度的水平导数；以及
使用每个所述海拔高度的速度场、所述径向速度的测量结果在三维轴上的投影、以及相应速度场的所述垂直速度的水平导数之间的几何关系，确定所述速度场的第二近似值，以及
呈现所述风流的速度场的所述第二近似值。
19.根据权利要求18所述的感测方法，其中，通过执行一组迭代直到满足终止条件来确定所述速度场的所述第一近似值，其中，每次迭代包括以下步骤：
用操作参数的当前值执行所述风流的CFD的模拟，以减小所述成本函数的值；
确定利用所述操作参数的当前值针对所述CFD确定的所述成本函数的值对所述操作参数的值的变化的灵敏度；以及
根据所确定的灵敏度来更新所述操作参数。
20.一种非暂时性计算机可读存储介质，该非暂时性计算机可读存储介质上实现有用于执行以下方法的由处理器能够执行的程序，所述方法用于从在每个海拔高度处的径向速度的测量结果的集合来确定在一组不同海拔高度处的风流的速度场，所述方法包括以下操作：
接受针对每个海拔高度的在视线点处的径向速度的测量结果的集合；
通过使用减小误差的加权组合的成本函数的操作参数来模拟所述风流的计算流体动力学CFD，确定在每个所述海拔高度处的所述速度场的第一近似值，每个误差对应于所述海拔高度之一，并且包括在相应海拔高度处在视线点处的测量速度与通过所述CFD针对所述相应海拔高度利用所述操作参数的当前值所模拟的在所述视线点处的模拟速度之间的差，其中，所述误差包括对应于第一海拔高度的第一误差和对应于第二海拔高度的第二误差，其中，在所述误差的加权组合中的所述第一误差的权重不同于在所述误差的加权组合中的所述第二误差的权重；
从所述速度场的所述第一近似值，确定在每个所述海拔高度处的垂直速度的水平导数；以及
使用每个所述海拔高度的速度场、所述径向速度的测量结果在三维轴上的投影、以及相应速度场的所述垂直速度的水平导数之间的几何关系，确定所述速度场的第二近似值，以及
呈现所述风流的速度场的所述第二近似值。
CN201880091302.0A 2018-03-20 2018-10-10 风流感测系统和用于确定风流的速度场的方法 Active CN111868533B (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US15/926,864 2018-03-20
US15/926,864 US10795054B2 (en) 2018-03-20 2018-03-20 System and method for sensing wind flow passing over complex terrain
PCT/JP2018/038566 WO2019181032A1 (en) 2018-03-20 2018-10-10 Wind flow sensing system and method for determining velocity fields of wind flow

Publications (2)

Publication Number Publication Date
CN111868533A true CN111868533A (zh) 2020-10-30
CN111868533B CN111868533B (zh) 2022-05-10

Family

ID=64277760

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201880091302.0A Active CN111868533B (zh) 2018-03-20 2018-10-10 风流感测系统和用于确定风流的速度场的方法

Country Status (5)

Country Link
US (1) US10795054B2 (zh)
EP (1) EP3769091B1 (zh)
JP (1) JP6896180B2 (zh)
CN (1) CN111868533B (zh)
WO (1) WO2019181032A1 (zh)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110598362B (zh) * 2019-09-30 2021-03-16 西安交通大学 充分发展湍流条件下的大涡模拟计算速度边界条件生成方法
FR3103562B1 (fr) * 2019-11-25 2022-06-24 Leosphere Procédé de traitement de données télémétriques pour l’estimation d’une vitesse du vent.
US11442076B2 (en) 2020-04-02 2022-09-13 Mitsubishi Electric Research Laboratories, Inc. System and method for wind flow turbulence measurement by LiDAR in a complex terrain
US11294063B2 (en) 2020-04-02 2022-04-05 Mitsubishi Electric Research Laboratories, Inc. System and method for fast wind flow measurement by LiDAR in a complex terrain
FR3111196B1 (fr) 2020-06-08 2022-06-10 Ifp Energies Now Procédé de détermination des composantes de la vitesse du vent au moyen d’un capteur de télédétection par Laser
CN111651896B (zh) * 2020-06-18 2021-08-24 浙江理工大学 一种基于实际风速及实际复杂地形的风电场流场计算方法
FR3111991B1 (fr) 2020-06-29 2022-06-17 Ifp Energies Now Procédé de détermination des composantes de la vitesse du vent au moyen d’un capteur de télédétection par Laser et au moyen d’une cohérence temporelle
CN112326999A (zh) * 2020-10-09 2021-02-05 北京观详光电技术有限公司 一种风速测量的方法、设备和终端
CN113177325B (zh) * 2021-05-17 2024-03-08 中国华能集团有限公司 复杂地形下校正标准k-ε模型可调参数方法、装置及存储介质
WO2022269655A1 (ja) * 2021-06-21 2022-12-29 三菱電機株式会社 風速予測装置、風速予測方法及びレーダ装置
KR102500311B1 (ko) * 2021-11-30 2023-02-14 제주대학교 산학협력단 풍속 데이터 기반의 기계학습 알고리즘을 활용한 타워 쉐도우 보정 방법 및 시스템
CN115436921B (zh) * 2022-11-10 2023-03-21 中国海洋大学 大气风场影响下的激光雷达飞机尾涡环量校正方法

Citations (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4735503A (en) * 1985-06-05 1988-04-05 Deutsche Forschungs- Und Versuchsanstalt Fur Luft-Und Raumfahrt Method for determining the direction and speed of wind in the atmosphere
JPH08166438A (ja) * 1994-12-15 1996-06-25 Kaijo Corp 風速測定用ドップラーレーダー
US20070256491A1 (en) * 2006-05-03 2007-11-08 The Boeing Company Accoustic profiler for wind, temperature, and turbulence
CN101389967A (zh) * 2004-03-26 2009-03-18 里索国家实验室 确定风轮机经受的风速和风向的方法和装置
CN102087130A (zh) * 2010-11-19 2011-06-08 北京工业大学 基于cfd技术的多声路超声流量计弯管安装声路的优化方法
CN102663251A (zh) * 2012-04-09 2012-09-12 华北电力大学 基于计算流体力学模型的风电场功率物理预测方法
US20130282285A1 (en) * 2011-01-05 2013-10-24 Leosphere Method and device for determining the movements of a fluid from remote measurements of radial velocities
JP2014238836A (ja) * 2013-06-10 2014-12-18 富士通株式会社 計算流体動力学モデルの準備のための方法及び装置
CN104268429A (zh) * 2014-10-15 2015-01-07 湖北大学 基于星载sar的近海岸海上风能资源遥感方法及系统
CN104331621A (zh) * 2014-11-05 2015-02-04 中国大唐集团新能源股份有限公司 一种风资源计算方法
CN104699936A (zh) * 2014-08-18 2015-06-10 沈阳工业大学 基于cfd短期风速预测风电场的扇区管理方法
DE102014216368A1 (de) * 2014-08-18 2016-02-18 Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. Mehrarmiges lidar-system
CN105608326A (zh) * 2015-12-31 2016-05-25 长沙理工大学 一种山区复杂地形风场大涡模拟入口边界条件输入方法
CN105675901A (zh) * 2014-12-05 2016-06-15 空中客车运营简化股份公司 用于估计飞行器空速的方法和设备
CN105792738A (zh) * 2013-12-04 2016-07-20 皇家飞利浦有限公司 用于改善功能狭窄分析的局部ffr估计和可视化
CN106644372A (zh) * 2016-12-28 2017-05-10 北京金风科创风电设备有限公司 检测风电机组的流体气动数据的方法和设备
CN107003411A (zh) * 2014-12-12 2017-08-01 三菱电机株式会社 激光雷达装置
CN107014451A (zh) * 2017-05-03 2017-08-04 东南大学 基于广义回归神经网络推测超声波流量传感器系数的方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8378518B2 (en) * 2009-03-26 2013-02-19 Terra Telesis, Inc. Wind power generator system, apparatus, and methods
CN102845088A (zh) 2010-04-16 2012-12-26 京瓷株式会社 无线通信系统、高功率基站、低功率基站及通信控制方法
DK2877741T3 (da) * 2012-07-27 2019-05-20 Univ Texas Tech System System og fremgangsmåde til at evaluere vindstrømningsfelter ved anvendelse af fjernregistreringsindretninger

Patent Citations (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4735503A (en) * 1985-06-05 1988-04-05 Deutsche Forschungs- Und Versuchsanstalt Fur Luft-Und Raumfahrt Method for determining the direction and speed of wind in the atmosphere
JPH08166438A (ja) * 1994-12-15 1996-06-25 Kaijo Corp 風速測定用ドップラーレーダー
CN101389967A (zh) * 2004-03-26 2009-03-18 里索国家实验室 确定风轮机经受的风速和风向的方法和装置
US20070256491A1 (en) * 2006-05-03 2007-11-08 The Boeing Company Accoustic profiler for wind, temperature, and turbulence
CN102087130A (zh) * 2010-11-19 2011-06-08 北京工业大学 基于cfd技术的多声路超声流量计弯管安装声路的优化方法
US20130282285A1 (en) * 2011-01-05 2013-10-24 Leosphere Method and device for determining the movements of a fluid from remote measurements of radial velocities
CN102663251A (zh) * 2012-04-09 2012-09-12 华北电力大学 基于计算流体力学模型的风电场功率物理预测方法
JP2014238836A (ja) * 2013-06-10 2014-12-18 富士通株式会社 計算流体動力学モデルの準備のための方法及び装置
CN105792738A (zh) * 2013-12-04 2016-07-20 皇家飞利浦有限公司 用于改善功能狭窄分析的局部ffr估计和可视化
DE102014216368A1 (de) * 2014-08-18 2016-02-18 Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. Mehrarmiges lidar-system
CN104699936A (zh) * 2014-08-18 2015-06-10 沈阳工业大学 基于cfd短期风速预测风电场的扇区管理方法
CN104268429A (zh) * 2014-10-15 2015-01-07 湖北大学 基于星载sar的近海岸海上风能资源遥感方法及系统
CN104331621A (zh) * 2014-11-05 2015-02-04 中国大唐集团新能源股份有限公司 一种风资源计算方法
CN105675901A (zh) * 2014-12-05 2016-06-15 空中客车运营简化股份公司 用于估计飞行器空速的方法和设备
CN107003411A (zh) * 2014-12-12 2017-08-01 三菱电机株式会社 激光雷达装置
CN105608326A (zh) * 2015-12-31 2016-05-25 长沙理工大学 一种山区复杂地形风场大涡模拟入口边界条件输入方法
CN106644372A (zh) * 2016-12-28 2017-05-10 北京金风科创风电设备有限公司 检测风电机组的流体气动数据的方法和设备
CN107014451A (zh) * 2017-05-03 2017-08-04 东南大学 基于广义回归神经网络推测超声波流量传感器系数的方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
HYUN-GOO KIM: "Correction of LiDAR measurement error in complex terrain by CFD: case study of the Yangyang pumped storage plant", 《WIND ENGINEERING》 *
KUMAR, GANESH: "Design Criteria and CFD analysis of a Small Wind Turbine Installed at Banihal Top in Snow Bound High Altitude Area: A Case Study", 《PROCEEDINGS OF 2014 1ST INTERNATIONAL CONFERENCE ON NON CONVENTIONAL ENERGY (ICONCE 2014)》 *
吕振峰等: "复杂地形对风速分布影响的数值模拟研究", 《能源工程》 *
张向旭等: "山区峡谷大跨度拱桥桥址区风参数数值模拟研究", 《铁道科学与工程学报》 *

Also Published As

Publication number Publication date
US10795054B2 (en) 2020-10-06
JP6896180B2 (ja) 2021-06-30
CN111868533B (zh) 2022-05-10
JP2021509718A (ja) 2021-04-01
US20190293836A1 (en) 2019-09-26
EP3769091A1 (en) 2021-01-27
WO2019181032A1 (en) 2019-09-26
EP3769091B1 (en) 2023-03-08

Similar Documents

Publication Publication Date Title
CN111868533B (zh) 风流感测系统和用于确定风流的速度场的方法
JP7403690B2 (ja) 複雑地形においてＬｉＤＡＲで高速の風の流れを測定するためのシステムおよび方法
Norton et al. Why atmospheric modeling is good science
Lundquist et al. An immersed boundary method enabling large-eddy simulations of flow over complex terrain in the WRF model
KR101914550B1 (ko) 레이더의 표적 위치 추적 방법
JP5961188B2 (ja) 流体の移動の半径方向速度の遠隔測定から該流体の運動を決定する方法及び装置
CN110389391A (zh) 一种基于空间域的重磁位场解析延拓方法
US20120120230A1 (en) Apparatus and Method for Small Scale Wind Mapping
JP7351029B2 (ja) 複雑地形においてＬｉＤＡＲで風の流れの乱流を測定するためのシステムおよび方法
Zaron et al. Bottom topography mapping via nonlinear data assimilation
Mukherjee et al. Accuracy of Cartosat-1 DEM and its derived attribute at multiple scale representation
CN108709956A (zh) 基于落球定位信息的大气参数测量方法和设备
CN106323271A (zh) 基于特征奇异值的航天器相对姿态测量矢量选取方法
Verkhoglyadova et al. Estimation of winds from GPS radio occultations
Chandrasekar et al. Assimilation of multi-channel radiances in mesoscale models with an ensemble technique to improve track forecasts of Tropical cyclones
JP7470413B2 (ja) 情報処理装置、情報処理方法、及びプログラム
Le Dimet et al. Variational data assimilation: Optimization and optimal control
Plu et al. A cyclogenesis evolving into two distinct scenarios and its implications for short-term ensemble forecasting
CN113656918A (zh) 一种应用于成品高架库场景的四旋翼仿真测试方法
Lin et al. Modeling and analysis of light detection and ranging point cloud error in vibration state of airborne platform
Rogers et al. Methodology Improvements for Three-Dimensional UAV-Based Travel-Time Acoustic Atmospheric Tomography
Yang et al. High-fidelity LiDAR Simulation System Based on Real Pointcloud Features
Crasto et al. Accuracy quantification of a CFD model in complex site by cross-checking wind measurements
Ollinaho Feasibility of assimilating ASCAT surface winds into a Limited Area Model
Wang et al. An Improved Ensemble-Based Data Assimilation Method with Adaptive Observation Error Adjustment for Wildland Fire Spread Prediction

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