CN112287296B - 一种基于双波段闪烁仪的地表水热通量测算方法 - Google Patents
一种基于双波段闪烁仪的地表水热通量测算方法 Download PDFInfo
- Publication number
- CN112287296B CN112287296B CN202011088336.XA CN202011088336A CN112287296B CN 112287296 B CN112287296 B CN 112287296B CN 202011088336 A CN202011088336 A CN 202011088336A CN 112287296 B CN112287296 B CN 112287296B
- Authority
- CN
- China
- Prior art keywords
- heat flux
- humidity
- calculating
- mws
- las
- 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
- 230000004907 flux Effects 0.000 title claims abstract description 97
- 238000000034 method Methods 0.000 title claims abstract description 75
- 239000002352 surface water Substances 0.000 title claims abstract description 40
- 238000004364 calculation method Methods 0.000 claims abstract description 24
- 230000008569 process Effects 0.000 claims abstract description 15
- 238000007781 pre-processing Methods 0.000 claims abstract description 11
- 238000012216 screening Methods 0.000 claims abstract description 8
- 230000008571 general function Effects 0.000 claims abstract description 5
- 230000006870 function Effects 0.000 claims description 25
- 238000001914 filtration Methods 0.000 claims description 12
- 230000003287 optical effect Effects 0.000 claims description 8
- 238000006073 displacement reaction Methods 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 6
- 238000004379 similarity theory Methods 0.000 claims description 6
- 230000003746 surface roughness Effects 0.000 claims description 6
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 6
- 230000005855 radiation Effects 0.000 claims description 4
- 230000001133 acceleration Effects 0.000 claims description 3
- 230000002452 interceptive effect Effects 0.000 claims description 3
- 238000001556 precipitation Methods 0.000 claims description 3
- 229920006395 saturated elastomer Polymers 0.000 claims description 3
- 230000014509 gene expression Effects 0.000 claims description 2
- 230000009977 dual effect Effects 0.000 description 6
- 239000002689 soil Substances 0.000 description 4
- 241000722941 Achillea Species 0.000 description 3
- 241000196324 Embryophyta Species 0.000 description 2
- 206010057040 Temperature intolerance Diseases 0.000 description 2
- JJWKPURADFRFRB-UHFFFAOYSA-N carbonyl sulfide Chemical compound O=C=S JJWKPURADFRFRB-UHFFFAOYSA-N 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000001704 evaporation Methods 0.000 description 2
- 230000008020 evaporation Effects 0.000 description 2
- 238000009499 grossing Methods 0.000 description 2
- 230000008543 heat sensitivity Effects 0.000 description 2
- 230000006698 induction Effects 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 208000003643 Callosities Diseases 0.000 description 1
- 244000025254 Cannabis sativa Species 0.000 description 1
- 206010020649 Hyperkeratosis Diseases 0.000 description 1
- 240000008042 Zea mays Species 0.000 description 1
- 235000005824 Zea mays ssp. parviglumis Nutrition 0.000 description 1
- 235000002017 Zea mays subsp mays Nutrition 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- SGRYPYWGNKJSDL-UHFFFAOYSA-N amlexanox Chemical compound NC1=C(C(O)=O)C=C2C(=O)C3=CC(C(C)C)=CC=C3OC2=N1 SGRYPYWGNKJSDL-UHFFFAOYSA-N 0.000 description 1
- 229960003731 amlexanox Drugs 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 235000005822 corn Nutrition 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 230000002262 irrigation Effects 0.000 description 1
- 238000003973 irrigation Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 230000001681 protective effect Effects 0.000 description 1
- 230000000087 stabilizing effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
- 238000009834 vaporization Methods 0.000 description 1
- 230000008016 vaporization Effects 0.000 description 1
- 238000005303 weighing Methods 0.000 description 1
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)计算:
作为本发明的一个优选实施例,通过式(4)~(6)计算空气折射指数的结构参数:
作为本发明的一个优选实施例,所述数据筛选,剔除降水时刻的数据、剔除饱和数据、剔除信号强度较弱的数据、剔除弱湍流数据。
式(10)~(11)中,AT,LAS、Aq,LAS分别表示LAS的温度和湿度折射指数的系数,AT,MWS、Aq,MWS分别表示MWS的温度和湿度折射指数的系数。
作为本发明的一个优选实施例,通过双波段交互BC方法计算气象结构参数,实时计算温度和湿度折射指数的系数RTq,通过求解空气折射指数结构参数和气象结构参数的关系式:
求解公式(12)~(14)得到:
式(15)~(17)中,AT,LAS、Aq,LAS分别表示LAS的温度和湿度折射指数的系数,AT,MWS、Aq,MWS分别表示MWS的温度和湿度折射指数的系数;
其中,LAS与MWS用同一公式计算温度和湿度折射指数的系数AT、Aq,计算公式如下:
式(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函数作为大气稳定度普适函数。
作为本发明的一个优选实施例,根据莫宁-奥布霍夫相似理论,温度结构参数、湿度结构参数与特征温度、特征湿度有如下关系:
式(20)~(21)中,z是有效高度,单位为m;d0是零平面位移高度,单位为m;T*是特征温度,单位为K;q*是特征湿度,为百分数;fT、fq为普适函数;LMO为奥布霍夫长度。
作为本发明的一个优选实施例,奥布霍夫长度和摩擦风速计算公式如下:
式(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:
式(26)和(27)中,为平均空气密度,单位为kg/m3;为平均比湿,单位为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
表2
图3为以阿柔站为例的地表水热通量测算方法详细流程图。
如图1和图3所示,所述测算过程如下:
步骤S1,基于双波段闪烁仪获取高频数据。
本实施例中,所述双波段闪烁仪采用OMS(型号:BLS900&RPG-MWSC-160) 双波段闪烁仪,来获取站点观测区域的高频数据。
步骤S2,将所述高频数据预处理后,计算光强方差。
本步骤中,所述高频数据的预处理,包括对高频数据去野点、和\或高频滤波处理、和\或低频滤波处理。
优选地,本步骤中采用数据平均绝对偏差(Mean Absolute Deviation,MAD) 方法去除野点值,具体的,所述MAD方法,遍历所有数据点,当一个窗口一个点值超过绝对偏差中位数的10倍,该点被识别为一个野点值,并被剔除,直至没有数据点被识别为野点值,这个过程才会停止,完成去野点环节。
本步骤中,完成高频数据的预处理后,通过式(1)~(3)计算光强方差
式(1)-(3)中,分别表示LAS、MWS、OMS的光强方差;ILAS、IMWS分别表示LAS、MWS的光强信号;分别表示LAS、 MWS闪烁信号在平滑窗口的平均值,上划线表示均值。优选地,LAS信号的平滑窗口为12s,MWS信号的平滑窗口为67s。上述公式计算的过程,即高频滤波处理后,求取的光强方差再通过低频滤波因子进行处理。
上述光强方差的计算过程,通过高频滤波和低频滤波移除了噪声和水汽吸收的影响。其中,所述高频滤波通过公式(1)~(3)实现,低频滤波通过 Matlab2018a中smooth库实现。
步骤S3,根据所述光强方差,计算空气折射指数的结构参数。
本步骤中,所述计算空气折射指数的结构参数,首先获取站点的A系列参数,包括AT,LAS、Aq,LAS、AT,MWS、Aq,MWS,根据A系列参数计算空气折射指数的结构参数。
其中,Gint函数在不同的站点是不同的,其主要的影响因子是近红外闪烁仪和MWS的波长、孔径直径、近红外闪烁仪与微波闪烁仪的间距和光径路线长度。Gint函数值计算公式如下:
式(7)-(16)计中,DtLAS、DrLAS分别表示LAS发射端和接收端的孔径直径;DtMWS、DrMWS分别表示MWS发射端和接收端的孔径直径;d(x)表示大孔径闪烁仪与微波闪烁仪的间距;N为光径路线长度;k为空间波数;kLAS、kMWS、 kOMS分别表示LAS、MWS、OMS空间波数(k=2π/λ);Jo、J1都表示贝塞尔函数,下标是阶数;x为在光径路线上的位置。
步骤S4,根据所述空气折射指数的结构参数,计算气象结构参数。
进行数据筛选,剔除降水时刻的数据、饱和数据、信号强度较弱的数据和弱湍流数据。
本步骤中,所述进行大气稳定度判断,选取温度和湿度折射指数的系数Rtq和净辐射Rn共同判断大气稳定状况,可以更好地反映大气稳定度。
进一步地,在TW方法中,空气折射指数结构参数和气象结构参数的公式如下所示:
进一步地,在BC方法中,空气折射指数结构参数和气象结构参数的公式如下所示:
根据公式(22)~(24)求解得到:
式(25)~(27)中,AT,LAS、Aq,LAS分别表示LAS的温度和湿度折射指数的系数,AT,MWS、Aq,MWS分别表示MWS的温度和湿度折射指数的系数,A系列参数与空气温度、湿度相关,需要实时计算,LAS与MWS用同一公式计算温度和湿度折射指数的系数,计算公式如下:
式(28)~(29)中,P表示大气压(Pa);T表示空气温度(K);q为空气比湿;R表示湿空气的空气常数;Rd、RV分别表示干空气和水蒸汽的气体常数; bt1、bt2、bq2的单位都是K Pa-1,与光束的波长有关;上划线表示均值。
步骤S5,根据所述气象结构,计算地表水热通量。
本步骤中,所述地表水热通量包括感热通量和潜热通量。优选地,根据莫宁-奥布霍夫相似理论,进行大气稳定度判断,选取合适的大气稳定度普适函数,通过迭代计算摩擦风速u*、特征温度T*、特征湿度q*,进而计算感热通量和潜热通量。
优选地,本步骤中,选取温度和湿度折射指数的系数Rtq和净辐射Rn共同判断大气稳定状况,可以更好地反映大气稳定度。选取Li函数作为稳定度函数进行计算,不稳定条件为:
稳定条件为:
其中,根据莫宁-奥布霍夫相似理论,温度结构参数、湿度结构参数与特征温度、特征湿度有如下关系:
式(29)~(31)中,z是有效高度(m);d0是零平面位移高度(m);T*温度特征量(K);q*温度特征量(%);fT、fq为普适函数;LMO为奥布霍夫长度。
奥布霍夫长度(LMO)和摩擦风速(u*)计算公式如下:
式(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),公式如下:
利用此发明处理阿柔站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 (10)
1.一种基于双波段闪烁仪的地表水热通量测算方法,其特征在于,所述方法包括如下步骤:
步骤S1,基于双波段闪烁仪获取高频数据;
步骤S2,将所述高频数据预处理后,计算光强方差;
步骤S3,根据所述光强方差,计算空气折射指数的结构参数;
步骤S4,根据所述空气折射指数的结构参数,计算气象结构参数;
步骤S5,根据所述气象结构参数,计算地表水热通量;
其中,所述双波段闪烁仪高频数据的预处理,包括对高频数据去野点、和\或高频滤波处理、和\或低频滤波处理;
所述高频数据去野点,采用数据平均绝对偏差方法去除野点值,遍历所有数据点,当一个窗口一个点值超过绝对偏差中位数的10倍,该点被识别为一个野点值,并被剔除,直至没有数据点被识别为野点值,过程停止,完成去野点环节;
所述步骤S2中计算光强方差,具体通过式(1)~(3)计算:
式(1)~(3)中,分别表示大孔径闪烁仪LAS、微波闪烁仪MWS、双波段闪烁仪OMS的光强方差;ILAS、IMWS分别表示LAS、MWS的光强信号;分别表示LAS、MWS闪烁信号在平滑窗口的平均值,上划线表示均值;
所述的地表水热通量测算方法,通过式(4)~(6)计算空气折射指数的结构参数:
其中,Gint函数在不同的站点是不同的,其主要的影响因子是近红外闪烁仪和MWS的波长、孔径直径、近红外闪烁仪与微波闪烁仪的间距和光径路线长度,Gint函数值计算公式如下:
上述式计中,DtLAS、DrLAS分别表示LAS发射端和接收端的孔径直径;DtMWS、DrMWS分别表示MWS发射端和接收端的孔径直径;d(x)表示大孔径闪烁仪与微波闪烁仪的间距;N为光径路线长度;k为空间波数;kLAS、kMWS、kOMS分别表示LAS、MWS、OMS空间波数,k=2π/λ;Jo、J1都表示贝塞尔函数,下标是阶数;x为在光径路线上的位置。
3.根据权利要求2所述的地表水热通量测算方法,其特征在于,进行数据筛选:剔除降水时刻的数据、剔除饱和数据、剔除信号强度较弱的数据、剔除弱湍流数据。
5.根据权利要求4所述的地表水热通量测算方法,其特征在于,通过双波段交互BC方法计算气象结构参数,实时计算温度和湿度折射指数的系数RTq,通过求解空气折射指数结构参数和气象结构参数的关系式(12)~(14):
得到计算气象结构参数、温度和湿度折射指数的系数RTq:
式(15)~(17)中,AT,LAs、Aq,LAS分别表示LAS的温度和湿度折射指数的系数,AT,MWS、Aq,MWS分别表示MWS的温度和湿度折射指数的系数;
其中,LAS与MWS用同一公式计算温度和湿度折射指数的系数AT、Aq,计算公式如下:
式(18)~(19)中,P表示大气压,单位为Pa;T表示空气温度,单位为K;q为空气比湿,单位为g/g;R表示湿空气的空气常数,单位为J/(mol·K);Rd、RV分别表示干空气和水蒸汽的气体常数,单位为J/(mol·K);bt1、bt2、bq2的单位都是K Pa-1,与光束的波长有关;上划线表示均值。
6.根据权利要求5所述的地表水热通量测算方法,其特征在于,所述步骤S5计算地表水热通量,根据莫宁-奥布霍夫相似理论,进行大气稳定度判断,选取合适的大气稳定度普适函数,通过迭代计算摩擦风速、特征温度、特征湿度,进而计算感热通量和潜热通量。
7.根据权利要求6所述的地表水热通量测算方法,其特征在于,选取湿度折射指数的系数RTq和净辐射Rn共同判断大气稳定状况,选取Li函数作为大气稳定度普适函数。
9.根据权利要求8所述的地表水热通量测算方法,其特征在于,奥布霍夫长度和摩擦风速计算公式如下:
式(22)和(23)中,g为重力加速度,单位为m/s2;u*为摩擦风速,单位为m/s;κ为vonKarman常数;Z0为地表粗糙度,单位为m;Zu为风速的测量高度,单位为m;T*是特征温度,T表示空气温度,单位为K,上划线为均值;
除了Gint函数、温度和湿度结构参数的系数、有效高度等参数局地化计算,还包括地表粗糙度、零平面位移高度的局地化计算,计算公式如下:
d0≈2/3hveg……………………………………(24)Z0≈0.13hveg……………………………………(25)
式(24)和(25)中,hveg为植被株高,单位为m。
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 CN112287296A (zh) | 2021-01-29 |
CN112287296B true 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) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113156542B (zh) * | 2021-03-12 | 2022-05-24 | 中国气象局气象探测中心 | 一种基于闪烁法测量区域蒸散总水量的方法 |
Citations (6)
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 | 北京师范大学 | 一种大尺度水热通量观测系统 |
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 | 中国科学院东北地理与农业生态研究所 | 一种湿润地区日蒸散量的估算方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
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 |
-
2020
- 2020-10-13 CN CN202011088336.XA patent/CN112287296B/zh active Active
Patent Citations (6)
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 | 北京师范大学 | 一种大尺度水热通量观测系统 |
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)
Title |
---|
A critical revision of the estimation of the latent heat flux from two-wavelength scintillometry;H. C. Warddeng;《Royal Meteorological Society》;1912-1922 * |
双波长交互法测算华北人工林平均水热通量的应用分析;张功等;《中国农业气象》;第第39卷卷(第第6期期);380-389 * |
Also Published As
Publication number | Publication date |
---|---|
CN112287296A (zh) | 2021-01-29 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Hartogensis et al. | Derivation of an effective height for scintillometers: La Poza experiment in Northwest Mexico | |
Nichols et al. | Evaluation of the evaporative fraction for parameterization of the surface energy balance | |
Weaver | Temperature and humidity flux-variance relations determined by one-dimensional eddy correlation | |
Allen et al. | Evapotranspiration information reporting: I. Factors governing measurement accuracy | |
Kustas et al. | Estimation of the soil heat flux/net radiation ratio from spectral data | |
Gavilan et al. | Measuring versus estimating net radiation and soil heat flux: Impact on Penman–Monteith reference ET estimates in semiarid regions | |
CN102455282B (zh) | 测量土壤含水量的方法 | |
Huang et al. | Modeling evapotranspiration for cucumber plants based on the Shuttleworth-Wallace model in a Venlo-type greenhouse | |
Marek et al. | Post-processing techniques for reducing errors in weighing lysimeter evapotranspiration (ET) datasets | |
Biscoe et al. | Barley and its environment. I. Theory and practice | |
Gultepe et al. | Roundhouse (RND) mountain top research site: Measurements and uncertainties for winter alpine weather conditions | |
Lakshmi et al. | The influence of the land surface on hydrometeorology and ecology: new advances from modeling and satellite remote sensing | |
Doran et al. | The Boardman regional flux experiment | |
Baker et al. | Influence of Stand Geometry on Light Interception and Net Photosynthesis in Cotton 1 | |
CN111079256A (zh) | 一种基于遥感的灌溉水有效利用系数测算方法 | |
CN112287296B (zh) | 一种基于双波段闪烁仪的地表水热通量测算方法 | |
CN111160680A (zh) | 一种基于信息同化融合的农业干旱评估方法 | |
Liu et al. | Micrometeorological methods to determine evapotranspiration | |
Zhang et al. | Diurnal and seasonal variations of surface albedo in a spring wheat field of arid lands of Northwestern China | |
Nyolei et al. | Evapotranspiration simulation from a sparsely vegetated agricultural field in a semi-arid agro-ecosystem using Penman-Monteith models | |
Graf et al. | Validation of a minimum microclimate disturbance chamber for net ecosystem flux measurements | |
Zhao et al. | Using infrared thermal imaging technology to estimate the transpiration rate of citrus trees and evaluate plant water status | |
Helfter et al. | A noninvasive optical system for the measurement of xylem and phloem sap flow in woody plants of small stem size | |
Cain et al. | Spatially averaged sensible heat fluxes measured over barley | |
CN110929653A (zh) | 一种基于遥感的灌溉水有效利用系数测算方法 |
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 |