ISSN 1000-3665 CN 11-2202/P
  • 中文核心期刊
  • GeoRef收录期刊
  • Scopus 收录期刊
  • 中国科技核心期刊
  • DOAJ 收录期刊
  • CSCD(核心库)来源期刊
  • 《WJCI 报告》收录期刊
欢迎扫码关注“i环境微平台”

承德市土壤重金属空间结构与分布特征

安永龙, 万利勤, 李霞, 殷志强, 卫晓峰, 何泽新, 贾凤超

安永龙, 万利勤, 李霞, 殷志强, 卫晓峰, 何泽新, 贾凤超. 承德市土壤重金属空间结构与分布特征[J]. 水文地质工程地质, 2020, 47(6): 119-131. DOI: 10.16030/j.cnki.issn.1000-3665.202007081
引用本文: 安永龙, 万利勤, 李霞, 殷志强, 卫晓峰, 何泽新, 贾凤超. 承德市土壤重金属空间结构与分布特征[J]. 水文地质工程地质, 2020, 47(6): 119-131. DOI: 10.16030/j.cnki.issn.1000-3665.202007081
AN Yonglong, WAN Liqin, LI Xia, YIN Zhiqiang, WEI Xiaofeng, HE Zexin, JIA Fengchao. Spatial structure and distribution characteristics of heavy metals in the soil in Chengde[J]. Hydrogeology & Engineering Geology, 2020, 47(6): 119-131. DOI: 10.16030/j.cnki.issn.1000-3665.202007081
Citation: AN Yonglong, WAN Liqin, LI Xia, YIN Zhiqiang, WEI Xiaofeng, HE Zexin, JIA Fengchao. Spatial structure and distribution characteristics of heavy metals in the soil in Chengde[J]. Hydrogeology & Engineering Geology, 2020, 47(6): 119-131. DOI: 10.16030/j.cnki.issn.1000-3665.202007081

承德市土壤重金属空间结构与分布特征

基金项目: 

中国地质调查局地质调查项目 DD20190822

详细信息
    作者简介:

    安永龙(1988-),男,硕士,工程师,从事地质学、地球化学相关研究。E-mail:aylzfj@163.com

    通讯作者:

    万利勤(1973-),女,博士,高级工程师,从事水文地质、地球化学相关研究。E-mail:wanlq@cigem.cn

  • 中图分类号: P632+.1

Spatial structure and distribution characteristics of heavy metals in the soil in Chengde

  • 摘要:

    基于地统计学和GIS相结合的方法,对承德全域表层土壤重金属As、Cd、Cr、Cu、Hg、Ni、Pb、Zn、微量元素Se正态分布特征、主导分布趋势及相互作用规律进行了分析,确定了不同元素最适宜的地统计插值模型并厘定出其空间分布规律。结果表明:As、Cd、Cr、Cu、Hg、Ni、Pb、Zn的质量含量平均值分别为8.28,0.200,60.85,24.37,0.034,27.76,26.65,77.10 mg/kg,Cd、Cu、Hg和Pb变异系数分别为385%、143%、350%、118%,分异性强。Zn含量均值受土壤类型影响显著,Cr、Cu、Ni含量均值则受土地利用类型影响显著。经过不同趋势阶数元素插值误差的综合对比,确定As、Cr、Pb、Zn、Ni、Se适宜选择无趋势参数,Hg和Cu适宜选择一阶趋势参数,而Cd适宜选择二阶趋势参数。As的理论模型为指数模型,主要受到结构性因素的影响;Cd、Cr、Cu、Hg、Ni、Pb、Zn、Se的理论模型为线性模型,主要受到随机性因素的影响。通过普通克里格插值图 可见区内9种元素具有北低南高的特点,中部地区形成了一条较宽的Pb高值带,与Cd相似。按照含量分布特点,土壤中Cr和Ni、Cu和Hg、Zn和Pb、Se和Cd之间的高值空间展布区具有相似性且来源相同,仅As具有个性,分析结果与传统统计学结果数据保持部分一致性。

    Abstract:

    Based on the combination of geo-statistics and GIS, this paper analyzed the normal distribution characteristics, dominant distribution trends and interaction regulations of heavy metals including As, Cd, Cr, Cu, Hg, Ni, Pb, Zn, and trace elements Se in the topsoil in whole area of Chengde. Both the most suitable geo-statistics interpolation model and spatial distribution law for these elements were determined. The results showed that the average contents of As, Cd, Cr, Cu, Hg, Ni and Pb, were respectively 8.28,0.200, 60.85, 24.37, 0.034, 27.76, 26.65, 77.1 mg/kg. The coefficients of variation of Cd, Cu, Hg, and Pb were respectively 385%, 143%, 350%, and 118%, which indicated the level of dispersion was great. The average value of Zn contents was significantly affected by soil types. While the average value of Cr, Cu, Ni contents was significantly affected by land use types. After comprehensively comparing the interpolation errors of order elements of different trends, the following results were obtained that non-trend parameters were suitable for As, Cr, Pb, Zn, Ni, and Se, first-order trend parameters were suitable for both Hg and Cu, while second-order trend parameters were suitable for Cd. The theoretical model of As was an exponential model, which was mainly affected by structural factors. The theoretical models of Cd, Cr, Cu, Hg, Ni, Pb, Zn, and Se were linear models, which were mainly affected by random factors. According to the Kriging method, it could be seen that the 9 elements in the area presented a characteristic of low level in the north and high in the south. And a relatively wide Pb high-value band was formed in the central area, similar to Cd. According to the characteristics of contents distribution, the high-value spatial distribution areas between Cr and Ni, Cu and Hg, Zn and Pb, Se and Cd were similar, and the sources were also the same, and only As was significantly different from them. The analysis results were partially consistent with the traditional statistical results.

  • 土壤是成土母岩经过风化、剥蚀、搬运、沉积等地质作用后叠加成土作用形成的自然物质体。在成土的每个阶段,都可能会受到外界环境(气候、海拔、地貌、人类活动等[1-2])的作用,因此在区域化土壤的物理化学性质和空间序列中形成了特有的空间变异性。传统的统计学方法往往忽略了土壤样本的固有属性数据和空间属性数据之间的交互式影响[3],导致元素的异质性和空间依赖性无法很好地表达[4]。目前已有不少国内外学者针对不同区域时空尺度,利用地统计学结合GIS的方法解决土壤元素空间变异性问题。如崔萌等[5]运用地统计学和主成分分析法对北京市大桃主产区土壤5种重金属含量和空间结构特征进行了研究,结果显示Pb、As、Cr、Hg含量受结构性因素影响,而Cd含量受结构性因素和随机性因素共同影响;汪璇等[6]运用地统计学方法对三峡库区表层土壤部分微量元素进行了空间变异特征分析,结果表明7种微量元素为中等自相关性,2种为弱自相关性,且4个方向上变异程度不明显;Roger等[7]通过环境因子作为数据对土壤养分指数完成克里格插值,表明克里格插值与线性回归共同分析可以一定程度提高土壤养分的预测精度。

    承德作为北京的“菜篮子”基地,土壤安全直接关系着首都人民的健康,在这样严峻的背景下摸清该地区土壤重金属本底及空间分布状况具有重要意义[8-9]。已有部分学者在区内按不同流域开展了土壤重金属调查研究工作,如孙厚云等[10]对承德市滦河流域12种重金属地球化学基线值进行厘定,发现表层土壤Cd、V、Ti等基线值高于河北省背景值,且Cd、Pb、Cu的累积程度较大;王梦雨等[11]对承德市柳河流域农田土壤重金属进行调查,发现Cd、Zn、Ni存在不同程度积累现象。本次研究基于地统计学和GIS技术,对承德市全域范围内表层土壤重金属As、Cd、Cr、Cu、Hg、Ni、Pb、Zn等的空间结构及分布特征进行解析,并对土壤重金属与Se之间的相关性进行了探索研究,旨在掌握承德地区表层土壤重金属空间分布状况及其特征,并进一步判断重金属与Se空间上是否存在一定关联,为承德地区着力打造健康安全的生态特色农果产业链,支撑“承德样板”[12]落地,开展国土空间规划和用途管控工作提供基础支撑。

    承德市处于华北和东北两个地区的过渡地带,南邻北京,东靠辽宁,北接内蒙古高原,地理位置十分重要。地势北高南低,北部为坝上高原地貌,中南部为山地、丘陵地貌,地貌展布具有独特的地理优势,为特色农果产品和道地药材种植提供了有利的地域条件,杏仁、山楂、板栗、金莲花、黄芩等都是承德的名优特产。承德市河流分属三个流域四大水系,即滦河流域、三河流域、辽河流域,年产水量37.6×108 m3,是京津唐的重要供水源地。区内属温带大陆性季风型半湿润山地气候,四季分明,雨热同期。年平均温度5.6 ℃,无霜期60~180 d,年降水量451~850 mm。土壤类型以褐土、棕壤、潮土、灰色森林土、粗骨土、栗钙土等为主。

    承德矿产资源分布较多,以铁矿资源最为丰富,是我国仅次于攀枝花的第二大钒钛磁铁矿资源基地[13]。铁矿和金矿集中分布于双滦区和滦平县北部以及宽城满族自治县周边,银矿、钼矿、铜矿等分布较为分散。

    本次研究利用了2017年承德地区土地质量地球化学调查数据,采样点分布情况见图 1。土壤样品采集过程中采用GPS定位结合地形图定点方式,充分考虑土壤类型、土地利用方式。以土地利用类型为采样的基本单元,土壤采样深度为0~20 cm,并以“S”或“N”形采样,去除杂草、砾石、虫壳、粪便等杂物,一个组合样由3~5个子样点等组分混合而成,子样点与采样单元内采集的土壤类型一致。采集样品的原始重量不低于1 kg,过20目筛后样品重量不少于500 g。尽量选择在平稳耕地、农用菜地、林果地、山坡下侧土层较厚等地带采样,避开明显点状污染或表土已被破坏的干扰地段,保证了样品的典型性和代表性。

    图  1  研究区土壤样点分布
    Figure  1.  Distribution of soil sampling sites in the study area

    测试的指标为类重金属As,重金属元素Cd、Cr、Cu、Hg、Ni、Pb、Zn,微量元素Se。采用氢化物发生原子荧光仪测定As和Se,采用ICP-OES (PE,USA)测定Cd、Cr、Cu、Hg、Ni、Pb、Zn。分析测试中准确度和精密度采用国家一级土壤标准物质(GBW07349)控制,加10%空白样与平行样控制,且合格率符合规范要求,指标的加标回收率均符合国家标准。

    应用Excel 2016进行土壤数据基本参数统计分析,应用SPSS 19.0完成土壤数据正态分布性检验、相关系数分析、聚类分析、因子分析等。为了避免在利用统计学和地统计学方法时可能出现的比例效应[14],所分析的数据必须全部服从正态分布或近似正态分布,否则需要将数据进行对数变换或者Box-Cox变换。聚类分析方法采用最近邻元素,选择区间Pearson相关性作为度量标准。因子分析采用最大四次方值法,抽取基于特征值大于1,最大收敛性迭代次数为25次。

    目前国内外判断数据是否服从正态分布的方法较多[15],主要与所使用的统计学软件类型(SAS、STATA、SPSS等)和样本量的多少有关。使用SPSS软件时,规则为探索性检验中基于理论假定的概率图、Q-Q图[16]、P-P图等不限定样本量;偏度峰度检验法适用于样本量N>200时;样本量7<N≤2 000时,以夏皮罗-威尔克(Shapiro-Wilk)法(也称W检验)为准;样本量N>2 000时,以Lilliefors检验(K-S检验[17])为准。基于处理软件和样本数量的限定,本文选用偏度峰度检验法和W 检验法,结果表明,偏度值和峰度值远大于1,W 检验的显著性均小于0.01;对9项元素进行对数变换后偏度值和峰度值基本在1附近,W 检验的显著性均大于0.05,符合正态分布的要求,满足地统计学分析的假设条件,可以进行半变异函数的计算。

    应用ArcGIS 10.2中Geostatistics analysis模块的普通克里格法结合半方差模型,完成最优内插过程。

    应用GS+9.0软件完成9项元素半变异函数的计算和多种理论模型的拟合。半变异函数作为描述随机场和随机过程空间相关性的重要统计量,是分析土壤指标空间变异结构性和随机性的有效性方法[18]。该函数如下:

    γ(h)=12N(h)i=1N(h)[Z(xi)Z(xi+h)]2 (1)

    式中:γ(h)——半变异函数值;

    N(h)——样点对的数目;

    h——步长,即两分隔样点的距离;

    Z(xi)—Z(xi)在空间位置xi上的土壤观测值;

    Z(xi+h)—Z(xi)在空间位置xi+h上的土壤观测值。

    在土壤地球化学工作中,通常区域化变量往往在不同的方向上会表现出不同的变异性,即便是在同一方向上也存在不同尺度的多维度的变异性。为了精准反映区域化变量的主要空间变异性,需要根据半变异函数的决定系数(R2)和残差(RSS)按照土壤元素所对应的最优理论模型进行半变异函数的拟合。本次所涉及的半方差函数模型主要有以下两种:

    指数模型(Exponential):

    \[γ(h)=C0+C[1exp(hα)],h>0γ(h)=0,h=0\] (2)

    线性模型(Linear):

    γ(h)=C0+Ch/α,h>0 (3)

    式中:C0——块金常数(间距为0时的半方差,通常由随机因素引起的变异);

    C——拱高(通常由系统因素引起的变异);

    α——变程(半方差达到基台值的样本间距,体现随机变量在空间上的自相关程度)。

    研究区内共采集土壤表层样品401件,表层土壤重金属元素是以统计学为基本依据,通过系统性计算各指标的最小值、最大值、中位数、平均值、标准差等统计特征值,来描述表层土壤重金属元素的含量区间和分布规律,见表 1。表层土壤重金属元素的中位数与平均值相差较大,表明其中心趋向分布可能被异常值影响而使其呈非标准正态分布。将各指标平均值与河北省背景值[19]对比,发现土壤中Cd、Cu和Pb的平均含量均高于河北省背景值,其中以Cd最高,是河北省背景值的2倍;As、Cr、Hg、Ni和Zn的平均含量均与河北省背景值基本一致。

    表  1  壤重金属的描述性统计
    Table  1.  Descriptive statistics of Soil heavy metals
    下载: 导出CSV 
    | 显示表格

    大量研究表明,变异系数(Cv)不仅可以反映总体土壤样本中各样点元素含量的平均变异程度,也可反映土壤元素的离散程度[20]。变异系数值越大,元素的空间分布越不均匀,离散程度越高。参照本次测试数据特点,将变异程度划分为三种类型:Cv<70%为均匀分布、70%≤Cv<100%为中强分异(中高起伏)、Cv≥100%为强分异(很大起伏)。依据该分类标准,Cd、Cu、Hg和Pb的变异系数分别为385%、143%、350%、118%,表明表层土壤中Cd、Cu、Hg和Pb属于强分异,其中以Cd和Hg的分异程度最强,表明局部地区可能受到人为活动影响,存在点源污染;As、Cr、Ni的变异系数分别为96%、93%、86%,属于中强分异;Zn的变异系数为64%,分异程度最弱。

    在置信区间为95%内,不同土壤类型中单因素方差分析结果为F=2.913、p=0.006,判断土壤类型对土壤Zn含量均值具有显著性影响(表 2);不同土地利用类型中单因素方差分析结果表明,Cr的F=3.408、p=0.003,Cu的F=3.097、p=0.006,Ni的F=2.355、p=0.030,判断土地利用类型对土壤Cr、Cu、Ni含量均值具有显著性影响(表 3)。

    表  2  9种元素在不同土壤类型中的均值特征及单方差分析结果
    Table  2.  Mean characteristics of 9 elements in different soils and the results of single variance analysis
    下载: 导出CSV 
    | 显示表格
    表  3  9种元素在不同土地利用类型的均值特征及单方差分析结果
    Table  3.  The mean characteristics of 9 elements in different land use types and the results of single variance analysis
    下载: 导出CSV 
    | 显示表格

    因此按照不同土壤类型,对Zn含量均值进行统计分析,发现潮土>栗钙土>粗骨土>褐土>棕壤>草甸土>灰色森林土>山地草甸土。按照不同土地利用类型,对Cr、Cu、Ni含量均值进行统计分析,发现Cu在内陆滩涂中含量高,Cr和Ni在建设用地中含量高,可能与人为活动影响有关。

    对土壤重金属元素和Se进行相关性分析(spearman分析),相关系数R越大,元素之间的关联程度越高;相关系数R越小,元素之间的关联程度越低或同源性越差。由各元素间的相关系数表可知(表 4),具有显著性相关关系的约占总关系数的25%,极显著相关关系约占总关系数的47.2%,其中Cr和Ni、Cd和Pb、Cd和Zn、Pb和Zn、Cu和Hg呈极显著的正相关,相关系数分别为0.832,0.698,0.739,0.731,0.635。Se与As、Hg、Pb、Zn呈极显著的正相关,与Cu和Ni呈显著的正相关,而与Cd和Cr之间的相关性较弱。Cr与As和Pb之间呈一定的负相关关系。由此推断Cr和Ni、Cd和Pb和Zn、Cu和Hg、Se和As之间可能具有相同的来源。

    表  4  9个土壤元素的相关性分析
    Table  4.  Statistical correlation analysis among 9 kinds of soil elements
    下载: 导出CSV 
    | 显示表格

    通过对研究区表层土壤中9种元素数据进行聚类分析,树状图结果表明,Cr和Ni最先聚合为一类,随后Cd、Pb、Zn聚合为第一类,Cu和Hg聚合为第二类,而As和Se聚合为第三类,如图 2所示。

    图  2  土壤元素树状图
    Figure  2.  Dendrogram of the soil elements

    在一定程度上通常采用特征值来表征主成分影响力度大小,本次提取特征值大于1的主成分。由主成分分析结果可知(表 5),表层土壤中9种元素可被划分为四个主成分,第一个主成分(F1)的方差贡献率为34.50%,占四种主成分中较大比例,可以作为衡量该区重金属状况的一个综合性指标。第二个主成分(F2)的方差贡献率为20.55%,第三个主成分(F3)的方差贡献率为14.25%,第四个主成分(F4)的方差贡献率为12.58%,这四个成分的累积方差贡献率为81.87%,具备原始数据的大部分信息。

    表  5  土壤元素的主成分分析结果
    Table  5.  Principal component analysis results of soil elements
    下载: 导出CSV 
    | 显示表格

    为了更加精准确定各主成分所包含的元素信息,经过方差最大正交旋转后,对9种元素的四种主成分上因子载荷进行统计,如图 3。结果发现,Cd、Pb、Zn在第一主成分中为高载荷,而在第二主成分、第三主成分、第四主成分中均为低载荷;Cr和Ni在第二主成分中为高载荷,而在第一主成分、第三主成分、第四主成分中均为低载荷;Cu和Hg在第三主成分中为高载荷,而在第一主成分、第二主成分、第四主成分中均为低载荷;As和Se在第四主成分中为高载荷,而在第一主成分、第二主成分、第三主成分中均为低载荷,见表 6

    图  3  表层土壤元素主成分载荷
    Figure  3.  Factor loading analysis of elements in surficiai soils
    表  6  主成分分析旋转后的成分载荷矩阵
    Table  6.  Rotated component matrix of PCA (principal component analysis)
    下载: 导出CSV 
    | 显示表格

    综上所述,传统统计学中主成分分析与相关性分析和聚类分析的结果在一定程度上保持一致,即Cr和Ni,Cd、Pb、Zn三者,Cu和Hg,As和Se来源可能相同。

    表层土壤一般容易受到外界因素(气候条件、地形地貌、人为活动等[21])的共同影响,因此区域性的表层土壤空间分布具有显著的趋势性和异向性。目前可通过Arcgis 10.2中趋势分析生成的三维趋势面曲线形状以及其变化的陡峭程度来判断空间数据的总体变化趋势。一般分为三种类型,即无(没有趋势效应)、一阶(区域化变量沿一定方向呈直线变化)和二阶(区域化变量沿一定方向呈多项式变化)[22],通过异向性轴向自动搜索模块,对各元素数据进行处理,从而确定各向异性的趋势效应特征。以研究区As、Cu、Cd的空间分布趋势效应分析图为例,分别代表了无、一阶、二阶趋势效应(图 4),图中东西方向由x轴表示,南北方向由y轴表示,每个样点实测值的大小由z轴表示;正北方向投影面上浅绿色曲线表示东西向的全局性趋势效应变化,正东向投影面上紫色曲线表示南北向的全局性趋势效应变化。从图上可以看出,研究区土壤中Cd表现为在东西方向上先减小后增大的趋势,Cu表现为在东西方向和南北方向上呈直线分布,而As表现为无阶效应。

    图  4  土壤元素趋势分析图
    Figure  4.  Trend analysis of soil elements

    在Kringing插值过程中,选用不同的趋势类型会导致插值过程中产生的误差不同,因此需要综合分析各类参数特点,最大程度将各类误差降为最低,进而使半方差函数模型及其参数最适合。本次研究参照以下标准判断,平均误差(mean error,ME)的绝对值最接近0,均方根误差(root-mean-square error,RMSE)与平均标准误差(average standard error,ASE)应尽量相等且值越小越优,若前者小于后者,表明低估了预测值,反之则说明高估了预测值。标准化均方根误差(root-mean-square standardized error,RMSSE)值最接近1,如果标准化均方根误差小于1,说明高估了预测值,反之为低估预测值[23],标准化平均误差(mean standardized error,MSE)的绝对值应最接近0。由表 7可知,As、Cr、Pb、Zn、Ni、Se的无趋势预测与一阶趋势预测和二阶趋势预测相比,ME的绝对值更加接近0,ASE和RMSE最接近,MSE的绝对值更加接近0,RMSSE的值最接近1,因此Kringing插值时应选择无趋势更加精准,同理Hg和Cu选择一阶,Cd选择二阶。

    表  7  不同趋势阶数插值误差比较
    Table  7.  Comparisions of semi-variogram parameters with different trends
    下载: 导出CSV 
    | 显示表格

    半变异函数(Semivariogram)是在区域化变量满足二阶平稳和本征假设的前提下,按照一定的抽样间隔获得样本方差的数学期望值,是地统计学中分析空间格局变量时最主要的应用工具。通常应用半变异函数及其参数可以对区域化变量的分布进行结构性和随机性的判断,进而结合克里格法完成对无样品空白区的插值和预测[24]。有三个相应的参数可直接反应空间变异性程度,块金值(C0)是在取样距离为0时,半变异函数为一个非0定值,反映的是最小抽样尺度以下由仪器测量误差等随机因素引起的微变异;变程(A0)也称半变异函数达到基台值时空间最大间隔距离,反映了变量空间自相关范围;块金系数C0/(C+C0)是块金值与基台值的比值,反应了空间异质性。而空间变异性程度主要由随机性变异和结构性变异大小所决定,结构性因素会增加变量的空间相关性,随机性因素则会降低变量的空间相关性,当块金系数小于25%时,表明区域变量的空间相关性强烈,此时空间变异主要由地形地貌、土壤类型、土地利用方式等结构性因素控制[25-26];若块金系数大于75%,则表明区域变量几乎不具有或具有较弱的空间相关性,其空间变异程度受随机性因子控制比重较高;当块金系数介于25%~75%,表明区域变量的空间相关性强度居中,其空间变异由随机性因素与结构性因素共同控制[27]。土壤表层重金属元素空间变异性研究的关键是半变异函数模型的确定,模型拟合的效果可按照R2值最大和RSS值最小的原则检验[28]

    表 8可知,9种土壤元素的R2值除了Pb为0.176,其余均在0.334~0.773之间,说明上述重金属元素和Se的理论半变异与实验半变异模型较为吻合,拟合效果趋优。这9种元素中有8项元素的变程远远超过了8 000 m,最小的变程也在3 000 m以上,因此可对其进行Kriging插值。

    表  8  土壤元素含量变异函数理论模型及其相关参数
    Table  8.  Soil elements content and variogram model parameters
    下载: 导出CSV 
    | 显示表格

    As的半方差函数属于指数模型,块金系数为10.8%,变程为3 000 m,因此As具强烈的空间相关性,反映出研究区表层土壤中As主要受到地形地貌、成土母质、地质背景等结构性因素的影响,这与Chen等[29]和何厅厅等[30]结论保持一致,他们认为As的空间分布主控因素为土壤属性和地质背景,受人为活动影响较小;Cr、Cd、Hg、Pb、Zn、Cu、Ni、Se的块金系数分别为83.7%、83.0%、87.4%、100%、81.7%、89.6%、85.7%、86.9%,且半方差函数均属于线性模型(图 5)。

    图  5  土壤元素半方差函数图
    Figure  5.  Semivariogram of the soil elements

    这8项元素具较弱的空间相关性,表明研究区表层土壤中这8项元素主要受到人为活动等随机性因素的影响,结论与Facchinelli等的[31]一致,其中Cu、Zn和Pb的来源与人类活动较为密切。不存在中等强度的空间变异性情况,这可能也与本次研究区边缘地带样品点较为稀疏有一定关系。

    GIS在处理海量数据和图形视觉化表达方面具有很突出的特点,选取模型插值成图后可以直观地反映土壤元素的空间变异性。空间插值法能够使用一定量的样本量对区域范围内空白样本进行预测,从而分析整个区域范围的相关特征。空间插值的方法较多,主要分为确定性插值法(反距离加权插值法、径向基函数插值法、局部多项式插值等[32-33])和地统计学插值法(克里格法/协同克里格法、面插值法、经验贝叶斯克里金法等[34-35])。克里格法又称空间局部插值法,是在一定区域内对区域化变量进行无偏最优估计的一种方法,可分为泛克里格插值、简单克里格插值等,由于满足内蕴假设等条件[36],本文采用地统计学中常用的普通克里格空间插值法[37]

    选取最优趋势参数和半方差函数模型进行普通克里格空间插值,得到研究区土壤元素的空间分布图(图 6)。参照《土地质量地球化学评价规范》[38]和河北省土壤元素背景值调查[19], As含量范围为2.9~29.9 mg/kg,整体含量值偏低;Ni含量范围为12.7~61.9 mg/kg,Pb含量范围为16.4~38.6 mg/kg,Se含量范围为0.055~0.401 mg/kg,整体含量值处于中等;Cd含量范围为0.082~0.742 mg/kg,Cr含量范围为28.0~228.4 mg/kg,Cu含量范围为11.3~142.4 mg/kg,Hg含量范围0.013~0.126 mg/kg,Zn含量范围为43.5~199.4 mg/kg,整体含量值偏高。从区域空间分布来看,承德市北部地区9种元素含量整体偏低,而南部地区含量整体偏高,体现出坝上高原地区整体受人为活动影响较低的特点。Cr和Ni空间分布具有很强的一致性,以丰宁县和隆化县为界,低值区主要分布于承德地区的北部,南部鹰手营子矿区和宽城满族自治县附近有低值区的零星分布;Cu和Hg保持了一定的空间分布一致性,高值区主要分布于承德南部地区,而Hg在丰宁县北部地区亦有高值区;Pb在中部形成了一条高值带,这又与Cd分布特征相似,这与区内铁矿、金矿等金属矿床的集中分布状况具有一定相关性。除去Pb中部的高值带,Zn和Pb二者在承德南部以及中部地区的高值区近似分布, Ni的高值区较好地表现出由南向北逐渐递变的过度特性,可能与区内金属矿产分布状况有一定关系;As低值区主要分布于承德中部隆化县-滦平县地区,围场满族蒙古族自治县以北、兴隆县以西地区。虽然上文统计学分析并未提及Se和Cd具有某种同源关系,但是二者高值区都分布在承德南部地区,西北部地区亦呈现局部高值,滦平县-双滦区、宽城满族自治县达到最高值,而西南部边缘地带和北部地区呈现低值区,也表现出弱空间异构性,这可能与Cd和Se都具有亲硫性和亲生物性,土壤中Cd、Se经常共生在一起有关,因此为富硒土地资源的安全开发利用方面带来了严峻挑战。

    图  6  土壤元素的空间分布
    Figure  6.  Spatial distribution of soil elements

    综上所述,单凭某一种方法去判断元素之间的同源关系是不准确的,需要在传统统计学结论的基础上进一步参考变异函数模型并叠加克里格插值结果进行优化。综合分析认为,研究区Cr和Ni、Cu和Hg、Zn和Pb、Cd和Se之间具有很好的空间分布相关性,来源具有一致性。建议应在现有土地利用现状属性基础上,充分考虑土壤重金属元素空间分布格局要素,在有矿业活动且土壤呈现重金属高值地区(如滦平县-双滦区北部山区)进行周期性土壤监测。同时,进一步开展区内土壤Cd和Se两者定量化控制因素的厘定及模型机理的研究,从而实现农作物富硒降镉的目标。

    (1)传统统计学分析中Cd、Cu、Hg和Pb属于强分异,变异系数分别为385%、143%、350%、118%,As、Cr和Ni属于中强分异,变异系数分别为96%、93%、86%,Zn属于弱分异,变异系数为64%。Zn含量均值受土壤类型影响显著,Cr、Cu、Ni含量均值则受土地利用类型影响显著。结合相关性分析、聚类分析和主成分结果,可知Cr和Ni、Cd和Pb和Zn、Cu和Hg来源相同。

    (2)通过异向性轴向自动搜索功能结合不同趋势阶数插值误差综合对比,确定了As、Cr、Pb、Zn、Ni、Se适宜选择无趋势,Hg和Cu适宜选择一阶,而Cd适宜选择二阶。

    (3)地统计学分析结果表明,As采用指数模型可较好地拟合,Cd、Cr、Cu、Hg、Ni、Pb、Zn、Se采用线性模型可较好地拟合。其中As的块金效应小于25%,因此具强烈的空间相关性,受结构性因素影响较大;Cd、Cr、Cu、Hg、Ni、Pb、Zn、Se的块金效应大于75%,具较弱的空间相关性,受随机性因素影响较大。

    (4)普通克里金插值图直观地反映了承德地区土壤重金属的空间结构特征。研究区土壤的9种重金属元素都呈现明显的北低南高趋势,Cr和Ni、Cu和Hg、Zn和Pb、Se和Cd之间高值区分布具有很强的空间一致性,Hg在丰宁县北部地区亦有高值区,Zn和Pb高值区都分布于承德南部有向中部地区发展的趋势, Pb在中部地区形成了一条较宽的高值带,而其余分布特征与Cd相似,Ni的高值区较好地表现出由南向北逐渐递变的过度特性,在一定程度上较好地印证了地统计学所得出的规律。

  • 图  1   研究区土壤样点分布

    Figure  1.   Distribution of soil sampling sites in the study area

    图  2   土壤元素树状图

    Figure  2.   Dendrogram of the soil elements

    图  3   表层土壤元素主成分载荷

    Figure  3.   Factor loading analysis of elements in surficiai soils

    图  4   土壤元素趋势分析图

    Figure  4.   Trend analysis of soil elements

    图  5   土壤元素半方差函数图

    Figure  5.   Semivariogram of the soil elements

    图  6   土壤元素的空间分布

    Figure  6.   Spatial distribution of soil elements

    表  1   壤重金属的描述性统计

    Table  1   Descriptive statistics of Soil heavy metals

    下载: 导出CSV

    表  2   9种元素在不同土壤类型中的均值特征及单方差分析结果

    Table  2   Mean characteristics of 9 elements in different soils and the results of single variance analysis

    下载: 导出CSV

    表  3   9种元素在不同土地利用类型的均值特征及单方差分析结果

    Table  3   The mean characteristics of 9 elements in different land use types and the results of single variance analysis

    下载: 导出CSV

    表  4   9个土壤元素的相关性分析

    Table  4   Statistical correlation analysis among 9 kinds of soil elements

    下载: 导出CSV

    表  5   土壤元素的主成分分析结果

    Table  5   Principal component analysis results of soil elements

    下载: 导出CSV

    表  6   主成分分析旋转后的成分载荷矩阵

    Table  6   Rotated component matrix of PCA (principal component analysis)

    下载: 导出CSV

    表  7   不同趋势阶数插值误差比较

    Table  7   Comparisions of semi-variogram parameters with different trends

    下载: 导出CSV

    表  8   土壤元素含量变异函数理论模型及其相关参数

    Table  8   Soil elements content and variogram model parameters

    下载: 导出CSV
  • [1]

    TATENO R, TAKEDA H. Forest structure and tree species distribution in relation to topography-mediated heterogencity of soil nitrogen and light at the forest floor[J]. Ecological Research, 2003, 18(5): 559-571.

    [2] 张志坚, 刘苑秋, 吴春生, 等. 基于地统计学和GIS的江西省森林土壤养分空间分布特征[J]. 水土保持研究, 2018, 25(1): 38-46.

    ZHANG Z J, LIU Y Q, WU C S, et al. Spatial distribution characteristics of forest soil nutrients in Jiangxi Province based on geostatistics and GIS[J]. Research of Soil and Water Conservation, 2018, 25(1): 38-46. (in Chinese)

    [3] 甘国娟, 刘伟, 邱亚群, 等. 湘中某冶炼区农田土壤重金属污染及生态风险评价[J]. 环境化学, 2013, 32(1): 132-138.

    GAN G J, LIU W, QIU Y Q, et al. Heavy metal pollution and ecological risk assessment of the paddy soils in a smelting area in Central Hunan[J]. Environmental Chemistry, 2013, 32(1): 132-138. (in Chinese)

    [4]

    ZHAO Y C, SHI X Z, YU D S, et al. Soil organic carbon density in Hebei Province, China: estimates and uncertainty[J]. Pedosphere, 2005, 15: 293-300.

    [5] 崔萌, 孙向阳, 李素艳, 等. 北京市桃主产区土壤重金属空间结构特征及来源[J]. 福建农林大学学报(自然科学版), 2019, 48(2): 238-243.

    CUI M, SUN X Y, LI S Y, et al. Spatial structure characteristics and origins of soil heavy metals in the main producing area for peach in Beijing[J]. Journal of Fujian Agriculture and Forestry University(Natural Science Edition), 2019, 48(2): 238-243. (in Chinese)

    [6] 汪璇, 王成秋, 唐将, 等. 基于地统计学和GIS的三峡库区土壤微量营养元素空间变异性研究[J]. 土壤通报, 2009, 40(2): 359-365.

    WANG X, WANG C Q, TANG J, et al. Geostatistics and GIS-based spatial distribution of microelements in the Three Gorges Reservoir area[J]. Chinese Journal of Soil Science, 2009, 40(2): 359-365. (in Chinese)

    [7]

    ROGER A, LIBOHOVA Z, ROSSIER N, et al. Spatial variability of soil phosphorus in the Fribourg canton, Switzerland[J]. Geoderma, 2014, 217: 26-36.

    [8] 安永龙, 黄勇, 刘清俊, 等. 北京城区表层土壤多元素分布特征及重金属元素污染评价[J]. 地质通报, 2016, 35(12): 2111-2120.

    AN Y L, HUANG Y, LIU Q J, et al. The distribution of surface soil elements and the pollution assessment of heavy metal elements in Beijing[J]. Geological Bulletin of China, 2016, 35(12): 2111-2120. (in Chinese)

    [9] 邢洪连, 郭华明, 王轶, 等. 河北保定市安新-清苑县土壤重金属形态分布及风险评估[J]. 水文地质工程地质, 2016, 43(2): 140-146.

    XING H L, GUO H M, WANG Y, et al. Fraction distribution and risk assessment of soil heavy metals in Anxin-Qingyuan County in Baoding of Hebei[J]. Hydrogeology & Engineering Geology, 2016, 43(2): 140-146. (in Chinese)

    [10] 孙厚云, 卫晓锋, 甘凤伟, 等. 承德市滦河流域土壤重金属地球化学基线厘定及其累积特征[J]. 环境科学, 2019, 40(8): 3753-3763.

    SUN H Y, WEI X F, GAN F W, et al. Determination of heavy metal geochemical baseline values and its accumulation in soils of the Luanhe river basin, Chengde[J]. Environmental Science, 2019, 40(8): 3753-3763. (in Chinese)

    [11] 王梦雨, 张静. 承德市柳河流域农田土壤重金属含量调查与评价[J]. 湖北农业科学, 2018, 57(1): 27-28.

    WANG M Y, ZHANG J. Survey and assessment on heavy metals in farmland soil along Liuhe river in Chengde city[J]. Hubei Agricultural Sciences, 2018, 57(1): 27-28. (in Chinese)

    [12] 殷志强, 卫晓锋, 刘文波, 等. 承德自然资源综合地质调查工程进展与主要成果[J]. 中国地质调查, 2020, 7(3): 1-12.

    YIN Z Q, WEI X F, LIU W B, et al. Progresses and main achievements of comprehensive geological survey project of natural resources in Chengde[J]. Geological Survey of China, 2020, 7(3): 1-12. (in Chinese)

    [13]

    YU C J, LI H Q, JIA X P, et al. Improving resource utilization efficiency in China’s mineral resource-based cities: a case study of Chengde, Hebei province[J]. Resources Conservation and Recycling, 2015, 94: 1-10.

    [14] 李小曼, 刘勤, 徐梦洁, 等. 苏南村镇土壤重金属空间变异性研究[J]. 土壤通报, 2016, 47(1): 179-185.

    LI X M, LIU Q, XU M J, et al. Spatial variability of heavy metal contents in towns of southern Jiangsu Province[J]. Chinese Journal of Soil Science, 2016, 47(1): 179-185. (in Chinese)

    [15] 马兴华, 张晋昕. 数值变量正态性检验常用方法的对比[J]. 循证医学, 2014, 14(2): 123-128.

    MA X H, ZHANG J X. The comparison among the common normality tests for numerical variables[J]. The Journal of Evidence-Based Medicine, 2014, 14(2): 123-128. (in Chinese)

    [16] 陶吉兴, 傅伟军, 姜培坤, 等. 基于Moran’s I和地统计学的浙江森林土壤有机碳空间分布研究[J]. 南京林业大学学报(自然科学版), 2014, 38(5): 97-101.

    TAO J X, FU W J, JIANG P K, et al. Using Morans I and geostatistics to analyze the spatial distribution of organic carbon in forest soil of Zhejiang province[J]. Journal of Nanjing Forestry University(Natural Science Edition), 2014, 38(5): 97-101. (in Chinese)

    [17] 宋旭, 高灯州, 曾从盛, 等. 基于GIS地统计学的海坛岛农田养分变异研究[J]. 实验室科学, 2017, 20(1): 25-28.

    SONG X, GAO D Z, ZENG C S, et al. Spatial variation of soil nutrient of the Haitan island based on GIS and Geo-statistic[J]. Laboratory science, 2017, 20(1): 25-28. (in Chinese)

    [18]

    CAMBARDELLA C A, MOORMAN T B, NOVAK J M, et al. Field-scale variability of soil properties in central Iowa soils[J]. Soil Science Society of America Journal, 1994, 58(5): 1501-1511.

    [19] 中国环境监测总站. 中国土壤元素背景值[M]. 北京: 中国环境科学出版社, 1990.

    Environmental Monitoring of China. Background values of soil elements in China[M]. Beijing: China Environment Science Press, 1990. (in Chinese)

    [20] 钟巧, 王勇辉, 焦黎. 夏尔希里地区土壤重金属含量特征及空间变异分析[J]. 水土保持研究, 2016, 23(3): 360-365.

    ZHONG Q, WANG Y H, JIAO L. Characteristics and spatial variability of the analysis of soil heavy metals in Xiaerxili area[J]. Research of Soil and Water Conservation, 2016, 23(3): 360-365. (in Chinese)

    [21] 周稀, 邓欧平, 潘洪旭, 等. 基于GIS的西河流域土壤氮素空间变异特征及影响因素研究[J]. 西南农业学报, 2016, 29(4): 896-902.

    ZHOU X, DENG O P, PAN H X, et al. Spatial distribution characteristics and influence factors analysis of soil nitrogen in west river valley based on GIS[J]. Southwest China Journal of Agricultural Sciences, 2016, 29(4): 896-902. (in Chinese)

    [22]

    ROBINSON T P, METTERNICHT G. Testing the performance of spatial interpolation techniques for mapping soil properties[J]. Computers and Electronics in Agriculture, 2006, 50(2): 97-108.

    [23] 王瑞, 何中青, 丁建方, 等. 洪泽湖农场土壤碱解氮含量的地统计学和GIS分析[J]. 安徽农业科学, 2011, 39(31): 19122-19126.

    WANG R, HE Z Q, DING J F, et al. Geostatistical and GIS analysis on soil Alkali-hydrolyzed nitrogen in Hongze lake farm[J]. Journal of Anhui Agriculture Sciences, 2011, 39(31): 19122-19126. (in Chinese)

    [24] 安永龙, 杜子图 , 黄勇. 基于地统计学和GIS技术的北京市大兴区礼贤镇土壤养分空间变异性研究[J]. 现代地质, 2018, 32(6): 1311-1321.

    AN Y L, DU Z T, HUANG Y. Spatial variation analysis of soil nutrients in Lixian town of Daxing district in Beijing based on geostatistics and GIS[J]. Geoscience, 2018, 32(6): 1311-1321. (in Chinese)

    [25]

    CHIEN Y J, LEE D Y, GUO H Y, et al. Geostatistical analysis of soil properties of mid-west Taiwan soils[J]. Soil Science, 1997, 162(4): 291-298.

    [26] 曹祥会, 龙怀玉, 周脚根, 等. 河北省表层土壤有机碳和全氮空间变异特征性及影响因子分析[J]. 植物营养与肥料学报, 2016, 22(4): 937-948.

    CAO X H, LONG H Y, ZHOU J G, et al. Analysis of spatial variability and influencing factors of topsoil organic carbon and total nitrogen in Hebei Province[J]. Plant Nutrition and Fertilizer Science, 2016, 22(4): 937-948. (in Chinese)

    [27]

    LI J, MIN Q W, LI W H, et al. Spatial variability analysis of soil nutrients based on GIS and geostatistics: A case study of Yisa Township, Yunnan, China[J]. Journal of Resources and Ecology, 2014, 5(4): 348-355.

    [28]

    WANG Y Q, SHAO M A, LIU Z P, et al. Regional spatial pattern of deep soil water content and its influencing factors[J]. Hydrological Sciences Journal, 2012, 57(2): 265-281.

    [29]

    CHEN M, MA L Q, HOOGEWEG C G, et al. Arsenic background concentrations in Florida, USA surface soils: determination and interpretation[J]. Environmental Forensics, 2001, 2(2): 117-126.

    [30] 何厅厅, 赵艳玲, 王亚云, 等. 矿区农田土壤重金属的空间变异: 以陕西省某金矿区为例田[J]. 贵州农业科学, 2013, 41(10): 190-193.

    HE T T, ZHAO Y L, WANG Y Y, et al. Spatial variability of farmland soil heavy metal in the mining area: taking a gold mine in Shaanxi as a case[J]. Guizhou Agricultural Sciences, 2013, 41(10): 190-193. (in Chinese)

    [31]

    FACCHINELLI A, SACCHI E, MALLEN L. Multivariate statistical and GIS-based approach to identify heavy metal sources in soils[J]. Environmental Pollution(Barking, Essex), 2001, 114(3): 313-324.

    [32] 王兴, 刘莹, 王春晖, 等. 海洋盐度分布的插值方法应用与对比研究[J]. 海洋通报, 2016, 35(3): 324-330.

    WANG X, LIU Y, WANG C H, et al. Study on the application and comparison of interpolation methods for the marine salinity distribution[J]. Marine Science Bulletin, 2016, 35(3): 324-330. (in Chinese)

    [33] 张优, 王娟, 张杰, 等. GIS与地统计学的土壤水分空间插值方法[J]. 四川师范大学学报(自然科学版), 2019, 42(5): 703-710.

    ZHANG Y, WANG J, ZHANG J, et al. Study on interpolation method of soil moisture based on GIS and statistical models[J]. Journal of Sichuan Normal University(Natural Science), 2019, 42(5): 703-710. (in Chinese)

    [34] 史文娇, 岳天祥, 石晓丽, 等. 土壤连续属性空间插值方法及其精度的研究进展[J]. 自然资源学报, 2012, 27(1): 163-175.

    SHI W J, YUE T X, SHI X L, et al. Research progress in soil property interpolators and their accuracy[J]. Journal of Natural Resources, 2012, 27(1): 163-175. (in Chinese)

    [35] 陈思萱, 邹滨, 汤景文. 空间插值方法对土壤重金属污染格局识别的影响[J]. 测绘科学, 2015, 40(1): 63-67.

    CHEN S X, ZOU B, TANG J W. Impact of spatial interpolation methods on identifying structure of heavy metal polluted soil[J]. Science of Surveying and Mapping, 2015, 40(1): 63-67. (in Chinese)

    [36] 王艳妮, 谢金梅, 郭祥. ArcGIS中的地统计克里格插值法及其应用[J]. 软件导刊, 2008, 7(12): 36-38.

    WANG Y N, XIE J M, GUO X. Application of Geostatistical Interpolation Method in Arcgis.[J]. Software Guide, 2008, 7(12): 36-38. (in Chinese)

    [37] 杨全合, 安永龙. 基于地统计学和GIS的通州区于家务乡土壤肥力综合评价[J]. 西南农业学报, 2019, 32(4): 882-891.

    YANG Q H, AN Y L. Comprehensive evaluation of soil fertility in Yujiawu town of Tongzhou district using geostatistics and GIS southwest China[J]. Southwest China Journal of Agricultural Sciences, 2019, 32(4): 882-891. (in Chinese)

    [38] 中国人民共和国国土资源部. 土地质量地球化学评价规范: DZ/T 0295-2016[S]. 北京: 中国标准出版社, 2016.

    Ministry of Land and Resources of the People’S Republic of China.Specification of land quality geochemical assessment:DZ/T 0295—2016 [S].Beijing:Standards Press of China,2016.(in Chinese)

  • 期刊类型引用(15)

    1. 孙厚云,马峰,陈自然,朱喜,卫晓锋. 高强度运输活动影响下承德钒钛磁铁矿集区土壤重金属生态风险与来源解析. 环境科学. 2025(04): 2486-2500 . 百度学术
    2. 向峰,李中志,熊熙,李斌勇. 粒子群局部优化的反距离权重插值算法. 计算机应用. 2023(02): 385-390 . 百度学术
    3. 刘琼,杨楠,高萌萌,李小磊,王轶,郅二铨,张旭航. 面向国土空间规划的农业生产适宜性评价——以辽中平原为例. 中国农业资源与区划. 2023(02): 183-193 . 百度学术
    4. 安永龙,殷秀兰,金爱芳,李文娟,鲁青原. 河北张家口市某种植区土壤营养元素生态化学计量与空间分异特征. 地质通报. 2023(Z1): 443-459 . 百度学术
    5. 黄丹,黄勇,安永龙,冯辉,孙朝,李欢. 北运河流域(北京段)表层土壤多环芳烃空间分布特征及来源解析. 水文地质工程地质. 2023(03): 159-171 . 本站查看
    6. 安永龙,殷秀兰,李文娟,金爱芳,鲁青原. 张家口市万全区某种植区土壤重金属污染评价与来源分析. 环境科学. 2023(06): 3544-3561 . 百度学术
    7. 史帅航,白甲林,余洋. 西南地区某矿产集采区土壤重金属迁移规律及生态风险评价. 金属矿山. 2022(02): 194-200 . 百度学术
    8. 邹友琴,李群,李卓悦,刘耀驰,李勇丽,章萍. Fe-Mn二元氧化物强化活性铝硅酸盐矿物对Tl(Ⅰ)的去除及其机制. 吉林大学学报(地球科学版). 2022(03): 992-1003 . 百度学术
    9. 王春光,彭景颂,李婉莹,李健,张双腾,余洋. 黄河中游重要生产煤矿山土壤重金属特征分析与评价. 矿业安全与环保. 2022(05): 124-130 . 百度学术
    10. 邵海,殷志强,王轶,邢博,彭令,王瑞丰. 河北坝上高原如意河流域风积沙厚度空间展布预测方法. 地质通报. 2022(12): 2138-2145 . 百度学术
    11. 崔潇,周妍如,刘孝阳,白中科. 平朔露天煤矿复垦区不同地质层组岩土质量综合评价. 水文地质工程地质. 2021(02): 164-173 . 本站查看
    12. 董燕,孙璐,李海涛,张作辰,张源,李刚,郭小彪. 雄安新区土壤重金属和砷元素空间分布特征及源解析. 水文地质工程地质. 2021(03): 172-181 . 本站查看
    13. 陈自然,孙厚云,卫晓锋,黄行凯,李多杰,张晓敏,郭颖超. 承德黑山钒钛磁铁矿矿集区土壤重金属空间结构特征与生态风险评价. 矿产勘查. 2021(04): 1019-1029 . 百度学术
    14. 孙厚云,卫晓锋,张晓敏,贾凤超,李多杰,刘卫,李健,陈自然. 河北承德中部伊逊河红旗地区土壤生源要素空间分布格局及其影响因素. 矿产勘查. 2021(04): 1008-1018 . 百度学术
    15. 杨树明,余小芬,邹炳礼,解燕,刘加红,王瑞宝,吕亚琼,蔡永占,张素华,邱学礼. 曲靖植烟土壤pH和主要养分空间变异特征及其影响因素. 土壤. 2021(06): 1299-1308 . 百度学术

    其他类型引用(5)

图(6)  /  表(8)
计量
  • 文章访问数:  431
  • HTML全文浏览量:  171
  • PDF下载量:  349
  • 被引次数: 20
出版历程
  • 收稿日期:  2020-07-30
  • 修回日期:  2020-09-03
  • 网络出版日期:  2020-12-09
  • 刊出日期:  2020-10-31

目录

/

返回文章
返回