Optimization of integrated observation station layout for terrestrial surface natural resources
GAO Bingbo,1, WANG Jinfeng,2, HU Maogui2, XU Chengdong2, LIU Huilan3, ZHOU Chenghu2通讯作者:
收稿日期:2020-09-4修回日期:2020-10-7网络出版日期:2020-10-25
基金资助: |
Received:2020-09-4Revised:2020-10-7Online:2020-10-25
作者简介 About authors
高秉博,男,宁夏海原人,副教授,研究方向为空间统计方法及其在资源环境中的应用。E-mail:
摘要
关键词:
Abstract
Keywords:
PDF (11219KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex收藏本文
本文引用格式
高秉博, 王劲峰, 胡茂桂, 徐成东, 刘慧兰, 周成虎. 中国陆表自然资源综合观测台站布点优化. 资源科学[J], 2020, 42(10): 1911-1920 doi:10.18402/resci.2020.10.08
GAO Bingbo, WANG Jinfeng, HU Maogui, XU Chengdong, LIU Huilan, ZHOU Chenghu.
1 引言
人类社会的可持续发展需要合理并可持续地利用自然资源。联合国规划署将自然资源定义为“在一定时间、地点的条件下能够产生经济价值的、以提高人类当前和将来福利的自然环境因素的总称”[1]。自然资源包括海洋自然资源与陆地自然资源,其中陆地自然资源与人类关系更加密切[2]。陆地自然资源主要包括土地资源、水资源、气候资源、生物资源和矿产资源。其中土地资源、水资源(地表水)、气候资源、生物资源(四种资源合起来简称“土水气生)等陆地表层资源易于观测,且与人类生产生活直接相关[3]。陆表的“土水气生”自然资源具有空间分布不均的特点,不同区域的资源种类及其数量具有较大差异,需要因地制宜合理利用。不合理的自然资源利用方式不仅影响自然资源的可持续性,还会引起严重的生态环境问题[4]。中国是人口大国,尤其要注重对自然资源的合理开发与利 用[5]。目前中国由于不合理的自然资源利用方式已经造成了严重的资源环境问题,如东北地区土壤有机质流失、西北地区荒漠化、华北地区地下水超采、南方地区土壤重金属污染与酸化等。科学合理的资源管理政策需要建立在长期综合的观测数据基础之上[6]。自然资源具有位置属性、数量属性、质量属性、生态经济及相关属性[7]。自然资源综合观测网络不仅需要能摸清自然资源各项属性指标的总体情况,更需要及时捕获自然资源各项属性指标的空间分布与动态变化,为自然资源高效监管和科学利用提供数据基础。目前,国际上已经建立了若干全球、大洲或者国家尺度的地球系统观测网络,如国际长期生态学研究网络(International Long-Term Ecological Research, ILTER)、全球地球关键带观测实验研究网络(Critical Zone Observatories, CZO)、美国国家生态系统观测网络(National Ecological Observatory Network, NEON)、英国环境变化网络(Environmental Change Network,ECN)、澳大利亚陆地生态研究网络(Terrestrial Ecosystem Research Network,TERN)等。国内也建立了气象、地震、环境、水文、生态等全国性的观测网络,但是它们均为专项观测,在服务于自然资源管理方面,观测要素不够全面。如气象观测台站主要观测大气温度、湿度、风速、太阳辐射等气象指标;地震观测台站主要观测地磁、重力、水位、水温等地震相关指标;环境观测网分为水环境、土壤环境、大气环境3类观测网络,主要观测水土气中污染物含量等指标;水文观测站主要观测水温、水位、流速、含沙量等水文指标,均不能全部涵盖陆表“土水气生”资源的各项属性指标。同时,由于各种专项观测网络在规划设计阶段缺乏有效的协同和统一的建设目标,其监测目的、指标与布设方法差异较大,目前尚未能联合形成一个高效的自然资源监测网络。中国生态系统研究网络(CERN)针对农田、森林、草地、荒漠、沼泽等生态系统制订了生物、土壤、水分和大气观测指标,但主要以生态系统为观测对象,未能包含针对冰川、冻土等陆表自然资源观测指标。而在陆表自然资源综合观测中,需要涵盖土地资源、地表水资源、气候资源和生物资源(主要为植物资源)的各项指标。另外,虽然国内针对一些典型流域和生态系统还建立了综合观测网络,如黑河流域地表过程综合观测网络等,但是它们覆盖的空间范围有限,不能反映全国自然资源状况和趋势。因此需要构建服务于自然资源管理的综合观测网络。
陆表自然资源位置、数量、质量、生态经济及相关属性的综合观测涉及多种类型的指标,需要结合空天地一体化的手段开展观测[8,9,10]。本文所选取的观测对象为陆地表层“土水气生”4类自然资源,主要以地面台站观测指标为对象,研究如何优化布设地面观测台站,以较少的代表性观测台站获得对观测指标在全国的整体状况、空间分布及其变化趋势的精确推断。综合观测台站的观测对象借鉴CERN的指标配置方式,为“土水气生”中不同的自然资源子类型分别设置不同的观测指标,在观测台站建设时,依据台站布设位置的主要自然资源类型配置观测指标。
基于空间统计学的空间抽样与统计推断是具有严格数学基础的科学方法,已广泛应用于社会经济、资源环境、土地利用和公共健康等领域的调查监测[11,12,13]。空间抽样方法包括随机抽样、便利抽样和目的性抽样3类[14]:①随机抽样方法给予抽样空间中的每个候选单元一个固定的入样概率,并基于入样概率进行样本选择,包括简单随机抽样及增加了限制条件的各种演化版本,如分层随机抽样、系统抽样、整群抽样等;②便利抽样优先考虑抽样成本,在易于获取样品的地方布样,如沿道路布设样点;③目的性抽样在一个具体且直接的目标指导下进行样点布设,目标函数法是其中最主要的一种方法。依据空间抽样与统计推断“三位一体”理论,空间抽样优化目标函数需要依据统计推断目标和总体的性质进行针对性地设计[15,16]。自然资源综合观测不仅需要了解全国各类陆表自然资源的总量,更需要了解各类自然资源的空间分布。因此局部空间预测是陆表自然资源观测的统计推断目标。目的性抽样在局部空间预测中有更高的效率,因此自然资源综合观测台站设计宜采用目的性抽样方法。目前,常见的空间抽样优化目标函数包括地理空间分布、特征空间分布和预测误差估计。在没有先验知识与辅助数据支撑时,样点的空间均匀分布是常见优化目标,如平均最短距离最小(Minimization of the Mean of the Shortest Distances,MMSD)、加权平均最短距离(Weighted Mean of the Shortest Distances,WMSD)等方法[17,18]。当具有与观测目标变量相关的多个辅助变量时,可以采用特征空间分布作为优化目标,如拉丁超立方体(Latin Hypercube Sampling,LHS)、均衡设计(Equal Range Stratification,ER)等方法[19,20,21]。也有同时优化地理空间分布和特征空间分布的空间抽样方法,如空间条件拉丁超立方体(Spatial Conditioned Latin Hypercube Sampling,ScLHS)[22]。当采样设计时能事先掌握目标变量的空间变异函数时,一般采用统计推断方法的预测误差估计作为优化目标函数,如平均克里金误差最小(Minimization of Average Kriging Variance,MAKV)[23]。在陆表自然资源综合观测台站布点优化设计中,统计推断方法的选择应充分考虑总体的性质。中国陆表自然资源的分布在大尺度上具有明显的空间分异特性,同时在小尺度上存在较高的空间自相关性[24]。因此,本文采用分层异质性表面无偏估计方法体系进行全国自然资源观测台站的布点设计[25,26]。
本文针对中国陆表自然资源综合观测台站布设优化需求,基于空间抽样与统计推断“三位一体”理论,探索陆表自然资源综合观测台站布设的理论方法,试图回答观测对象是什么、需要多少个台站以及台站布设在什么位置3个问题,用于指导实际的综合观测台站布设,服务于国家自然资源综合观测与科学管护。
2 数据来源与研究方法
2.1 数据来源
为了反映全国自然资源的分异性,本文采用郑度等[24]绘制的“中国生态地理区域图”作为自然资源分异性区划图(图1)。该区域图的一级生态地理分区包含49个区(以下统称“生态地理区”):首先,按照温度、植被和土壤等自然要素的组合划分出11个温度带;其次,在温度带内部,按照干湿程度引起的生态地域差异,划分为湿润、半湿润、半干旱和干旱4类干湿地区;最后,在干湿区中,按照生态地貌类型组合差异引起的生态地理特征差异,划分得到49个生态地理区。该区域图综合考虑了气温、水分、地貌、植被、土壤等因素,能够较好地刻画自然资源的空间分异情况[24]。区域图中的分区之间,自然资源赋存具有显著差异。图1
新窗口打开|下载原图ZIP|生成PPT图12009—2018年全国陆表平均NDVI分布图
Figure 1Distribution of average Normalized Difference Vegetation Index (NDVI) values of China, 2009-2018
多年的植被情况也能够较好地表征陆表自然资源赋存差异。卫星遥感数据产品归一化植被指数(Normalized Difference Vegetation Index, NDVI),能够利用不同地表对红外和近红外反射的差异反映地表植被覆被情况。NDVI值介于-1到1之间(注:MODIS的NDVI产品为了计算方便,将原始NDVI值放大10000倍),负值表示地面覆盖为水或者雪等,对可见光具有较高反射率;0表示地表覆被为岩石或裸土等;正值表示地表有植被覆盖,且NDVI值随植被覆盖度增大而增大。不仅水体、林地、草原、荒漠等不同陆表自然资源具有不同的NDVI取值,不同地形地貌还具有不同的NDVI空间变异程度,因此本文采用全国陆表近10年(2009—2018)的MODIS平均NDVI,反映自然资源在各生态地理区内的空间变异特性和在各生态地理区之间的空间关联性,采用从平均NDVI中提取的空间协方差参数指导观测台站布设。
2.2 陆表自然资源综合观测对象及其空间分布表征
陆表自然资源综合观测台站的观测对象主要包括土地资源、地表水资源、气候资源和生物资源(主要为植物资源)。其中土地资源的子类型主要包括土壤、地形和利用方式;地表水资源主要包括河流、沼泽、湖泊、冰川和永久性积雪;气候资源主要包括大气水、光能和热能以及风能;生物资源主要包括森林、草场、作物、荒漠等。虽然陆表自然资源在大尺度上的空间分异性和在小尺度上的空间相关性已被认知,但目前尚缺乏定量表征其空间分布特征的手段。“中国生态地理区域图”中的生态地理区综合考虑了温度、干湿程度、植被、土壤和地形地貌因素,能够直接表征土壤、地形、大气水、光能和热能、风能、森林、草场、作物和荒漠等子类型的空间分异特征。而NDVI能够定量反映植被覆盖程度,能够直接反映生物资源的空间分异性和相关性;同时由于水体对可见光具有高反射率,其NDVI为负值,也可用NDVI直接表征地表水资源的空间分布特征。对于气候资源,由于不同气候条件会产生不同的自然植被类型,同时也会影响人工植被类型以及农作物熟制等,使得不同气候条件下年均植被覆盖度具有差异,因此可以用NDVI间接反映气候资源的空间分异性和相关性。同样,对于土地资源,不同的地形形成的垂直地带性一般通过植被进行表征,而土地利用方式则直接决定地表是否有植被覆盖以及覆盖植被的类型、覆盖度、生长期等属性,因此可以用常年NDVI间接反映不同地形和土地利用现状的空间差异。对自然资源综合观测所涉及观测对象空间分布表征依据的说明详见表1。Table 1
表1
表1陆表自然资源空间分布表征依据
Table 1
陆表自然资源类型 | 子类型 | 空间分布表征变量 | 表征依据 |
---|---|---|---|
土地资源 | 土壤 | 生态地理区 | 生态地理划分区中考虑了全国土壤差异,能反映全国土壤的分异性 |
地形 | (1)生态地理区 (2)多年NDVI | (1)生态地理划分区中考虑了大尺度的地形差异 (2)中小尺度地形形成垂直地带性,影响植被分布,因此NDVI能间接表征中小尺度地形的空间分布特征 | |
利用方式 | 多年NDVI | 不同土地利用方式具有不同的植被覆盖度,如裸土、建筑用地、农田、草原和林地等,年均NDVI具有较大的差异,因此NDVI能间接表征土地利用方式的差异 | |
地表水资源 | 河流 | 多年NDVI | 对可见光高反射,NDVI为负值,可用NDVI直接表征 |
沼泽 | |||
湖泊 | |||
冰川 | |||
永久性积雪 | |||
气候资源 | 大气水 | (1)生态地理区 (2)多年NDVI | (1)生态地理区考虑了降水条件,刻画了大尺度的大气水的空间分异性 (2)不同大气水条件影响自然植被类型,形成如森林、草原和荒漠等差别,影响人工栽培植被的作物类型、熟制等,均在年均NDVI上有差异,因此多年NDVI能够间接表征大气水的空间分布特征 |
光能和热能 | (1)生态地理区 (2)多年NDVI | (1)生态地理区考虑了光能和热能条件,刻画了大尺度的光能和热能的空间分异性 (2)不同的光能热能影响自然植被与人工植被类型,在年均NDVI上有差异,因此可以用多年NDVI间接表征光能和热能的空间分布 | |
风能 | 生态地理区 | 生态地理区包含温度带和地形地貌,能过表征风能资源的空间分异 | |
生物资源(植物) | 森林 | (1)生态地理区 (2)多年NDVI | (1)生态地理区划分反映了植被差异 (2)植被覆盖度越大,NDVI值越高;多年NVDI能直接表征森林空间分布特征 |
草场 | (1)生态地理区 (2)多年NDVI | (1)生态地理区划分反映了植被差异 (2)植被覆盖度越大,V值越高;多年NVDI能直接表征草场空间分布特征 | |
作物 | (1)生态地理区 (2)多年NDVI | (1)生态地理区划分反映了植被差异 (2)植被覆盖度越大,NDVI值越高;多年NVDI能直接表征作物空间分布特征 | |
荒漠 | (1)生态地理区 (2)多年NDVI | (1)生态地理区划分了不同的荒漠区 (2)植被覆盖度越低,NDVI值越低;多年NVDI能直接表征荒漠空间分布特征 |
新窗口打开|下载CSV
2.3 综合观测台站布点设计
2.3.1 自然资源空间分层异质性探测为了验证采用全国陆表近10年平均NDVI反映陆表自然资源空间分布特征的合理性,使用地理探测方法量化分析NDVI与“生态地理区域”之间的关联性。地理探测器是一种基于空间分层异质性度量连续变量与类别变量之间相关关系的统计指标,可以量化计算线性和非线性相关关系,并提供显著性检验[27,28,29]。空间异质性指数计算方法如下:
式中:q为空间异质性指数;N为研究区所有NDVI像素数;
式中:L同公式(1);n为样本量;
2.3.2 综合观测台站布设优化目标函数设计
陆表自然资源空间部分同时存在空间相关性和异质性,因此采用能够同时考虑空间相关性和异质性的分层异质性表面无偏最优估计方法进行统计推断。自然资源综合观测除了需要估计各类资源的全国全局均值以外,还需要进行空间预测推断局部自然资源状况,即空间插值。为此,以分层异质性表面无偏最优点估计(Point Mean of the Surface with stratified Nonhomogeneity, P-MSN)方法估计结果的平均误差方差作为优化目标函数。自然资源观测网络优化目标函数如公式(3)所示:
式中:
式中:
式中:
2.3.3 综合观测台站布设优化求解
采用空间模拟退火算法(Spatial Simulated Annealing, SSA)求解优化目标。空间模拟退火算法是一种全局优化方法,最初起源于对热力学中退火过程的模拟:固体物质中的缓慢降温处理以使整体能量处于最低状态[30]。它在产生新解时,在当前解中任意选定的样点周围以随机的半径和角度进行搜索,并用得到的点替换当前解中的点。该方法简单易行,能有效避免陷入局部最优的问题,在空间采样优化中被广泛应用。算法具体步骤如下:
(1)制定采样空间,设定初始温度T0,迭代终止条件如最低温度、连续拒绝次数、总迭代次数等。
(2)生成初始样本
(3)扰动
式中:t为当前温度;exp()为以自然常数e为底的指数函数。
(4)判断是否达到迭代终止条件,如达到转步骤(5);如未达到,进行降温并转步骤(3)继续迭代执行。
(5)输出
同时,为了进行对比,采用普通克里金(Ordinary Kriging, OK)估计的平均误差方差最小为优化目标函数进行观测台站优化布设。由于OK方法没有考虑空间分异性,其优化结果不能保证每个生态地理区都有观测点。因此将基于OK方法的优化方案进行调整,调整方法分为两个步骤:①筛选出没有布设观测的m个生态地理区,在区内随机选择一个点,加入观测网;②找到观测点最多的m个分区,分别从中随机选取一个台站并去除,保证观测台站数目不变。基于OK方法的平均估计误差方差的计算方法如下:
式中:N同式(1);
对于基于P-MSN和OK的两种优化方法分别设置观测台站数为100、300、500、800、1000、3000和50000,通过优化算法生成7个布设方案。
3 结果与分析
3.1 自然资源空间分层异质性探测
使用地理探测器计算植被指数在生态地理区域中的空间分异性指数,检验其显著性。空间分异性指数高达0.813,显著性检验结果为极显著(0.000),说明生态地理区域划分能够精确地刻画植被指数的空间分异,反映自然资源的空间分异格局。从图1也能看出,不同的生态地理区之间的植被指数差异较为明显,说明生态地理区能够较好地反映出10年的植被指数空间分异情况。因此,采用近10年的植被指数代表自然资源空间分布特征具有充分的合理性。3.2 自然资源综合观测台站布点设计结果
采用基于P-MSN和OK的两种优化方法分别优化布设了100、300、500、800、1000、3000、5000个自然资源综合观测台站。针对每一套方案,分别计算台站布设后对全国陆表单元进行插值所得平均误差方差,其中基于P-MSN优化方法的平均估计误差方差按照公式(3)计算,基于OK优化方法的平均估计误差方差的计算方法如公式(7)。两种方法所得布设方案的估计误差方差随观测台站数变化曲线如图2所示。在相同的观测台站数量下,基于P-MSN的台站优化布设方法比基于OK的台站优化布设方法具有更小的误差方差,而且台站数量越少时,基于P-MSN的台站优化布设方法的优势越明显。同时,由图2可知,随着台站数量的增加,误差方差逐渐降低,并且降低的速率不断变缓。当台站数量超过1000时,增加台站数量对于误差方差的改进不大。图2
新窗口打开|下载原图ZIP|生成PPT图2估计误差方差随观测台站数量变化曲线图
Figure 2Estimated error variance as measured by the number of stations
图3分别展示了台站数量为100、300、500、800、1000时两种优化方法的台站布设方案。从图中可以看出基于OK的台站优化布设方法在全国布设的观测台站(图3b、3d、3f、3h、3j)较为均匀,而基于P-MSN的台站优化方法的布设结果(图3a、3c、3e、3g、3i)在不同的生态地理区密集程度有较大的差别。在西藏东高山深谷、东喜马拉雅南翼、琼雷山地丘陵、青东祁连山地、台湾南部低山等生态地理分区中,由于地形变化较大导致自然资源变化较为剧烈,台站布设较为密集。在塔里木盆地与吐鲁番盆地、阿拉善及河西走廊、昆仑高山高原、羌塘高原湖盆、贵州高原、四川盆地等生态地理分区中,由于自然资源较为单一且空间变异较小,台站布设较为稀疏。而在淮南与长江中下游、江南丘陵、江南与南岭山地、东北东部山地等生态地理分区中,由于自然资源的复杂性和其空间变异处于中等水平,因此观测台站布设密度居中。同时,从图3中可以看出,基于P-MSN的优化方法在不同生态地理区边界处的台站布设较为合理,由于P-MSN方法考虑了层间的空间相关性,对于相邻的两个生态地理区,在边界处既没有出现观测站点的过分密集,也没有出现都未覆盖的空白区域。
图3
新窗口打开|下载原图ZIP|生成PPT图3观测网布设方案
Figure 3A scheme of monitoring network layout
两种方法所得的台站布设方案对省份、地市和区县的覆盖情况如表2所示。两者的行政区覆盖度非常接近。当布设100个台站时,基于P-MSN的优化方法能覆盖79.41%的省份,基于OK的优化方法能覆盖82.35%的省份。当布设5000个台站时,基于P-MSN的优化方法能覆盖除澳门以外的其他33个省份,基于OK的优化方法能覆盖除香港和澳门以外的其他32个省份。当选择1000个台站时,使用基于P-MSN的优化方法,能够覆盖94.12%的省份、83.53%的地市和23.37%的区县。基于OK的优化方法能够覆94.12%的省份、84.68%的地市和23.74%的区县。总体来看,当观测台站数量小于1000时,基于P-MSN的优化方法在区县覆盖方面更有优势,而基于OK的优化方法在地市和省份覆盖方面更好。当台站数量大于1000时,两种方法对三级行政区的覆盖程度接近。因此,基于P-MSN的优化方法在取得较小的估计误差方差时并未降低对省市县行政区的覆盖,能够同时兼顾对行政区与自然资源空间变化情况的代表性。
Table 2
表2
表2观测台站覆盖区域
Table 2
台站数量 | 覆盖区县比例/% | 覆盖地市比例/% | 覆盖省份比例/% | |||
---|---|---|---|---|---|---|
P-MSN | OK | P-MSN | OK | P-MSN | OK | |
100 | 3.32 | 3.28 | 20.81 | 21.97 | 79.41 | 82.35 |
300 | 9.14 | 8.87 | 46.53 | 47.98 | 82.35 | 85.29 |
500 | 13.99 | 13.51 | 65.90 | 64.45 | 91.18 | 88.24 |
800 | 20.01 | 19.47 | 77.46 | 76.59 | 94.12 | 91.18 |
1000 | 23.37 | 23.74 | 83.53 | 84.68 | 94.12 | 94.12 |
3000 | 48.02 | 48.09 | 95.09 | 96.82 | 94.12 | 94.12 |
5000 | 61.09 | 62.00 | 97.40 | 97.11 | 97.06 | 94.12 |
新窗口打开|下载CSV
4 结论与讨论
4.1 结论
本文依据空间抽样与统计推断“三位一体”理论进行自然资源综合观测台站布点设计,通过分析自然资源的空间变异特征,采用中国生态地理区域图刻画自然资源区域分异,并采用近10年NDVI进一步表征分区内的自然资源空间变异和分区之间相关关系;在此基础上选用能够同时考虑空间相关性和异质性的分层异质性表面无偏最优估计方法构建优化目标函数,并使用空间模拟退火算法完成台站布点设计。主要研究结论如下:(1)多年平均NDVI能够较好地表征陆表自然资源的空间分布特征。
(2)分层异质性表面无偏最优估计方法能够同时考虑自然资源的空间相关性和分异性,基于该方法的优化布设方案在自然资源变化较大的区域具有较密的台站,在自然资源变化平缓的区域具有较稀疏的台站,与基于克里金的方法相比,能够以同样的观测台站数产生更小的误差。
(3)台站数量与误差方差曲线图可以作为确定台站数量的依据,依据该曲线,本文推荐1000点作为台站数量。
4.2 讨论
自然资源种类较多,观测指标庞杂,观测手段各异。本文仅以“土水气生”陆地表层自然资源为观测对象,针对地面台站观测方式,基于全国自然资源的空间分异性和相关性,采用分层异质性表面无偏估计方法体系和空间模拟退火算法进行全国自然资源观测台站的优化布设的理论方法研究。距离实际指导自然资源观测台站布设还有一定的距离。在实际的自然资源观测中,需要进一步考虑其他类型自然资源的观测。同时还需要与遥感、非固定点等其他观测手段结合,互相补充和印证,形成自然资源观测网。另外,需要在观测网的优化设计中将中国已有自然资源专项观测台站作为固定点,通过增加设备与观测指标形成综合观测能力,避免重复建设。参考文献 原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子
[J]. ,
[本文引用: 1]
[J].
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J].
[本文引用: 1]
[J/OL]. , (
URL [本文引用: 1]
[J/OL].
URL [本文引用: 1]
[J]. ,
[本文引用: 1]
[J].
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J].
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J].
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J].
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J].
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]
[M]. ,
[本文引用: 1]
[EB/OL]. (
URL [本文引用: 1]
[EB/OL]. (
URL [本文引用: 1]
[J]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]
[M]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]
[本文引用: 3]
[本文引用: 3]
[J]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J].
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]
[J]. ,
[本文引用: 1]