2. 中国科学院地理科学与资源研究所陆地表层格局与 模拟院重点实验室, 北京 100101
2. Key Laboratory of Land Surface Pattern and Simulation, Institute of Geographical Sciences and Natural Resources, Chinese Academy of Sciences, Beijing 100101, China
随着人口增长和经济发展,全球水资源供需矛盾十分突出,在干旱和半干旱地区水资源短缺已成为制约经济社会发展的瓶颈。在一定的自然地理背景和时空尺度下,人类满足自身需求以及维持整个社会经济活动而对其赖以生存和发展的水资源和水生态环境产生的影响和冲击,简称水资源压力。水资源压力的大小受到自然条件、人口规模、生活质量、经济总量、经济结构和技术条件等多重因素的影响[1]。1992年Falkenmark[2]定义人均水资源量为水资源压力指数。联合国以世界粮农组织用水资源开发程度来衡量水资源压力。不少学者[3-7]对水资源压力进行研究,但存在很多不足。第一是忽略了生态用水压力以及动态变化:以前的研究主要从人口和经济的角度分析区域水资源的承载能力,而对生态环境需水与水资源利用的动态关系重视不足,不能客观地反映水资源是否能够维护生态环境安全的需求。第二,从研究区域而言,主要集中在东中部及城市,西北内陆地区很少涉及。第三是研究尺度过大,主要集中在全国和省级尺度上,缺乏对县域水资源压力的研究。流域内县际气候条件、土壤质地、植被结构等自然要素以及经济社会状况都存在差异,特别是上中下游之间。在流域水资源配置和县域水资源管理中迫切需要掌握每个县来自自然和人类的水资源压力,从而因县制宜地制定水资源管理、产业结构调整和有效地指导流域水资源的配置和管理的政策。而目前的区域水资源可持续利用评价方法和指标体系注重多层次、多指标, 县级尺度的相关数据很难获得;同时实际操作存在诸多不便。因此,以前的水资源压力评价方法并不适用于县级尺度。
本文基于SWAT模型,结合统计方法,构建集中反映水资源压力的3个指标评估黑河流域11个县1995—2014年的水资源压力及其分压力的空间格局和时间序列上的动态。先利用SWAT模型获得县域的可利用水资源量和生态用水量,然后结合统计数据进行进一步研究。以期进一步完善水资源压力的理论,探索县域水资源压力研究的新方法, 为黑河流域水资源配置和县域水资源管理提供科学的指导。
1 方法与数据 1.1 研究区域黑河全长821 km, 是中国西北地区第二大内陆河, 发源于南部祁连山区, 流域东与石羊河流域相邻, 西与疏勒河流域相接, 北至内蒙古自治区额济纳旗境内的居延海, 流域范围介于98°~102°E, 37°50′~42°40′N之间, 包括肃州、肃南、祁连、山丹、民乐、临泽、嘉峪关、高台、甘州、金塔和额济纳旗等11个县级行政单位,总面积达14.29万km2。该流域气候干燥,降水稀少而集中, 多大风, 日照充足, 太阳辐射强烈,昼夜温差大。流域地势南高北低,按海拔高度和自然地理特点分为上游祁连山地、中游走廊平原和下游阿拉善高原3个地貌类型区。水资源短缺是该流域的限制性因素,自从20世纪50年代起流域水环境不断恶化,严重影响了流域生态系统的健康、社会经济的发展和人类的生产生活;同时研究区也是流域水资源管理的实验区[8-9]。
1.2 SWAT模型SWAT(soil and water assent tool)模型是美国农业部(USDA)农业研究中心(ARS)1993年为美国水文模型(HUMUS)项目开发的大、中尺度流域环境模型, 已经在美国18个主要流域以及其他国家得到广泛应用[10-12]。SWAT对水文循环过程的模拟包括产流、坡面汇流和河道汇流3部分[13]。径流模拟计算沿用SRRWB模型的水量平衡公式的方法, 计算每个水文子单元水量平衡公式如下
$ \Delta SW = P - Q - Ea - S + GGW, $ | (1) |
式中: ΔSW、P、Q、Ea、S和GGW分别为土壤含水量、降水、地表径流、实际蒸散发、深层下渗和浅层回归流。S可分为产流的下渗(PERC)和汇流的河道下渗(TLOSS),地表径流采用改进的SCS模型的径流曲线数方法计算。土壤表层到根区的剖面分成多层, 各层分别计算土壤水的出入渗、侧流和蒸发。模型采用土壤蓄水演算技术(storage routing technology)计算出土壤每层的下渗量;地下径流侧流每一层的侧流量采用动力蓄水模型(kinematic storage model)计算;潜在蒸散发(PET)提供Penman-Monteith,Priestley-Taylor和Hargreaves方法3种可选估算方案;实际蒸散发(Ea)在计算完PET后计算。植被蒸散发是PET、土壤根区深度和植被叶面积指数(LAI)的函数, 土壤蒸散发是土壤深度、土壤水分和PET的函数;融雪径流的计算基于度日因子方法;汇流采用马斯京根(Muskingum)方法进行演算。通过汇流最终到达主河道的水流为产水量(water yield),产水量是区域各种水分量综合作用后最后的水量结果[14]。
研究中用决定系数(R2)、纳西系数(NS)、标准差比率(RSR)和百分比偏差(precent bias, PBIAS)衡量模型校准期和验证期的模拟效果[15-16]。决定系数衡量模拟值和实测值的拟合效果,范围为0~1,值越大表示模拟效果越好。一般来说,0.65<NS<0.75, 0.50<RSR<0.60且10%<PBIAS<15%,表明效果好;0.50<NS<0.65,0.60<RSR<0.70且15%<PBIAS<25%,模拟结果满足要求[15]。
1.3 数据来源模型需要的数据包括土地利用数据、土壤数据、气象数据和DEM。2008年土地利用数据来自于中国科学院环境数据中心(http://www.resdc.cn/)。日气象数据下载于中国气象数据共享中心(http://cdc.cma.gov.cn/)。土壤数据利用FAO的世界和谐土壤数据(harmonized world soil dasebase,HWSD);在土壤属性数据库构建时利用SPAW(soil plant atmosphere water)软件。分辨率为30 m的DEM数据来自于中国地理数据云(http://www.gscloud.cn/)。模型所需的径流实测数据来自于黑河流域水资源管理局。流域基础地理数据来自于寒区旱区数据中心(http://westdc.westgis.ac.cn/)。
县域可利用水资源总量和生态用水可以通过模型模拟,然后统计获得。生产用水和生活用水需要利用统计数据。但长时间序列县级生产和生活用水数据无法直接获得。而黑河流域生产用水主要包括农业生产用水和工业生产用水。根据青海省、甘肃省、内蒙古自治区年鉴的农业用水量与农业灌溉面积获得每年灌溉定额的基础上,乘以县域的农业灌溉面积可获得县域农业生产用水;根据青海省、甘肃省、内蒙古自治区年鉴工业生产总值和工业用水总量获得工业生产用水定额的基础上,乘以县域工业生产总值可获得县域工业用水量。生活用水分为城市生活用水和农村生活用水,同样利用年鉴的数据获得整个流域的城市生活用水定额和农村生活用水定额,然后分别乘以县域城市人口和农村人口数,进而可获得县域生活用水。
1.4 水资源压力指数构建水资源总压力(TP)分为源于自然系统的压力和人类系统的压力。自然系统的水资源压力主要是生态用水压力(EP)。人类系统的水资源压力主要可分解为生产用水压力(PP)和生活用水压力(DP)。可用公式表示为
$ TP = EP + PP + DP + \mu , $ | (2) |
式中:
$ EP = \frac{{EW}}{{WA}}, $ | (3) |
$ PP = \frac{{PW}}{{WA}}, $ | (4) |
$ DP = \frac{{DW}}{{WA}}. $ | (5) |
式中:WA为可用水资源总量;EW是生态用水量;PW是生产用水量;DW是生活用水量;μ是其他的水资源压力。
1.5 Mann-Kendall检测县域水资源压力和分压力在时间维度的趋势变化利用Mann-Kendall(MK)检测进行研究。MK检验是世界气象组织推荐并广泛使用的非参数检验方法。MK是基于秩的非参数方法, 由于其不要求所分析数据服从某一概率分布, 同时也不受个别异常值的干扰, 能够客观地表征样本序列的整体变化趋势, 因而被广泛应用于气候参数和水文序列的分析中。
2 结果与分析 2.1 模型校准和验证利用黑河流域的扎马什克、莺落峡、新地、祁连、肃南、高崖,正义峡、狼心山等水文站的月实测径流数据对模型进行率定和验证。模型率定期是1995—2007年,验证期为2008—2013年。狼心山水文站由于数据限制,率定期为2004—2009年,验证期为2010—2013年。通过统计R2都在0.7以上,0.65<NS<0.75;0.50<RSR<0.60且10%<PBIAS<15%。这表明模拟结果良好,完全可以用于进一步的研究。
2.2 县域水资源压力的空间格局从表 1可看出甘州的水资源总压力最高,达10.56;这主要归因于生产用水压力巨大,达9.588 9。甘州的生态用水压力也比较大,达0.957 0;但生活用水压力相对其他水资源压力分量要小得多。水资源压力排在第2位的是额济纳,水资源总压力为2.439 2。水资源总压力高的原因和甘州类似,主要得益于生产用水压力(PP=1.967 8)过高;与此同时生态用水压力也较高,但生活用水压力相当小。甘州、额济纳和山丹的水资源总压力超过1,表明用水量超过正常状态下自然系统可以提供的水资源总量,即用水量已经超过水资源的承载能力。甘州区位于河西走廊的中部,是张掖市的经济和行政中心。工业规模的扩大对水资源产生重大压力。但是甘州区的水资源压力主要受农业生产用水的控制,而不是工业用水。农业用水量巨大一方面是水资源利用效率不高,采用漫灌的方式,2014年甘州区的农业灌溉定额高达104 m3/ha[17];另一方面是农业用地面积不断扩大[18],1995—2014年间农业灌溉面积增加1.675万ha。额济纳也是因为同样的原因导致水资源压力过大。这些水资源赤字是通过过度开采地下水来填平的,1995—1999年地下水位缓慢下降;2000—2008年显著下降;井灌区平均下降0.15~0.30 m/a[17]。
山丹、临泽和民乐的水资源总压力次于额济纳,分别为1.285 5、0.996 9和0.947 5;但原因与甘州、额济纳不同,水资源压力主要来自于生态用水,生态用水压力分别达到0.957 0、0.643 2和0.538 3。其他水资源总压力超过0.5的行政区从高到低依次是祁连、肃南和高台。
生态用水压力最大的行政区是山丹县,其次是甘州、祁连、肃南、临泽和民乐;这些区域的生态用水压力已经超过0.5。这主要和植被状况相关,山丹是黑河中游植被状况最好的区域[19],2016年森林覆盖率达到29.33%。这些植被需要耗费较多的水来维持生态系统功能的发挥。
而生产用水量最大的是下游的额济纳;主要是农业用水量大所导致。其次是黑河中游张掖市所属的县,其中甘州区最大。各行政区的生活用水压力小于生态和生产用水压力;其分布格局是黑河上游行政区(祁连和肃南)<下游行政区(金塔和额济纳)<中游行政区(民乐、山丹、甘州、临泽、高台、肃州、嘉峪关)。张掖市集中了黑河全流域95%的耕地、91%的人口和89%的国内生产总值,是流域内水资源的主要利用区[20]。因此水资源的压力不仅来源于农业用水量,还来自于工业生产用水量和生活用水。众多的人口导致生活用水压力也要比其他区域大。而经济的繁荣和广阔的耕地面积导致其为水资源高压区。
从构成总压力的主分量来看,祁连、肃南、民乐、山丹、临泽、高台、肃州、嘉峪关和金塔属于生态用水压力主导型;而甘州和额济纳属于生产用水压力主导型;生态用水对各行政单位的压力都很小,最大值出现在甘州,也不到0.1。
从区域上看,黑河上游区域水资源压力主要来自于生态用水;生产用水压力和生活用水压力很小。中游区域的生产用水压力和生活用水压力明显变大。下游区域的生态用水压力和生产用水压力要远远大于生活用水压力。
2.3 县域水资源压力的变化趋势 2.3.1 生态用水压力从图 1可以看出,祁连、肃南、高台、山丹和额济纳的生态用水压力呈现下降趋势;其中肃南的下降速率最大,年平均下降幅度为0.004 9。这说明在这5个行政区的生态用水量有所减少。而黑河流域其余的县级行政单位生态用水压力在1995—2014年呈现出上升趋势,其中最大上升速率出现在甘州,年平均下降幅度为0.006 6。这表明这6个县级行政区生态用水量在1995—2014年间增加,其中甘州增加速度最快。总体上看,各县级行政单位的变化程度比较微弱,除肃南、民乐和甘州外,变化趋势不显著。
Download:
|
|
生态用水量受到气温、植被状况和水分的影响。伴随着气温的上升,蒸发能力会增强,在保证充足水分的情况下,植被的蒸腾将会增加[21-22]。如果植被面积减少,生态用水量则会减少。同时生态用水量也会受到风速、湿度等各种气象因子的影响[23-25],伴随着风速的增加,生态用水量也会减小。因此,生态用水压力是多种因素综合作用的结果,在相互作用中存在着相互抵消的现象。在全球变暖的背景下,祁连、肃南、高台、山丹的最终结果是生态用水量有一定程度的减少;而其他行政区域则为增加。相关研究表明祁连和肃南地区的风速在过去30年有增加趋势;而其他行政区域与中国政府一系列的植树造林和生态工程有着密切的联系,特别是山丹和甘州。
2.3.2 生产用水压力从图 2可以看出,上游(黑河祁连)以及张掖市的肃南、甘州、山丹、民乐等县级行政区呈现出下降趋势。其中肃南、山丹、甘州、民乐的下降趋势显著,下降幅度分别是0.001 9,0.015 5, 0.410 4和0.006 4/a。临泽、高台、肃州、嘉峪关、金塔和额济纳等6个行政区表现上升趋势,其中肃州、嘉峪关、金塔和额济纳的变化趋势显著。上升幅度最大的是额济纳,平均每年上升0.082 2。总体上讲,黑河流域各行政区的生产用水压力有增有减,变化趋势比较显著;其中甘州区下降最大,额济纳增加最多。
Download:
|
|
生产用水包括农业生产用水和工业生产用水,而农业生产用水起着决定性的作用。根据年鉴统计,1995—2014年祁连县的农业灌溉面积有下降趋势,从而导致生产用水减少。但中下游地区农业用水量在不断地增加;一方面是由于农业灌溉面积不断扩张;同时与气温的升高也有一定的关系[26-27]。而工业生产用水量进一步增加了中游区域的水资源压力。在过去的20多年里,中国政府推行一系列的农业节水技术和节水政策,一部分区域取得了一定的成效,但区域间存在着差异。这可能是中游有些县级行政区域生产用水量减少的原因。
2.3.3 生活用水压力从图 3可以看出,祁连、民乐、山丹、临泽、高台、肃州、嘉峪关、金塔和额济纳都表现出显著的上升趋势,其中临泽县的上升幅度最大,平均为0.002 6/a。但是肃南和甘州表现出下降趋势。总体上看,生活用水分量对水资源的压力贡献量小,变化趋势尽管显著,但是变化幅度要比生态和生产用水幅度小得多。
Download:
|
|
生活用水量对黑河流域的水资源压力最小。生活用水量的变化与人口变化和城市化进程有着密切的联系。人口密集、增长速度快的区域生活用水量的增加速率更快,很显然,黑河流域中游行政区生活用水的压力变化幅度要大于上游和下游行政区。同时城市化进程也影响着生活用水压力的变化,因为城市人均生活用水量要远远高于农村人均生活用水量。农村人口转化为女城市人口则会导致生活用水量的增加。中游行政区域生活用水压力指数的增加也归因于城市化。而肃南和甘州的生活用水量压力呈现出下降趋势可能和居民节水意识的提高以及水价有关[28-29]。
3 结论与讨论 3.1 结论本文基于SWAT水文模型和Man-Kendall测试法,结合1995—2014经济社会统计数据对黑河流域县域的水资源压力进行分析,主要结论如下:
1) 从水资源压力的县域空间差异来看,黑河流域的甘州区水资源总压力最大,其次是额济纳旗。它们的用水量已经超过自然系统所能提供的水资源量,出现水量赤字。而这种水量赤字则通过开采地下水来填平。水资源总压力主要来自于生产用水压力,生态用水压力也有较大的贡献。山丹县也出现水量赤字,但水资源压力主要来自于生态用水量。临泽和民乐的水资源压力也很高,接近可用水资源的极限。山丹的生态用水压力最大;甘州的生产用水压力和生活用水压力最大。
2) 从时间序列的变化趋势来看,祁连、肃南、高台、山丹和额济纳旗的生态用水压力呈现下降趋势;其中肃南的下降速率最大;其余的县级行政单位生态用水压力在1995—2014年呈现出上升趋势,其中最大上升速率出现在甘州。这表明这6个县级行政区生态用水量在1995—2014年间增加,其中甘州增加速度最快。黑河流域各行政区的生产用水压力有增有减,变化趋势比较显著。祁连、肃南、高台、山丹和额济纳旗的生态用水压力呈现下降趋势;其中肃南的下降速率最大。祁连、肃南、甘州、山丹、民乐等县级行政区呈现下降趋势。临泽、高台、肃州、嘉峪关、金塔和额济纳旗表现上升趋势。甘州下降速率最大,额济纳旗增加速率最大。生活用水压力除肃南和甘州表现下降趋势外,其他所有行政区表现显著性地增加。这种县域变化趋势的差异是县的气候、经济、人口、农业、植被、政策、城市化进程等因素综合作用和权衡的结果。
3) 综合利用水文模型、统计方法,结合社会经济统计数据的,构建反映水资源主要矛盾的水压力,能够在有限数据的情况下有效地评价干旱和半干旱地区县域的水资源利用状况。
3.2 讨论黑河流域地处偏僻的西北,由于各种原因,长时间序列的数据很难获得,而县级尺度的数据更是有限。因此,复杂指标的水资源评价体系在该区域不适用。弥补数据限制和构建简单但能反映水资源收支的评价体系是黑河流域县域水资源评价的难点。本研究基于半分布式的水文模型——SWAT获得长时间序列的县域可利用水资源量和生态用水量,结合社会经济统计数据构建水资源压力值体系评价黑河流域县域水资源利用状况,可为数据有限的干旱和半干旱区县域水资源评价提供方法借鉴。
我们的实验表明流域内县域的水资源利用结构和变化趋势存在很大差异。这是因为不同的县有其具体的县情,在经济、人口、农业、植被、政策、城市化进程方面存在着差异。因此,水资源管理者和决策者需要根据县情采取不同的水资源管理政策。例如甘州迫切需要优化产业结构、提高生产中的水资源利用率;而山丹县生态用水压力高,需要优化造林的树种结构,以防止森林消耗过多水资源过多地挤占生产和生活用水。
[1] |
中国科学院可持续发展战略研究组. 2007中国可持续发展战略报告:水治理与创新[M]. 北京: 科学出版社, 2007.
|
[2] |
Falkenmark M. The massive water scarcity now threatening Africa:why isn't it being addressed?[J]. Ambio, 1989, 18(2):112–118.
|
[3] |
朱法君, 邬扬明. 浙江省各地市水资源压力指数评价[J]. 长江科学院院报, 2010, 27(9):5–28.
|
[4] |
廖乐, 吴宜进, 毕旭. 湖北省各主要地市水资源压力指数评价[J]. 环境保护科学, 2012, 38(3):82–94.
|
[5] |
张莉莎. 北京水资源压力量化及趋势分析[J]. 经济视角, 2013(1):1–6.
|
[6] |
郭卫华, 周永章. 基于PLSR的中山市水资源压力演变特征与趋势[J]. 水资源保护, 2014, 30(1):23–26.
|
[7] |
靳美娟. 陕西省各地市水资源压力指数评价[J]. 河南科学, 2014, 32(2):279–282.
|
[8] |
Liu W, Feng Q. Water resource environment of the Heihe River Basin, China[J]. J Exp Bot, 2003, 54:47–47.
DOI:10.1093/jxb/erg015 |
[9] |
Cheng G D, Li X, Zhao W Z, et al. Integrated study of the water-ecosystem-economy in the Heihe River Basin[J]. Natl Sci Rev, 2014, 1(3):413–428.
DOI:10.1093/nsr/nwu017 |
[10] |
Ghaffari G, Keesstra S, Ghodousi J, et al. SWAT-simulated hydrological impact of land-use change in the Zanjanrood Basin, Northwest Iran[J]. Hydrol Process, 2010, 24(7):892–903.
DOI:10.1002/hyp.v24:7 |
[11] |
Park J Y, Yu Y S, Hwang S J, et al. SWAT modeling of best management practices for Chungju dam watershed in South Korea under future climate change scenarios[J]. Paddy Water Environ, 2014, 12(1):65–75.
|
[12] |
Fukunaga D C, Cecilio R A, Zanetti S S, et al. Application of the SWAT hydrologic model to a tropical watershed at Brazil[J]. Catena, 2015, 125:206–213.
DOI:10.1016/j.catena.2014.10.032 |
[13] |
Baker T J, Miller S N. Using the soil and water assessment tool (SWAT) to assess land use impact on water resources in an East African watershed[J]. Journal of Hydrology, 2013, 486:100–111.
DOI:10.1016/j.jhydrol.2013.01.041 |
[14] |
卢晓宁, 韩建宁, 熊东红, 等. 基于SWAT模型的忠县虾子岭流域地表径流特征浅析[J]. 长江科学院院报, 2010, 27(11):15–19.
DOI:10.3969/j.issn.1001-5485.2010.11.004 |
[15] |
Moriasi D N, Arnold J G, Van Liew M W, et al. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations[J]. Transactions to Asabe, 2007, 50(3):885–900.
DOI:10.13031/2013.23153 |
[16] |
Troin M, Caya D. Evaluating the SWAT's snow hydrology over a Northern Quebec watershed[J]. Hydrol Process, 2014, 28(4):1858–1873.
DOI:10.1002/hyp.9730 |
[17] |
黄钰, 焦国军, 孔霞. 甘州区井灌区地下水资源分析研究[J]. 河南科技, 2015, 555(1):92–93.
|
[18] |
毛彦成, 张勃, 张华. 绿洲土地利用/覆盖变化的社会经济与自然驱动力分析:以张掖市甘州区为例[J]. 干旱区资源与环境, 2007, 21(2):90–93.
|
[19] |
田静, 苏洪波, 陈少辉, 等. 黑河中游绿洲荒漠化的时空变化遥感分析[J]. 资源科学, 2011, 33(2):347–355.
|
[20] |
蔡国英, 徐中民. 黑河流域中游地区国民经济用水投入产出分析:以张掖市为例[J]. 冰川冻土, 2013, 35(3):770–775.
|
[21] |
Zhang Y Q, Leuning R, Chiew F H S, et al. Decadal trends in evaporation from global energy and water balances[J]. J Hydrometeorol, 2012, 13(1):379–391.
DOI:10.1175/JHM-D-11-012.1 |
[22] |
Liu Y, Ren L, Yang X, et al. Effects of precipitation and potential evaporation on actual evapotranspiration over the Laohahe basin, northern China[J]. P Int Ass Hydrol Sci, 2015, 371:173–179.
|
[23] |
Roderick M L, Farquhar G D. The cause of decreased pan evaporation over the past 50 years[J]. Science, 2002, 298(5597):1410–1411.
|
[24] |
Liu B, Ma Z G, Xu J J, et al. Comparison of pan evaporation and actual evaporation estimated by land surface model in Xinjiang from 1960 to 2005[J]. J Geogr Sci, 2009, 19(4):502–512.
DOI:10.1007/s11442-009-0502-5 |
[25] |
McJannet D L, Webster I T, Cook F J. An area-dependent wind function for estimating open water evaporation using land-based meteorological data[J]. Environ Modell Softw, 2012, 31:76–83.
DOI:10.1016/j.envsoft.2011.11.017 |
[26] |
Han S J, Xu D, Wang S L. Decreasing potential evaporation trends in China from 1956 to 2005:Accelerated in regions with significant agricultural influence[J]. agricultural and Forest Meteorology, 2012, 154:44–56.
|
[27] |
Meijninger W M L, Jarmain C. Satellite-based annual evaporation estimates of invasive alien plant species and native vegetation in South Africa[J]. Water Sa, 2014, 40:95–107.
DOI:10.4314/wsa.v40i1.12 |
[28] |
贾国宁, 黄平. 居民用水阶梯式水价及其节水效果测算模型研究[J]. 自然资源学报, 2013, 28(10):1788–1793.
DOI:10.11849/zrzyxb.2013.10.013 |
[29] |
李琳琳, 陈晓光. 水价对中国北方城市居民用水需求影响分析[J]. 水利经济, 2007, 25(2):41–42.
DOI:10.3880/j.issn.1003-9511.2007.02.013 |