近30年河北大海陀自然保护区山地草甸植被(NDVI)变化及其对气候的响应
2018-05-16余琦殷邢韶华武占军李永霞任志和刘永胜
宋 超,余琦殷,邢韶华,*,武占军,李永霞,任志和,刘永胜
1 北京林业大学自然保护区学院,北京 100083 2 国家林业局经济发展研究中心,北京 100714 3 河北大海陀自然保护区管理处,赤城 075599
全球气候变暖引起了全球平均气温上升、降雨格局的改变,这些环境因子的变化又以不同程度不同方式影响着地上植物群落结构、物候变化以及地表和土壤微环境[1]。植被是陆地生态系统的主体之一,其生长状况与其所处环境密切相关,如气候、水质和地形等环境因子。气候变化改变了植物的生长环境,从而影响到了植被的生长状况[2]。国内外学者已开展了大量关于植被动态变化的研究,其研究方法主要为以下两种：一、传统野外调查与观测;二、遥感植被监测。传统野外调查方法有利于获取多种植被类型的详细生长数据,其研究内容较为丰富,但调查所需时间长且费力,造成经济成本较高,且地理范围和时间序列有限,不同调查人员收集的数据类型和可靠性存在差异,因此传统野外调查方法很难应用于监测大区域尺度上的植被动态变化[3]。自20世纪70年代,卫星遥感技术的迅速发展,为解决上述问题提供了良好的契机[3- 4]。遥感影像具有覆盖面积大,实时性、周期性、准确可靠性等优点,能够从景观尺度上来定量监测土地利用的变化,为解决当前资源与环境之间的问题提供了有力的技术支撑,已被广泛应用于大区域尺度的植被动态变化研究中[5]。国内外学者对植被与气候变化之间的关系研究较多,时间尺度为目前国际常用的时间序列,主要基于归一化差值植被指数(NDVI)数据研究不同植被类型对气候变化的响应[6- 7]。
草甸是以多年生、中生草本植物为主体的群落类型[8]。在草甸植被的研究中,国内学者通常以NDVI为表征,探讨草甸植被的变化趋势及其对气候变化的响应机制,以掌握不同地区草甸植被的演变趋势[9- 11],研究区域主要在我国西部典型草原草甸分布区。河北大海陀自然保护区山地草甸与我国西部草甸地区相比每年6—8月平均降水量和平均气温相对较高,在这种水分和温度相对充足的情况下,草甸植被的变化及其对气候变化的响应如何尚没有太多研究与报道。本研究以草甸植被的NDVI为特征指标,研究河北大海陀自然保护区山地草甸植被变化以及其对气候变化的响应,为进一步探讨我国东部地区山地草甸植被变化其对气候变化的响应机制提供支撑。
1 研究区域
河北大海陀国家级自然保护区地处河北省赤城县西南部,其地理坐标为115°42′57″—115°57′0″E, 40°32′14″—40°41′40″N,总面积12634hm2[12],主峰为海陀山,海拔2241m[13],是河北赤城县与北京市的分界岭[6]。
河北大海陀自然保护区地处暖温带的中山地区,属于东亚大陆性季风气候区,受地理位置影响,该区形成的地带性植被为暖温带落叶阔叶林。山地草甸分布于2000m以上,主要成分为多年生杂草,主要包括薹草属(Carex)、早熟禾属(Poa)、罂粟属(Papaver)、风毛菊属(Saussurea)、金莲花属(Trollius)、银莲花属(Anemone)、马先蒿属(Pedicularis)、梅花草属(Parnassia)等[6,14]。
2 数据来源与研究方法
2.1 数据来源
数据主要包括大海陀自然保护区近30年的TM遥感影像数据;赤城县1980—2015年的气候数据,包括年平均气温、7月份平均气温、1月份平均气温、年降水量、6—8月平均降水量等多个气候指标;大海陀自然保护区及周边地区的地形数据等[6]。
自然保护区遥感影像资料从地理空间数据云(www.gscloud.cn)中下载获得,使用LANDSAT系列数据,该系列数据每16天循环扫描一次[6]。具体选择的遥感数据见表1。
地形数据来地理空间数据云的GDEMDEM- 30M分辨率数字高程数据[6]。
近36年的温度和降水数据均采用国家气象局发布的数据,温度和降水数据取自1980—2015年逐年的1km2栅格数据[6]。
表1 遥感数据基本信息
2.2 数据处理
2.2.1 归一化差值植被指数的计算
归一化差值植被指数(Normalized Difference Vegetation Index,NDVI)是反映植被覆盖的一个重要指数,其时间序列的变化对应着植被的生长和变化[15]。NDVI是遥感估算植被覆盖度研究中最常用的植被指数,为植物生物生长状态及植被空间分布密度的最佳指示因子,与植被分布密度呈线性相关,长期用作监测植被变化情况[16]。NDVI对土壤背景的变化较为敏感,检测植被覆盖度的幅度较宽,当植被覆盖度小于15%时,土壤背景与植被可以区分开;当植被覆盖度为25%—80%时,NDVI值随植被覆盖度的增大呈线性增加;当植被覆盖度大于80%时,NDVI的检测能力逐步下降[6,17]。
根据实地调查结果,草甸主要位于大海陀自然保护区海拔2000m以上的地区,通过GDEMDEM- 30M 分辨率数字高程数据提取该范围矢量文件,并裁剪计算平均NDVI值[6]。
(1)
式中,R指红色波段,NIR指红外波段,其中在Landsat4- 5 TM分别为Band3、Band4;在Landsat 8 OLI_TRIS分别为Band4、Band5。由于季节对NDVI指数的影响,各年度选用相同的月份用于NDVI指数计算,根据遥感影像云量、草甸长势各月份间长势情况,选取8月底至9月下旬的遥感数据,共计12a[6]。
统计大海陀自然保护区海拔2000m以上地区的平均NDVI值,采用Office Excel与SPSS 20.0统计分析软件进行数据分析,然后利用滑动平均法对数据进行处理,得到NDVI值的变化趋势。NDVI取值范围为-1—1,其中NDVI≤0说明无植被覆盖,NDVI>0时有植被覆盖,且值越大,植被覆盖度越高,植被生长状况越好[6]。
2.2.2 气候变化数据整理与分析
运用ArcGIS将收集到的温度和降水数据转换成面数据,采用Office Excel与SPSS 20.0统计分析软件进行单年区域加权平均和多年滑动平均处理,得出温度和降水量的变化趋势;运用基于秩的非参数Mann-Kendall统计检验方法(以下简称M-K检验法)计算温度和降水量变化趋势是否明显。
(1)滑动平均法
滑动平均法能够精确地表示测量结果,抑制随机误差产生的影响。滑动平均法即是对非平稳的数据,在适当的小区间上视为接近平稳,从而计算局部的平均,以减小随机起伏误差,以此类推,沿全长N个数据逐一在小区间上进行不断的局部平局,即可得出较平滑的结果[18]。滑动平均法常被用于气候数据的分析,根据前人研究经验[19- 20],本文选取5年滑动平均法[6]。
滑动平均法为沿全长N个数据,不断逐一滑动取m个相邻数据进行加权平均运算,得出平滑数据,其一般算式为：
(2)
运用统计学方法一元线性方程定量描述气候因子的回归变化趋势,即以年代为时间因子,计算保护区的年均温、1月平均气温、7月平均气温、年降水量、6—8月平均降水量的滑动平均数据的线性趋势。以年份为自变量,要素(气温、降水量)分别为因变量,建立一元回归方程对要素和年份数据进行线性回归分析,计算直线的拟合斜率b和相关系数R值。用xi为第i年要素(气温、降水量),用Ti表示xi所对应的时间,建立一元线性回归模型：
xi=a+bti(i=1,2,…,n)
(3)
(4)
(5)
式中,a为截距,b为斜率(变化趋势),a和b运用最小二乘法估计。当斜率b>0时,说明该要素变化随年度呈上升趋势;若b<0时,说明该要素变化随年度呈下降趋势[6]。
(2)M-K检验法
该方法常用于水质、径流量、温度、降水等时间序列变化趋势的显著性检验,主要是通过计算统计量τ,方差σt2和标准化变量M,来判断序列趋势是否显著[21- 24]：
(6)
(7)
(8)
式中,S为序列中数值间i
2.2.3 影响草甸变化的因素分析
地理要素之间做相关分析主要目的是揭示各要素之间的密切程度,而这种密切程度的测定,则要通过计算它们之间的相关系数并检验来进行[6]。
(1)相关性分析
相关分析是考量两个连续变量之间的相关关系,实际对任何类型的变量,都可以运用相关的指标分析两者之间的相关关系,它能以数字方式准确地描述变量间的线性相关程度[6]。
Spearman秩相关系数是基于各因子间秩的大小,属于非参数检验,其表达式为：
(9)
式中,rs(i,j)为Spearman秩相关系数;n为因子总年数;dk=(xik-xjk),xik和xjk分别为因子i和因子j在k年份中的重要值[6]。
(2)偏相关分析
偏相关系数是在对其他变量影响进行控制的条件下,衡量多个变量中某两个变量之间的线性相关程度的指标。通过偏相关系数与相关系数的比较,来确定两个变量之间的内在线性联系会更真实、更可靠[25]。偏相关分析的计算方法是：先计算3个因子间的相关系数,然后通过3个简单相关系数来计算偏相关系数。其公式为：
(10)
式中，rij(k)就是在控制了k因子影响所计算的因子i和因子j之间的偏相关系数;rij、rik、rjk分别是i、j、k3个因子之间的两两简单相关系数[6]。
3 结果分析
3.1 草甸植被NDVI变化
图1 大海陀自然保护区草甸NDVI指数Fig.1 The NDVI index of Dahaituo Nature Reserve meadow
讨论NDVI变化特征时,生长值的峰值作为常用的特征量之一,能很好地反映植被生长发育程度和未来发展趋势等[26],故选取8月底至9月下旬植被生长状况较好时段的遥感影像进行分析[6](表1)。
由大海陀自然保护区1987—2015年间选取的12年的草甸植被NDVI指数变化趋势图(图1)可以看出,草甸植被NDVI在1987—2004年呈现小幅持续增加的趋势,在2004—2015年间呈下降趋势,NDVI值在2004年达到最高值0.58,随后逐渐下降,在2013年达到了最低值0.23[6]。
3.2 区域气候变化
3.2.1 降水量变化
为了获取研究期(1980—2015年)内大海陀自然保护区的降水量变化,采用统计学方法计算了该时期的年降水量及6—8月平均降水量,结果如图2、图3所示[6]。
图2 1980—2015年大海陀自然保护区6—8月平均降水量变化趋势Fig.2 The change trend of average precipitation in June to August in Dahaituo Nature Reserve meadow during 1980 to 2015
图3 1980—2015年大海陀自然保护区年降水量变化趋势Fig.3 The change trend of annual precipitation in Dahaituo Nature Reserve meadow during 1980 to 2015
从图2、图3中可以看出,大海陀自然保护区草甸区域近36年的年降水量变化整体表现为波动循环的趋势,且总体略有上升,其间伴随有小幅波动。其中年降水量2007年最低,为619.48mm;2000年最高,达1356.35mm;6—8月平均降水量2007年最低,为102.46mm;1997年最高,达230.79mm[6]。
从年代变化来看,1980—1991年降水量平均为915.04mm,1992—2003年降水量平均为950.92mm,2004—2015年降水量平均为912.42mm,长期变化的趋势并不明显[6]。
经M-K检验法检验,M(年降水量)=0.74<1.94,M(6—8月平均降水量)=1.20<1.94,均未达到显著水平,也说明了大海陀自然保护区年降水量和6—8月平均降水量在近36年的时间序列上变化不显著[6]。
3.2.2 气温变化
(1)年均温变化趋势
1980—2015年大海陀自然保护区草甸气温变化趋势如图4所示。近36年该地区的平均气温为6.47℃,年平均增温速率为0.55℃/10a。从年代变化来看,1980—1991年均温为5.92℃,1992—2003年均温为6.14℃,2004—2015年均温为7.34℃;其中虽然在1980—1991年间有阶段性下降的趋势,但从总体来看仍呈阶段性上升的趋势,期间的最高年均温为7.95℃,出现在2007年;最低年均温为5.02℃,出现在1980年[6]。
经M-K检验法检验,M(年均温)=4.54>2.58,到达了极显著水平(P<0.01),说明了该地区的气温上升趋势极显著,与全球变暖的趋势较为一致[6]。
图4 1980—2015年大海陀自然保护区年均温变化趋势Fig.4 The change trend of average annual temperature in Dahaituo Nature Reserve meadow during 1980 to 2015
(2) 1月、7月平均气温变化趋势
图5、图6为大海陀自然保护区草甸1月、7月的平均气温变化趋势。该地区1月、7月的月平均气温均呈现出了上升的趋势,1月平均气温的上升趋势明显高于7月平均气温的上升趋势,1月的气温增温速率为0.76℃/10a;7月气温增温速率为0.27℃/10a[6]。
经M-K检验法检验,M(1月均温)=2.10>1.94,M(7月均温)=0.65<1.94,其中1月平均气温的变化到达了显著水平(P<0.05),说明了该地区的气温上升趋势主要由以1月为主的冬季升高而引起[6]。
图5 1980—2015年大海陀自然保护区1月份平均气温变化趋势Fig.5 The change trend of average temperature in January in Dahaituo Nature Reserve meadow during 1980 to 2015
图6 1980—2015年大海陀自然保护区7月份平均气温变化趋势Fig.6 The change trend of average temperature in July in Dahaituo Nature Reserve meadow during 1980 to 2015
3.3 草甸植被NDVI与气候变化的关系
3.3.1 草甸植被NDVI与气候变化的相关分析
分别将草甸植被NDVI值与降水量、温度等环境因子作相关分析,结果见表2。
表2 草甸植被NDVI与环境因子的相关性分析
*P<0.05;**P<0.01
在Spearman秩相关分析中可以看出,草甸植被NDVI与温度均呈负相关,且与年均温及1月均温达到了显著相关水平(P<0.05);且相比之下,1月均温的显著性指数达到了-0.769。说明当年的年均温越低,草甸植被NDVI指数越高,草甸植被生长状况越好,且主要受1月均温的影响,受7月均温的影响较小[6]。
草甸植被NDVI与年降水量、6—8月平均降水量均呈正相关,但均不显著。
3.3.2 草甸植被NDVI与气候变化的偏相关分析
分别将代表草甸植被NDVI与降水量、年均温等环境因子的变化作偏相关分析,分别控制变量年均温和年降水量,对草甸植被NDVI变化进行分析(表3),可以看出,草甸植被NDVI与年均温的偏相关系数为-0.738,通过P<0.01的显著性水平检验,且负相关程度大于未作偏相关分析的相关系数,说明排除了年降水量影响后,年均温与草甸植被NDVI的相关程度更加明显;草甸植被NDVI与年降水量变化的偏相关系数为-0.249,未通过显著性水平检验,说明草甸植被NDVI与同期降水量属于不显著的负相关[6]。
表3 草甸植被NDVI与环境因子偏相关分析
4 结论与讨论
4.1 结论
(1)大海陀自然保护区草甸植被NDVI值呈先上升后下降趋势,在2004年达到最高值0.58,在2013年降到最低值0.23。
(2)大海陀自然保护区草甸区域的降水量、6—8月平均降水量呈波动循坏趋势,总体略有上升,但均未达到显著水平。该地区年均降水量为915.04mm,在2000年年降水量达到最高值1356.35mm,在2007年年降水量降到最低值201.46mm。
通过Spearman秩相关分析,草甸植被NDVI与年降水量、6—8月平均降水量呈正相关,但均不显著。通过偏相关分析,草甸植被NDVI与年降水量变化的偏相关系数为-0.249,未通过显著性检验,为不显著负相关关系。
(3)大海陀自然保护区草甸区域的年均温、1月平均气温、7月平均气温均呈上升趋势,其中年均温达到极显著水平,1月平均气温达到显著水平。该地区年均温为6.47℃,在2007年,年均温达到最高值7.95℃,在1980年年均温达到最低值5.02℃。
通过Spearman秩相关分析,草甸植被NDVI与温度均呈负相关,且与年均温及1月平均气温为显著负相关。说明当年的年均温越低,草甸植被NDVI指数越高,草甸植被生长状况越好,且主要受1月均温的影响,受7月均温的影响较小。通过偏相关分析,草甸植被NDVI与年均温的偏相关系数为-0.738,通过显著性检验,说明排除了年降水量影响后,年均温与草甸植被NDVI的相关程度更加明显。
4.2 讨论
在山地自然生态系统中,地形因子能影响温度和降水的再分配,一般而言,随着海拔高度的上升,气温逐渐降低,降水量和相对湿度呈增加趋势。植被生长在低海拔地区主要受到水分条件的限制,降水量的增加将促进植被的生长,在高海拔区域,植被生长会受到热量条件的限制,温度上升,或者温度上升和降水量增加共同改善了植被的生长条件,有利于植被生长[27- 28]。根据以上研究中对草甸所在区域海拔的划分,本研究中大海陀自然保护区草甸分布区属于低海拔地区,近30年来,降水没有显著变化,草甸植被NDVI与年降水量也没有显著的相关性,但却与年均温表现为负相关性,在降水量不变的情况下,仅温度升高,反而对草甸植被的生长具有抑制作用,即年均温升高,NDVI指数降低。有学者在对我国东北多年冻土区植被 NDVI的研究后认为,伴随着气温的显著升高和降水量下降,植被 NDVI 显著下降,降水量是影响植被 NDVI 的主要因子[15],这与本研究得出的结果类似[6]。
姚玉璧等人在对甘肃玛曲县亚高山草甸研究后,发现山地草甸类草地植被从返青到分蘖期,气温对草地牧草产量形成正效应,拔节到抽穗期牧草产量对温度变化敏感,开花到籽实成熟期,热量对牧草产量的影响由正效应转向负效应,在籽实成熟期牧草产量形成对气温变化进入第二个敏感期;除籽实成熟期降水量对牧草产量形成为负相关外,其余时段降水量对牧草产量形成均为正效应[29],本研究选取8月底至9月下旬草甸遥感影像资料,根据姚玉璧的草甸植被生长期分类,大海陀自然保护区草甸在8月底至9月下旬为籽实期至黄枯期,温度对牧草产量影响为负效应,与本研究得出的结论一致[6]。
在对我国东部草甸所在地区的植被覆盖变化研究中,张建亮等人对长白山地区的的研究中采用了7—8月份NDVI,发现植被变化受气温的影响程度均强于降雨,温度的波动对植被的影响更为敏感[30- 31]。黄晓霞等人对小五台亚高山草甸的研究中,从景观格局的角度研究发现,制约草甸群落分布的最为重要的条件是热量状况,其次是水分条件的影响[32]。这些与本研究得出的结论一致。
大海陀自然保护区草甸的植被变化除了受降雨量和温度等气象因子影响外,与旅游活动等人为因素有一定关系。随着当地政府对保护区的保护力度加大,草甸植被逐步改善,大海陀自然保护区于2003年被评为国家级自然保护区。同时,当地政府加大了保护区的投资开发力度,大海陀逐渐成为北京周边旅游的热门地区之一,草甸成为游人活动的集中区域[33]。但由于地处相对偏僻,社会经济发展比较落后,对旅游的管理落后,缺乏科学规划和专业人员,可能成为导致草甸植被生长状况出现逐步下降的因素之一。
参考文献(References):
[1] 宋飘, 张乃莉, 马克平, 郭继勋. 全球气候变暖对凋落物分解的影响. 生态学报, 2014, 34(6): 1327- 1339.
[2] 杨波. 秦巴山间地区植被覆盖变化及其对气候变化的响应[D]. 兰州: 西北师范大学, 2014.
[3] Pettorelli N, Vik J O, Mysterud A, Gaillard J M, Tucker C J, Stenseth N C. Using the satellite-derived NDVI to assess ecological responses to environmental change. Trends in Ecology & Evolution, 2005, 20(9): 503- 510.
[4] Kerr J T, Ostrovsky M. From space to species: ecological applications for remote sensing. Trends in Ecology & Evolution, 2003, 18(6): 299- 305.
[5] 潘竟虎, 刘菊玲, 王建. 基于遥感与GIS的江河源区土地利用动态变化研究. 干旱区地理, 2004, 27(3): 419- 425.
[6] 宋超. 河北大海陀自然保护区山地草甸植被变化及影响因素研究[D].北京: 北京林业大学,2016.
[7] 蒋智. 经济发展对植被覆盖变化的影响分析——以重庆市北碚区为例[D]. 重庆: 西南大学, 2010.
[8] 吴征镒. 中国植被. 北京: 科学出版社, 1980.
[9] 李辉霞, 刘国华, 傅伯杰. 基于NDVI的三江源地区植被生长对气候变化和人类活动的响应研究. 生态学报, 2011, 31(19): 5495- 5504.
[10] 张琪, 袁秀亮, 陈曦, 罗格平, 李龙辉. 1982—2012年中亚植被变化及其对气候变化的响应. 植物生态学报, 2016, 40(1): 13- 23.
[11] 顾润源, 周伟灿, 白美兰, 李喜仓, 邸瑞琦, 杨晶. 气候变化对内蒙古草原典型植物物候的影响. 生态学报, 2012, 32(3): 767- 776.
[12] 中华人民共和国环境保护部. 河北大海陀等28处国家级自然保护区的面积和范围. (2013-07- 17) [2014-04- 20]. http://www.mep.gov.cn/gkml/hbb/bh/201307/t20130722_256012.htm.
[13] 白凤鸣. 海陀纪事. 北京: 化学工业出版社, 2011.
[14] 蒋瑞海, 王国荣. 海陀大观. 石家庄: 河北科学技术出版社, 2006.
[15] 毛德华, 王宗明, 宋开山, 刘殿伟, 张柏, 张素梅, 罗玲, 张春华. 东北多年冻土区植被NDVI变化及其对气候变化和土地覆被变化的响应. 中国环境科学, 2011, 31(2): 283- 292.
[16] 李苗苗. 植被覆盖度的遥感估算方法研究[D]. 北京: 中国科学院研究生院(遥感应用研究所), 2003.
[17] 胡良军, 邵明安. 论水土流失研究中的植被覆盖度量指标. 西北林学院学报, 2001, 16(1): 40- 43.
[18] 裴益轩, 郭民. 滑动平均法的基本原理及应用. 火炮发射与控制学报, 2001, (1): 21- 23.
[19] 袁宏霞. 内蒙古东部农牧交错区近60年的气候变化研究[D]. 呼和浩特: 内蒙古师范大学, 2015.
[20] 阿依仙木古丽·阿布来提. 新疆喀什市近30年气候变化及其对农业影响的研究[D]. 乌鲁木齐: 新疆师范大学, 2013.
[21] 徐丽梅, 郭英, 刘敏, 胡乔利, 陈福军, 杨树聪, 沈彦俊. 1957年至2008年海河流域气温变化趋势和突变分析. 资源科学, 2011, 33(5): 995- 1001.
[22] 张峰, 雷晓辉, 蒋云钟, 白建锋. 珠江上游流域近55年气温变化特征分析. 水资源与水工程学报, 2012, 23(2): 20- 25.
[23] 贺伟, 布仁仓, 熊在平, 胡远满. 1961—2005年东北地区气温和降水变化趋势. 生态学报, 2013, 33(2): 519- 531.
[24] 章诞武, 丛振涛, 倪广恒. 基于中国气象资料的趋势检验方法对比分析. 水科学进展, 2013, 24(4): 490- 496.
[25] 张远东, 张笑鹤, 刘世荣. 西南地区不同植被类型归一化植被指数与气候因子的相关分析. 应用生态学报, 2011, 22(2): 323- 330.
[26] 李娜. 1999- 2006年石羊河流域植被对气候变化的响应研究[D]. 兰州: 兰州大学, 2010.
[27] Zeng B, Yang T B. Impacts of climate warming on vegetation in Qaidam area from 1990 to 2003. Environmental Monitoring and Assessment, 2008, 144(1/3): 403- 417.
[28] 邓少福. 祁连山气候变化对植被的影响研究(2000- 2011)[D]. 兰州: 兰州大学, 2013.
[29] 姚玉璧, 张秀云, 段永良. 气候变化对亚高山草甸类草地牧草生长发育的影响. 资源科学, 2008, 30(12): 1839- 1845.
[30] 南颖, 刘志锋, 董叶辉, 李秀霞, 吉喆. 2000—2008年长白山地区植被覆盖变化对气候的响应研究. 地理科学, 2010, 30(6): 921- 928.
[31] 张建亮, 刘方正, 崔国发. 长白山国家级自然保护区植被时空变化及其驱动因子. 生态学报, 2016, 36(12): 3525- 3536.
[32] 黄晓霞, 江源, 刘全儒, 黄秋如. 小五台亚高山草甸与生境关系分析. 植物生态学报, 2007, 31(3): 437- 444.
[33] 刘鹏. 河北大海陀国家级自然保护区生态旅游开发研究[D]. 石家庄: 河北师范大学, 2009.
