2016, Vol. 35
2. 中国气象局北京城市气象研究所, 北京 100089
空气动力学粗糙度长度Z0作为影响陆气之间能量、水分交换的动力因子，不仅能反映下垫面特征，而且能有效刻画陆气之间物质、能量及水分交换过程强弱、大小，是众多包含动量通量、显热通量和潜热通量的“陆面-大气”模式中一个非常重要的参数(杨耀先等，2014)。Z0大小不但取决于下垫面地表粗糙性质，如粗糙元种类、大小、形状、高度、密度、植被盖度、排列方式以及运动与否等，而且还与流经地表的流体的性质有密切关系，如风速、风向等(张强和吕世华，2003；石雪峰等，2006；胡张保和俞炳丰，2008；王玉玉等，2014)。由于下垫面的复杂非均匀性特征，Z0在区域范围内存在明显空间差异，而这种不均匀性对网格地表湍流通量的估算有重要影响(陈斌等，2010)，但由于下垫面尤其是城市下垫面的建筑物粗糙元特征在空间区域上难以准确获取，这给在区域范围内估算含有植被和城市的复杂下垫面Z0的估算带来了很大困难。
目前常用的Z0计算方法分为两类：一类是气象观测法，最常用的为廓线法，即根据通量塔或气象塔的风速资料代入到对数风速廓线方程中应用最小二乘方法求解，或者采用质量守恒定律、中心压力法、阻力法等方法求解；还有一种常用方法是利用单一高度上快速响应的超声风温仪观测的风、温湍流资料，根据方差随稳定度变化的规律进行计算，这种方法又称方差法。由于气象观测方法一般基于少量的单一通量塔(站)或气象塔，估算结果只能代表塔(站)周围下垫面Z0状况，难以在区域空间范围上开展应用。另一类是形态学方法，即根据粗糙元的几何形状、分布密度等计算Z0，目前应用的比较广泛的有根据粗糙元高度、粗糙元顶部面积、粗糙元迎风面积计算。由于粗糙元的几何形态有可能在不同空间大小内通过卫星遥感或GIS技术获得，这为在区域范围上估算复杂下垫面的Z0带来了可能。
由于粗糙元性质的差异，当前对Z0的研究主要针对植被下垫面和城市下垫面。国内对植被下垫面空气动力学参数的研究较多，且大部分采用了气象观测法来估算Z0。覃文汉(1994)应用中心压力法对山东禹城小麦、水稻、大豆及玉米作物的Z0研究发现，Z0不仅随作物群体结构变化，而且还受大气层结稳定度和群体内外湍流交换强弱的影响，呈现出明显的日变化规律。刘和平等(1997)、钟中和韩士杰(2002)、赵晓松等(2004)利用气象观测法分别估算了长白山阔叶林的Z0，结果表明Z0存在着较明显的季节变化，在生长季Z0较小而非生长季恰好相反。何奇瑾等(2007)对芦苇湿地和蔡福等(2013)对玉米农田的研究表明Z0的季节变化受风速、株高和叶面积指数影响。张学仕等(2010)对江苏句容次生林研究表明Z0随叶面积指数LAI增大而增大，当LAI增大到一定值时，Z0随着LAI的增大而变化不大。近年来国外除了气象观测法，形态学方法也得到较为广泛应用。Raupach(1994)提出了用植被冠层高度和植被叶面积指数估算植被Z0的简单形态学方法。Schaudt et al(2000)考虑到Raupach(1994)方法的实验数据主要是固态物质和茂密森林，不能代表稀疏植被，同时参考Lindroth(1993)方法对Raupach方法进行了补充，分针叶、阔叶、混交林三种类型来计算这三种类型的Z0。Jasinski et al(1999，2005)在Raupach(1994)方法基础上提出了基于四种植被类型Z0的计算参数，并利用MODIS(中分辨率成像光谱仪)卫星数据估算了美国南部大平原的归一化粗糙度长度Z0/Zh(其中Zh为植被高度)。Borak et al(2005)则根据这四个类型的参数推导了美国南部平原植被区域相对植被高度的时间序列Z0/Zh。杨阿强等(2011)则以MODIS卫星产品数据作为数据源，采用形态学方法估算了中国东部2001-2003年植被地区时间序列的Z0/Zh，但大部分情况下由于植被高度Zh的难以获取，Z0还不能直接计算得到。李沁怡等(2014)利用2010年高分辨率全球土地覆盖卫星观测数据估算了中国区域的Z0，但城市区域的Z0仅为一个固定值，无法反映城市区域内Z0的变化。
对于城市下垫面，粗糙元的分布非常复杂，空间分布上存在很大的非均匀性，如何确定城市下垫面粗糙度长度Z0，目前尚无最佳方案(Grimmond et al，1998)。王金星等(2002)和胡张宝等(2008)分别概述了城市空气动力学参数的多种计算方法，主要有最小二乘拟合迭代法、牛顿迭代法、温度方差法(TVM)、Martano法(单层超声风速仪法)、陈家宜法(超声风温仪法)、粗糙元形态法等。高志球等(2002)、陆龙骅等(2002)分别利用北京325 m气象塔观测资料估算了气象塔附近的下垫面Z0，并发现城市下垫面的Z0也随季节变化；梅凡民等(2002)利用风洞试验分析了粗糙元几何参数对Z0的影响；李倩等(2003)、徐阳阳等(2009)则分析了城市化发展对Z0的影响。张强和吕世华(2003)对城市下垫面的Z0研究指出，城市建筑高度、建筑覆盖率是Z0最主要的影响因子，并给出了一个以建筑高度和建筑覆盖率为主要变量来估算Z0的简单公式，但并没有对估算结果通过实例进行验证分析。韩素芹等(2008)利用天津255 m气象塔资料估算了气象塔各方位Z0，并分析了城市化对Z0的影响。Grimmond and Oke(1999)引入多种形态学动力计算方法来计算城市Z0，并对各种方法进行了比较评估，认为Raupach(1994)的方法同其他方法相比在计算Z0时效果比较理想，但仍同实际观测数据有较大差别，并指出如何利用城市地理信息和卫星资料，结合边界层观测实验资料来确定城市形态和粗糙元空间分布特征是今后的研究趋势。
目前，国内对植被和城市下垫面Z0的估算仍以气象观测法为主，而且针对的都是小范围或单个点的估算结果，很少有既含植被又含城市的大范围地区Z0的空间分布估算，这是气象观测法难以做到的。同时在城市地区形态学方法应用上，大部分是针对的实测或模拟的单点建筑群开展的相关研究，也很少有大范围区域城市地区Z0的估算研究。而遥感技术在大范围内获取的植被地表参数信息和高精度的城市GIS地理基础信息资料为区域尺度上获取空间区域的空气动力学参数Z0提供了可能。本研究将以北京地区为例，采用卫星数据和GIS技术开展区域尺度上空气动力学Z0的估算方法研究，对获取大范围非均匀复杂下垫面的Z0具有重要参考意义。
2 研究区域和方法 2.1 研究区域和思路以北京市为研究区域。北京市位于华北平原北端，在39°28′N-41°05′N，115°25′E-117°30′E之间，属于典型的暖温带半湿润大陆性季风气候，春秋短而冬夏长，常年平均气温10 ~ 12 ℃，年降水量在470 ~ 600 mm之间。通过北京市2012年主要植被覆盖类型(图 1)可以看出，下垫面类型比较复杂，北京西、北和东北为山区森林地带，主要植被类型为阔叶林、针叶林和灌丛，东南是海拔在20 ~ 60 m之间的北京小平原，主要下垫面类型为农田、水体、绿地、裸地和城镇用地，其中城市建成区主要位于北京市六环内。
|
图 1 2012年北京市主要植被覆盖类型 Figure 1 Main vegetation types of Beijing city in 2012 |
由于植被和城市为两种差异很大的地表粗糙元，因此把研究区划分为北京六环外以植被为主的北京主要植被地区和六环内以建成区为主的城市地区，分别来估算下垫面空气动力学参数粗糙度长度Z0。另外由于北京地区没有实测的区域Z0值分布，而且以往研究对植被地区归一化的粗糙度长度Z0/Zh较多，因此本文将利用卫星数据和GIS技术分别来估算主要植被地区的Z0/Zh和Z0以及城市地区的Z0，对这些估算结果和以前相关研究成果进行对比分析验证，用来评估本文方法模型的可靠性与合理性。
2.2 植被粗糙度长度估算Raupach(1994)提出了一种简单的估算各向同性的植被粗糙元Z0的形态学模型方法，可以用下式表达：
| $\frac{{{Z}_{d}}}{{{Z}_{h}}}=1.0-\frac{1.0-exp-\left[ {{\left( {{C}_{d}}\cdot 2{{\lambda }_{F}} \right)}^{0.5}} \right]}{{{\left( {{C}_{d}}\cdot 2{{\lambda }_{F}} \right)}^{0.5}}}=\left( \frac{\beta \cdot \Lambda }{2+\beta \cdot \Lambda } \right)\left[ 1.0-\alpha {{\gamma }^{-1}}{{\Lambda }^{-1/2}} \right],$ | (1) |
| $~\frac{{{Z}_{0}}}{{{Z}_{h}}}=\left( 1.0-\frac{{{Z}_{d}}}{{{Z}_{h}}} \right)exp\left( -k\frac{{{U}_{h}}}{{{u}_{*}}}+{{\psi }_{h}} \right),$ | (2) |
| $~\frac{{{U}_{h}}}{{{u}_{*}}}=\gamma ={{\left( {{C}_{s}}+{{C}_{R}}{{\lambda }_{F}} \right)}^{-0.5}}\cdot exp\left( \frac{c\cdot {{\lambda }_{F}}\cdot \gamma }{2} \right),$ | (3) |
| $\frac{{{U}_{*}}}{{{u}_{h}}}={{\left( \frac{{{u}_{*}}}{{{U}_{h}}} \right)}_{max}}={{\gamma }_{max}}^{-1}~,当{{\lambda }_{F}}＞{{\lambda }_{Fmax}}或\Lambda ＞{{\Lambda }_{max}}$ | (4) |
| $\alpha ={{C}_{d}}{{\left( \frac{2b}{{{Z}_{h}}} \right)}^{1/2}}~,\beta =\frac{{{C}_{R}}}{{{C}_{S}}},\Lambda =2{{\lambda }_{F}}$ | (5) |
式中：Zd为零平面位移高度(单位：m)；Z0为粗糙度长度(单位：m)；Zh为植被粗糙元高度(单位：m)；Zd/Zh为归一化的零平面位移高度；Z0/Zh为归一化的粗糙度长度；λF为单位地表面积上植被粗糙元迎风面积比；Λ为植被冠层面积指数(对各向同性植被粗糙元有Λ=2λF)；Λmax为Z0/Zh最大值时对应的植被冠层面积指数；ψh为粗糙层影响函数；Cs为地表阻力系数；CR为粗糙元阻力系数；c为粗糙元遮阴系数；b为粗糙元冠层宽度(单位：m)；Cd为经验系数；Uh为风速(单位：m·s-1)；u*为摩阻速度(或剪切速度)(单位：m·s-1)；k =0. 4为卡曼常数。
在估算Z0时，Zh和Λ为模型输入变量，而CR、Cs(或β)、α、c、γmax、ψh和Λmax等7个参数为模型参数，不同植被类型具有不同的形态模型参数。Jasinski et al(2005)在估算美国南部大平原的归一化粗糙度长度Z0/Zh时，以Raupach(1994)方法为基础，结合实地测量数据，根据MODIS卫星的IGBP地表覆盖类型产品建立了森林、灌丛、草地、农田的Raupach模型的各种参数。因此，结合研究区域的植被覆盖类型(图 1)，本研究采用Jasinski et al(2005)确定的各种模型参数来估算研究区域植被地区的Z0。模型参数确定后，根据式(1)和式(2)，估算Z0就只剩下模型输入参数中植被冠层面积指数Λ和植被冠层高度Zh。
2.2.1 植被冠层面积指数Borak et al(2005)在计算美国南部植被地区的Z0/Zh和杨阿强等(2011)在估算中国东部植被地区的Z0/Zh时，均采用了Zeng et al(2002)定义的植被Λ计算方法，即：
| $~\left[ \Lambda \cong {{I}_{g}}+{{I}_{s}} \right],$ | (6) |
| ${{I}_{s}}^{n}=max\left\{ \left[ \alpha {{I}_{s}}^{n-1}+max{{I}_{g}}^{n-1}-{{I}_{g}}^{n},0 \right],{{I}_{S}}{{,}_{min}} \right\}~,$ | (7) |
式中：Ig为绿色叶面积指数，即叶面积指数LAI；Is为干枯叶面积指数，包含树的茎干面积指数和枯萎凋落的叶面积指数；n代表月份；Isn代表第n月份的绿色叶面积指数；Is，min为不同植被覆盖类型的月最小干枯叶面积指数；1-α代表每月枯叶凋落率。不同植被覆盖类型Is，min和α有不同，Zeng et al(2002)定义了MODIS IGBP不同植被覆盖类型的Is，min和α参数表，根据该类型参数表和研究区域植被覆盖类型(图 1)可获取北京地区林地类型的Is，min和α分别为1. 0和0. 5，农田类型的Is，min和α分别为0. 1和0，裸地类型的Is，min和α分别为0. 5和1. 0。
绿色叶面积指数Ig可直接从遥感反演的叶面积指数LAI直接获取。在本研究中直接利用EOS/Terra 卫星的MODIS叶面积指数产品MOD15A2(第5版本)，可从美国国家航空与航天局网站(http://reverb.echo.nasa.gov)免费下载获取。MOD15A2是由8 天合成的1 km 分辨率的HDF格式文件，具体包括LAI、光合有效辐射分量(FPAR)、质量评价等数据层。作为L4级产品，MOD15A2已根据参数文件中的参数，对图像进行了几何纠正；其次，还根据质量评价数据集，对LAI进行了基本质量评估，LAI数据总体质量较好。在这里采用平均值法可获取北京地区2012年1、4、7、10月份的叶面积指数LAI。对MOD15A2北京区域的LAI产品质量分析发现，平原地区LAI值为零的范围过大，接近北京平原城市用地面积的2倍。这很可能是由于MOD15A2所采用的1 km土地覆盖类型中城市用地面积过大导致(在该产品中均定义了城市用地的LAI均默认为零)，因此需对这部分非城市用地LAI进行补缺修正。补缺修正方法采用王修信等(2014)建立的北京林地LAI与归一化植被指数NDVI的多项式回归关系式：
| $\begin{align} & LAI=59.\text{ }408NDV{{I}^{2}}-49.\text{ }469NDVI+12.\text{ }070\text{ } \\ & {{R}^{2}}=0.\text{ }583 \\ \end{align}# | (8) |
在这里分别选用2012年1、4、7、10月北京地区MODIS/Terra1B 晴空遥感影像第一、二波段计算归一化植被指数NDVI，根据式(8)就可获取北京补缺范围用地的叶面积指数LAI，合并前面的MOD15A2 LAI产品后即可得到北京主要植被地区的绿色叶面积指数Ig，再根据式(6)和(7)就可估算植被Λ。通过估算的北京地区2012年1、4、7、10月的植被Λ(图 2)可以看出，夏季7月份Λ均处于一年四季中的高值，山区森林一般在2. 0以上，局部可达6. 0，平原农田一般为1. 0 ~ 3. 0；其次为秋季10月份，春季4月和冬季1月Λ都很低，仅春季部分农田地区Λ高于冬季，平原地区Λ一般在1. 0以下，山区森林大部分为1. 0 ~ 2. 0。
|
图 2 2012年1(a)、4(b)、7(c)和10(d)月北京地区植被冠层面积指数估算结果 Figure 2 Vegetation canopy area index of Beijing in January(a)，April(b)，July(c)and October(d)in 2012 |
北京地区森林植被主要分布于山区，平原地区则主要以农田植被为主，也分布着一些以林地为主的绿化植被。森林植被高度随季节变化较小，而平原区农田植被则由于生育物候期的变化，植株高度变化明显。因此针对北京地区植被高度的估算在这里不考虑森林植被高度的季节变化，仅考虑平原农田植被高度的季节变化。森林地区的植被高度Zh可以通过卫星雷达反演获取(周志强等，2011)。在本研究中采用美国卫星Geoscience Laser Altimeter System(GLAS，地球科学激光测高系统)反演的全球1 km空间分辨率的森林植被高度数据，可以从美国国家航空与航天局地球观测数据信息中心网站(http://webmap.ornl.gov/wcsdown)免费下载获取。GLAS数据已经在全球的许多区域成功地用于森林地区植被结构参数的提取。在这里采用Simard et al(2011)制作的2005年全球1 km森林植被高度图中的北京区域数据，该反演结果与对崔国发等(2000)、宋秀瑜(2010)研究区域森林植被高度空间分布趋势具有较好一致性，但是对比发现GLAS反演的森林高度数值偏高明显，本文结合GLAS反演的森林高度与崔国发等(2000)、宋秀瑜(2010)确定的阔叶林、针叶林和灌丛等植被类型平均高度进行了北京森林植被高度订正。北京平原地区农作物由于其固有的生理物候特征，植被高度在不同季节基本上是固定的，可由当地农业气象观测站获取。从北京地区夏季植被高度估算结果(图 3)可知，北京地区夏季的山区森林植被高度主要在4 ~ 15 m之间，其中5 m以上的林区占到近80%。平原地区夏季农田以玉米为主，植被高度一般在0. 2 ~ 2. 0 m，局部绿化林地高度在2. 0 ~ 6. 0 m之间。
|
图 3 北京地区夏季植被高度估算结果分布 Figure 3 Distribution of plant height estimation |
求得植被Λ与植被高度Zh后，根据式(1)和(2)就可估算出零平面位移高度Zd和粗糙度长度Z0。
2.3 城市地区粗糙度长度估算在本研究中采用Grimmond et al(1999)根据Raupach(1994)模型建立的城市形态学模型,可由下式表达：
| $~\frac{{{Z}_{d}}}{{{Z}_{h}}}=1.0-\frac{1.0-exp\left[ -{{\left( 7.5\times 2\times {{\lambda }_{F}} \right)}^{0.5}} \right]}{{{\left( 7.5{{\lambda }_{F}} \right)}^{0.5}}},$ | (9) |
| $~\frac{{{Z}_{0}}}{{{Z}_{h}}}=\left( 1.0-\frac{{{Z}_{d}}}{{{Z}_{h}}} \right)exp\left( -0.4\times \frac{{{U}_{H}}}{{{u}_{*}}}+0.193 \right),$ | (10) |
| $~\frac{u*}{{{U}_{h}}}=min\left[ {{\left( 0.003+0.3\times {{\lambda }_{F}} \right)}^{0.5}},0.3 \right],$ | (11) |
其中：λF为单位地表面积上城市建筑迎风面积比，也称建筑截面积指数。建筑截面积指数是一个涉及建筑形态多个因子的复杂变量，与城市建筑覆盖率λP有密切关系。在Grimmond and Oke(1999)的研究中，采用了复杂的实地测量建筑形态因子来估算λF，本研究中由于无实地测量建筑形态因子，采用了一种简单方法来估算λF。对Grimmond and Oke(1999)研究结果发现大多数城市地区归一化的粗糙度长度Z0/Zh峰值对应的建筑覆盖率λP介于0. 3 ~ 0. 4之间，集中在0. 35左右，因此本研究采用城市建筑覆盖率0. 35为峰值来反推建筑截面积指数λF与城市建筑覆盖率λP的关系，最终得到λF=0.8λp较为理想，从而可以模拟出归一化粗糙度长度Z0/Zh与归一化的零平面位移高度Zd/Zh随建筑覆盖率λP的变化曲线，并估算出在不同建筑覆盖率不同建筑物高度下的Z0，如表 1所示，归一化粗糙度长度Z0/Zh在λP=0. 35时最高，约为0. 1328。在大多数不同建筑覆盖率情形下，低矮楼群(≤6 m)一般Z0 ＜1. 0，中高楼群(6 m＜Zh≤15 m)一般 在1. 0 ~ 2. 0 m之间，高楼群(Zh≥24 m)一般Z0 ＞2. 0，这与Grimmond and Oke(1999)研究中的低矮楼群(5 ~ 8 m)的Z0为0. 3 ~ 0. 8 m、中高度楼群(7 ~ 14 m)的Z0为0. 7 ~ 1. 5 m、高楼群(11 ~ 20 m)的Z0为0. 8 ~ 1. 5 m、很高楼群(＞20 m)的Z0 ＞2. 0比较接近，表明该估算结果具有较好的合理性。
|
表 1 不同建筑覆盖率下的Z0/Zh与不同建筑物高度下的Z0估算结果 Table 1 Z0/Zh of different area percent of building and Z0 of different building height |
根据上述模型，城市地区的Z0可由建筑覆盖率和建筑高度两个因子简单估算得到。北京城市地区的建筑覆盖率和建筑高度由2009年北京六环区域1∶2000基础地理信息中的建筑物信息图层中获取，该图层包含建筑物斑块位置和建筑物楼层数信息，利用GIS空间分析技术可估算不同分辨率网格单元的建筑覆盖率和建筑物楼层数，其中建筑覆盖率为该空间单元内的建筑物面积所占百分比，建筑物楼层数为该空间单元内所有网格建筑物楼层数平均值。根据中国《民用建筑设计通则》(GB50252-2005)，建筑物按其高度和层数可分为低层(1 ~ 3层)、多层(4 ~ 6层)、中高层(7 ~ 9层)、高层(10 ~ 29层)和超高层(30层以上)建筑。
从2009年北京六环内城市地区的100 m网格分辨率的建筑覆盖率和建筑物楼层数(图 4)可以看出，在100 m分辨率网格下，五环区域内建筑物面积百分比一般在10 ~ 80%，其中在二环内城中心和城东南部处于建筑覆盖率高值区，建筑覆盖率达40%以上，但同时二环区域内城中心也是一个建筑高度低值区，这是由于该区域为北京旧城保护区，分布着大面积的平房、四合院等低矮旧建筑，绿化面积少，因此建筑覆盖率高而建筑楼层数处于低值，而城东南部则分布着大面积的建材城和旧民居建筑，建筑覆盖率也较高。二环到四环区域为一建筑楼层数高值区，其中朝阳、海淀、丰台10层以上的高层建筑较多，城区北部的大型居民社区天通苑、回龙观和东南部的亦庄经济开发区以及通州、大兴、房山等城中心也处于建筑楼层数数高值区，一般在4层以上。
|
图 4 2009年北京城市地区建筑面积覆盖率(a)和建筑楼层数(b) Figure 4 Percentage of area(a)and number of stories(b)for buildings in Beijing urban area in 2009 |
因此，利用建筑物楼层数图，按每楼层约3 m进行建筑物高度估算，可得到北京城市地区的建筑高度图，结合建筑覆盖率图，根据式(9)~(11)，即可求取北京城市地区的粗糙度长度Z0。
3 结果分析 3.1 北京主要植被地区归一化糙度长度北京主要植被地区估算的2012年归一化糙度长度Z0/Zh结果(图 5)显示，全年植被Z0/Zh在0. 01 ~ 0. 13之间。结合图 1的植被覆盖类型可知，山区森林Z0/Zh大部分为0. 06 ~ 0. 10，局部超过0. 10；其中冬季非生长季1月份森林的Z0/Zh为0. 08 ~ 0. 10，夏季生长季7月份大部分为0. 06 ~ 0. 08，这与周艳莲等(2006)研究得出的夏季长白山森林和千烟洲森林的Z0/Zh为0. 06 ~ 0. 09比较接近。夏季森林地区的Z0/Zh明显小于冬季，这是由于夏季生长季森林的处于高值，而过于密集的树叶使风不能从缝隙穿过，直接从树冠顶部流过，反而降低了森林对气流的阻力，因此在冬季反而出现Z0/Zh高值、夏季出现低值情形。平原农田的Z0/Zh变化较大，冬季非生长季1月份农田的很小，Z0/Zh处于低值，为0. 01 ~ 0. 08，到春季4月份平原农田以处于生长盛期的冬小麦为主，明显增加，Z0/Zh增加到0. 06 ~ 0. 13，这与周艳莲等(2006)研究得出的小麦在0. 06 ~ 0. 08比较接近；到了夏季生长季7月份(平原农田作物以玉米为主)和秋季10月份平原作物均处于较高值，Z0/Zh大部分为0. 10 ~ 0. 13，这与蔡福等(2013)得出的玉米作物在夏季生长盛期的Z0/Zh为0. 10 ~ 0. 14接近。另外，北京主要植被地区2012年归一化糙度长度Z0/Zh的季节变化符合Borak et al(2005)和杨阿强等(2011)对森林植被和农田植被 的季节变化描述，即森林冬季Z0/Zh大于夏季而平原农田冬季Z0/Zh小于夏季。
|
图 5 2012年1(a)、4(b)、7(c)和10(d)月北京主要植被地区估算的Z0/Zh Figure 5 Z0/Zh in Beijing main vegetation area in January(a)，April(b)，July(c)and October(d)in 2012 Z0/Zh |
通过北京主要植被地区估算的2012年粗糙度长度Z0结果(图 6)可知，全年植被Z0在0 ~ 1. 3 m之间，存在着明显的季节变化。结合图 1的植被覆盖类型可知，山区森林的Z0大部分在0. 6 ~ 1. 3 m之间，怀柔中部与北部、延庆西北和东北部、门头沟西部和房山西北部的Z0位于高值区，这与这些地区森林植被高度较高有直接相关。另外山区森林冬季1月份和春季4月份的Z0明显大于夏季Z0，原因与前面描述的冬季森林地区的 明显大于冬季类似。平原地区的Z0大部分在0 ~ 0. 2 m之间，局部林地区域的Z0为0. 2 ~ 0. 4 m。由于北京平原地区主要是低矮农作物、休闲裸地和水体等，这与Wieringa(1992)把水体、裸地、冰面、平坦草地、低矮作物等的Z0值定义为0. 0002 ~ 0. 25 m之间接近。结合图 2的植被冠层面积指数可知，冬季1月份由于裸地面积较大，作物冠层面积指数很低，也很低，大部分在0. 1 m以下；春季4月份由于作物冠层面积指数增加，Z0值也有增加；夏季7月份作物最大，大部分Z0在0. 1 ~ 0. 2 m之间；秋季10月份减小，相应地Z0也有所降低，大部分降至0. 1 m以下。
|
图 6 2012年1(a)、4(b)、7(c)和10(d)月北京植被地区估算的Z0(单位：m) Figure 6 Z0 in Beijing vegetation area in January(a)，April(b)，July(c)and October(d)in 2012. Unit：m |
通过北京市城区估算的地表粗糙度长度Z0(图 7)可知，主城区的Z0大部分在0. 1 ~ 8. 0 m之间(图 7a)，空间变化差异很大，且明显大于周边平原农田地区的Z0，其中Z0的高值区(＞ 0. 5)大部分集中于北京五环区域内。结合图 4可知，Z0在1. 0 m以上的地区主要分布于北京四环区域内及朝阳东扩CBD地区、石景山地区、望京地区、昌平回龙观、天通苑、大兴城区、通州城区和亦庄经济开发区，这些都是建筑楼群密度大、楼层高的区域，而在二环区域内存在一个Z0相对低值区，大部分Z0在0. 5 m以下，这主要是这个区域为北京旧城历史文化区，以大量的低矮房屋为主，高楼很少，Z0明显低于周边其他地区。
|
图 7 2009年北京主城区(a)和大气所铁塔周边地区(b)估算的Z0 Figure 7 Z0 of main urban area(a)and around the atmospheric tower(b)of Beijing urban area in 2009 |
北京城市动力粗糙度Z0的空间分布研究非常少，只有对北京城市大气物理所附近的325 m气象铁塔(简称大气所铁塔)的动力粗糙度研究较多，因此在这里开展大气所铁塔周边地区的动力粗糙度Z0估算结果(图 7b)与建筑群覆盖率和高度之间关系分析，并与前人研究成果进行对比验证分析。
大气所铁塔位于北京市内德胜门以北约217 m的中国科学院大气物理研究所分部院内，距离北三环路约1 km，经纬度为(39°58′29. 20″N，116°22′20. 00″E)，海拔49 m。该塔建于1974年，1979年8月正式投入使用。目前铁塔周围建起了多座高层建筑，主要建于20世纪90年代后，已从当初城郊下垫面变为现在的城市粗糙复杂下垫面，2001年后城市下垫面变化很小。铁塔现位于三环路与四环路之间，塔东面200 m处为南北走向的八达岭高速公路，北面50 m处为东西走向的北土城西路，南面100 m耸立一排高度为80 m左右的冠城园住宅区。该塔目前是国际上研究城市边界层和监测城市低层大气污染垂直分布的一个难得的理想观测平台，在城市复杂下垫面边界层及城市大气污染研究中将发挥出越来越重要的作用。
(1) 粗糙度长度与建筑覆盖率和建筑高度关系分析：从大气所铁塔周边地区的Z0与相应的建筑覆盖率和建筑高度散点分布和模拟的趋势线(图 8)可以看出，Z0与建筑覆盖率呈单峰曲线变化，随建筑覆盖率升高，Z0有增加趋势，但建筑覆盖率增加到一定值，Z0开始下降，这符合前面建立的城市形态学模型估算Z0原理。Z0与建筑高度则呈线性增加变化趋势，表明建筑高度越高，Z0值越大，这与韩素芹等(2008)利用天津255 m铁塔资料分析得出的“粗糙度与城市建筑群高度之间存在线性关系”的结论基本一致。
|
图 8 大气所铁塔周边地区Z0随建筑覆盖率(a)和建筑高度(b)变化的散点分布和趋势线 Figure 8 Scatter plot and trend curve of Z0 with building percent(a)and building height(b)around the atmospheric tower |
(2) 验证分析：表 2为不同学者对大气所铁塔附近Z0的估算结果，根据Liu and Sun(2010)对南京城市地区Z0的代表范围约为500 m的研究结果，本文也以大气所铁塔附近500 m估算的Z0和表 2中的前人研究成果进行比较，即根据大气所铁塔经纬度位置在图像上选定该像元位置，以该像元为中心，对比分析周围5个像元(像元空间分辨率为100 m)距离范围内的 。由表 2可知，虽然各位学者得到的结果差异比较大，但综合各种结果，2003年以后大气所铁塔附近Z0应该在4. 0 m以上，而且由于2001年后城市下垫面变化比较小，Z0的变化也不会有明显变化。由图 7可知，在铁塔周围500 m范围内，西南方向为大面积的高层楼群，Z0在4. 0 ~ 9. 0 m之间，西北方向也有大面积的高层楼群，Z0一般在2. 0 m以上，局部达6. 0 m以上，东南方向既有明显＜0. 1 m的Z0低值区，也有较多4. 0 ~ 9. 0 m的Z0高值区，而东北方向大部为0. 1 ~4. 0 m的Z0中低值区。这与徐阳阳等(2009)得到的2003年铁塔Z0为4. 59 m和李倩等(1999)得出的1999年铁塔Z0西南方向为5. 31 m、西北为2. 75 m、东南为2. 99 m、东北为2. 46 m比较接近，表明本文对城市地区估算的Z0结果比较可信。
|
|
表 2 不同学者对大气所铁塔附近Z0的估计 Table 2 Estimation of Z0 of the iron town by different researchers |
以遥感和GIS技术为基础，以形态学模型为基础，开展了空间区域尺度上空气动力学粗糙度长度Z0的估算方法研究，建立了一种不需考虑复杂的气象观测数据、仅考虑卫星资料和基础地理信息资料相对简便的估算空间区域尺度下垫面Z0方法，并以北京为例分别估算了北京空间区域内主要植被地区和城市地区的Z0，与相关研究结果进行了比较和验证分析，结果显示：2012年北京主要植被地区Z0存在着明显的季节变化，山区森林的Z0一般为0. 6 ~ 1. 3 m，冬季高于夏季；平原农田的Z0一般为0 ~ 0. 4 m，夏季高于冬季，这与前人对类似植被Z0的季节变化特征研究结果具有一致性，该结果表明利用卫星资料可有效估算植被地区粗糙度长度Z0，并能反映其时空变化特征。北京城市地区Z0一般为0. 1 ~ 8. 0 m，空间变化明显大于植被地区，其中Z0高值区(在1. 0 m以上)主要分布于北京四环区域内及四环外大型居民社区、经济开发区和郊区县城中心，而在二环区域内存在一个Z0低值区，显示北京城市下垫面存在空间高度非均匀性；结合对大气所铁塔周围Z0值的比较验证分析表明，本文建立的基于遥感与GIS的技术方法可有效估算空间区域精细化的城市地区粗糙度长度Z0，这为未来估算城市复杂下垫面空气动力学粗糙度长度Z0提供一种有效途径。
本文虽然基于卫星数据与GIS地理信息采用形态学法分别估算了研究区域主要植被地区和城市地区的空气动力粗糙度Z0，但估算的结果还存在诸多不确定性：首先是Z0的地面验证非常困难，至今为止没有测定或估算Z0的统一标准，而通常采用的气象观测法估算的Z0不但与或城市建筑覆盖率有直接关系，还与随时空变化的风速、风向有密切关系，且不同方法计算的结果不一致性和代表周边范围的有效距离也难以确定，这给基于网格空间单元的Z0地面验证带来很大困难，未来开展形态学方法和气象观测法的结合研究是提高地表Z0精度的有效方法；第二是对研究区划分为植被地区和城镇地区分别来估算Z0存在不确定性，一方面是山区植被地区没有考虑地形影响，估算结果可能存在着较大偏差；另一方面城市地区也包含一定面积的植被，这部分植被对城市下垫面Z0的影响并没有考虑，如何考虑地形影响、构建复杂城市下垫面植被与建筑混合的Z0估算模型是未来需要研究的内容；第三是基于卫星资料的植被参数叶面积指数和植被高度反演存在误差，尤其是较粗分辨率的森林植被高度的遥感反演还存在诸多不确定性，而植被高度是下垫面Z0估算的关键因子，未来利用多源卫星资料开展更高空间分辨率的植被高度遥感反演有助于植被地区Z0估算精度的提高；第四是建筑参数的估算对Z0很重要，而且不同参数的计算方法也影响着Z0。本文中建筑群高度是按照该区域内所有建筑物平均高度来估算的，实际上建筑群内的高差大小不一，差别很大，很可能存在平均高度相等的高度较为均匀的建筑群和高差很大的建筑群，但这两者实际对Z0的影响可能存在很大不同，如何更合理地估算该区域建筑高度对Z0的影响还需要进一步研究。还有建筑覆盖率随估算的空间尺度(空间分辨率)变化而变化，从而造成Z0估算的变化。未来建立中国各大城市不同分辨率的建筑信息数据库是研究城市地区Z0的重要措施，美国已于2007年建立覆盖44个主要城市地区空间分辨率分别为250 m和1 km的第二版本全国建筑统计信息共享数据库(NBSD2，the second generation compilation of the National Building Statistics Database)(Burian et al，2007)，里面包含了涵盖建筑覆盖率、建筑高度、建筑截面积指数、高宽比、开阔度等等多个建筑形态因子信息，通过这些信息可以更好地改进城市地表Z0估算模型，从而对城市地表Z0的观测、估算和验证都具有重要意义。
尽管当前利用形态学模型来估算Z0还存在诸多不确定性，但正如Grimmond and Oke(1999)所指出的那样，在开展城市下垫面空气动力学参数研究方面，如何利用城市地理信息和卫星资料，结合边界层观测实验资料进行研究将是今后的一个重要研究趋势。随着高分辨率卫星资料和基础地理信息资料的应用普及和各种参数遥感反演精度的提高，具有空间差异分布的下垫面空气动力学参数的估算精度将逐步会得到改善。
致谢 感谢中国科学院遥感与数字地球研究所杨阿强博士给予的无私帮助。| Borak J S, Jasinski M F, Crago R. 2005. Time series vegetation aerodynamic roughness fields estimated from MODIS observations[J]. Agricultural and Forest Meteorology, 135 (1-4): 252–268. DOI:10.1016/j.agrformet.2005.12.006 | |
| Burian S J,Augustus N,Jeyachandran I,et al. 2007. National building statistics database,version 2[R]. Los Alamos National Laboratory unclassified report. | |
| Grimmond C S B, King T S, Roth M, et al. 1998. Aerodynamically roughness of urban areas:Derived from wind observations[J]. Boun-Layer Meteor, 89 (1): 1–24. DOI:10.1023/A:1001525622213 | |
| Grimmond C S B, Oke T R. 1999. Aerodynamic properties of urban areas derived from analysis of surface form[J]. J Appl Meteor, 38 (9): 1262–1292. DOI:10.1175/1520-0450(1999)038<1262:APOUAD>2.0.CO;2 | |
| Lindroth A. 1993. Aerodynamic and canopy resistance of short-rotation forest in relation to leaf area index and climate[J]. Boun-Layer Meteor, 66 (3): 265–279. DOI:10.1007/BF00705478 | |
| Liu G, Sun J N. 2010. Impact of surface variations on the momentum flux above the urban canopy[J]. Theor Appl Climatol, 101 (3-4): 411–419. DOI:10.1007/s00704-009-0219-5 | |
| Jasinski M F, Crago R. 1999. Estimation of vegetation aerodynamic roughness of natural regions using frontal area density determined from satellite imagery[J]. Agricultural and Forest Meteorology, 94 (1): 65–77. DOI:10.1016/S0168-1923(98)00129-4 | |
| Jasinski M F, Borak J, Crago R. 2005. Bulk surface momentum parameters for satellite-derived vegetation fields[J]. Agricultural and Forest Meteorology, 133 (1-4): 55–68. DOI:10.1016/j.agrformet.2005.07.017 | |
| Raupach M R. 1994. Simplified expressions for vegetation roughness length and zero-plane displacement as a function of canopy height and area index[J]. Bound-Layer Meteor, 71 (1-2): 211–216. DOI:10.1007/BF00709229 | |
| Schaudt K J, Dickinsonb R E. 2000. An approach to deriving roughness length and zero-plane displacement height from satellite data,prototyped with BOREAS data[J]. Agricultural and Forest Meteorology, 104 (2): 143–155. DOI:10.1016/S0168-1923(00)00153-2 | |
| Simard M, Pinto N, Fisher J B, et al. 2011. Mapping forest canopy height globally with spaceborne lidar[J]. J Geophys Res, 116 : G04021. DOI:10.1029/2011JG001708 | |
| Wieringa J. 1992. Updating the davenport roughness classification[J]. Journal of Wind Engineering and Industrial Aerodynamics, 41 (1): 357–368. | |
| Zeng X B, Shaikh M, Dai Y J, et al. 2002. Coupling of the common land model to the NCAR community climate model[J]. J Climate, 15 (14): 1832–1854. DOI:10.1175/1520-0442(2002)015<1832:COTCLM>2.0.CO;2 | |
| 蔡福, 周广胜, 明惠青. 2013. 玉米农田空气动力学参数动态及其与影响因子的关系[J]. 生态学报 , 33 (17): 5339–5352. DOI:10.5846/stxb Cai Fu, Zhou Gguangsheng, Ming Hhuiqing, et al. 2013. Dynamics of aerodynamic parameters over a rainfed maize agroecosystem and their relationships with controlling factors[J]. Acta Ecologica Sinica, 33 (17): 5339–5352. DOI:10.5846/stxb | |
| 陈斌, 徐祥德, 丁裕国, 等. 2010. 地表粗糙度非均匀性对模式湍流通量计算的影响[J]. 高原气象 , 29 (2): 340–348. Chen Bin, Xu Xiangde, Ding Yuguo, et al. 2010. The impact of heterogeneity of land surface roughness length on estimation of turbulent flux in model[J]. Plateau Meteor, 29 (2): 340–348. | |
| 崔国发, 成克武, 申国珍, 等. 2000. 北京喇叭沟门林区森林植被现状及分类[J]. 北京林业大学学报 , 22 (4): 46–51. Cui Guofa, Cheng Kewu, et al. 2000. Status and classification of vegetation in Labagoumen forest region in Beijing[J]. Journal of Beijing Forestry University, 22 (4): 46–51. | |
| 高志球, 卞林根, 逯昌贵, 等. 2002. 城市下垫面空气动力学参数的估算[J]. 气象应用报 , 13 (special): 26–33. Gao Zhiqiu, Bian Lingen, Lu Changgui, et al. 2002. Estimation of aerodynamic parameters in urban areas[J]. J Appl Meteor Sci, 13 (special): 26–33. | |
| 韩素芹, 刘彬贤, 解以扬. 2008. 利用255 m铁塔研究城市化对地面粗糙度的影响[J]. 气象 , 34 (1): 54–58. Han Suqin, Liu Binxian, Xie Yiyang. 2008. Research of urbanization on surface roughness with 255 m meteorological tower[J]. Meteor Mon, 34 (1): 54–58. | |
| 何奇瑾, 周广胜, 周莉, 等. 2007. 盘锦芦苇湿地空气动力学参数动态特征及其影响因素分析[J]. 气象与环境学报 , 23 (4): 7–12. He Qijin, Zhou Guangsheng, Zhou Li, et al. 2007. Dynamic characteristics of aerodynamic parameters in Panjin reed wetland and their controlling factors[J]. J Meteor Environ, 23 (4): 7–12. | |
| 胡张保, 俞炳丰. 2008. 城市下垫面空气动力学参数确定方法综述[J]. 气象与环境学报 , 24 (5): 55–60. Hu Zhangbao, Yu Bingfeng. 2008. Review on methods calculating aerodynamic parameters over urban underlying surface[J]. J Meteor Environ, 24 (5): 55–60. | |
| 李倩, 刘辉志, 胡非, 等. 2003. 城市下垫面空气动力学参数的确定[J]. 气候与环境研究 , 8 (4): 443–450. Li Qian, Liu Huizhi, Hu Fei, et al. 2003. The determination of the aerodynamical parameters over urban land surface[J]. Climatic Environ Res, 8 (4): 443–450. | |
| 李沁怡, 蔡旭晖, 宋宇. 2014. 中国高分辨率地表粗糙度分布研究[J]. 高原气象 , 33 (2): 474–482. DOI:10.7522/j.issn.1000-0534.2012.00191 Li Qinyi, Cai Xuhui, Song Yu. 2014. Research of the distribution of national scale surface roughness length with high resolution in China[J]. Plateau Meteor, 33 (2): 474–482. DOI:10.7522/j.issn.1000-0534.2012.00191 | |
| 刘和平, 朱廷曜. 1997. 森林冠层空气动力学参数的确定[J]. 北京大学学报(自然科学版) , 33 (4): 522–528. Liu Heping, Zhu Tingyao. 1997. Determination of aerodynamic parameters of Changbai Mountain forest[J]. Acta Scientiarum Universitatis Pekinenesis, 33 (4): 522–528. | |
| 陆龙骅, 卞林根, 程彦杰, 等. 2002. 冬季北京城市近地层的气象特征[J]. 应用气象学报 , 13 (special): 34–42. Lu Longhua, Bian Lingen, Cheng Yanjie, et al. 2002. Meteorological characteristics of the ground layer in Beijing in winter[J]. J Appl Meteor Sci, 13 (special): 34–42. | |
| 梅凡民, 蒋缠文, 江姗姗, 等. 2002. 粗糙元几何参数的交互作用对床面空气动力学粗糙度的影响[J]. 中国沙漠 , 32 (6): 1534–1541. Mei Fanming, Jiang Chanwen, Jiang Shanshan, et al. 2002. The integrated effects of roughness elements geometric parameters on aerodynamic roughness length[J]. J Desert Res, 32 (6): 1534–1541. | |
| 覃文汉. 1994. 应用压力中心法确定农田空气动力参数[J]. 气象学报 , 52 (1): 99–106. Qin Wenhan. 1994. Aerodynamic parameters of crop canopies estimated with a center-of-pressure technique[J]. Acta Meteor Sinica, 52 (1): 99–106. | |
| 石雪峰, 夏建新, 吉祖稳, 等. 2006. 空气动力学粗糙度与植被特征关系的研究进展[J]. 中央民族大学学报(自然科学版) , 15 (3): 218–225. Shi Xuefeng, Xia Jjianxin, Ji Zuwen, et al. 2006. A review:the study in relationship between aerodynamic roughness length and the characteristics of vegetation[J]. Journal of the CUN(Natural Sciences Edition), 15 (3): 218–225. | |
| 宋秀瑜. 2010. 密云水库上游水源涵养林主要乔木结构特征研究-以潮关西沟林场为例[J]. 中国水土保持 (4): 54–57. Song Xiu Yu. 2010. Study of main arbor construction characteristic of water conservation forest of Miyun reservoir upstream-a case of Chaoguanxigou tree farm[J]. China Soil and Water Conservation (4): 54–57. | |
| 徐阳阳, 刘树华, 胡非, 等. 2009. 北京城市化发展对大气边界层特性的影响[J]. 大气科学 , 33 (4): 859–867. Xu Yangyang, Liu Shuhua, Hu Fei, et al. 2009. Influence of Beijing urbanization on the characteristics of atmospheric boundary layer[J]. Chinese J Atmos Sci, 33 (4): 859–867. | |
| 杨阿强, 孙国清, 卢立新, 等. 2011. 基于MODIS 资料的中国东部时间序列空气动力学粗糙度和零平面位移高度估算[J]. 气象科学 , 31 (4): 516–524. Yang Aqiang, Sun Guoqing, Lu Lixin, et al. 2011. Deriving aerodynamic roughness length and zero-plane displacement height from MODIS product for Eastern China[J]. J Meteor Sci, 31 (4): 516–524. | |
| 杨耀先, 李茂善, 胡泽勇, 等. 2014. 藏北高原高寒草甸地表粗糙度对地气通量的影响[J]. 高原气象 , 33 (3): 626–636. DOI:10.7522/j.issn.1000-0534.2013.00199 Yang Yaoxian, Li Maoshan, Hu Zeyong, et al. 2014. Influence of surface roughness on surface-air fluxes in alpine meadow over the Northern Qinghai-Xizang Plateau[J]. Plateau Meteor, 33 (3): 626–636. DOI:10.7522/j.issn.1000-0534.2013.00199 | |
| 张强, 吕世华. 2003. 城市表面粗糙度长度的确定[J]. 高原气象 , 22 (1): 24–32. Zhang Qiang, Lü Shihua. 2003. The determination of roughness length over city surface[J]. Plateau Meteor, 22 (1): 24–32. | |
| 张学仕, 蒋琰, 薛建辉, 等. 2010. 江苏句容下蜀次生栎林的空气动力学参数研究[J]. 南京林业大学学报(自然科学版) , 34 (6): 61–65. Zhang Xueshi, Jiang Yan, Xue Jianhui, et al. 2010. Aerodynamic parameters of secondary oak forest in Xiashu Jurong county of Jiangsu province[J]. Journal of Nanjing Forestry University(Natural Science Edition), 34 (6): 61–65. | |
| 钟中, 韩士杰. 2002. 长白山阔叶红松林冠层空气动力学参数的计算[J]. 南京大学学报(自然科学) , 38 (4): 565–571. Zhong Zhong, Han Sshijie. 2002. The calculation of aerodynamic parameters of Korea Pine canopy in Changbai Mountains[J]. Journal of Nanjing University(Natural Sciences Edition), 38 (4): 565–571. | |
| 赵晓松, 关德新, 吴家兵, 等. 2004. 长白山阔叶红松林的零平面位移和粗糙度[J]. 生态学杂志 , 23 (5): 84–88. Zhao Xiaosong, Guan Dexin, Wu Jiabin, et al. 2004. Zero-plane displacement and roughness length of the mixed forest of broad-leaved and Korean-pine in Changbai Mountain[J]. Chinese Journal of Ecology, 23 (5): 84–88. | |
| 周艳莲, 孙晓敏, 朱治林, 等. 2006. 几种不同下垫面地表粗糙度动态变化及其对通量机理模型模拟的影响[J]. 中国科学D辑 , 36 (suppl): 244–254. Zhou Yanlian, Sun Xiaomin, Zhu Zhilin, et al. 2006. Dynamic change of some different surface roughness length and influence on simulation of flux mechanism model[J]. Science in China(Series D), 36 (suppl): 244–254. | |
| 周志强, 岳彩荣, 徐天蜀, 等. 2012. 森林高度遥感估测研究综述[J]. 现代农业科技 (2): 198–199. Zhou Zhiqiang, Yue Cairong, Xu Tianshu, et al. 2012. A review of estimation of remote sensing for forest height[J]. Modern Agriculture Science and Technology (2): 198–199. | |
| 王金星, 卞林根, 高志球, 等. 2002. 城市边界层湍流和下垫面空气动力学参数观测研究[J]. 气象科技 , 30 (2): 65–72. Wang Jinxing, Bian Lingen, Gao Zhiqiu, et al. 2002. Observation Study of urban boundary layer turbulence and aerodynamic parameters of underlying surface[J]. Meteor Sci Technol, 30 (2): 65–72. | |
| 王修信, 孙涛, 朱启疆, 等. 2014. 林地叶面积指数遥感估算方法适用分析[J]. 生态学报 , 34 (16): 4612–4619. Wang Xiuxin, Sun Tao, Zhu Qijiang, et al. 2014. Assessment of different methods for estimating forest leaf area index from remote sensing data[J]. Acta Ecologica Sinica, 34 (16): 4612–4619. | |
| 王玉玉, 姚济敏, 韩海东, 等. 2013. 科其喀尔冰川表碛区空气动力学粗糙度分析[J]. 高原气象 , 33 (3): 762–768. DOI:10.7522/j.issn.1000-0534.2013.00140 Wang Yuyu, Yao Jimin, Han Haidong, et al. 2013. Analysis of aerodynamic roughness of the debris-covered Keqicar Glacier[J]. Plateau Meteor, 33 (3): 762–768. DOI:10.7522/j.issn.1000-0534.2013.00140 |
2. Institute of Beijing Urban Meteorology, China Meteorological Administration, Beijing 100089, China