综述

面积平均通量与光闪烁方法

  • 王介民
展开
  • 中国科学院西北生态环境资源研究院, 甘肃 兰州 730000

王介民 (1937 -), 男, 山西万荣人, 研究员, 主要从事大气物理、 大气遥感、 陆面过程研究. E-mail:

收稿日期: 2021-05-10

  修回日期: 2021-09-02

  网络出版日期: 2021-12-28

基金资助

高分辨率对地观测系统国家重大专项(21-Y20B01-9001-19/22)

国家自然科学基金项目(42101411)

Area Averaged Fluxes and Scintillometry

  • Jiemin WANG
Expand
  • Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou 730000, Gansu, China

Received date: 2021-05-10

  Revised date: 2021-09-02

  Online published: 2021-12-28

本文亮点

本文是对近20年新发展的双波段光闪烁方法的一个综述。陆面过程模式常常是基于局地或斑块尺度上的观测发展的, 其与大气模式较大网格尺度的不匹配, 显然会对后者的效能带来很大影响。如何扩展野外实验站点的代表性尺度已成为当前更好地了解陆面过程, 进而改善陆面过程模式与大气模式网格尺度匹配问题的关键。可用的面积平均通量观测方法, 包括以涡动相关方法为主的多点微气象观测、 飞机观测、 卫星和地面遥感等5种。其中, “光闪烁方法”是当前最为可行的、 可以大到10 km尺度的感热通量和潜热通量观测方法, 特别是, 它可以应用于复杂下垫面包括山谷地区和城市等。光闪烁方法的理论涉及电磁波传输和大气湍流。文章从折射指数、 结构参数、 湍流谱等基本概念开始, 对由对数光强方差计算折射指数的结构参数等基本公式的推导, 由光程权重函数、 空间谱权重函数、 光强的时间序列谱等对闪烁仪主要工作尺度的了解, 由折射指数结构参数计算温度、 湿度的结构参数的方法, 以及利用近地层相似理论计算感热通量和潜热通量等光闪烁方法的理论、 公式和计算步骤等做了较系统的阐述。进而, 在介绍双波段闪烁仪通量足迹函数之后, 对光闪烁方法与涡动相关方法从特征、 优势与缺点三方面做了比较; 指出结合使用涡动相关和光闪烁两种方法, 可以更好地进行面积平均通量分析, 进而用于模式的发展和检验, 以及更好的流域尺度的能量和水循环研究。最后, 对光闪烁方法的应用, 从较均匀下垫面、 复杂下垫面、 城市地区、 遥感模式的地面“真值”及在大气模式中的应用等几方面做了介绍; 特别是, 结合黑河流域阿柔和大满两站2020年部分资料的分析, 彰显了双波段闪烁仪对较大尺度感热、 潜热通量观测的明显优势。但是, 有关方法特别是观测水汽通量的微波闪烁仪的研制究竟为时尚短; 相关硬件、 软件、 资料处理方法等, 许多地方都还需要研究改进。相对于涡动相关方法, 光闪烁方法的理论和数据处理都更为复杂; 有关台站和资料分析人员, 需要更好的物理和微气象学基础。国内闪烁仪在城市地区的应用, 至今仍开展较少; 青藏高原和内地一些湖面蒸发研究中的难点, 也有望在应用光闪烁方法的过程中有所突破。基于光闪烁法等多种面积平均通量观测, 对陆面过程研究做尺度扩展, 并借以推动中—大尺度大气模式的发展, 更是我们殷切期盼的。

本文引用格式

王介民 . 面积平均通量与光闪烁方法[J]. 高原气象, 2021 , 40(6) : 1377 -1393 . DOI: 10.7522/j.issn.1000-0534.2021.zk017

Highlights

This is a review of the Optical-Microwave Scintillometer (OMS) system newly developed in last two decades, which can measure area averaged sensible and latent heat fluxes over a scale of 1 -10 km, especially over heterogeneous surfaces such as cross a valley or over urban areas.Among the methods of area averaged flux measurements, such as eddy-covariance based multi-point observations, air craft observations, satellite, and surface remote sensing etc., scintillometry is probably the most feasible technique in getting areal fluxes up to 10 kilometers.The basic theory of scintillometry includes electromagnetic wave propagation, atmosphere turbulence, and micrometeorology, which are more sophisticated than that of the popular eddy covariance (EC) system.Based on the introduction of concepts such as refractive index, structure parameter and turbulence spectra etc., basic scintillometer theories and equations are presented briefly, including: (1) calculation of structure parameter of refractive index via the variance of received light log-intensities; (2) the understanding of the working scale of scintillometry via the light-path weighting function, the spatial spectral weighting function, and the temporal spectral characteristics; (3) the derivation of structure parameters of temperature and humidity via the structure parameters of refractivity; (4) the calculation of fluxes by using the typical functions of Monin-Obukhov similarity; (5) the footprint analysis of scintillometry.A comparison between scintillometry and EC are presented in three aspects: ‘Characteristics’, ‘Advantages’ and ‘Weakness’.It is clear that a combined use of EC & Scintillometry can provide better area averaged fluxes, and, refined flux aggregation schemes.Then, the application of scintillometry is introduced for rather homogeneous surfaces, complicated surfaces, urban areas, and, the ‘ground truth’ for remote sensing and the application of areal averaged fluxes in atmospheric models.The example utilizations of the OMS systems in the Arou alpine-meadow station and the Zhangye oasis station, of the Heihe River basin, clearly show the advantages of scintillometry over the EC in the measurements of larger scale evapotranspiration.Nowadays there are hundreds flux stations operating over various climate regions and surface states of the world.Upscaling of the spatial representativeness of these stations becomes a key point to better understanding the land surface processes, and improving the spatial matching between land surface models and meso/large scale atmospheric models.Comparatively, the time of development of scintillometry, particularly the microwave scintillometers in measuring water vapor fluxes, are still short.Further improvements of relevant hardware, data sampling, and data processing software etc.are still needed.

1 陆面过程模式与大气模式网格尺度的匹配问题

陆面过程模式(Land Surface Model, LSM), 从1969年Manabe将简单的“水桶模式”包含到全球大气环流模式(Atmospheric General Circulation Model, AGCM)起, 50年间, 经过加入植被物理过程的第二代模式, 加入植被光合作用的第三代模式等, 最新的模式已经包含土壤-植被-大气系统中几乎所有重要的物理、 化学和生物过程(如CLM5, Lawrence et al, 2019)。陆面过程模式的发展为当代数值天气预报模式(Numerical Weather Prediction, NWP)和区域/全球气候模式(RCM/GCM)性能的提高发挥了巨大作用(Sellers et al, 1997a)。
陆面过程模式的核心是依据地形、 植被、 土壤等参数, 在大气模式提供的辐射、 风、 温、 湿、 压、 降水等的驱动下, 为大气模式提供动量、 感热、 潜热和水汽、 CO2等通量及下垫面的温、 湿要素等。其发展和检验, 通常是基于局地或斑块尺度上的观测与分析进行的: 如土壤和植被参数一般是在植株到田块尺度(10-2~102 m)上得到的, 各种通量的观测, 其代表性尺度一般也只有数百米到1 km。而大气模式的网格尺度, 包括NWP和GCM等, 则为10~100 km或更大; 而且, 通常采取“大叶”概念, 即设定模式网格尺度上的陆面特性是均匀的, 并直接采用局地尺度上得到的参数或参数化方案运行模式(Sellers et al, 1997b)。
两者尺度的不匹配显然会对大气模式的效能带来很大影响。自然下垫面, 在模式网格尺度上, 通常是非均匀的, 包括地形起伏、 植被多样性、 土壤湿度变化等等。20世纪80年代后期开始, 探讨次网格尺度陆面非均匀性的参数化方法, 成为一个非常活跃的研究领域(Avissar and Pielke, 1989Giorgi and Avissar, 1997); 文献中多种方法纷至沓来, 如“有效参数方法”, “Mosaic/Tile(斑块)方法”, “统计-动力学方法”等(de Vrese et al, 2016)。对不同方法的合理性、 可行性以及最终模拟结果的分析, 是一个非常复杂的问题。由此带来的模式的边界条件、 输入参数和参数化、 模式的检验和评价等问题, 随着陆面过程模式和气候系统模式不断地复杂化, 愈来愈受到极大的关注(Dickinson et al, 2006)。
对复杂陆面过程的不同简化, 对上述网格尺度上陆面非均匀特性的不同处理, 形成了多个陆面过程模式。为了解各模式的优势与不足, 20世纪90年代起启动了几个国际性的陆面过程模式比较计划, 如PILPS(Henderson-Sellers et al, 1995, 参加的有26个陆面过程模式), ALMIP(Boone et al, 2009, 参加的有14个陆面过程模式)等。这些计划, 总的说是成功的, 如增进了对各个模式结构的了解, 包括其结果的不确定性范围和需要改进之处等。但较为遗憾的是, 模式比较的结果, 不同模式, 即便对同一地区同一组气象驱动数据的模拟, 也有明显的差别; 而且没有哪一个在所有方面都比其他更好(Lu et al, 2020)。其原因, 除不同模式对一些过程(如土壤水文、 生物物理过程等)的处理方法, 及模式结构、 参数、 算法等的不同外, 主要还源于次网格尺度非均匀性处理上的困难。
为了更好地了解陆面过程, 包括解决上述陆面过程模式与NWP/GCM模式网格尺度的匹配问题, 1980年代中起, 在“全球变化”(包括碳循环)研究的促动下, 国际上开展了一系列大型“陆面过程实验”(Sellers et al, 1997a), 如20世纪80年代的HAPEX-MOBILHY, FIFE, HEIFE到90年代的EFEDA, BOREAS, AMAZON, GAME-Tibet等(王介民, 1999), 以及本世纪的HiWATER(Li et al, 2013), TIPEX-III(Zhao et al, 2018), LITFASS(Beyrich et al, 2012)等。在全球不同区域进行的这些实验, 不仅采取各种最先进的地面观测手段(如涡动相关方法的普遍应用), 而且充分利用快速发展的各种卫星遥感及相关模式, 获取不同尺度的地表辐射和植被、 土壤含水量等陆面参数, 以及地气相互作用研究中最关心的动量、 能量、 水汽和几种温室气体通量。这些实验, 大多十分关注陆面参数的尺度扩展及其与大气模式网格尺度的耦合, 包括次网格尺度非均匀性的参数化方法, 模式网格尺度上面积平均通量的“聚合(integration)”方法, 以至较大尺度上面积平均通量的直接观测等(Sellers et al, 1997b; Mahrt et al, 2001)。

2 面积平均通量的观测方法

陆面过程模式与大气模式网格尺度的不匹配, 网格尺度上非均匀性参数化的诸多困难, 推动了基于观测获取面积平均通量的实验研究的发展。一些推算地表参数(包括各通量)的遥感模型, 特别是一些应用于大气模式的全球性的遥感产品(如MODIS产品, https: //modis.gsfc.nasa.gov/data/ dataprod/index.php#atmosphere), 其地表分辨率都在1~5 km或更大; 其检验也需要相应尺度的地面观测。
在地表能量收支各分量中, 潜热通量常常是仅次于净辐射的最重要组分。与之相关的地表蒸散量(Evapotranspiration, ET), 包括植被蒸腾和土壤、 水体以及植被冠层截留蒸发等, 不仅同时链接地表的能量收支和水分收支, 也在植物光合作用等生态过程中起重要作用。ET又涉及从叶片到植株、 田块、 景观、 区域, 以至全球等多个尺度。其模式网格或遥感像元尺度的面积平均值, 是许多模式和应用研究最关心的参数之一。欧洲多国的EVA-GRIPS等计划(Mengelkamp et al, 2006), 便是力图通过各种观测与模式分析, 研究非均匀下垫面大气模式网格及遥感像元尺度上(所谓中γ尺度, 2~20 km)的ET参数化问题。
当前, 面积平均通量的观测主要有如下几种方法:
(1) 基于多点微气象观测, 以涡动相关系统(EC)为主, 包含辐射、 风温湿梯度、 土壤温湿度以及植被参数等, 通过土地利用面积加权平均等聚合方法得到较大尺度的面积平均。
(2) 基于飞机, 特别是携带快速传感器的涡动相关系统的无人机或有人机, 观测边界层风、 温、 湿梯度及较大区域的感热、 潜热和CO2通量。
(3) 基于探空或地面遥感系统(风廓线雷达、 声雷达、 微波辐射仪、 激光雷达等)探测边界层风、 温、 湿廓线及混和层高度, 利用局地相似理论或边界层水热收支方程等计算感热潜热通量。
(4) 利用各种空间/频谱分辨率的卫星遥感影像, 结合有关模式推算不同尺度的地表土地利用、 叶面积指数、 地表反照率、 地表温度、 土壤含水量等参数, 以及感热、 潜热通量和ET值。
(5) 利用光闪烁方法即近20~30年新发展的大孔径闪烁仪和微波闪烁仪观测0.5~10 km尺度上的感热和潜热通量。
这5种方法, 都是当前可行的; 但各有其优势和不足。
多点微气象观测, 主要是基于多点涡动相关观测的方法, 已在LITFASS(Beyrich et al, 2012)和HiWATER(Li et al, 2013)等几个实验中应用。如HiWATER曾在张掖绿洲5.5×5.5 km2尺度上建了17个EC站进行各有关通量与辐射、 土壤、 植被等参数的观测。如此大规模的配置, 一般很难实现。站点的增加不仅带来设备、 经费、 人力等的困难, EC等先进系统的维护、 资料处理、 质量保证等更需要一个专业知识和技术水平都很高的工作集体才能完成。而且, 由多个点到面, 如果下垫面较复杂, 动力和热力的不均匀会形成“内边界层”甚至局地环流, 简单的面积加权平均等聚合方法显然会带来较大的不确定性。
飞机观测, 国内外都有应用。特别是无人机低空遥感, 有一些明显的优势; 但也有费用较高和空域管理等限制。携带涡动相关通量直测系统的飞机至今国外很少, 国内更是绝无仅有; 可行性依然不高。
利用GPS探空及地面的激光、 微波、 声波遥感等的边界层探测, 对了解边界层结构及其发展, 必不可少; 但用于计算水热通量, 理论和实践上都有较大的不确定性。
各种卫星遥感, 发展迅猛; 其应用将愈来愈普遍和深入。但当前涉及通量(特别是ET)计算的数十种模型, 包括应用广泛的SEBAL(Bastiaanssen et al, 1998), SEBS(Su, 2002)和Two-source model(双源模式)(Kustas and Norman, 1999)等, 其理论框架和模式参数等仍有许多可改进之处; 其推算的卫星像元尺度各通量, 也需要相应尺度的地面观测检验(这也是上述EVA-GRIPS等计划的目标); 而且, 作为模式依据的卫星过境“瞬时”像元, 时间尺度上与检验其结果的30 min EC通量的差别较大, 也带来一些不可忽视的问题。
相对于以上四种方法, “光闪烁方法”应当是当前最为可行的、 可以大到10 km尺度的面积平均通量观测方法。各种不同的光闪烁仪现已应用于全球多个地区特别是地形起伏等复杂下垫面的地气相互作用观测研究, 并已在一些台站通过约20 多年的业务应用考验(Beyrich et al, 2013)。各种有关杂志论文的发表近20 年来逐年呈指数式增长。数以百计的研究已表明, 光学波段的大孔径闪烁仪可以较准确地观测一些复杂下垫面上的面积平均感热通量; 小孔径激光闪烁仪可以同时估算田块尺度上的动量通量和感热通量; 红外-微波双波段闪烁仪系统可以在一些复杂下垫面同时观测感热通量和潜热通量。对十分复杂的城市地区, 全球数十个大中城市利用光闪烁仪的观测实验表明, 一些其他方法很难确定的下垫面特征参数, 由光闪烁仪的应用取得了很有意义的结果(Ward et al, 20152017)。
以下将对光闪烁方法的原理、 特点及应用等做较详细的介绍。

3 光闪烁方法(Scintillometry

3.1 概述

光闪烁现象, 若干世纪之前人们就注意到了; 一个众所周知的例子就是地面看到的星光闪烁。早期的有关研究便是星光闪烁造成的图像模糊对天文观测的影响; 近代的研究则更多地关注其对数字通讯, 卫星导航, 以及激光系统等的干扰。上世纪中后期人们更认识到, 对近地层大气的光闪烁分析, 还可应用于气象、 水文、 农业、 环境等科学领域。
一个恒定强度光源的光通过一定大气路径后, 接收端测到的光强脉动, 即所谓的“光闪烁”, 主要由空气折射指数的起伏引起; 后者则主要取决于光路上不同湍涡的温度和水汽浓度的变化。对此, 太阳光谱段与微波谱段有明显不同。在可见光和近红外谱段(如闪烁仪常用的波长0.6~1 μm), 折射指数主要与温度相关; 而在微波谱段(如常用的波长1~10 mm), 湿度脉动则起主要作用。靠近地面的动量通量和感热通量驱动大气湍流的发展; 加上地表的蒸发蒸腾, 近地层大气中不同湍涡的温度湿度起伏造成折射指数的脉动, 后者, 与大气湍流输送强度直接相关。从而, 通过光闪烁强度的观测, 就可以得到近地层的湍流通量, 包括动量通量、 感热通量、 潜热通量等(Wesely, 1976)。
随着微波闪烁仪近十多年来的迅速发展, 与普通大孔径闪烁仪并用的“红外—微波”双波长闪烁仪业已推出; 其观测结果, 在一些简单的下垫面上, 与涡动相关仪所测的感热通量和潜热通量有非常好的一致性。这对较复杂下垫面1~10 km尺度面积平均的水热交换研究及模式检验等提供了重要契机(Wang, 2018)。

3.2 特征尺度参数及主要类型

如上所述, 闪烁仪包括发射端与接收端。发射端发射一束波长为
λ
、 强度恒定的光(电磁波), 距离
L
以外的接收端接收到的则是经过湍流大气散射的强度起伏的光; 光强的起伏主要由光路上空气折射指数的变化引起。图1为最新发展的微波闪烁仪(MW, 波长1~10 mm)与大孔径闪烁仪(IR, 波长0.8~0.9 μm)并用的双波段系统示意图。左端发射, 右端接收, 两条光线交叉。对任一波段, 如接收到的光强为
I
, 按湍流大气中的光传播理论(Tatarskii, 1961), 在所谓的“非饱和”情况下, 对数光强的方差
σlnI2
应与该波段空气折射指数的结构参数
Cn2
成正比。还可计算两波段对数光强的协方差, 推导相应的折射指数交叉结构参数
Cn1n2
。由此, 即可推算温度及湿度的结构参数(
CT2
Cq2
), 进而计算感热、 潜热通量(
H
LvE
)(详参后)。作为此过程的重要条件, 是否“非饱和”, 还与闪烁仪的另一个尺度参数——孔径
D
相关。特别对近红外波段, 大气湍流越强, 光程越长, 越容易出现“饱和”。但加大孔径
D
, 就可以避免此问题; 是为“大孔径闪烁仪”的由来。
图1 双波段闪烁仪系统示意图

MW为微波段, IR为近红外波段。发射端的恒定光强通过湍流大气后, 接收端收到的是脉动光强, 由之可了解近地层大气的湍流强度并进而推算感热和水汽通量

但是接收端接收到的光, 除由发射端直线传输来的一部分外, 还可能通过其他许多路径到达。按经典光学的惠更斯-菲涅尔定律, 影响接收光强的主要路线是在所谓的第一菲涅尔区内。后者是一个椭球体, 其最大直径
F=λL
, 是又一个非常重要的尺度参数。实际上, 孔径
D
与第一菲涅尔区尺度
F
的相对大小决定了影响测量的主要涡旋尺度。
所以, 光闪烁仪的特征尺度参数主要是波长
λ
, 孔径
D
, 光程长度
L
和第一菲涅尔区尺度
F
这四个。据此, 可将现有光闪烁仪分为三种, 其特征参数及主要性能如下(表1):
表1 几种闪烁仪的基本尺度参数与可测参数

Table 1

Basic scale and measurable parameters of major types of scintillometer
类型
λ
D
L
F=λL
可测参数
SAS0.65 μm2.5 mm50~250 m≈1 cm
l0,Cn2ε, CT2u*,H

LAS

XLAS

≈0.9 μm≈15 cm≈30 cm0.5~5 km1~10 km

≈4 cm

≈10 cm

Cn,las2CT2(+u*)H
MWS(与LAS并用)1~10 mm30~40 cm1~10 km≈1~4 m
Cn,las2,Cn,mws2
CT2,Cq2+u*H,LvE
(1) 小孔径激光闪烁仪(SAS)。当前最常用的即所谓“位移光束激光闪烁仪(DBLS)”, 工作波长约0.65 μm。其主要影响的涡旋尺度在湍流谱的耗散区—惯性区一段(0.005~0.1 m), 可观测湍流的内尺度
l0
并由之推算湍流动能耗散率
ε
和折射指数结构参数
Cn2
; 由后二者可推算摩擦速度
u*
(或动量通量τ)及感热通量
H
。还可由双光束的相关性推算与光程垂直的横向水平风速
vc
(2) 大孔径闪烁仪(LAS)和超大孔径闪烁仪(XLAS)。一般, 其发射和接收两端的孔径
D>2F
; 工作波长在0.85~0.9 μm近红外波段上。近20多年来, LAS已得到广泛应用, 其孔径约为15 cm, 光程在5 km以内; XLAS的孔径约为30 cm, 光程可达10 km。两者的主要影响涡旋尺度都与其孔径大小相当。单个闪烁仪只能观测与其光程尺度相应的面积平均感热通量; 对双发射光束闪烁仪, 如Scintec的BLS900, 还能观测横向风速
vc
(3) 微波闪烁仪(MWS)。一般工作于毫米波段(波长1~10 mm)。由于其技术较为复杂, 2000年代中才研制成功, 近几年才有商业产品提供(德国, RPG-MWSC-160)。微波的传播受空气湿度影响较大, 且微波闪烁仪必须与近红外波段的大孔径闪烁仪同时使用, 组成所谓的“光学-微波”双波长闪烁仪系统(OMS), 才能得到一些必要的参数并进而推算感热通量和潜热通量(蒸发蒸腾量)。

3.3 基本理论与数据分析要点

相对于大气探测的一般方法, 特别是涡动相关方法, 光闪烁法涉及的理论明显较为复杂; 主要包括3个方面: (1)电磁波传播; (2)大气湍流; (3)微气象学。
限于篇幅, 本文只能就实际应用中特别是光闪烁资料处理中必需掌握的一些原理和基本方程, 作简要介绍。光闪烁法的基础, 是20世纪60~70年代建立的“电磁波在湍流大气中的传播理论”。有兴趣的读者可参看Tatarskii(1961)的开创性著作(有中译本), Monin and Yaglom(1975)的专著, Andreas(1990)编辑的文集, 以及近30年来发表的大量论文。
由于小孔径激光闪烁仪观测尺度较小, 大孔径闪烁仪只能测感热通量, 以下的介绍将着重于近10多年新发展的有巨大应用前景的“红外—微波双波段闪烁仪”的有关方面。 先介绍几个基本概念:

3.3.1 折射指数, 结构参数, 湍流谱

(1) 折射指数

前曾强调, 光闪烁主要由空气折射指数的起伏引起。空气的折射指数(n)除与光波长有关外, 主要取决于大气的温度、 湿度和气压等要素。这是一个非常接近于1 的数(如对红外光, 标准大气下, n≈1.00027)。如定义折射率(Refractivity)
nλ=(n-1)106
, 按Owens(1967)的实验及Ward et al(2013)的写法:
nλ=m1+(m2-m1)RvqRpT
式中:
p
为气压(单位: Pa);
T
为温度(单位: K);
q
为比湿(单位: kg·kg-1);
R=Rd+q(Rv-Rd)
Rd
(=287.05)和
Rv
(=461.53)分别为干空气和水汽的气体常数;
m1
m2
为与波长有关的实验参数(具体表式参Ward et al, 2013)。
式(1), 对光学波段(
0.36 μm<λ<3 μm
), 湿度的影响很小。
nλ=nopt=0.7761+7.52×10-3λ2pT
对微波段, 湿度有明显影响。
λ>3 mm
下, 有较确切的实验结果:
nλ=nmw=0.776+3750T-0.056RvqRpT

(2) 结构函数与结构参数

大气湍流场包含各种不同尺度的湍涡及无规则的随机运动和变化。Kolmogorov(1941)在经典的“湍流级串”理论基础上进一步假设: 湍流的小尺度结构在统计上是均匀、 各向同性的, 且已与大尺度的结构无关; 在这个尺度上, 湍能没有生成也没有耗散, 湍流靠“惯性”运动并以量值为湍能耗散率的能量, 由较大湍涡向较小湍涡下传(直至小于某个尺度之后, 分子粘性起主要作用, 湍能耗散为热)。这一区域, 后来称为“惯性区”; 其上下边界分别定义为湍流的“外尺度”(
L0
)和“内尺度”(
l0
)。为重点研究惯性区较小尺度涡旋的统计特征, Kolmogorov创造性地引入“结构函数”概念: 对空间两点
r1
r+r1
的湍流场某要素, 如风速分量
u
, 其结构函数定义为,
Duu(r)ur+r1-ur12¯
即借助要素空间两点之差, 排除尺度大于
r=r
的湍涡的影响。
在惯性区内,
Duu
仅是湍能耗散率
ε
(量纲
m2·s-2
)和空间距离
r
(量纲m)的函数。由量纲分析易于导得,
Duur=C(εr)2/3Cu2r2/3,  (l0rL0)
式中: 常数
C
需由实验得到; 导出量
Cu2
称为风速
u
的“结构常数”或“结构参数”:
Cu2ur+r1-ur12¯r2/3,  (l0rL0)  
类似式(6), 可计算闪烁仪理论常用到的对于温度、 湿度和折射指数的结构参数
CT2
Cq2
Cn2
式(5), 即由Kolmogorov于1941年得到的
Duurr2/3
的关系, 是Kolmogorov湍流模型的基本规律之一; 后来称为“2/3定律”。
同理, 还有两个要素如温度和湿度的交叉结构函数
DTqr
和交叉结构参数
CTq
CTqDTqrr2/3=Tr1+r-Tr1qr1+r-qr1¯r2/3

(3) 湍流谱

这是一个更常用的概念。光闪烁法最常用的三维空间谱, 就是湍流动能或标量方差在不同尺度湍涡间的分布。常见的谱图, 即谱密度
Φ
与空间波数
κ
κ=2π/r
r
为湍涡尺度)的关系, 可由量纲分析得到并经实验验证; 不论对风速分量或
T
q
n
等标量, 都有如图2所示的形式。
图2 三维空间湍流谱示意图

横坐标为波数

κ
κ=2π/r
r
为湍涡尺度)。惯性区, 在双对数坐标图中为一段斜率为
-11/3
的直线, 其两端分别为湍流外尺度
L0
和内尺度
l0

我们仍然最关心“惯性区”的情况, 因为闪烁仪不论红外或微波, 都主要工作在这个尺度上(下详)。由图2可知, 惯性区的谱在双对数坐标图中为一段斜率为
-11/3
的直线。其左右边界, 即湍流的外尺度
L0
和内尺度
l0
的数值, 随观测高度及大气状态变化, 且因气象要素(风速、 温度或气压等)的不同而不同。多数情况下,
L0
随湍流强度和观测高度增加而增加, 为101~102 m;
l0
则与湍能耗散率及分子粘滞率有关, 为10-2~10-3m。
类似以上相关函数的“2/3律”, 三维谱与波数的“-11/3”关系也可由量纲分析得到(如对风速
u
, 有
ϕu(κ)κ-11/3
)。后者也是Kolmogorov湍流模型的基本规律之一, 称为“-11/3律”(注意: 对一维空间谱, 为“-5/3律”)。
结构函数与谱密度, 其实是从不同侧面对Kolmogorov湍流特性的描述。按维纳-辛钦定理, 相关函数与谱互为傅里叶变换; 如对风速分量
u
, 可由之导得如下关系[Tatarskii, 1961, (3.23)式]:
Duur=8π01-sin κrκrϕuκκ2dκ
Duur
做微分, 即可得到:
Φuκ=0.033Cu2κ-11/3
。类似, 可以得到光闪烁研究最关心的折射指数在惯性区段的三维空间谱:
Φnκ=0.033Cn2κ-11/3
与以上单要素的功率谱分析类似, 还可以得到二要素‘协谱’的表式。对本文着重分析的双波段闪烁仪, 红外与微波二波段的折射指数(如分别为
n1, n2
)的协谱有与式(9)类似的形式(Ludi et al, 2005):
Φn1n2κ=0.033Cn1n2κ-11/3
类似式(7)的定义,
Cn1n2
为二波段折射指数的交叉结构参数。

3.3.2 光闪烁基本公式: 由对数光强方差计算折射指数的结构参数

光(电磁波)在湍流大气中的传播研究, 特别是大气湍流引起光散射并造成接收端“光闪烁”的定量分析, 需要从Maxwell方程组出发。Tatarskii(1961)基于观测实验和Kolmogorov的大气湍流理论, 用平缓扰动法(不是直接对电场而是先对电场求对数再做扰动)求解Maxwell波动方程, 得到了接收端对数光强方差与传播距离及折射指数湍流强度之间的关系。按当前常用的闪烁仪参数, 针对点光源(小孔径闪烁仪), Tatarskii导得的光闪烁最基本公式为:
σlnI2=16π2kλ20Ldx0dκκΦnκsinκ2xL-x2kλL2
式中:
I
为接收光强;
σlnI2
为对数光强方差; 积分变量
x
为沿光程(
L
)距离,
κ
为波数;
λ
为闪烁仪工作波长,
kλ=2π/λ
为相应波数;
Φnκ
为折射指数三维谱。将式(9)代入式(11)并设
Cn2
为光路上平均值, 即得到由
σlnI2
计算
Cn2
的具体算式:
σlnI2=0.528π2kλ2Cn20Ldx0dκκ-8/3sinκ2xL-x2kλL2
积分即有,
σlnI2=0.496Cn2kλ7/6L11/6
这是Tatarskii得到的开创性的结果: 大气光闪烁强度和传播距离
L
成11/6次方增长的关系, 和折射指数湍流强度
Cn2
成正比增长的关系。这些在以后的实验中都被证明是正确的1
Tatarskii推导式(11)时的“平缓扰动法”, 隐含着大气状态为“弱湍流”的假设。当大气湍流较强时, 式(11)及据之导得的式(13), 即
Cn2
σlnI2
成正比的关系不再成立; 通常称之为出现“饱和”(而且光程
L
愈长, “饱和”愈易出现)。Wang et al(1978)提出加大闪烁仪孔径以避免“饱和”的设想, 并在积分方程(11)中加入了一个孔径平均因子, 进而积分得到了一个对式(13)的改进式:
σlnI2=0.893Cn2D-7/3L3
与Tatarskii的结果不同, 此处光闪烁强度和传播距离
L3
成正比; 即在一定的湍流强度下, 闪烁仪的光程明显可以更长些。这一结果, 推动了大孔径闪烁仪(LAS)的发明; 式(14)至今仍是LAS资料处理中最基本的公式之一。
加了孔径平均因子后的式(12)现在一般写为:
σlnI2=0.528π2kλ2Cn20Ldx0dκκ-8/3H2F2
式中:
H
函数即为原式中的
sin
函数(称波散射函数);
F
函数即为与孔径
D
有关的孔径平均效应函数。即:
H=sinκ2xL-x2kλL
F=2J10.5κDrx/L0.5κDrx/L2J10.5κDt1-x/L0.5κDt1-x/L
式中:
Dr
Dt
分别为接收端和发射端孔径, 一般
Dr=Dt=D
J1
为第一类一阶贝塞尔函数。
对微波闪烁仪, 式(15)~(17)仍然成立。但由于其波长长很多, 闪烁仪主要工作在第一菲涅尔区尺度上(参下), 式(15)的积分结果反而与小孔径闪烁仪的式(13)相似。微波闪烁仪其实是一种小孔径闪烁仪; 只是因孔径平均效应函数即不同站点的
λ
L
D
的影响, 式(13)中的常数0.496略有改变。如对黑河流域张掖大满站的MWS系统,
λ=1.86 mm
L=1854 m
D=0.30 m
式(15)的积分结果为:
σlnI22=0.429Cn22kλ27/6L11/6
式中:
σlnI2
Cn2
kλ
等参数加了下标“2”, 表明是针对微波闪烁仪的结果。
式(15)和(18)即可得到从观测所得对数光强方差计算折射指数结构参数的表式。

3.3.3 光程权重函数与空间谱权重函数——闪烁仪工作的主要尺度

上述闪烁仪基本方程式(15), 二重积分相互独立。先对波数
κ
或距离
x
积分, 得到如下两个式子:
σlnI2=0LPWFxdx
σlnI2=0SWFκdκ
式中:
PWFx
为沿路径各点的贡献;
SWFκ
为不同波数(或湍涡尺度)的贡献。
PWFx=0.528π2kλ2Cn20dκκ-8/3H2F2
SWFκ=0.528π2kλ2Cn20Ldxκ-8/3H2F2
归一化的
PWF
(令
xPWFx=1
)称为路径权重函数。其图像, 对红外波段的LAS和微波段的MWS如图3所示, 光程中部的权重最大, 到两端逐渐趋于0; MWS的比LAS的变化更平缓些。
图3 归一化的闪烁仪光程权重函数

以大满站配置为例:

L=1854 m
; 对LAS,
λ=880 nm
D=0.15 m
; 对MWS,
λ=1.86 mm
D=0.30 m

相应的,
SWF
称为空间谱权重函数, 其图像, 对LAS和MWS, 分别画在图4中。图4(b)为单对数坐标, 横坐标为湍涡尺度
r=2π/κ
, 由之更易于看到闪烁仪的主要工作尺度。以图所用大满站的配置为例, LAS的峰值波数
κ28.6
, 相应的湍涡尺度
r0.22 m
, 与LAS的孔径(0.15 m)相近; MWS的峰值波数
κ3.96
, 相应的湍涡尺度
r1.59 m
, 与MWS的第一菲涅尔区尺度(
F=λL=1.86 m
)相近。特别值得注意的是, 由图4(b)可见, 二者的尺度范围都很窄(SWF降为其0.1时, LAS的尺度为0.09~2 m, MWS的为0.85~11.6 m), 而且显然都在湍流谱的惯性区内。所以, 有的作者认为, 闪烁仪基本上是在单一尺度上工作的(van Kesteren, 2013)。
图4 闪烁仪的空间谱权重函数
SWF

图(b)横坐标为湍涡尺度

r=2π/κ
, 更易看清闪烁仪的工作尺度范围(仍以大满站配置为例, 参见图3)

3.3.4 对数光强的时间序列谱

上述
SWFκ
, 如式(20)所示, 可认为是对数光强的空间谱, 或“
κ
域谱”。与式(20)相应, 如果对对数光强的时间序列做傅里叶变换, 可得其在频域的分布:
σlnI2=0SlnIfdf
SlnIf
称为对数光强的时间序列谱(常常也简称时间谱),
f
为自然频率(Hz)。
基于Tatarskii点光源闪烁的基本公式, Clifford(1971)导得对数光强时序谱的理论表式。加上孔径平均因子((17)式)后, 类似(12)式的写法, 理论谱如下:
          SlnIf=0.528π2kλ2CnR20Ldx2πfvcdκκ-8/3×sin2κ2x(L-x)2kλLκvc2-2πf2F2
除类似空间谱
SWFκ
的分析可了解闪烁仪的湍流尺度响应特征外, 时间谱的分析主要用于:
(1) 检查由于水汽吸收(或电路漂移)等引起的低频端干扰, 据之对原始信号构建适当的数字高通滤波进行预处理。图5给出重庆虎头村站微波闪烁仪一个时次的对数光强时序谱例, 包括实际资料谱和理论谱(
vc
=1.35 m·s-1为观测值); 由水汽吸收引起的低频干扰明显可见, 特别在频率加权的半对数图(右图)中看得很清楚。对红外波段的大孔径闪烁仪, 亦有类似结果(但干扰的频率不同, 图略)。
图5 对数光强时序谱

重庆虎头村站微波闪烁仪的一个时次(2020年8月27日10:30,

vc=1.35 m·s-1
); 实际谱(红线)与理论谱(蓝线)的比较; 图中由水汽吸收引起的低频干扰明显可见

(2) 检查仪器硬件或安装缺点, 如塔体振动的影响、 交流电的干扰等。
(3) 由理论谱表式可见, 方程中多了一个参数: 横风(与光路垂直的)风速
vc
。这是一个对机场、 高速公路等安全影响较大的参数, 除3.2节提到其可由双光束闪烁仪的观测推算外, 还可如式(24)所示, 由时间谱的分析得到。

3.3.5 由折射指数结构参数计算温度、 湿度的结构参数

3.3.1节提到, 光闪烁主要由空气折射指数的起伏引起。而折射指数的起伏, 则主要由湍流大气的温、 湿度脉动引起; 气压脉动的影响相对很小。如令
n=n¯+n'
, 则考虑温度和湿度的相对变化, 应有:
n'=ATT'T¯+Aqq'q¯
式中: 系数
Ax=x¯n/x
(对
x=Tq
), 应为波长
λ
T
q
、 p的函数; 且如前述, 也可由实验确定(Ward et al, 2013)。
由上式, 以及定义
Cn2nr+r1-nr12¯/r2/3
, 折射指数的结构参数与温度、 湿度的结构参数应有如下关系(Andreas, 1988):
Cn2=AT2T¯2CT2+2ATAqT¯q¯CTq+Aq2q¯2Cq2
式中:
CTq
为温湿度交叉结构参数[式(7)]。
Cn2
分别由近红外和微波两个波段测量, 则仍不能据两个类似式(26)的式子求解三个未知数
CT2
Cq2
CTq
Hill et al(1988)提出一个“双波长方法(Two-wavelength Method)”, 额外假设温湿度交叉结构参数
CTq=±CT2Cq21/2
, 进而求解
CT2
Cq2
。此方法在处理OMS系统资料中虽已有不少的应用, 但有关假设, 意味着温湿度相关系数
rTq=CTqCT2Cq2=±1
即温湿度完全相关,
rTq
白天为
+1
, 夜间为
-1
。这常常与实际观测不符。特别在下垫面复杂情况下, 温湿度一般不是完全相关(如前述大满站的一段资料显示
rTq
为-0.4~0.8)。
有作者提出用
rTq
的实测值代入, 解决这一困难。但
rTq
的实测其实也不易。
为解决‘双波长方法’计算中的不确定性问题, Ludi et al(2005)提出了一个全新的“双波长相关法(Bichromatic correlation method)”。首先, 在OMS系统观测资料分析时, 不仅计算LAS和MWS各自的对数光强方差, 而且计算二者对数光强的协方差, 并通过下式计算二波段折射指数的交叉结构参数
Cn1n2
[参见式(10)]:
Coms=0.528π2kλ1kλ2Cn1n20Ldx0dκκ-8/3H1H2F1F2J0κd(x)
式中:
Coms=cov(lnI1,lnI2)
, 即为二波段对数光强的协方差; 下标1, 2分别表示LAS(红外波段)和MWS(微波波段); 函数
H
F
分别与式(16)和(17)相似(具体略去);
J0
为第一类0阶贝塞尔函数,
d(x)
为二光束距离
d
随光程
x
的变化(参图1), 后者可由发射和接收两端二传感器的距离计算。
这样, 类似式(26)就有如下三个式子以求解三个气象结构参数:
Cn12=AT12T¯2CT2+2AT1Aq1T¯q¯CTq+Aq12q¯2Cq2
Cn22=AT22T¯2CT2+2AT2Aq2T¯q¯CTq+Aq22q¯2Cq2
Cn1n2=AT1AT2T¯2CT2+2AT1Aq2+AT2Aq1T¯q¯CTq+Aq1Aq2q¯2Cq2
这是一个三元一次线性方程组, 可按标准方法求解。

3.3.6 利用近地层相似理论计算感热通量和潜热通量

通常, 大气近地层的感热通量
H=ρCpw'T'¯
, 水汽通量
E=ρw'q'¯
(其中
w'
为垂直风速脉动)。由之, 可引入近地层相似关系常用到的两个尺度参数:
θ*=-w'T'¯/u*
q*=-w'q'¯/u*
; 其中
u*=τ/ρ
τ
为动量通量(切应力)。
近地层中, 稳定度参数
ζ=zeff-d0Lob
式中:
zeff
为光程有效高度,
d0
为下垫面的零平面位移; 表征稳定度的奥布霍夫长度:
Lob=-T¯kgu*3w'T'¯+0.61Tw'q'¯
式中:
k=0.4
为卡曼常数;
g=9.8
为重力加速度。按莫宁-奥布霍夫相似理论, 无量纲化的结构参数
CT2
Cq2
应仅为稳定度参数
ζ
的函数, 即有如下相似关系:
CT2z2/3θ*2=fT(ζ)
Cq2z2/3q*2=fq(ζ)
其中亦有
z=zeff-d0
。早期的研究(如Wyngaard et al, 1971, 对Kansas实验的分析)已给出相似函数
fT(ζ)
fq(ζ)
的如下形式:
不稳定层结下(
ζ<0
):
fTζ=cT11-cT2ζ-2/3
fqζ=cq11-cq2ζ-2/3
稳定层结下(
ζ>0
):
fTζ=cT11+cT2ζ2/3
fqζ=cq11+cq2ζ2/3
各系数
cT
cq
需要由实验得到。表2列出两组系数: 一组是通常用于LAS的, 由Wyngaard et al(1971)基于Kansas实验得到, 后经Andreas(1988)修正卡曼常数后的结果; 一组是用于OMS的, 由Kooijmans and Hartogensis (2016)在总结十几个实验(包括双波段闪烁仪结合涡动相关通量观测)的基础上得到的结果, 应当是当前最可信的一组参数。
表2 由温度和湿度结构参数利用相似理论计算感热和水汽通量的相似函数参数

Table 2

<strong>Parameters of similarity function in calculating sensible and water vapor fluxes</strong>,<strong> based on temperature and humidity structure parameters</strong>,<strong> respectively</strong>
文献来源不稳定层结(
ζ<0
稳定层结(
ζ>0
cT1
cT2
cq1
cq2
cT1
cT2
cq1
cq2
Wyngaard et al (1971); Andreas (1988)4.96.14.92.2
Kooijmans & Hartogensis (2016)5.66.54.57.35.51.14.51.1
式(34)和(35)得到
θ*
q*
后, 即可由下二式分别计算感热通量
H
和水汽通量
E
或潜热通量
LvE
H=-ρCpu*θ*
LvE=-ρLvu*q*
摩擦速度
u*
需借助风速(
u
)及其观测高度(
zu
)和下垫面粗糙度(
z0
), 同样依靠相似关系, 在计算
θ*
q*
Lob
的同时, 用循环迭代法计算。

3.3.7 通量足迹函数(Footprint)分析

与涡动相关(EC)等其他通量观测方法一样, 对光闪烁仪如OMS的通量足迹(Footprint)函数分析是了解其通量空间代表性的必要步骤。通量观测的尺度扩展, 从点到面, 从较小尺度到更大的如大气模式网格尺度, 有关通量站点的Footprint分布是必不可少的计算依据。
某位置(
x,y,zeff
)观测到的通量
F(x,y,zeff)
, 应为其上风向地面源区各点[源强
F(x',y',0)
]的贡献累积, 可用下式计算:
Fx, y,zeff=-xFx',y',0ϕx-x',y-y',zeffdx'dy'
式中:
ϕx-x',y-y',zeff
即所谓Footprint函数, 与有效观测高度
zeff
及风速、 地面粗糙度、 大气稳定度等参数有关。
Footprint概念最早是20世纪60年代Pasquill在分析大气污染物的输送与扩散时提出的; 微气象上, 用于通量源区贡献分析的模式, 其实1990年前后才逐渐出现(Leclerc and Foken, 2014蔡旭晖, 2008)。各种通量足迹模式中, 最复杂的是求解大气动力学方程的“闭合模式”或“大涡模拟”; 其次是基于随机微分方程处理湍流场中粒子运动的“拉格朗日随机扩散模式”。两者理论上都较成熟, 后者还可以用于下垫面较复杂的一些情况; 但都需要大的计算量, 一般分析, 特别在处理长期观测数据时很难实施。当前普遍采用的是两个较便捷的方法:
(1) 基于一些简化假设求解平流-扩散方程, 利用所得解析解做足迹分析。常用的如Kormann-Meixner解析模式(Kormann and Meixner, 2001)。
(2) 基于拉格朗日随机扩散模式的大量模拟结果, 借助量纲分析, 得到一个较简单的参数化方法; 后者类似解析公式, 可利用实时观测资料做足迹分析。常用的如Kljun模式(Kljun et al, 2015)。
对闪烁仪, 其通量Footprint分析时, 将光程看做一条线状路径, 沿光程逐点利用单点足迹模型计算, 再将结果以光程权重函数(PWF)为权重叠加(有关结果可参图7)。
图6 阿柔高原草地站OMS系统与EC系统通量观测结果比较(2019年5月30日至6月14日)

Fig.6

An intercomparison of Sensible (H) and Latent (LE) heat fluxes observed by optical-microwave scintillometry (OMS) and eddy-covariance system (EC), from 30 May to 14 June 2019, over the Arou alpine meadow station with rather homogeneous surface states
图7 大满站两个系统OMS和EC的潜热通量比较(a, b)及二者的足迹分布(c, EC的通量源区明显比OMS的小很多)

2020年6月头几天(a)与最后几天(b), 由于植被覆盖变化巨大(净辐射也明显变大), 二者潜热通量的比较明显不同

3.4 光闪烁方法与涡动相关方法的比较

涡动相关(Eddy-Covariance, EC)方法已是当前观测地气间动量、 热量、 水分及二氧化碳等温室气体通量交换的主要方法, 广泛应用于各个陆面过程实验及几乎所有的通量台站(https: //fluxnet.fluxdata.org/about/)。涡动相关仪器主要由超声风速温度仪和气体分析仪组成, 通过观测风速与温度、 水汽、 二氧化碳浓度等的快速变化, 计算垂直风速脉动与其余各量的协方差而“直接”得到各相应的通量, 避免了其他通量观测方法的诸多假设条件, 因而具有最高的准确度, 并成为通量观测的“标准”方法。随着传感器、 数据采集与计算技术的发展, 有关仪器不仅可以直接用于野外环境, 不带来人为干扰, 而且可以长期连续工作, 具有令人称羡的稳定性和可靠性。
正如任何事物都有另一面一样, 涡动相关方法也有其局限性: (1)它需要捕获各种大小不同的涡旋对湍流通量的贡献。因而不仅要求足够快的采样频率(10~20 Hz), 而且要求足够长的取平均时间(30~60 min)。(2)它要求较好的大气平稳性(至少在取平均时间30~60 min内)和湍流发展较好的大气环境。(3)它要求所在的下垫面不是太复杂, 最好是平坦均匀的, 以减小平流输送的影响(或要求垂直风速的平均值
w¯0
)。(4)EC传感器及其支架会引起自然流场的畸变并影响风速及通量的观测精度, 如常用的CSAT超声仪要求对着主导风向安装(或需施加一定的“阴影修正”)。(5)更为重要的一点是, EC观测的代表性尺度较小。作为一种固定的“单点”观测, 对常见的3~6 m观测高度, EC的通量源区(或“足迹”)尺度一般只有数百米。这与常用的大气模式网格尺度有较大的差别。
本文主要介绍的光闪烁方法, 是一种“非直接”的通量观测方法, 需要借助相似理论由气象结构参数推算通量。半经验的相似函数, 带来通量的不确定性。但是, 它也克服了涡动相关方法的几个缺点: (1)光闪烁方法观测的是由电磁波的发射端到接收端数百米至10 km尺度上的面积(加权)平均通量, 其代表性尺度明显大于涡动相关方法。(2)其发射接收端支架引起的自然风场畸变对通量的影响可以忽略不计。(3)与涡动相关方法更为明显不同的是, 可见光、 近红外、 微波等几个波段的闪烁仪, 都工作在大气湍流谱的惯性区, 或惯性区到耗散区的一个尺度较小且范围很窄的涡旋尺度上。在闪烁仪较长的光程尺度上, 存在大量的如此对光闪烁敏感的涡旋, 故而在很短的取平均时间内(<1 min)即可得到统计上稳定的通量结果。众所周知, 当前许多遥感ET(Evapotranspiration)模式都是基于瞬间拍摄的遥感影像资料运行的; 闪烁仪的短时特征对遥感模式的结果验证尤其重要。(4)同样, 由于闪烁仪可以取很短的取平均时间, 故而也可以忽略类似EC那样对大气平稳性的要求。(5)由上述闪烁仪敏感涡旋尺度较小的特点, 带来一个主要优势, 即其可工作于较复杂的下垫面, 包括地形起伏甚至城市地区。近20年来, 已有数十个大中城市应用闪烁仪观测得到了大量极有价值的研究结果, 包括更可信的下垫面有效参数等。通量计算中, 闪烁仪对近地层相似理论在复杂下垫面的应用也没有出现比平坦均匀下垫面更突出的问题(Ward, 2017)。
表3 对两种方法的特征、 优点和局限分别做了简要总结。各通量站的观测研究表明, 结合使用涡动相关和光闪烁两种方法, 可以更好地进行面积平均通量分析, 进而用于模式的发展和检验, 以及更好的流域尺度的能量和水循环研究。
表3 通量观测的光闪烁方法与涡动相关方法比较

Table 3

<strong>A comparison of optical-microwave scintillometry </strong>(<strong>OMS</strong>)<strong> and eddy-covariance </strong>(<strong>EC</strong>)<strong> methods in flux observations</strong>
涡动相关方法(EC)光闪烁方法(OMS)
特征

观测风速、 温度、 标量浓度等要素的湍流脉动量, 由计算协方差直接得到各通量;

影响有关通量输送的湍涡包括10-2~102 m等极宽的尺度范围。

由观测光强脉动计算折射指数及温、 湿度结构参数, 借助近地层相似理论推算通量;

影响通量结果的涡旋主要在与F或D相当的尺度上, 即湍流惯性子区内一个波谱很窄的范围。

优点

通量精度高; 相关理论及数据处理方法较为简单;

当前公认的通量观测标准设备。

可很快(<1 min) 得到统计上稳定的通量值;

可应用于地形起伏或如城市等复杂下垫面;

地面代表性较大(1~50 km2)。

局限

要求大气状态平稳且湍流发展较好;

要求足够快的采样频率及足够长的取平均时间;

地面代表性较小(0.05~1 km2)。

通量计算依赖近地层相似关系, 故有较大的不确定性;

通量计算需要与光程尺度匹配的气象及地表参数, 其不确定性也影响通量计算结果;

硬软件技术及数据处理方法, 仍在发展中。

4 光闪烁方法的应用简介

光闪烁方法用于通量观测和陆面过程研究, 实际上始于1991年的EFEDA实验: 在西班牙半干旱区稀疏葡萄园上, 荷兰Wim Kohsiek研制的欧洲第一台大孔径闪烁仪, 与独立进行的涡动相关仪30 min通量观测惊人的一致性(相关性0.98, 散点分布斜率0.97; 参见文献de Bruin et al, 1995), 令实验参加者大开眼界。其代表性尺度大等优点, 促发了各类闪烁仪的研制、 生产和应用。
从EFEDA开始, 30年来, 全球不同地区已有数百个台站采用光闪烁仪观测通量。国内安装有大孔径闪烁仪的通量站, 据2013年统计已有约30个; 2015年以来, 国内新建的装有光学-微波闪烁仪的台站也已超过15个。本文篇幅有限; 有关应用的较详细介绍, 感兴趣的读者请参考文献de Bruin and Wang(2017)

4.1 较均匀下垫面的应用

地形相对平坦、 下垫面较为均匀的情况如草地、 农田、 荒漠等, 国外有新西兰、 欧洲几国、 日本、 墨西哥及美国新墨西哥州等数十个工作报导。国内的工作, 以北京师范大学刘绍民研究组为主, 从2004年的小汤山观测开始, 到其后几年的海河流域馆陶站、 黑河流域阿柔高原草地站等, 做了大量工作(卢俐等, 2009; Liu et al, 2011; Liu et al, 2018)。
2019年, 阿柔草地站安装了全新的双波段闪烁仪系统(OMS)。图6给出新系统运行初期一段通量观测资料与EC的比较。该站海拔3003 m; OMS光程长2390 m, 有效高度13.2 m; 涡动相关系统(EC, 包括CSAT3超声风速温度仪和LI-7500A气体分析仪)安装高度3.15 m, 大致在光程中部。图6所示的15天里, 有时有点小雨, 但基本上不影响通量特别是闪烁仪系统的观测(受影响的EC资料已剔除)。两套系统观测的通量, 包括感热和潜热, 总的来说一致性非常好。

4.2 复杂下垫面的应用

多数通量站的下垫面都是复杂的: 常见的如斑块状的各类农田, 农田与裸地、 林地、 草地等的混合, 地形起伏甚至山区, 各类森林等。有关报告很多, 可参Beyrich等关于德国LITFASS(农田、 森林、 水面等混合)的系列工作, Ezzahar et al(2007)关于墨西哥半干旱区橄榄果园的试验, Evans (2009) 在英国起伏地形农田的工作等。利用大孔径闪烁仪(LAS)估算水面蒸发量的工作也有一些(McJannet et al, 2013)。
国内2016年以来新建的OMS通量站, 包括林科院小浪底森林站(Zhang et al, 2021), 地化所普定陈旗站等, 都是较复杂的山区起伏地形; 西北院平凉站在陇东黄土高原白庙塬上, 西南大学虎头村站则在宽约600 m的山谷中。黑河地区的两个站, 除上述阿柔草地站外, 张掖大满站则建在包含农田、 果园、 防护林及村庄等的复杂下垫面上(图7)。
复杂下垫面上, 涡动相关系统的局地尺度的影响常常比较明显。从大满站两个系统OMS和EC的潜热通量比较及二者的足迹分布(图7)(感热通量一致性较好, 图略)可以看出, 两个时段, 2020年6月前几天[图7(a)], 当地主要作物玉米的苗高仅约20 cm, 土地基本上还是裸露的; 平均植被指数, EC小通量源区的比OMS大源区的约小30%。结果, EC所测潜热通量明显比OMS的小。6月后几天[图7(b)], 玉米株高已约1.5 m, 下垫面几乎被完全覆盖, OMS与EC所测潜热通量的差别就很小了(由于植被旺盛, 净辐射最大值也已从6月初的不到700 W·m-2变成高于800 W·m-2)。

4.3 城市地区的应用

城市地区, 除可能的地形变化外, 主要是街区建筑物、 路面、 广场、 绿地等等的混合, 可能是一种最复杂的下垫面。涡动相关方法在城市的应用, 只能得到一些很局地的信息, 对整个城市的陆面过程的了解几乎没有代表性可言。而光闪烁方法, 由于其可以大到10 km的空间跨度及内在的空间积分特性, 其观测结果的代表性明显要好很多; 近15 ~20年来, 在全球包括伦敦、 东京、 马塞尔、 罗兹、 马赛、 赫尔辛基等超过30个城市及城郊地区, 得到愈来愈广泛的应用。Ward (2017) 对闪烁仪在城市和复杂环境的应用有一个很好的综述, 可参。
国内城市对闪烁仪的应用至今很少, 当前看到的只有北京大学张宏升研究组对长江三角洲城市群的一个研究(Zhang et al, 2016)。

4.4 遥感模式的地面“真值”及在大气模式中的应用

以上提到, 光闪烁方法不仅在空间尺度上对一些遥感模型(如众多遥感ET模型)的结果, 特别是地面分辨率1 km或更大的遥感产品, 有较好的代表性; 其取平均时间可以短到1 min的特点, 也与普通遥感通常取卫星过境时刻影像做反演在时间尺度上匹配更好。近几十年来, 这方面的工作很多, 可参考Hemakumara et al(2003)Schuttemeyer et al (2007)Kleissl et al(2009)等文献。国内的相关研究至今仍相对较少。
大气模式研究中对光闪烁仪资料的应用, 包括陆面过程模式的检验, 中尺度模式的分析等。这方面仅有一些试验性的工作, 但也是本文较关注的, 可参考Beyrich and Mengelkamp (2006)Schuttemeyer et al (2008)Steeneveld et al (2011)Lee et al (2015) 等文献。

5 结语

基于湍流大气中电磁波的传播理论和微气象学建立的光闪烁方法, 近30年的发展特别是近10年来光学-微波双波段闪烁仪(OMS)的应用, 已经证明它是在1~10 km尺度、 特别是非均匀下垫面上观测感热、 潜热通量的有效方法, 具有当代最常用的涡动相关方法难以企及的多个优势。但是, 有关方法特别是观测水汽通量的微波闪烁仪的研制究竟为时尚短; 相关硬件、 软件、 资料处理方法等, 许多地方都还需要研究改进。如“时间序列谱”一节提到的水汽吸收的低频干扰问题等, 其处理方法如高通或带通滤波器的截断频率的选择, 至今没有公认的方法可循。国内近年新建的10多个OMS通量站, 也多因仪器长期工作的稳定性较差而受到困扰。
光闪烁方法牵涉的理论究竟较多, 资料的预处理如数据“掉包”、 “野点”等(可能与仪器的硬件及采集软件缺陷有关), 通量计算中如“双波长法”或“双波长相关法”的选择, 相似函数或参数的选择, 摩擦速度(
u*
)计算中代表性的风速和地表粗糙度的确定, 以及整体的数据质量判断等, 相对于涡动相关方法都更为复杂。有关操作人员和资料分析人员, 需要较好的物理基础和微气象学基础。
光闪烁方法在复杂下垫面特别是城市及市郊的应用, 国外已经很多, 也已由之得到了涡动相关方法不可能得到的城市陆面过程研究所需的多个有效参数。但国内的有关工作至今仍开展很少。青藏高原及内地多个项目涉及湖面蒸发观测, 一直有不少难点。光闪烁方法特别是双波段闪烁仪的发展也为此提供了契机。
国内基于光闪烁方法检验遥感模式包括遥感ET模式的工作至今也不多。用1~10 km尺度的闪烁仪观测对青藏高原和内地许多通量站的陆面过程研究做尺度扩展, 并借以推动国内中-大尺度大气模式的发展, 更是我们殷切期盼的。

王介民, 2021. 面积平均通量与光闪烁方法[J].高原气象, 40(6): 1377-1393.

WANG Jiemin, 2021. Area Averaged Fluxes and Scintillometry[J].Plateau Meteorology, 40(6): 1377-1393.

1 当Tatarskii的俄文原著在1959年出版时他还不到30岁。上式即原书中的(9.43)式, 仅将原对数光幅方差改为对数光强方差。

林科院小浪底、 中科院西北院平凉、 西南大学虎头村、 中科院地化所陈旗、 武汉大学屈家岭、 中科院西北院与北师大共建的阿柔和大满等各OMS站的站长和同事提供了相关资料。徐菲楠博士协助处理了图7有关资料及绘制footprint图。谨致衷心感谢!

AndreasE L1988.Atmospheric stability from scintillation measurements[J].Applied Optics, 27: 2241-2246.

AndreasE L1990.Selected papers on turbulence in a refractive medium[C]//SPIE Milestones Series, 25, 693pp.SPIE Optical Engineering Press, Bellingham, Wash., USA.

AvissarRPielkeR A1989.A parameterization of heterogeneous land surface for atmospheric numerical models and its impact on regional meteorology[J], Monthly Weather Review, 117: 2113-2136.

BastiaanssenWMenentiMFeddesRalet1998.A remote sensing surface energy balance algorithm for land (SEBAL).1.Formulation[J].Journal of Hydrology, 212: 198-212.

BeyrichFMengelkampH T2006.Evaporation over a heterogeneous land surface: EVA_GRIPS and the LITFASS-2003 experiment-an overview[J].Boundary-Layer Meteorology, 121: 5-32.

BeyrichFBangeJHatogensisOalet2012.Towards a validation of scintillometer measurements: The LITFASS-2009 experiment[J].Boundary-Layer Meteorology, 144: 83-112.

Beyrich F, Dereszynski P, van Kesteren B2013.Some aspects and results of scintillometer long-term operation.4th Workshop on Scintillometers and Applications[Z].Tübingen, Germany, 7-9 October, 2013.

BooneAde RosnayPBasalmoGalet2009.The AMMA Land Surface Model Intercomparison Project[J].Bulletin of the American Meteorological Society90(12): 1865-1880.

CliffordS F1971.Temporal-frequency spectra for a spherical wave propagating through atmospheric turbulence[J].Journal of the Optical Society of America61(10): 1285-1292.

de BruinHvan den HurkBKohsiekW1995.The Scintillation Method tested over a vineyard area[J].Boundary-Layer Meteorology, 76: 25-40.

de BruinHWangJ M2017.Scintillometry: A review[M/OL].[2021-02-15].

de VresePSchulzJHagemannS2016.On the representation of heterogeneity in land-surface-atmosphere coupling[J].Boundary-Layer Meteorology, 160: 157-183.

DickinsonR EOlesonK WBonanGalet2006.The Community Land Model and its climate statistics as a component of the Community Climate System Model[J], Journal of Climate, 19: 2302-2324.DOI: 10.1175/jcli3742.1.

EvansJ G2009.Long-path scintillometry over complex terrain to determine areal-averaged sensible and latent heat fluxes[D].Ph.D.thesis, The University of Reading, 176.

EzzaharJChehbouniAHoedjesJ C Balet2007.The use of the scintillation technique for monitoring seasonal water consumption of olive orchards in a semi-arid region[J].Agricultural Water Management89(3): 173-184.

GiorgiFAvissarR1997.Representation of heterogeneity effects in earth system modelling: Experience from land surface modelling[J].Reviews of Geophysics, 35: 413-437.

Henderson-SellersAPitmanA JLoveP Kalet1995.The Project for Intercomparison of Land Surface Parameterization Schemes (PILPS): Phase 2 and Phase 3[J].Bulletin of the American Meteorological Society74(4): 489-503.

HemakumaraH MChandrapalaLMoeneA F2003.Evapotranspiration fluxes over mixed vegetation areas measured from large aperture scintillometer[J].Agricultural Water Management, 58: 109-122.

HillR JBohlanderR ACliffordS Falet1988.Turbulence-induced millimeter-wave scintillation compared with micrometeorological measurements[J].IEEE Transactions on Geoscience and Remote Sensing26(3): 330-342.

KleisslJHongS HHendrickxJ M H2009.New Mexico scintillometer network supporting remote sensing and hydrologic and meteorological models[J].Bulletin of the American Meteorological Society90(2): 207-218.DOI: 10.1175/2008BAMS2480.1.

KljunNCalancaPRotachM Walet2015.The simple two-dimensional parameterization for Flux Footprint Prediction (FFP)[J].Geoscientific Model Development Discussions8(8): 6757-6808.DOI: 10.5194/gmdd-8-6757-2015.

KolmogorovA N1941.The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers[J].C.R.Acad.Sci.URSS, 30: 301-305.

KooijmansL M JHartogensisO K2016.Surface-layer similarity functions for dissipation rate and structure parameters of temperature and humidity based on eleven field experiments[J].Boundary-Layer Meteorology, 160: 501-527.

KormannRMeixnerF X2001.An analytical footprint model for non-neutral stratification[J].Boundary-Layer Meteorology, 99: 207-224.

KustasHNormanJ1999.Evaluation of soil and vegetation heat flux predictions using a simple two-source model with radiometric temperatures for partial canopy cover[J].Agricultural and Forest Meteorology94(1): 13-29.

LawrenceD MFisherR AKovenC Dalet2019.The Community Land Model Version 5: Description of new features, benchmarking, and impact of forcing uncertainty[J].J.Advances in Modeling earth Systems, 11, 4245-4287.

LeclercM YFokenT2014.Footprints in Micrometeorology and Ecology[M].Springer, Heidelberg, New York, Dordrecht, London, 239 pp.

LeeS HLeeJ HKimB Y2015.Estimation of turbulent sensible heat and momentum fluxes over a heterogeneous urban area using a large aperture scintillometer[J].Advances in Atmospheric Sciences, 32: 1-14.

LiXChenG DLiuS Malet2013.Heihe Watershed Allied Telemetry Experimental Research (HiWATER): Scientific objectives and experimental design[J].Bulletin of the American Meteorological Society94(8): 1145-1160.

LiuS M.Xu Z W, Wang W Z, et al, 2011.A comparison of eddy-covariance and large aperture scintillometer measurements with respect to the energy balance closure problem[J].Hydrology & Earth System Sciences15(4): 1291-1306.

LiuS MLiXXuZ Walet2018.The Heihe Integrated Observatory Network: A basin scale land surface processes observatory in China[J].Vadose Zone Journal, 17: 1-21.

LuHZhengD HYangKalet2020.Last-decade progress in understanding and modeling the land surface processes on the Tibetan Plateau[J].Hydrology and Earth System Sciences, 24: 5745-5758.

LudiABeyrichFMatzlerC2005, Determination of the turbulent temperature-humidity correlation from scintillometric measurements[J] Boundary-Layer Meteorology, 117: 525-550.

MahrtLVickersDSunJalet2001.Calculation of area-averaged fluxes: Application to BOREAS[J].Environmental Science, 40: 915-920.

McJannetDCookFMcGloinRalet2013.Long-term energy flux measurements over an irrigation water storage using scintillometry[J].Agricultural and Forest Meteorology, 168: 93-107.

MengelkampHBeyrichFHeinemannGalet2006.Evaporation over a heterogeneous land surface-The EVA-GRIPS project[J].Bulletin of the American Meteorological Society, 87: 775-786.

MoninAYaglomA1975.Statistical fluid mechanics: mechanics of turbulence.Volume 2[M].The MIT Press, Cambridge, Massachusetts, and London, England, 862.

OwensJ C1967.Optical refractive index of air: Dependence on pressure, temperature and composition[J].Applied Optics.6(1): 51-59.

SellersP JDickinsonR ERandallD Aalet1997a.Modeling the exchanges of energy, water, and carbon between continents and the atmosphere[J].Science275(5299): 502-509.

SellersP JHeiserhM D, Ha11 F G, et al, 1997b.The impact of using area-averaged land surface properties--topography, vegetation condition, soil wetness-in calculations of intermediate scale (approximately 10 km’) surface-atmosphere heat and moisture fluxes[J].Journal of Hydrology190: 269-361.

SchuttemeyerDSchillingsCMoeneA Falet2007.Satellite-based actual evapotranspiration over drying semiarid terrain in west Africa[J].Journal of Applied Meteorology & Climatology, 46: 97-111.

SchuttemeyerDMoeneA FHoltslagA A Malet2008.Evaluation of two land surface schemes used in terrains of increasing aridity in West Africa[J].Journal of Hydrometeorology.9(2): 173-193.

SteeneveldG JTolkL FMoeneA Falet2011.Confronting the WRF and RAMS mesoscale models with innovative observations in the Netherlands: evaluating the boundary layer heat budget[J].Journal of Geophysical Research Atmospheres, 116: D23114.

SuB2002.The Surface Energy Balance System (SEBS) for estimation of turbulent heat fluxes[J].Hydrology and Earth System Sciences, 6: 85-100.

TatarskiiV I1961.Wave propagation in a turbulent medium[M].McGraw-Hill Book Company Inc., New York, 285pp(中译本: 湍流大气中波的传播理论 [M].译者: 温景嵩等, 科学出版社, 1978, 232 pp).

van KesterenBHartogensisO KVan DintherDalet2013.Measuring H2O and CO2 fluxes at field scales with scintillometry: Part I-Introduction and validation of four methods[J].Agricultural and Forest Meteorology178/179: 75-87

WardH CEvansJ GHartogensisO Kalet2013.A critical revision of the estimation of the latent heat flux from two-wavelength scintillometry[J].Quarterly Journal of the Royal Meteorological Society, 139: 1912-1922.

WardH CEvansJ GGrimmondC S Balet2015.Infrared and millimetre-wave scintillometry in the suburban environment-Part 1: Structure parameters[J].Atmospheric Measurement Techniques, 8: 1385-1405.

WardH C2017.Scintillometry in urban and complex environments: A review[J].Measurement Science & Technology28(6): 064005.DOI: 10.1088/1361-6501/aa5e85.

Wang J M, 2018.Area-averaged Flux Measurements and Scintillometry[Z].3rd Intl.Potsdam GHG Flux Workshop, Nanjing, China, 22-25 October, 2018.

WangTOchsG RCliffordS F1978.A saturation-resistant optical scintillometer to measure Cn2[J].Journal of the Optical Society of America 68(3): 334-338.

WeselyM L1976.The combined effect of temperature and humidity fluctuations on refractive index[J].Journal of Applied Meteorology15(1): 43-49.

WyngaardJ CIzumiYCollinsS A1971.Behavior of the refractive-index-structure parameter near the ground[J].Journal of the Optical Society of America61(12): 1646-1650.

ZhangHZhangH SCaiX Halet2016.Contribution of Low-Frequency Motions to Sensible Heat Fluxes over Urban and Suburban Areas[J].Boundary-Layer Meteorology161(1): 183-201.

ZhaoPXuX DChenFalet2018.The third atmospheric scientific experiment for understanding the earth-atmosphere coupled system over the Tibetan Plateau and its effects[J].Bulletin of the American Meteorological Society99(4): 757-776.

蔡旭晖, 2008.湍流微气象观测的印痕分析方法及其应用拓展[J].大气科学32(1): 123-132.

卢俐, 刘绍民, 徐自为, 等, 2009.不同下垫面大孔径闪烁仪观测数据处理与分析[J].应用气象学报20(2): 171-178.

王介民, 1999.陆面过程实验和地气相互作用研究—从HEIFE到IMGRASS 和GAME-Tibet/ TIPEX[J].高原气象18(3): 280-294.

文章导航

/