1 引言
不同于森林、草地、农田等均一下垫面，城市下垫面具有其特殊的复杂性。由于城市人口不断增加，城市不断扩张，越来越多的人认识到城市碳收支、城市热岛、城市污染等城市特殊气候效应的重要性。许多科学家在城市密集区展开了通量观测实验，得到了以下一些结果：在美国巴尔的摩，城市森林板块对CO2通量的吸收使得整个通量贡献区内为CO2的汇，郊区生态系统对城市CO2排放有明显的补偿作用(Grimmond et al.,2002)；苏格兰爱丁堡市城市平均CO2通量为26 μmol m−2 s−1，其高峰期与人为活动和交通流高度有关(Nemitz et al.,2002)；在瑞士巴塞尔高密度建筑群街区，夏季白天于17:00 CO2通量可达到最高峰14 μmol m−2 s−1，对应交通高峰(Vogt et al.,2005)；Coutts et al.(2007)在澳大利亚两个人口密度和植被覆盖密度不同的区域进行了超过1年的全天候CO2通量观测实验，结果显示峰值排放出现在早晨和下午，冬季早晨达到最大排放，夏季CO2排放量低于动机，排放强度和模态与北半球相似；Matese et al.(2009)在意大利佛罗伦萨，Crawford et al.(2011)也在墨尔本发现CO2排放具有季节变化的现象。由于城市复杂的地表特性，通量的源和汇分布十分不均，导致直接测量以及定量的描述城市通量变得非常困难(Liu et al.,2011)。近些年，随着观测设备的不断发展，涡动协方差观测系统在城市CO2观测中得到广泛应用，但由于涡动协方差观测系统传感器都架设在一定的高度上，其测量结果只反映特定点下垫面或某部分下垫面的物理过程，使得由点到面或者由线到面的空间代表性就成为通量观测中不可忽视的问题(宫丽娟等，2009)。通量观测足迹和源区的提出正是为了解决这一问题。足迹函数(Footprint)又称源权重函数，是描述近地面层表面源或汇的空间分布和仪器观测通量值之间关系的函数(Schimid，1997)。确定某一观测点通量数据的源区无疑是进行通量观测资料的质量控制和有关研究的基础。尤其是在观测城市这样的复杂下垫面，由于城市源的非均匀性给观测和分析带来困难，需要我们事先对观测地点进行足迹分析，来验证观测资料是否涵盖了能代表城市的各个属性。
通量足迹模型经过了十几年的发展，已发表并广泛应用的有Horst-Weil(Horst and Weil，1992)、Kormann and Meixner(Koemann and Meixner，2001)和FSAM(Flux Source Area Model)(Schmid and Lloyd，1999)模型等。Schmid and Lloyd(1999)提出了主要源区尺幅的参数化公式即FSAM模型，它物理机制明确，使用简单，得到广泛应用。Leclerc et al.(1999)提出了适合足迹函数的大涡模拟。Kormann and Meixner(2001)，Kormann et al.(2003)提出机遇非中性层结条件下的足迹模型，其结果与拉格朗日模型进行了比较，二者模拟的源区面积范围大致相同。拉格朗日模型的最大点比Kormann模型较近。米娜等(2006)应用FSAM模型对中国通量观测网络中8个站点通量数据的空间代表性进行研究。赵晓松等(2005)年利用足迹模型分析了长白山阔叶红松林通量观测的足迹函数及源区分布。以上这些文章都是针对均质下垫面而言。对于城市这样的复杂下垫面，尤其作为研究城市通量的观测点，其周围的景观具有异质性，风向和风速的日变化，季节变化都有来自于气象塔周围不同方向，不同距离区域的贡献，在分析任何城市通量前，都需要对测量地点进行足迹分析，目前还没有完整的针对城市下垫面通量的足迹分析。北京325 m气象塔以高度高，分层细，全天候的优点成为本文研究城市通量足迹分析的首选。利用气象塔7层水汽—二氧化碳通量资料和超声风速仪观测数据，找寻不同观测高度、不同风向与足迹模型的关系。由于城市下垫面的复杂性，城市边界层中各参数较之其他下垫面有很大的不同。使得很多模式，如KM、FSAM模型，由于输入参数具有严格的局限性而不能适用。因此，本文利用北京325 m气象塔涡动观测资料，研究讨论一种适用于城市下垫面的足迹模式(如图 1所示)。
图 1Fig. 1
图 1 Footprint示意图
Fig. 1 Sketch of Footprint
2 Hsieh 模型介绍
本文选用Hsieh et al.(2000)在拉格朗日随机扩散模型和维度分析原理基础上建立的解析模型，即Hsieh模型(如图 1所示)。
Hsieh模型有3个重要的输入变量：Zm/Zo、Zm/L、u*(Zm=Z-d，Z为观测高度，单位为m；d为零平面位移，单位为m；Zo是空气动力学粗糙度，单位为m)。本文粗糙度计算方法采用的是最小二乘拟合迭代法，利用风速廓线拟合迭代的方法，能得到零平面位移d和Zo，计算公式如下：
拟合成二元方程的形式为
其中，u为风速，θ为位温，其他参数分别表示为
Z0为需求的空气动力学粗糙度，Zt为热力学粗糙度，L为奥布霍夫长度。Ψm、ΨA为稳定度订正函数，k为冯卡门常数(k=0.4)。u*为摩擦速度。将零平面位移d以一定范围的步长增加变化，拟合出风速和温度廓线，计算出相应的一系列空气动力学粗糙度，拟合相关系数。当拟合关系数最大的时候对应的空气动力学粗糙度即为所求，不同方法确定方法得到的Zo值不同，因此得到的贡献区分布也会有一定的差异 [具体可以参照周艳莲等(2007)]；L是奥布霍夫长度，单位为m；u*是摩擦速度，以及3个重要的输出变量：x、F/So、(x为上风向有效风浪距离；F/So为通量贡献权重系数，范围在0～1之间；F为通量，单位为g m−2 s−1；So为源强，单位为g m−2 s−1；是足迹函数值)。下面简述其计算方法。
拉格朗日随机扩散模型中用到了三个参数：x，Zm/Zo，Zm/L。其中x又是F、Zm、Zo组成的函数，因此Hsieh用Zm和Zo重新定义了一个尺度长度Zu(单位为m)：
图 3Fig. 3
图 2 7层涡动协方差观测系统(由CAST3三维超声风速仪和Li-7500a二氧化碳/水汽分析仪构成)
Fig. 2 Eddy covariance system of 7 levels(CAST3 and Li-7500a)
这样，拉格朗日模型改编后有了新的三个变量，将他们用于维度分析：x，L，Zu。用L作为主变量，缩减3个变量至2个，则可以得到与之间的关系：
Hsieh总结出了这样的经验关系式：
其中，k是冯卡门常数，取0.4，D、p是作者设定的常数参数，在不同稳定度条件下有不同值。按照大气稳定度分成三类：稳定(0.2＜Z/L)、不稳定(Z/L＜－0.04)、中性不包含弱稳定(|Z/L|＜0.04)。其中，不稳定条件下：D=0.28，p=0.59；中性条件下：D=0.97，p=1；稳定条件下：D=2.44，p=1.33。则周期通量贡献系数计算式为
足迹函数是周期通量贡献系数的导数：
将(10)式足迹函数求导后结果取0，可以算出函数的最大值对应的x的值，就可以找到最大通量贡献区的位置：
公式(9～11)就组成了Hsieh足迹分析模式的基本部分。
3 研究区域简介及观测实验概况
3.1 区域简介
中国科学院大气物理研究所气象塔建成于1979年8月，气象塔位于39°58′N，116°22′E，位于北京正北方向。海拔高度49 m，塔高325 m，塔上设有15层观测平台，高度分别为8、15、32、47、63、80、102、120、140、160、180、200、240、280以及320 m。每层平台都装有风速仪、风向仪和温度仪。为减小塔体阻挡气流对测风仪造成的影响，分别在西北和东南两个盛行风向各装一台测风仪。主要为北京空气污染以及大气边界层、大气湍流扩散等研究服务。在世界众多气象塔中，北京气象塔具有其独特的优势，它的高度是亚洲第一、世界第二，塔上的观测平台分布较细致，更重要的是它经历了北京的城市化进程，其所在地由20世纪90年代以前的市郊变成现在的市区三四环之间(刘辉志和洪忠祥，2002)。距离气象塔西边约0.3 km处有一个小型公园由树木、草坪和河流组成。塔的东边是青藏高速有大量的车辆每天由南向北行驶。另外一条很忙碌的北辰路在塔的北边。塔的北边和南边是密集的居民区，居民区内有小面积低矮的草丛(图 3h)。
图 3Fig. 3
图 3 2012年5月～2012年11月7层高度的风向、风速(m s−1)及其分布频率玫瑰图和325 m气象塔周围环境照片(单位：%)
Fig. 3 Wind direction，wind speed(m s−1)，and their distribution frequency and a google map of 325 m tower(in percentage)
3.2 实验概况
为满足中国科学院战略先导科技专项“应对气候变化的碳收支认证及相关问题”的课题“碳卫星验证系统与综合观测”的第一子课题“高塔验证系统和综合观测”的需要，2012年于气象塔上(分别位于8 m、16 m、47 m、80 m、140 m、220 m和280 m处)新架设了7层涡动协方差观测系统(见图 2)，由最基本的CAST3三维超声风速仪(Campbell Scientific，USA)和Li-7500a 二氧化碳/水汽分析仪构成(Li-cor，USA)。采样频率为10 Hz，数据由采集器(CR5000X Campbell Scientific，USA)自动采集并存储。并在每个仪器端口处直接连接出一根数据线连入电脑，这样做是为保证每一层数据采集下来的时候是同步的。
3.3 实验数据预处理
数据选自2012年5月～2012年11月共7个月，7层二氧化碳/水汽通量数据和超声风速仪数据。采用Edire涡动系统通量处理软件(王介民等，2007)对北京气象塔观测数据进行质量控制及评价，对Li-7500a观测资料进行野点去处并计算平均值；平面拟合法进行坐标旋转；水汽二氧化碳相对于垂直风速的时间滞后修正；计算半小时统计平均值；计算摩擦速度、地表粗糙度。对水汽通量和二氧化碳通量进行空气密度脉动修正后，计算奥布霍夫长度和大气稳定度。超声数据原始数据进行阈值、去野点计算半小时三维风速平均值、水平风速平均值和风向角(姜明等，2011)。
4 结果讨论与分析
4.1 观测点主风向分析
通量观测中，通量贡献区域通常分布在观测点盛行风向一侧(马迪等，2009)。关注观测点的通量变化时，首先要满足实际风向的需要，以便整体了解城市复杂系统通量的时空变异情况。图 3给出了以气象塔为中心，22.5°步长，7层高度上全年风向频率统计的玫瑰图，同时按照风速大小分类，统计出每个子集中风速发生的频率。为方便分析，右下角给出以325 m气象塔为中心，半径为1 km的卫星地图。图中白色是325 m气象塔位置所在。
在全年尺度上(由于篇幅有限，本文未分季节讨论)可见，平均风速随着高度增加而增加，较大的风速多出现在西北方向。风向上，呈现高低层主风向不一致的现象，并有位移区(当来流偏北)或空腔区(当来流偏南)。观测塔低层8 m，16 m主要风向以西北风和东南风为主。47 m、80 m呈现比较紊乱的状态。越往高处，转为偏角度南北风。这是由于97年后北京三环不断发展，出现不少高层建筑物，从而改变了气象塔周围流场(彭珍和胡非，2006)。就如今建筑物的分布状况而言，当来流是偏北风时，前方建筑物(冠城园，位于观测点南面)的抬升作用，气流出现分流绕流；而当来流是西南风的时候，气象塔所处的位置刚好处于空腔区，气流还来不及恢复到与来流方向一致，反映在风向玫瑰图上(图 3a至图 3g)，可以看到西南方向很少存在气流。另一方面，由于气象塔周围建筑物的不断修建，因此近地面层的风速扰动逐年加大，湍流结构趋于复杂，反映在风向玫瑰图上就是城市冠层(47 m、80 m)风向都很紊乱，看不出明显的主风。这些现象不同于其他下垫面(草原、森林、农田等)，是人为造成的城市特有复杂风场。
4.2 不同风向下通量源区的变化特征
得到观测点不同高度上的主风向后，为验证“通量贡献区域通常分布在观测点盛行风向一侧”这一观点，本文分别计算了3个高度上(低层8 m、冠层47 m、高层280 m)，以30°为步长，在不稳定、稳定两个大气条件下90%贡献区范围。如图 4所示。
图 4Fig. 4
图 4 不同稳定度条件下90% 通量贡献区范围随风向的变化图(左列为不稳定条件，右列为稳定条件；由上至下观测高度依次为8 m、47 m、280 m)
Fig. 4 Footprint [for 90% of effective fetch(x)] varying with wind directions under different stability for 8 m，47 m，280 m(from top to bottom)
结果显示，只有当大气处于稳定条件下时，通量贡献区范围大小才与主风向有着无显著差异的一致性。观测高度8 m处，主风向为西北和东南风，通量贡献区则在西北方向延伸至300 m、47 m和280 m也在主风向上有通量贡献区范围较大的现象。而当大气处于非稳定条件下，如左列图所示，通量贡献区范围呈现出近圆形。无论在主风向或者在非主风向上，范围大小近似相同。可以看到8 m、47 m、280 m处，通量来源于观测点附近以30、300、9000 m为半径的区域。
4.3 不同稳定度对通量源区变化特征的影响
从上图 4中已经大致可以看出，稳定条件下通量贡献区范围要大于不稳定条件下的范围。为了进一步验证此结论，分别在主风向下(8 m、16 m主风向为112.5°和292.5°；47 m主风向为22.5°、180°、310°；80 m主风向为0°、22.5°、202.5°、227.5°；140 m主风向为0°、22.5°、50°、202.5°、220°；200 m、280 m主风向为180°、202.5°、315°、337.5°)计算了7层观测高度，三个大气稳定度下的足迹函数。得到7层高度周期通量贡献系数F/So在不同稳定度条件下随主风向有效风浪区(x)的变化(图 5 a、b、c)，以及7层高度足迹函数(源权重函数)随主风向有效风浪区(x)的变化(图 5d、e、f)。7层高度不同大气稳定度下90%通量贡献区平均范围和最大源强到观测点距离记录于表 1。观测高度为47 m处的通量贡献区范围与宋涛于2011年得到的47 m结果无明显差异(Song and Wang，2012)。
图 5Fig. 5
图 5(a、b、c)周期通量贡献系数F/S0随有效风浪距离x变化关系；(d、e、f)足迹函数(源权重函数)f(x，Zm)随有效风浪距离x变化关系。(a、d)为不稳定条件；(b、e)为中性条件；(c、f)为稳定条件
Fig. 5(a，b，c)Contribution coefficient of flux(F/S0)vs. the main wind distance of effective fetch(x);(d，e，f)footprint function f(x，Zm)vs. the main wind distance or effective fetch(x).(a，d)are results of unstable stability;(b，e)are results of neutral stability;(c，f)are results of stable stability
从图 5和表 1中共同可以看出，同一高度 下，不同稳定度的F/So随x的变化不同。通常在不稳定大气环境下(白天)，大气垂直运动非常剧烈，物质向上输送很快，浓度变化迅速，很快能被传感器感应到。因此，传感器测得的通量来源于主风向上较近的地方。若大气处于稳定条件下(夜间)时，湍流活动相对较弱，大气输送物质的能力也较弱，物质扩散得缓慢，传感器感应到浓度变化的这种物质可能从很远的地方缓慢输送过来。所以，稳定大气条件下，通量贡献区距离较长。大气处于中性条件下时，介于稳定与不稳定两者之间(胡非，1995)。以47 m为例，不稳定条件下，90%通量贡献区在观测点附近306 m左右，源强中心在16 m，出现在大气所实验区内；稳定条件下，通量贡献区范围能够延伸至1500 m范围远处。此时1500 m包含了北辰西路主干道、八达岭高速路段、居民区、学校、医院、小型公园等。源强中心延伸至79 m，地理上的位置是位于气象塔北面北辰西路主干道、西北面和东南面公园以及南面密集高层居民区。从结果上，可以得到观测高度通量数据结果的代表性。8 m、16 m的观测高度上的观测结果并不能代表城市下垫面的通量贡献区，可能只包含了草地通量贡献，或者水泥路的通量贡献；47 m以上能够代表城市下垫面通量贡献区，280 m的通量贡献区范围很大，已经包含了郊区和城区的权重平均。
表1(Table 1)
表 1 7层高度不同大气稳定度下90%通量贡献区平均范围和最大源强到观测点距离
Table 1 90% of average effective fetch(x) and the distance from measure site to the biggest f(x，Zm)site under different stabilities for seven height levels
稳定度 高度/m 90%通量贡献区 范围/m 观测点至源强中心 距离/m 不稳定 8 35 5 16 133 7 47 306 16 80 1189 26 140 3586 60.4 200 6080 212 280 9725 541 中性 8 91 7 16 223 12 47 441 23 80 2039 45 140 7255 117 200 13148 263 280 22049 1201 稳定 8 132 18 16 1040 28 47 1517 79 80 5425 134 140 18110 396 200 33680 712 280 56800 3001
表 1 7层高度不同大气稳定度下90%通量贡献区平均范围和最大源强到观测点距离
Table 1 90% of average effective fetch(x) and the distance from measure site to the biggest f(x，Zm)site under different stabilities for seven height levels
4.4 不同高度对通量源区变化特征的影响
为验证不同高度下观测得到的数据是具有城市下垫面的代表性，将各个高度主风方向上的参数平均值输入模式分析(同上)。如图 5中所示，无论在大气处于何种稳定度条件下，通量贡献区范围都随着高度的增加而不断扩大。以不稳定条件为例，观测高度为8 m处计算得到的贡献区范围为35 m，47 m处增加到300 m左右，140 m处为3000 m是47 m观测数据计算得到的区域范围的10倍 远，280 m则延伸到10 km，此时包含了西北郊区一带，此时观测到的通量代表了城市和郊区的权重平均。源强中心也有8 m观测高度处的5 m增加到280 m观测高度的541 m，几乎是8 m处计算得到通量观测范围的100倍。
如果能够建立通量贡献区范围与观测高度的关系，那贡献区范围半径距离x就可以作为观测高度Z的函数。能计算出任意高度处的通量贡献区范围，还能预计更高处对应的贡献区。本文将7层高度不同大气稳定度下平均贡献区范围与观测高度建立关系，得到以下拟合关系式(如图 6所示)：
图 6Fig. 6
图 6 90%通量贡献区范围与观测高度关系
Fig. 6 Footprint(for 90% of effective fetch(x))varying with measure heights
值得注意的是8 m以下的部分，呈现出通量贡献区范围与高度呈反比的趋势。但由于缺乏观测数据，有待验证。以上这3个线性关系式能够大致代表整个北京城市区域通量贡献区范围与观测高度的关系。同时，一阶拟合关系是否是最精确的拟合关系这一问题还有待考证，可能这种拟合关系更适合用二阶或者非线性关系，需要在以后的研究中加密观测，是一个值得有待研究的问题。总的来说，目前这种一阶的关系式能够帮助我们摆脱观测实验的局限性，可以估算任意高度上的通量贡献区范围。此外，各层通量代表性的不同，将会导致观测值与实际值之间的差异。观测高度较低时候，通量贡献区域主要是作为二氧化碳汇的草地，导致通量值偏小。观测高度越高，观测的范围越远，观测值蔓延到城郊，城郊以农田为主，也会对二氧化碳有所吸收。这些造成观测值与实际值的差异在本文中暂且没有考虑。尽管通量贡献区模型简化了气流流动的物理过程，在数学上变得简单有一定的局限性，但是依然能够提供有价值的通量信息，可供参考。
5 结论
本文选取2012年5～11月，北京325 m气象塔7层涡动观测数据，计算了不同风向、不同高度、不同大气稳定度下通量足迹函数和贡献源区大小。得到以下结论：
(1)Hsieh模型能够很好的计算城市复杂下垫面的通量贡献源区。给定模式需要的参数，就能够很好的分析10%～90%的通量贡献区域以及最强源和观测点的距离。在精准度不要求到特别精确级别的情况下，该模式结果较好，能够运用到整个北京城市下垫面。
(2)低层(8 m、16 m)的主风向为西北风和东南风；冠层(47 m、80 m)的主风呈现比较紊乱状态，大风速主要集中在西北方向；高层(140 m、200 m、270 m)主风向以偏北偏南风为主。在大气稳定条件下时，通量贡献区范围大小与主风向有着无显著差异的一致性，而不稳定条件下与主风向无关。
(3)在同一高度上，无论哪个高度都出现稳定条件下通量贡献区的范围要大于不稳定条件下范围的现象，中性条件则介于两者之间。8 m，16 m高度上的观测结果不能完全代表城市下垫面通量贡献区，47 m以上能够代表城市下垫面通量贡献区，280 m则已经包含了郊区和城区的权重平均。
(4)在同一大气稳定度条件下，高度越高通量贡献区范围越大，90%通量贡献区范围与观测高度成线性关系：y=38.6x－697.5(不稳定)，y=87.2x－1911.5(中性)，y=233.5x－4833.9(稳定)。这种线性关系可以预测没有观测高度或者更高处的通量贡献区范围。
致谢 感谢中国科学院大气物理研究所李爱国高工在本文使用的观测资料方面给与的帮助。