CN112287296A - 一种基于双波段闪烁仪的地表水热通量测算方法 - Google Patents

一种基于双波段闪烁仪的地表水热通量测算方法 Download PDF

Info

Publication number
CN112287296A
CN112287296A CN202011088336.XA CN202011088336A CN112287296A CN 112287296 A CN112287296 A CN 112287296A CN 202011088336 A CN202011088336 A CN 202011088336A CN 112287296 A CN112287296 A CN 112287296A
Authority
CN
China
Prior art keywords
heat flux
calculating
humidity
mws
temperature
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
CN202011088336.XA
Other languages
English (en)
Other versions
CN112287296B (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.)
Beijing Normal University
Original Assignee
Beijing Normal University
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 Beijing Normal University filed Critical Beijing Normal University
Priority to CN202011088336.XA priority Critical patent/CN112287296B/zh
Publication of CN112287296A publication Critical patent/CN112287296A/zh
Application granted granted Critical
Publication of CN112287296B publication Critical patent/CN112287296B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01K—MEASURING TEMPERATURE; MEASURING QUANTITY OF HEAT; THERMALLY-SENSITIVE ELEMENTS NOT OTHERWISE PROVIDED FOR
    • G01K17/00—Measuring quantity of heat
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/21—Design, administration or maintenance of databases
    • G06F16/215—Improving data quality; Data cleansing, e.g. de-duplication, removing invalid entries or correcting typographical errors
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24—Querying
    • G06F16/245—Query processing
    • G06F16/2457—Query processing with adaptation to user needs

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Algebra (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Software Systems (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Quality & Reliability (AREA)
  • Computational Linguistics (AREA)
  • Investigating Or Analyzing Materials Using Thermal Means (AREA)

Abstract

本发明提供了一种基于双波段闪烁仪的地表水热通量测算方法，用以解决现有技术中计算大尺度地表水热通量过程中的不完整和不准确问题。所述测算方法包括：S1，基于双波段闪烁仪获取高频数据；S2，将所述高频数据预处理后，计算光强方差；S3，根据所述光强方差，计算空气折射指数的结构参数；S4，根据所述空气折射指数的结构参数，计算气象结构参数；S5，根据所述气象结构参数，计算地表水热通量。本发明在地表水热通量测算方法中，在高频数据预处理的基础上，进行局地化参数计算，数据筛选，并进行大气稳定度判断，选取合适的大气稳定度普适函数，计算过程完整，从而获取更加准确的感热通量和潜热通量。

Description

一种基于双波段闪烁仪的地表水热通量测算方法
技术领域
本发明属于大气和水文领域，具体涉及一种基于双波段闪烁仪的地表水热通量测算方法。
背景技术
在大气、水文领域中，地表水热通量(感热通量、潜热通量/蒸散发)是一个可以直观反映地表水热状况的重要参量，其中大尺度(一般1～5km)地面观测值可以用来验证和校正地表蒸散发遥感产品，同时，也可以对陆面过程模型、水文模型等进行参数率定和模拟结果验证。
目前，大尺度地表水热通量的观测方法主要包括涡动相关仪的观测矩阵、机载涡动相关仪以及闪烁仪，其中闪烁仪是观测公里级地表水热通量的最佳手段。
闪烁仪能够测量几百米至几公里的大尺度地表水热通量，其中，大孔径闪烁仪(Large Aperture Scintillometer,LAS)较早得到应用，但是大孔径闪烁仪的发射波段(近红外波段)对温度敏感，只能观测获取公里级的感热通量，潜热通量需通过能量平衡方程(余项法)获取。由于微波波段对于水汽更加敏感，由此微波闪烁仪(MicrowaveScintillometer,MWS)逐渐兴起，并与大孔径闪烁仪组合成为双波段闪烁仪(Optical-Microwave Scintillometer,OMS)，为直接、准确获取公里尺度的平均水热通量(感热和潜热通量)提供了新的观测手段，成为观测水热通量的重要仪器之一。
目前双波段闪烁仪通常在草地、农田或森林下垫面架设，通过相关理论和计算方法，获取感热通量和潜热通量，并与涡动相关系统观测结果进行对比。但在不同的下垫面条件下均发现双波段闪烁仪计算结果仍然存在许多的不确定性，其中，在草地下垫面观测的感热通量比涡动相关仪(Eddy Covariance System, EC)观测结果偏小16％，潜热通量比EC观测结果偏大4％；在森林下垫面观测的感热通量比EC观测结果整体偏大1％，潜热通量比EC观测结果整体偏大16％；在森林下垫面观测的感热通量比EC观测结果偏大5％，潜热通量比EC观测结果偏大28％。所以急需一种全面的、系统的、基于双波段闪烁仪的地表水热通量测算方法，特别是在高频数据的预处理与数据筛选、局地化参数计算、大气稳定度的判断与普适函数的选取等测算方面需要进一步完善，用于处理双波段闪烁仪数据，以便获取更加准确的感热通量和潜热通量。
发明内容
本发明实施例的目的是提出一种基于双波段闪烁仪的地表水热通量测算方法，是一种完整的、处理双波段闪烁仪数据的计算过程，在高频数据预处理和数据筛选的基础上，进行局地化参数计算，并进行大气稳定度判断和选取合适的大气稳定度普适函数，计算过程完整，所获得的地表水热通量比较准确。
为了实现上述目的，本发明实施例采用的技术方案如下：
一种基于双波段闪烁仪的地表水热通量测算方法，所述方法包括如下步骤：
步骤S1，基于双波段闪烁仪获取高频数据；
步骤S2，将所述高频数据预处理后，计算光强方差；
步骤S3，根据所述光强方差，计算空气折射指数的结构参数；
步骤S4，根据所述空气折射指数的结构参数，计算气象结构参数；
步骤S5，根据所述气象结构参数，计算地表水热通量。
作为本发明的一个优选实施例，所述双波段闪烁仪高频数据的预处理，包括对高频数据去野点、和\或高频滤波处理、和\或低频滤波处理。
作为本发明的一个优选实施例，所述高频数据去野点，采用数据平均绝对偏差方法去除野点值，遍历所有数据点，当一个窗口一个点值超过绝对偏差中位数的10倍，该点被识别为一个野点值，并被剔除，直至没有数据点被识别为野点值，过程停止，完成去野点环节。
作为本发明的一个优选实施例，所述步骤S2中计算光强方差，具体通过式 (1)～(3)计算：
Figure BDA0002721156500000021
Figure BDA0002721156500000022
Figure BDA0002721156500000023
式(1)～(3)中，
Figure BDA0002721156500000024
分别表示LAS、MWS、OMS的光强方差；ILAS、IMWS分别表示LAS、MWS的光强信号；
Figure BDA0002721156500000025
分别表示LAS、MWS闪烁信号在平滑窗口的平均值，上划线表示均值。
作为本发明的一个优选实施例，通过式(4)～(6)计算空气折射指数的结构参数：
Figure BDA0002721156500000031
Figure BDA0002721156500000032
Figure BDA0002721156500000033
式(4)～(6)中，
Figure BDA0002721156500000034
CoOMS分别表示LAS、MWS、OMS的空气折射指数结构参数；GintLAS、GintMWS、GintOMS分别表示LAS、MWS、 OMS的Gint函数。
作为本发明的一个优选实施例，计算气象结构参数，先进行数据筛选，再根据剩余数据计算温度结构参数
Figure BDA0002721156500000035
湿度结构参数
Figure BDA0002721156500000036
和温湿度结构参数
Figure BDA0002721156500000037
作为本发明的一个优选实施例，所述数据筛选，剔除降水时刻的数据、剔除饱和数据、剔除信号强度较弱的数据、剔除弱湍流数据。
作为本发明的一个优选实施例，通过双波长TW方法计算气象结构参数，设定温度和湿度折射指数的系数
Figure BDA0002721156500000038
为常数，通过公式(7)～(9)求解
Figure BDA0002721156500000039
和
Figure BDA00027211565000000310
Figure BDA00027211565000000311
Figure BDA00027211565000000312
Figure BDA00027211565000000313
设定
Figure BDA00027211565000000314
求解公式(7)～(9)得到：
Figure 1
Figure BDA0002721156500000041
Figure BDA0002721156500000042
式(10)～(11)中，AT,LAS、Aq,LAS分别表示LAS的温度和湿度折射指数的系数，AT,MWS、Aq,MWS分别表示MWS的温度和湿度折射指数的系数。
作为本发明的一个优选实施例，通过双波段交互BC方法计算气象结构参数，实时计算温度和湿度折射指数的系数RTq，通过求解空气折射指数结构参数和气象结构参数的关系式：
Figure BDA0002721156500000043
Figure BDA0002721156500000044
Figure BDA0002721156500000045
求解公式(12)～(14)得到：
Figure BDA0002721156500000046
Figure BDA0002721156500000051
Figure BDA0002721156500000052
式(15)～(17)中，AT,LAS、Aq,LAS分别表示LAS的温度和湿度折射指数的系数，AT,MWS、Aq,MWS分别表示MWS的温度和湿度折射指数的系数；
其中，LAS与MWS用同一公式计算温度和湿度折射指数的系数AT、Aq，计算公式如下：
Figure BDA0002721156500000053
Figure BDA0002721156500000054
式(18)～(19)中，P表示大气压，单位为Pa；T表示空气温度，单位为 K；q为空气比湿，单位为g/g；R表示湿空气的空气常数，单位为J/(mol·K)； Rd、RV分别表示干空气和水蒸汽的气体常数，单位为J/(mol·K)；bt1、bt2、bq2的单位都是K Pa-1，与光束的波长有关；上划线表示均值。
作为本发明的一个优选实施例，所述步骤S5计算地表水热通量，根据莫宁 -奥布霍夫相似理论，进行大气稳定度判断，选取合适的大气稳定度普适函数，通过迭代计算摩擦风速、特征温度、特征湿度，进而计算感热通量和潜热通量。
作为本发明的一个优选实施例，选取湿度折射指数的系数RTq和净辐射Rn 共同判断大气稳定状况，选取Li函数作为大气稳定度普适函数。
作为本发明的一个优选实施例，根据莫宁-奥布霍夫相似理论，温度结构参数、湿度结构参数与特征温度、特征湿度有如下关系：
Figure BDA0002721156500000055
Figure BDA0002721156500000056
式(20)～(21)中，z是有效高度，单位为m；d0是零平面位移高度，单位为m；T*是特征温度，单位为K；q*是特征湿度，为百分数；fT、fq为普适函数；LMO为奥布霍夫长度。
作为本发明的一个优选实施例，奥布霍夫长度和摩擦风速计算公式如下：
Figure BDA0002721156500000061
Figure BDA0002721156500000062
式(22)和(23)中，g为重力加速度，单位为m/s2；u*为摩擦风速，单位为m/s；k为vonKarman常数；Z0为地表粗糙度，单位为m；Zu为风速的测量高度，单位为m；上划线为均值；
除了Gint函数、温度和湿度结构参数的系数、有效高度等参数局地化计算，还包括地表粗糙度、零平面位移高度的局地化计算，计算公式如下：
d0≈2/3hveg……………………………………(24)
Z0≈0.13hveg……………………………………(25)
式(24)和(25)中，hveg为植被株高，单位为m。
作为本发明的一个优选实施例，通过迭代得出摩擦风速u*、特征温度T*、特征湿度q*，最终得到感热通量H和潜热通量LE：
Figure BDA0002721156500000063
Figure BDA0002721156500000064
式(26)和(27)中，
Figure BDA0002721156500000065
为平均空气密度，单位为kg/m3；
Figure BDA0002721156500000066
为平均比湿，单位为g/g；Cp为空气比热，单位为J/(kg·K)；Lv是蒸发潜热，单位为KJ kg-1； H为感热通量，单位为W/m2；LE为潜热通量，单位为W/m2。
本发明具有如下有益效果：
在本发明中，通过双波段闪烁仪高频数据预处理，剔除高频数据中的野点值和噪声，计算更准确的光强方差；通过计算站点的有效高度和Gint函数，计算更准确的空气折射指数的结构参数；通过计算实时的温度和湿度结构参数的系数，可获取更加准确的气象结构参数；通过计算地表粗糙度和零平面位移高度以及选取更加合适的函数作为大气稳定度普适函数，可获得更加准确的地表水热通量(感热通量和潜热通量)。进而为验证和校正地表蒸散发遥感产品提供数据支持，同时，也可以对陆面过程模型、水文模型等进行参数率定和模拟结果验证。具体地，基于本发明已经处理两种不同地表类型观测站的双波段闪烁仪数据，并且得到较好的感热通量和潜热通量结果。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例基于双波段闪烁仪的地表水热通量测算方法流程示意图；
图2为本发明实施例测算方法中观测站点的位置和仪器布设图：(a)黑河流域；(b)阿柔站观测仪器布设概况；(c)大满站观测仪器布设概况；
图3为本发明实施例以阿柔站为例的地表水热通量测算方法详细流程图；
图4为本发明实施例中阿柔站2019年6月2日～2019年6月5日OMS与 EC观测的感热通量和潜热通量比较图；
图5为本发明实施例中大满站2020年5月4日～2020年5月7日OMS与 EC观测的感热通量和潜热通量比较图。
具体实施方式
下面通过参考示范性实施例，对本发明技术问题、技术方案和优点进行详细阐明。以下所述示范性实施例仅用于解释本发明，而不能解释为对本发明的限制。本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非在这里进行定义，否则不会用理想化或过于正式的含义来解释。
下面通过两个具体的观测站点的测算过程为例，对本发明实施例作进一步详细说明。
实施例一
以图2中的阿柔站为例，详细说明一下本发明实施例中地表水热通量的测算过程。如图2所示，阿柔站(100°27′53″E,38°2′40″N)位于青海省祁连县阿柔乡的河谷区域，海拔3033m，属高原山地气候，下垫面为高寒草甸，比较平坦与均匀。2019年5月17日，阿柔站架设了德国微波闪烁仪(RPG-MWSC-160，德国)，与大孔径闪烁仪(德国BLS900)结合组成双波段闪烁仪。闪烁仪西南-东北方向架设，光径路线长度为2390m，闪烁仪发射和接收端安装高度均为15.6m，采样频率为1000Hz。在闪烁仪的路径中间配有一套气象要素梯度观测系统(28m塔)、一套涡动相关仪以及物候相机、宇宙射线土壤水分测定仪、称重式雨量计等，并在闪烁仪源区内安装了1套土壤水分传感器网络(8个节点)。表1中给出了阿柔站的站点信息以及参数概况，表2中给出了阿柔站双波段闪烁仪的局地化参数值。
表1
Figure BDA0002721156500000081
表2
Figure BDA0002721156500000091
图3为以阿柔站为例的地表水热通量测算方法详细流程图。
如图1和图3所示，所述测算过程如下：
步骤S1，基于双波段闪烁仪获取高频数据。
本实施例中，所述双波段闪烁仪采用OMS(型号：BLS900&RPG-MWSC-160) 双波段闪烁仪，来获取站点观测区域的高频数据。
步骤S2，将所述高频数据预处理后，计算光强方差。
本步骤中，所述高频数据的预处理，包括对高频数据去野点、和\或高频滤波处理、和\或低频滤波处理。
优选地，本步骤中采用数据平均绝对偏差(Mean Absolute Deviation,MAD) 方法去除野点值，具体的，所述MAD方法，遍历所有数据点，当一个窗口一个点值超过绝对偏差中位数的10倍，该点被识别为一个野点值，并被剔除，直至没有数据点被识别为野点值，这个过程才会停止，完成去野点环节。
本步骤中，完成高频数据的预处理后，通过式(1)～(3)计算光强方差
Figure BDA0002721156500000092
Figure BDA0002721156500000093
Figure BDA0002721156500000094
Figure BDA0002721156500000095
式(1)-(3)中，
Figure BDA0002721156500000096
分别表示LAS、MWS、OMS的光强方差；ILAS、IMWS分别表示LAS、MWS的光强信号；
Figure BDA0002721156500000097
分别表示LAS、 MWS闪烁信号在平滑窗口的平均值，上划线表示均值。优选地，LAS信号的平滑窗口为12s，MWS信号的平滑窗口为67s。上述公式计算的过程，即高频滤波处理后，求取的光强方差再通过低频滤波因子进行处理。
上述光强方差的计算过程，通过高频滤波和低频滤波移除了噪声和水汽吸收的影响。其中，所述高频滤波通过公式(1)～(3)实现，低频滤波通过 Matlab2018a中smooth库实现。
步骤S3，根据所述光强方差，计算空气折射指数的结构参数。
本步骤中，所述计算空气折射指数的结构参数，首先获取站点的A系列参数，包括AT,LAS、Aq,LAS、AT,MWS、Aq,MWS，根据A系列参数计算空气折射指数的结构参数。
对应于LAS、MWS、OMS的三个空气折射指数的结构参数为
Figure BDA0002721156500000101
CoOMS，通过式(4)～(6)计算空气折射指数的结构参数：
Figure BDA0002721156500000102
Figure BDA0002721156500000103
Figure BDA0002721156500000104
式(4)～(6)中，
Figure BDA0002721156500000105
CoOMS分别表示LAS、MWS、OMS的空气折射指数结构参数；GintLAS、GintMWS、GintOMS分别表示LAS、MWS、 OMS的Gint函数。
其中，Gint函数在不同的站点是不同的，其主要的影响因子是近红外闪烁仪和MWS的波长、孔径直径、近红外闪烁仪与微波闪烁仪的间距和光径路线长度。Gint函数值计算公式如下：
Figure BDA0002721156500000106
Figure BDA0002721156500000107
Figure BDA0002721156500000108
Figure BDA0002721156500000109
Figure BDA00027211565000001010
Figure BDA00027211565000001011
Figure BDA00027211565000001012
Figure BDA00027211565000001013
Figure BDA00027211565000001014
Figure BDA0002721156500000111
式(7)-(16)计中，DtLAS、DrLAS分别表示LAS发射端和接收端的孔径直径；DtMWS、DrMWS分别表示MWS发射端和接收端的孔径直径；d(x)表示大孔径闪烁仪与微波闪烁仪的间距；N为光径路线长度；k为空间波数；kLAS、kMWS、 kOMS分别表示LAS、MWS、OMS空间波数(k＝2π/λ)；Jo、J1都表示贝塞尔函数，下标是阶数；x为在光径路线上的位置。
步骤S4，根据所述空气折射指数的结构参数，计算气象结构参数。
所述计算气象结构参数，先进行数据筛选，再根据剩余数据计算温度结构参数
Figure BDA0002721156500000112
湿度结构参数
Figure BDA0002721156500000113
和温湿度结构参数
Figure BDA0002721156500000114
进行数据筛选，剔除降水时刻的数据、饱和数据、信号强度较弱的数据和弱湍流数据。
本步骤中，所述进行大气稳定度判断，选取温度和湿度折射指数的系数Rtq和净辐射Rn共同判断大气稳定状况，可以更好地反映大气稳定度。
所述温度结构参数(
Figure BDA0002721156500000115
)、湿度结构参数(
Figure BDA0002721156500000116
)、温湿度结构参数(
Figure BDA0002721156500000117
)，通过双波长方法(Two-wavelength,TW)或双波段交互方法 (Bichromatic-correlation,BC)得到。
其中，所述TW方法假设温度和湿度折射指数的系数(
Figure BDA0002721156500000118
)为稳定的常数，结合空气折射指数的结构参数获取
Figure BDA0002721156500000119
进一步地，在TW方法中，空气折射指数结构参数和气象结构参数的公式如下所示：
Figure BDA00027211565000001110
Figure BDA00027211565000001111
Figure BDA00027211565000001112
设定
Figure BDA00027211565000001113
求解公式(17)～(19)得到：
Figure 2
Figure BDA0002721156500000121
Figure BDA0002721156500000122
而所述BC方法，实时计算温度和湿度折射指数的系数
Figure BDA0002721156500000123
结合空气折射指数的结构参数获取
Figure BDA0002721156500000124
进一步地，在BC方法中，空气折射指数结构参数和气象结构参数的公式如下所示：
Figure BDA0002721156500000131
Figure BDA0002721156500000132
Figure BDA0002721156500000133
根据公式(22)～(24)求解得到：
Figure BDA0002721156500000134
Figure BDA0002721156500000135
Figure BDA0002721156500000136
式(25)～(27)中，AT,LAS、Aq,LAS分别表示LAS的温度和湿度折射指数的系数，AT,MWS、Aq,MWS分别表示MWS的温度和湿度折射指数的系数，A系列参数与空气温度、湿度相关，需要实时计算，LAS与MWS用同一公式计算温度和湿度折射指数的系数，计算公式如下：
Figure BDA0002721156500000137
Figure BDA0002721156500000138
式(28)～(29)中，P表示大气压(Pa)；T表示空气温度(K)；q为空气比湿；R表示湿空气的空气常数；Rd、RV分别表示干空气和水蒸汽的气体常数； bt1、bt2、bq2的单位都是K Pa-1，与光束的波长有关；上划线表示均值。
步骤S5，根据所述气象结构，计算地表水热通量。
本步骤中，所述地表水热通量包括感热通量和潜热通量。优选地，根据莫宁-奥布霍夫相似理论，进行大气稳定度判断，选取合适的大气稳定度普适函数，通过迭代计算摩擦风速u*、特征温度T*、特征湿度q*，进而计算感热通量和潜热通量。
优选地，本步骤中，选取温度和湿度折射指数的系数Rtq和净辐射Rn共同判断大气稳定状况，可以更好地反映大气稳定度。选取Li函数作为稳定度函数进行计算，不稳定条件为：
Figure BDA0002721156500000141
稳定条件为：
Figure BDA0002721156500000142
其中，根据莫宁-奥布霍夫相似理论，温度结构参数、湿度结构参数与特征温度、特征湿度有如下关系：
Figure BDA0002721156500000143
Figure BDA0002721156500000144
式(29)～(31)中，z是有效高度(m)；d0是零平面位移高度(m)；T*温度特征量(K)；q*温度特征量(％)；fT、fq为普适函数；LMO为奥布霍夫长度。
奥布霍夫长度(LMO)和摩擦风速(u*)计算公式如下：
Figure BDA0002721156500000145
Figure BDA0002721156500000146
式(32)和(33)中，g为重力加速度(m/s2)；u*为摩擦风速(m/s)；k为 von Karman常数；Z0为空气动力学粗糙度(m)；Zu为风速的测量高度(m)；上划线为均值。式中，除了Gint函数、温度和湿度结构参数的系数、有效高度等参数局地化计算，还包括地表粗糙度、零平面位移高度的局地化计算，计算公式如下：
d0≈2/3hveg……………………………………(34)
Z0≈0.13hveg……………………………………(35)
式(34)和(35)中，hveg为植被株高，单位为m。
通过公式(30)～(35)迭代得出特征温度(T*)、特征湿度(q*)、摩擦风速(u*)，从而计算感热通量(H)和潜热通量(LE)，公式如下：
Figure BDA0002721156500000147
Figure BDA0002721156500000151
式(36)和(37)中，
Figure BDA0002721156500000152
为平均空气密度(kg/m3)；
Figure BDA0002721156500000153
为平均比湿；Cp为空气比热(J/(kg·K))；Lv是蒸发潜热(KJ·kg-1)；H为感热通量(W/m2)；LE为潜热通量(W/m2)。
利用此发明处理阿柔站2019年6月2日～2019年6月4日OMS数据，得到相应的感热和潜热通量结果。图4为阿柔站2019年6月2日～2019年6月5 日OMS与EC观测的感热和潜热通量比较图。如图4所示，无论是双波长方法还是双波段交互方法，OMS获取的感热通量和潜热通量的变化趋势和量级与EC 观测值均一致。
实施例二
再以图2中的黑河流域中游的大满站为例，如图2所示，大满站(100°22′ 20″E,38°51′20″N)位于甘肃省张掖市大满灌区内，海拔1556m，属温带大陆性气候，下垫面以农田下垫面为主，包括制种玉米、温室大棚、村庄和防护林等。2019年7月15日，大满站架设双波段闪烁仪(大孔径闪烁仪型号BLS900 和微波闪烁仪型号RPG-MWSC-160结合组成OMS)，该站闪烁仪架设方向为西南-东北方向，光径长度1854m，发射端安装高度为15.6m，接收端安装高度为 25.6m，采样频率为1000Hz。在闪烁仪路径中间配有一套气象要素梯度观测系统(40m塔)，两套涡动相关仪以及物候相机、宇宙射线土壤水分测定仪，并在闪烁仪源区内各安装了1套叶面积指数(6个上节点、28个下节点)和土壤水分传感器网络(10个节点)等。表1中给出了大满站的站点信息以及参数概况，表2中给出了大满站双波段闪烁仪的局地化参数值。
以大满站为例的地表水热通量测算方法详细过程与阿柔站的测算过程相同 (图3)。利用此发明处理大满站2020年5月4日～2020年5月7日OMS数据，得到相应的感热和潜热通量结果。图5为大满站2020年5月4日～2020年5月 7日OMS与EC观测的感热通量和潜热通量比较图。如图5所示，在大满站，无论是双波长方法还是双波段交互方法，OMS获取的感热通量和潜热通量变化趋势和量级也与EC观测值比较一致。
总体来说，基于本发明，阿柔站OMS与EC观测的感热通量和潜热通量变化趋势和量级基本一致；大满站OMS与EC观测的感热通量和潜热通量变化趋势和量级也基本一致。以上结果表明，本发明实施例基于双波段闪烁仪的地表水热通量测算方法，可以很好地对不同地表类型的水热通量进行测算，所获得的测算结果较为准确。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同或相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于装置或系统实施例而言，由于其基本相似于方法实施例，所以描述得比较简单，相关之处参见方法实施例的部分说明即可。以上所描述的装置及系统实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。
以上所述是本发明的优选实施方式，应当指出，本发明并不受限于以上所公开的示范性实施例，说明书的实质仅仅是帮助相关领域技术人员综合理解本发明的具体细节。对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，在本发明揭露的技术范围做出的若干改进和润饰、可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

Claims (14)

1.一种基于双波段闪烁仪的地表水热通量测算方法，其特征在于，所述方法包括如下步骤：
步骤S1，基于双波段闪烁仪获取高频数据；
步骤S2，将所述高频数据预处理后，计算光强方差；
步骤S3，根据所述光强方差，计算空气折射指数的结构参数；
步骤S4，根据所述空气折射指数的结构参数，计算气象结构参数；
步骤S5，根据所述气象结构参数，计算地表水热通量。
2.根据权利要求1所述的地表水热通量测算方法，其特征在于，所述双波段闪烁仪高频数据的预处理，包括对高频数据去野点、和\或高频滤波处理、和\或低频滤波处理。
3.根据权利要求2所述的地表水热通量测算方法，其特征在于，所述高频数据去野点，采用数据平均绝对偏差方法去除野点值，遍历所有数据点，当一个窗口一个点值超过绝对偏差中位数的10倍，该点被识别为一个野点值，并被剔除，直至没有数据点被识别为野点值，过程停止，完成去野点环节。
4.根据权利要求1至3任一项所述的地表水热通量测算方法，其特征在于，所述步骤S2中计算光强方差，具体通过式(1)～(3)计算：
Figure FDA0002721156490000011
Figure FDA0002721156490000012
Figure FDA0002721156490000013
式(1)～(3)中，
Figure FDA0002721156490000014
分别表示大孔径闪烁仪LAS、微波闪烁仪MWS、双波段闪烁仪OMS的光强方差；ILAS、IMWS分别表示LAS、MWS的光强信号；
Figure FDA0002721156490000015
分别表示LAS、MWS闪烁信号在平滑窗口的平均值，上划线表示均值。
5.根据权利要求4所述的地表水热通量测算方法，其特征在于，通过式(4)～(6)计算空气折射指数的结构参数：
Figure FDA0002721156490000016
Figure FDA0002721156490000017
Figure FDA0002721156490000018
式(4)～(6)中，
Figure FDA0002721156490000021
CoOMS分别表示LAS、MWS、OMS的空气折射指数结构参数；GintLAS、GintMWS、GintOMS分别表示LAS、MWS、OMS的Gint函数。
6.根据权利要求5所述的地表水热通量测算方法，其特征在于，计算气象结构参数，先进行数据筛选，再根据剩余数据计算气象结构参数，包括温度结构参数
Figure FDA0002721156490000022
湿度结构参数
Figure FDA0002721156490000023
和温湿度结构参数CTq。
7.根据权利要求6所述的地表水热通量测算方法，其特征在于，进行数据筛选：剔除降水时刻的数据、剔除饱和数据、剔除信号强度较弱的数据、剔除弱湍流数据。
8.根据权利要求7所述的地表水热通量测算方法，其特征在于，通过双波长TW方法计算气象结构参数，设定温度和湿度折射指数的系数RTq为常数，通过公式(7)～(9)求解
Figure FDA0002721156490000024
和CTq：
Figure FDA0002721156490000025
Figure FDA0002721156490000026
Figure FDA0002721156490000027
设定RTq＝±0.8，求解公式(7)～(9)得到：
Figure FDA0002721156490000028
Figure FDA0002721156490000031
Figure FDA0002721156490000032
式(7)～(11)中，AT，LAS、Aq，LAS分别表示LAS的温度和湿度折射指数的系数，AT，MWS、Aq，MWS分别表示MWS的温度和湿度折射指数的系数。
9.根据权利要求8所述的地表水热通量测算方法，其特征在于，通过双波段交互BC方法计算气象结构参数，实时计算温度和湿度折射指数的系数RTq，通过求解空气折射指数结构参数和气象结构参数的关系式(12)～(14)：
Figure FDA0002721156490000033
Figure FDA0002721156490000034
Figure FDA0002721156490000035
得到计算气象结构参数、温度和湿度折射指数的系数RTq：
Figure FDA0002721156490000036
Figure FDA0002721156490000037
Figure FDA0002721156490000041
Figure FDA0002721156490000042
式(15)～(17)中，AT，LAS、Aq，LAS分别表示LAS的温度和湿度折射指数的系数，AT，MWS、Aq，MWS分别表示MWS的温度和湿度折射指数的系数；
其中，LAS与MWS用同一公式计算温度和湿度折射指数的系数AT、Aq，计算公式如下：
Figure FDA0002721156490000043
Figure FDA0002721156490000044
式(18)～(19)中，P表示大气压，单位为Pa；T表示空气温度，单位为K；q为空气比湿，单位为g/g；R表示湿空气的空气常数，单位为J/(mol·K)；Rd、RV分别表示干空气和水蒸汽的气体常数，单位为J/(mol·K)；bt1、bt2、bq2的单位都是KPa-1；上划线表示均值。
10.根据权利要求9所述的地表水热通量测算方法，其特征在于，所述步骤S5计算地表水热通量，根据莫宁-奥布霍夫相似理论，进行大气稳定度判断，选取合适的大气稳定度普适函数，通过迭代计算摩擦风速、特征温度、特征湿度，进而计算感热通量和潜热通量。
11.根据权利要求10所述的地表水热通量测算方法，其特征在于，选取湿度折射指数的系数RTq和净辐射Rn共同判断大气稳定状况，选取Li函数作为大气稳定度普适函数。
12.根据权利要求11所述的地表水热通量测算方法，其特征在于，根据莫宁-奥布霍夫相似理论，温度结构参数、湿度结构参数与特征温度、特征湿度有如下关系：
Figure FDA0002721156490000045
Figure FDA0002721156490000046
式(20)～(21)中，z是有效高度，单位为m；d0是零平面位移高度，单位为m；T*是特征温度，单位为K；q*是特征湿度，为百分数；fT、fq为普适函数；LMO为奥布霍夫长度。
13.根据权利要求12所述的地表水热通量测算方法，其特征在于，奥布霍夫长度和摩擦风速计算公式如下：
Figure FDA0002721156490000051
Figure FDA0002721156490000052
式(22)和(23)中，g为重力加速度，单位为m/s2；u*为摩擦风速，单位为m/s；k为vonKarman常数；Z0为地表粗糙度，单位为m；Zu为风速的测量高度，单位为m；上划线为均值；
除了Gint函数、温度和湿度结构参数的系数、有效高度等参数局地化计算，还包括地表粗糙度、零平面位移高度的局地化计算，计算公式如下：
d0≈2/3hveg..........................................(24)
Z0≈0.13hveg..........................................(25)
式(24)和(25)中，hveg为植被株高，单位为m。
14.根据权利要求13所述的地表水热通量测算方法，其特征在于，通过迭代得出摩擦风速u*、特征温度T*、特征湿度q*，最终得到感热通量H和潜热通量LE：
Figure FDA0002721156490000053
Figure FDA0002721156490000054
式(26)和(27)中，
Figure FDA0002721156490000055
为平均空气密度，单位为kg/m3；
Figure FDA0002721156490000056
为平均比湿，单位为g/g；Cp为空气比热，单位为J/(kg·K)；Lv是蒸发潜热，单位为KJ kg-1；H为感热通量，单位为W/m2；LE为潜热通量，单位为W/m2。
CN202011088336.XA 2020-10-13 2020-10-13 一种基于双波段闪烁仪的地表水热通量测算方法 Active CN112287296B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011088336.XA CN112287296B (zh) 2020-10-13 2020-10-13 一种基于双波段闪烁仪的地表水热通量测算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011088336.XA CN112287296B (zh) 2020-10-13 2020-10-13 一种基于双波段闪烁仪的地表水热通量测算方法

Publications (2)

Publication Number Publication Date
CN112287296A true CN112287296A (zh) 2021-01-29
CN112287296B CN112287296B (zh) 2023-05-26

Family

ID=74496101

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011088336.XA Active CN112287296B (zh) 2020-10-13 2020-10-13 一种基于双波段闪烁仪的地表水热通量测算方法

Country Status (1)

Country Link
CN (1) CN112287296B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113156542A (zh) * 2021-03-12 2021-07-23 中国气象局气象探测中心 一种基于闪烁法测量区域蒸散总水量的方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS5974212A (ja) * 1982-10-20 1984-04-26 Kawasaki Steel Corp 鉄損の少ない無方向性電磁鋼板の製造に供する溶鋼の取鍋精錬方法
WO2007122256A1 (fr) * 2006-04-26 2007-11-01 Commissariat A L'energie Atomique Procede de preparation d'une couche nanoporeuse de nanoparticules et couche ainsi obtenue
CN102192885A (zh) * 2010-03-08 2011-09-21 北京师范大学 一种大尺度水热通量观测系统
US20160176186A1 (en) * 2014-12-19 2016-06-23 Palo Alto Research Center Incorporated Variable data lithography system with embedded plasmonic fillers in a printing plate
CN107655564A (zh) * 2017-05-11 2018-02-02 南京邮电大学 一种基于智能终端的多种技术融合的室内外环境检测方法
WO2018138340A1 (fr) * 2017-01-30 2018-08-02 Extreme Wheather Expertises - Exwexs Procédé et dispositif de météorologie et produit programme d'ordinateur associé
CN111368258A (zh) * 2020-03-04 2020-07-03 中国科学院东北地理与农业生态研究所 一种湿润地区日蒸散量的估算方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS5974212A (ja) * 1982-10-20 1984-04-26 Kawasaki Steel Corp 鉄損の少ない無方向性電磁鋼板の製造に供する溶鋼の取鍋精錬方法
WO2007122256A1 (fr) * 2006-04-26 2007-11-01 Commissariat A L'energie Atomique Procede de preparation d'une couche nanoporeuse de nanoparticules et couche ainsi obtenue
CN102192885A (zh) * 2010-03-08 2011-09-21 北京师范大学 一种大尺度水热通量观测系统
US20160176186A1 (en) * 2014-12-19 2016-06-23 Palo Alto Research Center Incorporated Variable data lithography system with embedded plasmonic fillers in a printing plate
WO2018138340A1 (fr) * 2017-01-30 2018-08-02 Extreme Wheather Expertises - Exwexs Procédé et dispositif de météorologie et produit programme d'ordinateur associé
CN107655564A (zh) * 2017-05-11 2018-02-02 南京邮电大学 一种基于智能终端的多种技术融合的室内外环境检测方法
CN111368258A (zh) * 2020-03-04 2020-07-03 中国科学院东北地理与农业生态研究所 一种湿润地区日蒸散量的估算方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
H. C. WARDDENG: "A critical revision of the estimation of the latent heat flux from two-wavelength scintillometry", 《ROYAL METEOROLOGICAL SOCIETY》 *
张功等: "双波长交互法测算华北人工林平均水热通量的应用分析", 《中国农业气象》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113156542A (zh) * 2021-03-12 2021-07-23 中国气象局气象探测中心 一种基于闪烁法测量区域蒸散总水量的方法
CN113156542B (zh) * 2021-03-12 2022-05-24 中国气象局气象探测中心 一种基于闪烁法测量区域蒸散总水量的方法

Also Published As

Publication number Publication date
CN112287296B (zh) 2023-05-26

Similar Documents

Publication Publication Date Title
Hartogensis et al. Derivation of an effective height for scintillometers: La Poza experiment in Northwest Mexico
Meijninger et al. Determination of area-averaged water vapour fluxes with large aperture and radio wave scintillometers over a heterogeneous surface–Flevoland field experiment
Gultepe et al. Roundhouse (RND) mountain top research site: Measurements and uncertainties for winter alpine weather conditions
CN103293117A (zh) 一种微脉冲差分吸收激光雷达水汽时空分布反演方法
Doran et al. The Boardman regional flux experiment
Kohsiek et al. An extra large aperture scintillometer for long range applications
Yu et al. Supplement of the radiance-based method to validate satellite-derived land surface temperature products over heterogeneous land surfaces
CN104406715A (zh) 一种遥感估算地表感热/潜热通量的精度评价方法及系统
Kharaghani Propagation of refraction errors in trigonometric height traversing and geodetic levelling
Liu et al. Micrometeorological methods to determine evapotranspiration
Angevine et al. Entrainment results from the Flatland boundary layer experiments
CN111160680A (zh) 一种基于信息同化融合的农业干旱评估方法
Zheng et al. Comparison of sensible and latent heat fluxes from optical-microwave scintillometers and eddy covariance systems with respect to surface energy balance closure
Graf et al. Validation of a minimum microclimate disturbance chamber for net ecosystem flux measurements
CN112287296A (zh) 一种基于双波段闪烁仪的地表水热通量测算方法
JP4630092B2 (ja) 森林のｃｏ２フラックス測定法
Zhao et al. Using infrared thermal imaging technology to estimate the transpiration rate of citrus trees and evaluate plant water status
Wang et al. Modelling water and energy fluxes with an explicit representation of irrigation under mulch in a maize field
Odhiambo et al. Surface layer scintillometry for estimating the sensible heat flux component of the surface energy balance
JP5039195B2 (ja) 森林のｃｏ２フラックス測定法
Odhiambo et al. Sensible heat flux by surface layer scintillometry and eddy covariance over a mixed grassland community as affected by Bowen ratio and MOST formulations for unstable conditions
Butterworth et al. Characterizing energy balance closure over a heterogeneous ecosystem using multi-tower eddy covariance
De Bruin et al. Variance method to determine turbulent fluxes of momentum and sensible heat in the stable atmospheric surface layer
Evans Long-path scintillometry over complex terrain to determine areal-averaged sensible and latent heat fluxes
Krauss et al. Complex study of characteristics of a Cb cloud developing over the Arabian peninsula under conditions of high dew point deficit in the atmosphere. Part 2. Analysis of the meteosat data

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