中国农业科学 2005,38(6):1204-1212 Scientia Agricultura Sinica
洞庭湖平原区土壤全磷含量地统计学和GIS分析
路 鹏
1,2
, 彭佩钦, 宋变兰, 唐国勇, 邹 焱, 黄道友, 肖和艾, 吴金水
1111111,2
, 苏以荣
1
2
(1中国科学院亚热带农业生态研究所，长沙 410125；中国科学院水利部水土保持研究所黄土高原土壤侵蚀与旱地农业国家重点实验室，杨凌 712100)
摘要：以洞庭湖平原区典型景观单元为试点，利用GPS定位共取得651个耕层（0～20 cm）土壤样品。通过对数转化、稳健统计、域法处理和Box-Cox转化4种数据处理方法对土壤全磷进行了正态分布性检验。结果表明，Box-Cox转化成功地使数据集服从正态分布并消弱了异常值的影响，应用地统计学分析方法进行了实验变异函数的计算和最适合模型的拟合，得出土壤全磷最好的理论模型为球状模型，随后用普通克立格估值方法绘制了土壤全磷的空间分布图。在经过不同趋势阶数土壤全磷克立格（Kriging）插值误差的综合比较的基础上，结果表明趋势效应参数宜选取二阶。Kriging 估值标准差（KSD）被认为是内插象素值的标准差，并为提高全磷的制图精度提供了有用的信息。应用地统计学和地理信息系统（GIS）的概率克立格法研究了洞庭湖平原典型区的土壤全磷的空间分布并进行了风险性评价。结果表明，土壤全磷不同含量水平下的概率分布图对风险性评价、合理施肥和控制磷素的面源污染的管理实践将是十分有益的。
关键词：土壤全磷；地统计学；空间变异性；概率克立格；洞庭湖平原区
Geostatistical and GIS Analyses on Soil Total P in the Typical Area
of Dongting Lake Plain
LU Peng1,2, PENG Pei-qin1, SONG Bian-lan1, TANG Guo-yong1, ZOU Yan1, HUANG Dao-you1,
XIAO He-ai1, WU Jin-shui1,2, SU Yi-rong1
(1The Institute of Subtropical Agriculture, Chinese Academy of Sciences, Changsha 410125; 2 State Key Laboratory of
Soil Erosion and Dry Land Farming on the Loess Plateau, Institute of Soil and Water Conservation,
Chinese Academy of Sciences and Ministry of Water Resources, Yangling 712100)
Abstract: A typical landscape unit of Dongting Lake plain was selected as an experimental site. Approximate grid approaches were employed for the sampling scenario in 2004 with 651 Global Position System (GPS) established spots sampled in topsoil (0-20 cm). The purpose of this study is to evaluate various data processing methods, including logarithmic transformation, robust statistics, excluding outliers and Box-Cox transformation for evaluating soil total P content with normal distribution. The result showed that Box-Cox transformation was applied in order to achieve normality in the data set and to dampen the effect of outliers. Geostatistical analyses were carried out, including calculation of experimental variograms and model fitting. The best theoretical model for semivariogram of soil total P were spherical model. The ordinary kriging estimates of soil total P concentration were mapped. The integrative comparisons of semivariogram parameters with different trends of the kriging prediction errors of soil total P indicated that the 2-order trend effect was preferable. Kriging standard deviations (KSD) were regarded as the standard deviations of the interpolated pixel values and provided valuable information for which will increase the accuracy of the total P mapping. Spatial distribution and hazard assesment of soil total P in the typical area of Dongting Lake plain were investigated using geostatistics and geographic information system (GIS) techniques such as probability kriging on the basis of the software ArcGIS Desktop. Probability distribution of soil total P at different levels will be helpful to conduct hazard assessment, optimal fertilization and develop manage
收稿日期：2005-04-04
基金项目: 中国科学院知识创新工程重要方向项目(KZCX3-SW-426)、国家自然科学基金重点项目(40235057)、国家重点基础研究发展规划项目(2002
CB412503)
作者简介: 路 鹏 (1976-)，男，博士研究生，主要从事土壤资源信息管理与环境方面的研究。E-mail: [email protected]。苏以荣为通讯作者，Tel:
0731-4615222; E-mail: [email protected]
6期 路 鹏等：洞庭湖平原区土壤全磷含量地统计学和GIS分析 1205
ment practices to control the non-point sources of P pollution.
Key words: Soil total P; Geostatistics; Spatial variability; Probability kriging ; Dongting Lake plain area
洞庭湖区位于湖南省的北部，其地处东经111˚40′～113˚10′，北纬27˚55′～30˚22′。尤其洞庭湖平原区自然条件优越，土地肥沃，是中国南方有名的“鱼米之乡”。但是，自从20世纪80年代以来，忽视有机肥的施用，盲目大量滥施化肥，已使整个地区的耕地质量明显下降，对农业生产带来不利影响[1]。由于磷肥的施用量逐年增加，土壤全磷含量总体水平也在不断上升，超过了作物生长的需求并且普遍出现盈余现象。这一方面会导致施肥经济效益的显著下降，另一方面也会加大土壤磷素向水体流失的风险，并对流域水环境质量构成严重威胁。因此，磷肥的合理施用及其生态环境效应已受到人们的广泛关注[2]。
近年来，农业面源污染问题受到政府和研究者的不断重视[3]。磷素是农业面源污染和引发湖泊发生富营养化的主要因子之一。因此，致力于研究区域尺度，特别是景观尺度土壤磷素的空间分布特征就成了当前环境与农业领域研究热点之一。了解景观尺度土壤全磷的空间变异性，无论是对于提高作物产量，还是对于面源污染的发生和发展以及减少环境的潜在危害都具有十分重要的现实意义。
地统计学方法是一种最优的空间插值方法，它和GIS已经成为空间不确定性和风险性评价非常有用的研究工具[4]。自从20世纪70年代以来，成功地应用地统计学方法解决空间变异问题在土壤学和环境科学领域中已有较多的报道[2,3,5,6]。
目前人均耕地日趋减少，施肥模式和种植结构不太合理，土壤营养元素平衡失调，已使农业生产处于徘徊状态。因此，笔者在实地调查和采样分析的基础上，运用地统计学与GIS相结合的方法，对洞庭湖平原典型区土壤全磷的空间变异特征进行了分析，并用克立格（Kriging）插值方法绘制了土壤全磷含量分布图和风险评价图，从而指导农业生产、合理施肥，以期减少地表径流和渗漏对地下水、湖泊的污染所造成不利的影响，并为整个洞庭湖平原区农业生态环境管理提供理论基础。
1 研究区域概况与研究方法
1.1 研究区域概况
研究区域为洞庭湖平原区典型景观单元，其位于湖南省北部的沅江市，洞庭湖垸的大同乡人新与人和
村，面积约为3.7 km2。沅江市（东经112˚16′～112˚56′，北纬28˚42′～29˚11′），属中亚热带向北亚热带过渡的湿润季风湿润气候，年均气温16.5℃，年均降水量 1 313 mm，年平均相对湿度81%，降水集中且多为暴雨。地势平坦，主要地貌类型为冲积平原和湖区平原，土壤以湖积物和河流冲积物发育而成的水稻土和潮土为主，土层深厚。
1.2 土壤样品采集与化学分析
2004年3月，根据不同的土地利用方式（包括水田、旱地、新开垦的旱地），在研究区域内共采取表层（0～20 cm）土壤样品651个。用美国Garmin公司推出的eTrex Vista手持式GPS导航和定位，记录采样点的经纬度和高程，并详细记录采样点周围的景观信息。定位的样点经格式转换成ArcView能够接受的shp格式，再进行ArcInfo投影转换，生成以米（m）为单位的平面坐标（投影类型为高斯克吕格投影，中央经线为112˚00′00″E），最后生成用于地统计学分析的样点分布图（图1）。土壤样品是用网格内随机取样的方法采集，在取样点所在田块内随机采集15个点作为一个混合样，然后风干研磨，过0.25 mm筛用于样品的测定。土壤全磷采用HClO4-H2SO4 消煮，钼锑抗比色法测定。
图1 土壤样品空间分布图
Fig.1 The spatial distribution of soil samples
1.3 地统计学基本理论
关于地统计学的基本原理和方法，文献中已经有很好的解析和综述[7]，本文只作简要的说明。半变异函数是区域化变量空间自相关性的表征，它依赖于分隔两点的矢量h。假设区域化变量满足二阶平稳假设
1206 中 国 农 业 科 学 38卷
和本征假设，其半方差函数可用下式表示：
2
N(h)
γ(h)=∑Z(xi)−Z(xi+h)
i2N(h)=1
（K-S检验）是通过SPSS11.0完成的。实验半方差函数的计算、理论模型的拟合及Kriging插值和图形绘制是由软件ILWIS3.1，ArcGIS8.3共同完成的。将采样点的定位数据导入GIS软件中，获得样点分布的空间数据库，通过关键字段实现空间数据与土壤属性数据之间的联接，作为实验半方差函数的计算和Kriging插值法制图的数据源。
本文采用对数转化、稳健处理、域法处理（样本平均值X±3s）和Box-Cox转化[8] 4种数据处理方法和柯尔莫哥洛夫-斯蜜诺夫检验（K-S检验），计算了偏度、峰度和K-S P等参数，并对其是否服从正态分布进行了检验。用普通Kriging插值和概率Kriging插值方法，估值样点数最大为20，产生以5 m×5 m分辨率的栅格图形。插值的结果以ASCII数据的格式输出存储，然后生成GRID图形，再将图形用研究区域边界图形切割，最后生成研究区域相应的土壤特性含量分布图。
1
[]
式中，γ(h)为半方差函数；h为样点的空间间隔距离，称为步长；N(h)为距离等于h的点对数；Z(xi)为处于点xi处变量的实测值；Z(xi+h)为与点xi偏离h处变量的实测值。
不同形式的Kriging内插法是目前地统计学应用最广泛的最优内插法，它是根据半方差分析所提供的空间自相关程度的信息来进行插值，因此可以对未采样点给出最优无偏估计，而且能同时提供估计值的误差和精确度。也即它是利用已知点的数据去估计未知
点X0的数值，其实质是一个实行局部估计的加权平均值：
Z(Xo)=∑λiZ(Xi)
i=1
n
式中，Z(X0)是在未经观测点X0上的内插估计值，Z(Xi)是在点X0附近的若干点上获得的实测值。λi是考虑了半方差图中表示空间的权重，所以，Z(X0)值的估计应该是无偏的，因为：∑λi=1，估计偏差是最
i=1n
2 结果与分析
2.1 土壤全磷的描述性统计和正态分布性检验
全磷直方图展示了强的负偏斜效应（图2），进一步来说，这种分布和可能存在的异常值能导致在地统计分析和Kriging插值中产生计算结果的偏差。因此，检验数据的正态分布性是使用空间统计学克立格方法进行土壤特性空间分析的前提，只有当数据服从正态分布时，克立格方法才有效。
由于原始数据集的概率分布表现出强的偏斜效应，因此，在地统计和插值分析之前必须进行数据的转化。一些参量如偏度、峰度和正态性的柯尔莫哥洛夫-斯蜜诺夫检验的显著水平（K-S P）列于表1。由表1
可以得出，经过对数变化、域法处理都没有使全磷
小的，并可由下列方程求出：
σ
2
d
=b
T
λμ
式中，b是被估计点与其它点之间的半方差矩阵，μ为拉格朗日参数。 1.4 数据处理
正态分布性检验和柯尔莫哥洛夫-斯蜜诺夫检验
图2 土壤全磷含量的直方图展示了强的偏斜概率分布
Fig. 2 Histogram of total soil P concentrations (n=651) showing a heavily skewed
6期 路 鹏等：洞庭湖平原区土壤全磷含量地统计学和GIS分析 1207
表1 土壤耕层全磷的描述性统计
Table 1 Descriptive statistics of soil total P in topsoil (g·kg -1)
处理方法 Treatment method
样本数 Sample
元素 Element
平均值 Mean
标准差 S.D.
最小值 Min.
最大值 Max
偏度
峰度
Skewness Kurtosis
K-S P
原始数据Original data 651
对数转化Logarithmic transformation 651 稳健统计Robust statistics 651 域法处理Excluding outliers 651 Box-Cox 651 TP 0.99 0.13 0.48 1.34 -0.52 0.65 0.02
TP -0.01 0.06 -0.32 0.13 -1.08 2.19 0 TP 0.99 0.12 0.74 1.19 -0.38 -0.34 0.09 TP 0.99 0.13 0.60 1.34 -0.45 0.37 0.03 TP 0.01 0.13 -0.30 0.47 0.38 0.37 0.24
含量成正态分布。经Box-Cox变换后全磷通过了柯尔莫哥洛夫-斯蜜诺夫检验（K-S检验），这表明已经服从正态分布且满足地统计学分析要求。尽管稳健处理使得土壤全磷成正态分布，但效果远不如Box-Cox变化。因此，Box-Cox转化后的数据被用来随后的Kriging插值分析。
然而，表1的描述性统计分析只能说明土壤全磷样点的含量特征，但是不能完全反映整个景观单元的区域特征，也就是说不能反映土壤全磷含量的空间结构性和随机性。因此，需要用地统计学和GIS相结合的方法来弥补这一缺陷。
2.2 各向同性下土壤全磷含量的空间变异分析 根据不同空间位置上全磷含量的分析数据，计算
实际半方差值γ(h)并绘制变异函数曲线图，这是地统计学中空间变异分析的基本步骤，也是进行Kriging插值的前提。变异函数曲线图表示全磷含量的区域化变量在距离与方向上的所有成对点观测值之间的空间相关性。得出的半变异函数图的起伏特征、原点处性状、趋势走向、不同的方向差异性等形状特点提供了丰富的空间结构信息。
土壤全磷的各向同性下的变异函数（图3）展示了很好的空间结构，很好的符合球状模型。步长间隔为60 m时的块金方差为0.0078，块金方差很小表明采样密度能够充分揭示研究区域土壤全磷的空间结 构[8]。进一步来说，其变程为338 m，说明空间自相关尺度远比步长间隔60 m大。从另一个方面说明了当前的采样设计是恰当的，并且预示了好的空间结构将展示在空间插值图中。
表2中，全磷半方差函数模型的拟合度均在0.90以上，说明半方差函数模型的选取能很好地体现其空间变异性。块金值与基台值之比是反映区域化变量空间异质性程度的重要指标，也反映了空间变异的成分中区域因素（自然因素）和非区域因素（人为因素）谁占主导作用。从C0/C0+C的大小可以看出，影响全磷含量空间分布的总变异中，由非区域因素引起的变
图3 土壤全磷各向同性下变异函数图 Fig.3 Isotropic semivariogram of soil total P
异达到了45.35%，这说明非区域因素起的作用表现的比较明显，也与近20年来人们普遍重视磷肥的施用有关。因此，施肥等人为因素成为影响全磷含量空间分布的主要因子。全磷的C0/C0+C小于0.5，反映了洞庭湖平原区其在所研究的尺度上具有中等的空间自相关格局，目前虽然受到植物吸收、施肥等小尺度因素的影响，但是还没有达到破坏其原有空间格局的程度[9]。 2.3 各向异性下土壤全磷含量的空间变异分析
为了解自然过程对全磷空间变异性的影响，笔者通过半变异/协方差函数云（semivariogram/ covariance cloud）工具和异向性轴轴向自动搜寻分析
了其空间方向效应，并结合地块具有明显的条带状和对称性的特点，拟以NE25˚、NE115˚、NE70˚、NE160˚4个方向计算其空间变异函数参数值（表3）和它们的
各向异性比K(h)（图4）。
由表3可知，全磷在NE70˚、NE160˚两个方向上有相同的基台值，变程不同，具有几何异向性结构特点。在NE25˚、NE115˚两个方向上，有不同的基台值，变程也不同，具有带状异向性结构特点。从图4可知，虽然全磷的各向异性比（K(h)）在大部分尺度上均位于各向同性线之上，但是起伏不大，说明全磷有向各向同性发展的趋势。这也从另一个角度说明了人为因
1208 中 国 农 业 科 学 38卷
表2 土壤全磷的半方差模型拟合参数
Table 2 Parameters fitted by semivariogram model for soil total P
土壤属性 Soil properties 土壤全磷 soil total P
模型类型 Model S
块金值 (C0) Nugget 0.0078
基台值 (C0+C)
Sill 0.0172
块金值/基台值 C0/C0+C 0.4535
变程 (m) Range 338
R2 RSS 0.9000
0.000008
S. 球状模型 Spherical model 下同 The same as below
表3 4个方向上全磷的理论模型及相应的参数
Table 3 Theory models and corresponding parameters of soil total P in four directions
土壤属性 Soil nutrient 土壤全磷 Soil total P
方向 Direction
模型 Model
块金值 (C0)
基台值 (C0+C)
块金值/基台值 C0/C0+C
变程 Range (m)
280 320 360 260
25˚ S 0.007 0.0172 0.407 70˚ S 0.007 0.0175 0.400 115˚ S 0.009 0.0160 0.563 160˚ S 0.006 0.0175 0.343
二阶（二阶多项式变化）。图5是全磷含量的全局趋势效应分析示意图。图中，X轴表示正东方向，Y轴表示正北方向，Z轴表示各样点测定值的大小；左后投影面上曲线表示东－西向全局性的趋势效应变化情况，右后投影面上曲线表示南－北向全局性的趋势效应变化情况。
在均考虑各向异性的情况下，对分别选择无趋势、一阶趋势和二阶趋势效应参数，结合普通Kriging插值方法所造成的插值误差进行了比较（表4）。判断
图4 不同方向上土壤全磷的各向异性比
Fig.4 Anisotropic ratio of soil total P at different directions
半方差函数模型及其参数是否合适可按以下标准综合进行：平均误差（ME）的绝对值最接近于0；标准化平均误差（MSE）最接近于0；均方根误差（RMSE）越小越好；平均标准误差（ASE）与均方根误差（RMSE）最接近，如果ASE>RMSE则高估了预测值，反之如果ASE1，则低估了预测值。从表4、图5可以得出，洞庭湖平原典型区土壤全磷的趋势效应不太明显，但在两个方向均成二阶多项式变化。综合考虑分析，在以后进行普通Kriging插值时，趋势参数宜选取二阶。
素如施肥、耕作管理等措施消弱了其空间异质性程度。
因此，随后的Kriging插值分析是根据各向异性结构的变异函数理论模型，并采用块段空间局部估计的方法绘制的空间分布格局。 2.4 土壤全磷含量的趋势分析
运用有趋势分析和异向性轴轴向自动搜寻功能的ArcGIS8.3软件的地统计学模块，可以方便地获得全磷异向性分布特征参数，以及趋势效应特征。趋势效应一般分为0（没有趋势效应）、一阶（线性变化）、
表4 不同趋势阶数土壤全磷Kriging插值误差的比较
Table 4 The comparisons of semivariogram parameters with different trend of soil total P
土壤属性 Soil properties 土壤全磷 Soil total P
趋势效应 Trend effect 无 Non 一阶 1-order
二阶 2-orfer
预测误差 Prediction error
ME RMSE ASE MSE RMSSE -0.000014 -0.000016
0.1166 0.1166
0.1111 0.1110
-0.000154 -0.000178
1.046 1.046
0.000027 0.1166 0.1110 0.000145 1.046
6期 路 鹏等：洞庭湖平原区土壤全磷含量地统计学和GIS分析
1209
图5 土壤全磷的趋势效应 Fig.5 Surface trend of soil total P
2.5 土壤全磷含量的空间分布格局分析
空间异质性是空间插值研究的隐含前提，要素的非均匀空间分布才需要空间插值；空间相关性则是空间插值研究的基础，缺乏这种相关性，空间插值就成为了一种数学游戏。因此，基于经过Box-Cox转化的采样点的数据，结合普通Kriging插值分析方法，选取二阶趋势参数并考虑各向异性，随后用Box-Cox转化的逆变换过程获得了洞庭湖平原典型区的土壤全磷含量的空间分布格局图（图6）。用一个5 m×5 m分辨率的栅格，将研究区的面积分成了440列（NE115˚方向，2 200 m）和560行（NE25˚方向，2 800 m）。当然，Kriging插值的结果也受到变异函数模拟精度、样点的分布、邻近样点的选取数等的影响[10]。
图6 土壤全磷的空间分布格局
Fig.6 Spatial distribution pattern of soil total P
全磷含量的空间分布格局进一步表明土壤表层（0～20 cm）养分含量分布具有高度的空间异质性，
并决定了空间格局的存在。从插值图中可方便地了解到全磷含量斑块大小、形状及空间分布等具有显著的差异，但不论在什么方向上，养分含量由高到低的分布梯度的规律总是存在的。显然，在此基础上借鉴土壤第二次普查的养分分级标准，可以判断一定区域内养分流动与非点源污染发展趋势及潜力。因此，这种养分的空间分布特征与研究区土壤在不同空间位置上的各种物理、化学和生物过程有着重要的联系，更重要的是土壤与作物长期相互作用，如稻田生态系统养
分生物地球化学循环的结果。
从图6中可以看出，全磷的方向性效应并不太明显，空间分布总的趋势是高值斑块区呈现零星分布的特点，并没有表现出条带状的分布格局。这从另一个角度确认了其变异函数所揭示的方向不太明显的特征。全磷的高值斑块区是水田分布的区域；低值斑块区恰恰是旱地、新开垦旱地和一些鱼塘的点缀分布区。最上部的低值区为一条河堤被新开垦成的旱地，这是其表现相对低的主要原因。这也体现出了土地利用方式的不同造成了全磷含量空间分布上的差异。 2.6 土壤全磷的标准差分布图
Kriging方差是与变异函数的结构和样点的分布有关系的，并且能提供Kriging插值结果的置信水平。Kriging的标准差（KSD）是Kriging方差的平方根，被认为是内插象素值的标准差，也是样本数据的标准差和大范围取样插值结果内在不确定性的综合反映。因此，其值大小可用于评价Kriging插值精度，值越小表明Kriging插值结果越可靠。
图7是土壤全磷Kriging插值的标准差分布图，其体现了模型误差的水平，其估值标准差变化于0.1052～0.1415 g·kg-1，平均值为0.1119 g·kg-1。在采样点的附近，KSD值较小；远离采样点的区域，KSD值一般较大。围绕采样点周围的近似椭圆形的KSD值是由变异函数的椭圆形结构造成的。KSD值的大小与采样点密度存在着明显的对应关系，这是由于KSD值仅与采样点密度、采样点的布设结构及实验函数模型有关，而与采样点实测结果的大小没有直接关系[8]。
2.7 土壤全磷含量的风险性评价
为了进一步了解洞庭湖垸典型区全磷含量的空间分布特点，本文基于ArcGIS Desktop 8.3软件，采用概率克立格方法并结合Box-Cox转化的数据，选取二阶趋势参数并考虑各向异性，对不同含量水平下全磷的概率分布做了风险性评价。值得注意的是，Box-Cox转化也应当应用于临界值的设定中。
1210 中 国 农 业 科 学 38卷
1.08 km2，占整个研究区面积的28.9%。其次是0.20～0.32，面积为0.74 km2，占整个研究区面积的19.8%。低的全磷含量水平主要以小概率区间分布，这说明研究区范围内全磷的整体水平是较高的。全磷含量≤1.0 g·kg-1的平均概率为0.5109，其空间分布以0.2～0.4、0.5～0.65和0.65～0.8的概率为主（图8-b），面积分别为0.85、0.82和0.81 km2，占整个研究区面积的22.8%、21.9%和21.8%。从图8-b也可以看出，0.0～0.2和0.8～1.0这两个概率区间所占的面积最小，其余
图7 土壤全磷的Kriging标准差分布图
Fig.7 Kriging standard deviation (KSD) of soil total P
4个概率区间占的面积相当。全磷含量>1.05 g·kg-1的平均概率为0.3263(图8-c)，其空间分布以0.25～0.40较低的概率区间为主，面积为1.05 km2，占整个研究区面积的28.3%。这部分田块区域成为农业面源污染的严重区。因此，开展全磷的风险性评价研究，能从区域尺度上发现农田土壤养分管理过程中可能存在的问题，这对整个洞庭湖平原区农业生态环境管理提供理论基础，将有助于农业新一轮的大发展。
土壤全磷含量的风险性评价如图8所示，从整体上看，全磷含量≤0.9 g·kg-1的平均概率为0.2480。其概率区间以0.0～0.1的分布最广（图8-a），面积为
图8 土壤全磷不同含量水平下的概率分布
Fig.8 Probability distribution of soil total P at different levels
3 讨论
正如传统的统计学研究一样，在地统计学中也需要变量服从正态分布。即使没有严格的要求，但是太高的偏度、峰度及异常值将影响变异函数的结构和Kriging插值的结果。数据集中的异常值往往能使变异函数展示奇怪的现象，然而数据转化能消弱极值间的差异[11]。对数转化被广泛应用于呈右偏态的数据集。然而，在环境科学领域中并不一直遵循对数正态分布，Box-Cox转化则经常用于这一领域。本文中全磷的数据经过Box-Cox转化后很好的符合了正态分布，为后
面获得较好的变异函数结构和准确的插值图打下了基础。目前，研究出了不考虑数据分布类型新的处理数据的方法，如非参数统计方法或称无分布方法。近年来提出的非参数地质统计学理论及其主要方法——指示Kriging方法是一种更为精确和有效的估计方法。这需要另文研究，以比较它们之间对插值结果精度的影响。
笔者应用地统计学和GIS技术相结合的方法，研究了全磷的空间变异特性。由于受成土因素的影响，区域土壤性质的空间分布常呈现明显的趋势特征和异向性分布[12]。但由于过去常用的地统计学软件都无趋
6期 路 鹏等：洞庭湖平原区土壤全磷含量地统计学和GIS分析 1211
势分析和异向性轴轴向自动搜寻功能，土壤性质数据在变异函数和Kriging插值建模的参数分析上均存在困难，其空间分布规律不能很好地得到体现。运用ArcGIS8.3软件的地统计学模块的趋势分析、半变异协方差函数云工具和异向性轴轴向自动搜寻功能，可以方便地获得全磷异向性分布特征参数，和趋势效应特征，并确定了全磷的空间方向效应，可以准确地选定各向异性的方向，进而揭示全磷的空间异质性。从而克服了以往的确定各向异性的主观性和片面性，以及插值结果异向性显示不明显的缺陷。
在进行土壤养分的风险性评价时，分布概率图中临界值的设定是一个值得商榷的问题。本研究是借鉴土壤第二次普查的养分分级标准，结合研究区内养分含量的分布设定的。但是，临界值到底是多大的时候，才有可能造成农业的面源污染，才被认为是“危险的”？这与作物的生产力、土壤的性质和施肥的种类和方式有关[13]，准确的确定其值的大小，还需要进一步研究。
假设Kriging估值遵循正态分布，则超过某一临界值的概率能被计算出来。一些研究者[14]争论由于普通克立格方差不依赖于数据值，并且不能用来作为局部估值精度的衡量标准。然而，有研究者[9]持不同的观点，认为这个方法应用于这个研究领域提供了一个合理的、公正的超过某一临界值的指示概率图，而不是用来计算强烈的非线性的地统计学和模拟两者之一。
本文的研究区域约为3.7 km2，采样密度较大（每公顷约2个），这将得到较准确的插值图。本文的分析方法可以作为确定具体实施治理时田块水平的辅助参考，而对于分析评价农业非点源污染问题的研究尺度还不够大，但是为非点源污染的重点区和控制区的有效识别以及进行相关的评价提供了一种思路和方法。然而，其它平原农区上的研究[15]，采样密度较小，主要集中在较大尺度上研究土壤速效磷、钾养分的时空分布上。地统计学方法本身就是一种先进的空间分析方法，通过变异函数可以确定变量因子的空间变异程度及空间变异尺度。因为空间异质性与观测的尺度密切相关，它是尺度的函数。所以，将来在亚热带地区研究地块（field）、流域（watershed）和区域（region）不同尺度上土壤养分的空间变异特征，分析不同尺度上空间信息的相似性和相异性，这些将为空间信息的尺度外推提供有益的启示。
4 结论
本文应用地统计学和GIS相结合的方法，以洞庭湖平原区为例，研究了其典型景观单元全磷含量的空间变异特征和风险性评价。主要结论如下：
全磷的原始数据存在较强的偏斜效应，并可能存在异常值。经过Box-Cox转化成功地使数据集服从正态分布并消弱了异常值的负面影响。全磷的好的变异函数结构揭示了好的空间格局将展示在其分布图中，也说明了当前的采样密度和设计方案能恰当的揭示这样的空间格局。
在步长间隔60 m下，土壤全磷具有中等程度的空间相关性，其相关距离为338 m。全磷在NE70˚、NE160˚两个方向上具有几何异向性结构特点。在NE25˚、NE115˚两个方向上具有带状异向性结构特点。但其方向性效应并不太显著，其空间分布总的趋势是高值斑块区呈现零星分布的特点，并没有表现出条带状的分布格局。这从另一个角度确认了其变异函数所揭示的方向不太明显的特征。在经过不同趋势阶数土壤全磷Kriging插值误差的综合比较的基础上，表明趋势效应参数宜选取二阶。基于Kriging插值的概率分布图和克立格标准差图为土壤养分的风险性评价和制定决定性的管理策略提供了十分有用的信息。
References
[1]
李元沅, 李迪华, 刘劲凡, 周伟民. 洞庭湖区稻麻土壤生态系统养分平衡状况的研究. 湖南农学院学报, 1991, 17(增刊): 209-217. Li Y Y, Li D H, Liu J F, Zhou W M. Nutrient balance of soil ecological system in rice and ramie fields in Dongting Lake region. Journal of Hunan Agricultural College, 1991, 17(suppl.): 209-217. (in Chinese) [2]
刘付程, 史学正, 于东升, 潘贤章. 太湖流域典型地区土壤磷素含量的空间变异特征. 地理科学, 2003, 23(1): 77-81.
Liu F C, Shi X Z, Yu D S, Pan X Z. Characteristics of spatial variability of total phosphorus in soil of the typical area of Taihu Lake watershed. Scientia Geographica Sinica, 2003, 23(1): 77-81. (in Chinese) [3]
朱益玲, 刘洪斌, 江希流. 江津市紫色土中N、P养分元素区域空间变异性研究. 环境科学, 2004, 25(1): 138-143.
Zhu Y L, Liu H B, Jiang X L.Investigation of the spatial variability of nitrogen and phosphorus in purple soils in Jiangjin city, Sichuan, China. Environmental Science, 2004, 25(1): 138-143. (in Chinese) [4]
Goovaerts P. Geostatistical modeling of uncertainty in soil science.
1212
中 国 农 业 科 学 38卷
Geoderma, 2001, 103: 3-26. [5]
白由路, 金继运, 杨俐苹, 梁呜早. 基于GIS的土壤养分分区管理模型研究. 中国农业科学, 200l, 34(1): 46-50.
Bai Y L, Jin J Y, Yang L P, Liang M Z. Research on the subarea management model of soil nutrients by GIS. Scientia Agricultura Sinica, 200l, 34(1): 46-50. (in Chinese)
[6] 高祥照, 胡克林, 郭 焱, 李保国, 马韫韬, 杜 森, 王运华. 土壤养分与作物产量的空间变异特征与精确施肥. 中国农业科学, 2003, 35(6): 660-666.
Gao X Z, Hu K L, Guo Y, Li B G, Ma Y T, Du S, Wang Y H. Spatial variability of soil nutrients and crop yield and site-specific fertilizer management. Scientia Agricultura Sinica, 2003, 35(6): 660-666. (in Chinese)
[7] 王政权. 地统计学及在生态学中的应用. 北京: 科学出版社, 1999.
Wang Z Q. Geo-statistics and Its Application in Ecology. Beijing: Science Press, 1999. (in Chinese)
[8] McGrath D, Zhang C S, Carton O T. Geostatistical analyses and hazard assessment on soil lead in Silvermines area, Ireland. Environmental Pollution, 2004, 127: 239-248.
[9] 胡克林, 张凤荣, 吕贻忠, 王 茹, 徐 艳. 北京市大兴区土壤重金属含量的空间分布特征. 环境科学学报, 2004, 24(3): 463-468. Hu K L, Zhang F R, LÜ Y Z, Wang R, Xu Y. Spatial distribution of concentrations of soil heavy metals in Daxing county, Beijing. Acta Scientiae Circumstantiae, 2004, 24(3): 463-468. (in Chinese)
[10] Goovaerts P. Geostatistics in soil science: state-of-the-art and
perspectives. Geoderm, 1999, 89: 1-45.
[11] Gringarten E, Deutsch C V. Teacher’s aide: variogram interpretation
and modeling. Mathematical Geology, 2001, 33(4): 507-534. [12] 黄元仿, 周志宇, 苑小勇, 张红艳. 干旱荒漠区土壤有机质空间
变异特征. 生态学报, 2004, 24(12): 2 776-2 781.
Huang Y F, Zhou Z Y, Yuan X Y, Zhang H Y. Spatial variability of soil organic matter content in an arid desert area. Acta Ecologica Sinica, 2004, 24(12): 2 776-2 781. (in Chinese)
[13] 郭胜利, 党廷辉, 郝明德. 施肥对半干旱地区小麦产量、NO-
3-N
累积和水分平衡的影响. 中国农业科学, 2005, 38(4): 754-760. Guo S L, Dang T H, Hao M D. Effects of fertilization on wheat yield, NO-3-N accumulation and soil water content in semi-arid area of China. Scientia Agricultura Sinica, 2005, 38(4): 754-760. (in Chinese)
[14] Deutsch C V, Journel A G. GSLIB: Geostatistical Software Library
and Users Guide( second ed). New York: Oxford University Press, 1997.
[15] 齐 伟, 徐 艳, 张凤荣. 黄淮海平原农区县城土壤养分平衡评价
方法及其应用. 中国农业科学, 2004, 37(2): 238-243.
Qi W, Xu Y, Zhang F R. Study on evaluation methods of soil nutrient balance and application at county level in Huang-huai-hai plain. Scientia Agricultura Sinica, 2004, 37(2): 238-243. (in Chinese)
（责任编辑 李云霞）
中国农业科学 2005,38(6):1204-1212 Scientia Agricultura Sinica
洞庭湖平原区土壤全磷含量地统计学和GIS分析
路 鹏
1,2
, 彭佩钦, 宋变兰, 唐国勇, 邹 焱, 黄道友, 肖和艾, 吴金水
1111111,2
, 苏以荣
1
2
(1中国科学院亚热带农业生态研究所，长沙 410125；中国科学院水利部水土保持研究所黄土高原土壤侵蚀与旱地农业国家重点实验室，杨凌 712100)
摘要：以洞庭湖平原区典型景观单元为试点，利用GPS定位共取得651个耕层（0～20 cm）土壤样品。通过对数转化、稳健统计、域法处理和Box-Cox转化4种数据处理方法对土壤全磷进行了正态分布性检验。结果表明，Box-Cox转化成功地使数据集服从正态分布并消弱了异常值的影响，应用地统计学分析方法进行了实验变异函数的计算和最适合模型的拟合，得出土壤全磷最好的理论模型为球状模型，随后用普通克立格估值方法绘制了土壤全磷的空间分布图。在经过不同趋势阶数土壤全磷克立格（Kriging）插值误差的综合比较的基础上，结果表明趋势效应参数宜选取二阶。Kriging 估值标准差（KSD）被认为是内插象素值的标准差，并为提高全磷的制图精度提供了有用的信息。应用地统计学和地理信息系统（GIS）的概率克立格法研究了洞庭湖平原典型区的土壤全磷的空间分布并进行了风险性评价。结果表明，土壤全磷不同含量水平下的概率分布图对风险性评价、合理施肥和控制磷素的面源污染的管理实践将是十分有益的。
关键词：土壤全磷；地统计学；空间变异性；概率克立格；洞庭湖平原区
Geostatistical and GIS Analyses on Soil Total P in the Typical Area
of Dongting Lake Plain
LU Peng1,2, PENG Pei-qin1, SONG Bian-lan1, TANG Guo-yong1, ZOU Yan1, HUANG Dao-you1,
XIAO He-ai1, WU Jin-shui1,2, SU Yi-rong1
(1The Institute of Subtropical Agriculture, Chinese Academy of Sciences, Changsha 410125; 2 State Key Laboratory of
Soil Erosion and Dry Land Farming on the Loess Plateau, Institute of Soil and Water Conservation,
Chinese Academy of Sciences and Ministry of Water Resources, Yangling 712100)
Abstract: A typical landscape unit of Dongting Lake plain was selected as an experimental site. Approximate grid approaches were employed for the sampling scenario in 2004 with 651 Global Position System (GPS) established spots sampled in topsoil (0-20 cm). The purpose of this study is to evaluate various data processing methods, including logarithmic transformation, robust statistics, excluding outliers and Box-Cox transformation for evaluating soil total P content with normal distribution. The result showed that Box-Cox transformation was applied in order to achieve normality in the data set and to dampen the effect of outliers. Geostatistical analyses were carried out, including calculation of experimental variograms and model fitting. The best theoretical model for semivariogram of soil total P were spherical model. The ordinary kriging estimates of soil total P concentration were mapped. The integrative comparisons of semivariogram parameters with different trends of the kriging prediction errors of soil total P indicated that the 2-order trend effect was preferable. Kriging standard deviations (KSD) were regarded as the standard deviations of the interpolated pixel values and provided valuable information for which will increase the accuracy of the total P mapping. Spatial distribution and hazard assesment of soil total P in the typical area of Dongting Lake plain were investigated using geostatistics and geographic information system (GIS) techniques such as probability kriging on the basis of the software ArcGIS Desktop. Probability distribution of soil total P at different levels will be helpful to conduct hazard assessment, optimal fertilization and develop manage
收稿日期：2005-04-04
基金项目: 中国科学院知识创新工程重要方向项目(KZCX3-SW-426)、国家自然科学基金重点项目(40235057)、国家重点基础研究发展规划项目(2002
CB412503)
作者简介: 路 鹏 (1976-)，男，博士研究生，主要从事土壤资源信息管理与环境方面的研究。E-mail: [email protected]。苏以荣为通讯作者，Tel:
0731-4615222; E-mail: [email protected]
6期 路 鹏等：洞庭湖平原区土壤全磷含量地统计学和GIS分析 1205
ment practices to control the non-point sources of P pollution.
Key words: Soil total P; Geostatistics; Spatial variability; Probability kriging ; Dongting Lake plain area
洞庭湖区位于湖南省的北部，其地处东经111˚40′～113˚10′，北纬27˚55′～30˚22′。尤其洞庭湖平原区自然条件优越，土地肥沃，是中国南方有名的“鱼米之乡”。但是，自从20世纪80年代以来，忽视有机肥的施用，盲目大量滥施化肥，已使整个地区的耕地质量明显下降，对农业生产带来不利影响[1]。由于磷肥的施用量逐年增加，土壤全磷含量总体水平也在不断上升，超过了作物生长的需求并且普遍出现盈余现象。这一方面会导致施肥经济效益的显著下降，另一方面也会加大土壤磷素向水体流失的风险，并对流域水环境质量构成严重威胁。因此，磷肥的合理施用及其生态环境效应已受到人们的广泛关注[2]。
近年来，农业面源污染问题受到政府和研究者的不断重视[3]。磷素是农业面源污染和引发湖泊发生富营养化的主要因子之一。因此，致力于研究区域尺度，特别是景观尺度土壤磷素的空间分布特征就成了当前环境与农业领域研究热点之一。了解景观尺度土壤全磷的空间变异性，无论是对于提高作物产量，还是对于面源污染的发生和发展以及减少环境的潜在危害都具有十分重要的现实意义。
地统计学方法是一种最优的空间插值方法，它和GIS已经成为空间不确定性和风险性评价非常有用的研究工具[4]。自从20世纪70年代以来，成功地应用地统计学方法解决空间变异问题在土壤学和环境科学领域中已有较多的报道[2,3,5,6]。
目前人均耕地日趋减少，施肥模式和种植结构不太合理，土壤营养元素平衡失调，已使农业生产处于徘徊状态。因此，笔者在实地调查和采样分析的基础上，运用地统计学与GIS相结合的方法，对洞庭湖平原典型区土壤全磷的空间变异特征进行了分析，并用克立格（Kriging）插值方法绘制了土壤全磷含量分布图和风险评价图，从而指导农业生产、合理施肥，以期减少地表径流和渗漏对地下水、湖泊的污染所造成不利的影响，并为整个洞庭湖平原区农业生态环境管理提供理论基础。
1 研究区域概况与研究方法
1.1 研究区域概况
研究区域为洞庭湖平原区典型景观单元，其位于湖南省北部的沅江市，洞庭湖垸的大同乡人新与人和
村，面积约为3.7 km2。沅江市（东经112˚16′～112˚56′，北纬28˚42′～29˚11′），属中亚热带向北亚热带过渡的湿润季风湿润气候，年均气温16.5℃，年均降水量 1 313 mm，年平均相对湿度81%，降水集中且多为暴雨。地势平坦，主要地貌类型为冲积平原和湖区平原，土壤以湖积物和河流冲积物发育而成的水稻土和潮土为主，土层深厚。
1.2 土壤样品采集与化学分析
2004年3月，根据不同的土地利用方式（包括水田、旱地、新开垦的旱地），在研究区域内共采取表层（0～20 cm）土壤样品651个。用美国Garmin公司推出的eTrex Vista手持式GPS导航和定位，记录采样点的经纬度和高程，并详细记录采样点周围的景观信息。定位的样点经格式转换成ArcView能够接受的shp格式，再进行ArcInfo投影转换，生成以米（m）为单位的平面坐标（投影类型为高斯克吕格投影，中央经线为112˚00′00″E），最后生成用于地统计学分析的样点分布图（图1）。土壤样品是用网格内随机取样的方法采集，在取样点所在田块内随机采集15个点作为一个混合样，然后风干研磨，过0.25 mm筛用于样品的测定。土壤全磷采用HClO4-H2SO4 消煮，钼锑抗比色法测定。
图1 土壤样品空间分布图
Fig.1 The spatial distribution of soil samples
1.3 地统计学基本理论
关于地统计学的基本原理和方法，文献中已经有很好的解析和综述[7]，本文只作简要的说明。半变异函数是区域化变量空间自相关性的表征，它依赖于分隔两点的矢量h。假设区域化变量满足二阶平稳假设
1206 中 国 农 业 科 学 38卷
和本征假设，其半方差函数可用下式表示：
2
N(h)
γ(h)=∑Z(xi)−Z(xi+h)
i2N(h)=1
（K-S检验）是通过SPSS11.0完成的。实验半方差函数的计算、理论模型的拟合及Kriging插值和图形绘制是由软件ILWIS3.1，ArcGIS8.3共同完成的。将采样点的定位数据导入GIS软件中，获得样点分布的空间数据库，通过关键字段实现空间数据与土壤属性数据之间的联接，作为实验半方差函数的计算和Kriging插值法制图的数据源。
本文采用对数转化、稳健处理、域法处理（样本平均值X±3s）和Box-Cox转化[8] 4种数据处理方法和柯尔莫哥洛夫-斯蜜诺夫检验（K-S检验），计算了偏度、峰度和K-S P等参数，并对其是否服从正态分布进行了检验。用普通Kriging插值和概率Kriging插值方法，估值样点数最大为20，产生以5 m×5 m分辨率的栅格图形。插值的结果以ASCII数据的格式输出存储，然后生成GRID图形，再将图形用研究区域边界图形切割，最后生成研究区域相应的土壤特性含量分布图。
1
[]
式中，γ(h)为半方差函数；h为样点的空间间隔距离，称为步长；N(h)为距离等于h的点对数；Z(xi)为处于点xi处变量的实测值；Z(xi+h)为与点xi偏离h处变量的实测值。
不同形式的Kriging内插法是目前地统计学应用最广泛的最优内插法，它是根据半方差分析所提供的空间自相关程度的信息来进行插值，因此可以对未采样点给出最优无偏估计，而且能同时提供估计值的误差和精确度。也即它是利用已知点的数据去估计未知
点X0的数值，其实质是一个实行局部估计的加权平均值：
Z(Xo)=∑λiZ(Xi)
i=1
n
式中，Z(X0)是在未经观测点X0上的内插估计值，Z(Xi)是在点X0附近的若干点上获得的实测值。λi是考虑了半方差图中表示空间的权重，所以，Z(X0)值的估计应该是无偏的，因为：∑λi=1，估计偏差是最
i=1n
2 结果与分析
2.1 土壤全磷的描述性统计和正态分布性检验
全磷直方图展示了强的负偏斜效应（图2），进一步来说，这种分布和可能存在的异常值能导致在地统计分析和Kriging插值中产生计算结果的偏差。因此，检验数据的正态分布性是使用空间统计学克立格方法进行土壤特性空间分析的前提，只有当数据服从正态分布时，克立格方法才有效。
由于原始数据集的概率分布表现出强的偏斜效应，因此，在地统计和插值分析之前必须进行数据的转化。一些参量如偏度、峰度和正态性的柯尔莫哥洛夫-斯蜜诺夫检验的显著水平（K-S P）列于表1。由表1
可以得出，经过对数变化、域法处理都没有使全磷
小的，并可由下列方程求出：
σ
2
d
=b
T
λμ
式中，b是被估计点与其它点之间的半方差矩阵，μ为拉格朗日参数。 1.4 数据处理
正态分布性检验和柯尔莫哥洛夫-斯蜜诺夫检验
图2 土壤全磷含量的直方图展示了强的偏斜概率分布
Fig. 2 Histogram of total soil P concentrations (n=651) showing a heavily skewed
6期 路 鹏等：洞庭湖平原区土壤全磷含量地统计学和GIS分析 1207
表1 土壤耕层全磷的描述性统计
Table 1 Descriptive statistics of soil total P in topsoil (g·kg -1)
处理方法 Treatment method
样本数 Sample
元素 Element
平均值 Mean
标准差 S.D.
最小值 Min.
最大值 Max
偏度
峰度
Skewness Kurtosis
K-S P
原始数据Original data 651
对数转化Logarithmic transformation 651 稳健统计Robust statistics 651 域法处理Excluding outliers 651 Box-Cox 651 TP 0.99 0.13 0.48 1.34 -0.52 0.65 0.02
TP -0.01 0.06 -0.32 0.13 -1.08 2.19 0 TP 0.99 0.12 0.74 1.19 -0.38 -0.34 0.09 TP 0.99 0.13 0.60 1.34 -0.45 0.37 0.03 TP 0.01 0.13 -0.30 0.47 0.38 0.37 0.24
含量成正态分布。经Box-Cox变换后全磷通过了柯尔莫哥洛夫-斯蜜诺夫检验（K-S检验），这表明已经服从正态分布且满足地统计学分析要求。尽管稳健处理使得土壤全磷成正态分布，但效果远不如Box-Cox变化。因此，Box-Cox转化后的数据被用来随后的Kriging插值分析。
然而，表1的描述性统计分析只能说明土壤全磷样点的含量特征，但是不能完全反映整个景观单元的区域特征，也就是说不能反映土壤全磷含量的空间结构性和随机性。因此，需要用地统计学和GIS相结合的方法来弥补这一缺陷。
2.2 各向同性下土壤全磷含量的空间变异分析 根据不同空间位置上全磷含量的分析数据，计算
实际半方差值γ(h)并绘制变异函数曲线图，这是地统计学中空间变异分析的基本步骤，也是进行Kriging插值的前提。变异函数曲线图表示全磷含量的区域化变量在距离与方向上的所有成对点观测值之间的空间相关性。得出的半变异函数图的起伏特征、原点处性状、趋势走向、不同的方向差异性等形状特点提供了丰富的空间结构信息。
土壤全磷的各向同性下的变异函数（图3）展示了很好的空间结构，很好的符合球状模型。步长间隔为60 m时的块金方差为0.0078，块金方差很小表明采样密度能够充分揭示研究区域土壤全磷的空间结 构[8]。进一步来说，其变程为338 m，说明空间自相关尺度远比步长间隔60 m大。从另一个方面说明了当前的采样设计是恰当的，并且预示了好的空间结构将展示在空间插值图中。
表2中，全磷半方差函数模型的拟合度均在0.90以上，说明半方差函数模型的选取能很好地体现其空间变异性。块金值与基台值之比是反映区域化变量空间异质性程度的重要指标，也反映了空间变异的成分中区域因素（自然因素）和非区域因素（人为因素）谁占主导作用。从C0/C0+C的大小可以看出，影响全磷含量空间分布的总变异中，由非区域因素引起的变
图3 土壤全磷各向同性下变异函数图 Fig.3 Isotropic semivariogram of soil total P
异达到了45.35%，这说明非区域因素起的作用表现的比较明显，也与近20年来人们普遍重视磷肥的施用有关。因此，施肥等人为因素成为影响全磷含量空间分布的主要因子。全磷的C0/C0+C小于0.5，反映了洞庭湖平原区其在所研究的尺度上具有中等的空间自相关格局，目前虽然受到植物吸收、施肥等小尺度因素的影响，但是还没有达到破坏其原有空间格局的程度[9]。 2.3 各向异性下土壤全磷含量的空间变异分析
为了解自然过程对全磷空间变异性的影响，笔者通过半变异/协方差函数云（semivariogram/ covariance cloud）工具和异向性轴轴向自动搜寻分析
了其空间方向效应，并结合地块具有明显的条带状和对称性的特点，拟以NE25˚、NE115˚、NE70˚、NE160˚4个方向计算其空间变异函数参数值（表3）和它们的
各向异性比K(h)（图4）。
由表3可知，全磷在NE70˚、NE160˚两个方向上有相同的基台值，变程不同，具有几何异向性结构特点。在NE25˚、NE115˚两个方向上，有不同的基台值，变程也不同，具有带状异向性结构特点。从图4可知，虽然全磷的各向异性比（K(h)）在大部分尺度上均位于各向同性线之上，但是起伏不大，说明全磷有向各向同性发展的趋势。这也从另一个角度说明了人为因
1208 中 国 农 业 科 学 38卷
表2 土壤全磷的半方差模型拟合参数
Table 2 Parameters fitted by semivariogram model for soil total P
土壤属性 Soil properties 土壤全磷 soil total P
模型类型 Model S
块金值 (C0) Nugget 0.0078
基台值 (C0+C)
Sill 0.0172
块金值/基台值 C0/C0+C 0.4535
变程 (m) Range 338
R2 RSS 0.9000
0.000008
S. 球状模型 Spherical model 下同 The same as below
表3 4个方向上全磷的理论模型及相应的参数
Table 3 Theory models and corresponding parameters of soil total P in four directions
土壤属性 Soil nutrient 土壤全磷 Soil total P
方向 Direction
模型 Model
块金值 (C0)
基台值 (C0+C)
块金值/基台值 C0/C0+C
变程 Range (m)
280 320 360 260
25˚ S 0.007 0.0172 0.407 70˚ S 0.007 0.0175 0.400 115˚ S 0.009 0.0160 0.563 160˚ S 0.006 0.0175 0.343
二阶（二阶多项式变化）。图5是全磷含量的全局趋势效应分析示意图。图中，X轴表示正东方向，Y轴表示正北方向，Z轴表示各样点测定值的大小；左后投影面上曲线表示东－西向全局性的趋势效应变化情况，右后投影面上曲线表示南－北向全局性的趋势效应变化情况。
在均考虑各向异性的情况下，对分别选择无趋势、一阶趋势和二阶趋势效应参数，结合普通Kriging插值方法所造成的插值误差进行了比较（表4）。判断
图4 不同方向上土壤全磷的各向异性比
Fig.4 Anisotropic ratio of soil total P at different directions
半方差函数模型及其参数是否合适可按以下标准综合进行：平均误差（ME）的绝对值最接近于0；标准化平均误差（MSE）最接近于0；均方根误差（RMSE）越小越好；平均标准误差（ASE）与均方根误差（RMSE）最接近，如果ASE>RMSE则高估了预测值，反之如果ASE1，则低估了预测值。从表4、图5可以得出，洞庭湖平原典型区土壤全磷的趋势效应不太明显，但在两个方向均成二阶多项式变化。综合考虑分析，在以后进行普通Kriging插值时，趋势参数宜选取二阶。
素如施肥、耕作管理等措施消弱了其空间异质性程度。
因此，随后的Kriging插值分析是根据各向异性结构的变异函数理论模型，并采用块段空间局部估计的方法绘制的空间分布格局。 2.4 土壤全磷含量的趋势分析
运用有趋势分析和异向性轴轴向自动搜寻功能的ArcGIS8.3软件的地统计学模块，可以方便地获得全磷异向性分布特征参数，以及趋势效应特征。趋势效应一般分为0（没有趋势效应）、一阶（线性变化）、
表4 不同趋势阶数土壤全磷Kriging插值误差的比较
Table 4 The comparisons of semivariogram parameters with different trend of soil total P
土壤属性 Soil properties 土壤全磷 Soil total P
趋势效应 Trend effect 无 Non 一阶 1-order
二阶 2-orfer
预测误差 Prediction error
ME RMSE ASE MSE RMSSE -0.000014 -0.000016
0.1166 0.1166
0.1111 0.1110
-0.000154 -0.000178
1.046 1.046
0.000027 0.1166 0.1110 0.000145 1.046
6期 路 鹏等：洞庭湖平原区土壤全磷含量地统计学和GIS分析
1209
图5 土壤全磷的趋势效应 Fig.5 Surface trend of soil total P
2.5 土壤全磷含量的空间分布格局分析
空间异质性是空间插值研究的隐含前提，要素的非均匀空间分布才需要空间插值；空间相关性则是空间插值研究的基础，缺乏这种相关性，空间插值就成为了一种数学游戏。因此，基于经过Box-Cox转化的采样点的数据，结合普通Kriging插值分析方法，选取二阶趋势参数并考虑各向异性，随后用Box-Cox转化的逆变换过程获得了洞庭湖平原典型区的土壤全磷含量的空间分布格局图（图6）。用一个5 m×5 m分辨率的栅格，将研究区的面积分成了440列（NE115˚方向，2 200 m）和560行（NE25˚方向，2 800 m）。当然，Kriging插值的结果也受到变异函数模拟精度、样点的分布、邻近样点的选取数等的影响[10]。
图6 土壤全磷的空间分布格局
Fig.6 Spatial distribution pattern of soil total P
全磷含量的空间分布格局进一步表明土壤表层（0～20 cm）养分含量分布具有高度的空间异质性，
并决定了空间格局的存在。从插值图中可方便地了解到全磷含量斑块大小、形状及空间分布等具有显著的差异，但不论在什么方向上，养分含量由高到低的分布梯度的规律总是存在的。显然，在此基础上借鉴土壤第二次普查的养分分级标准，可以判断一定区域内养分流动与非点源污染发展趋势及潜力。因此，这种养分的空间分布特征与研究区土壤在不同空间位置上的各种物理、化学和生物过程有着重要的联系，更重要的是土壤与作物长期相互作用，如稻田生态系统养
分生物地球化学循环的结果。
从图6中可以看出，全磷的方向性效应并不太明显，空间分布总的趋势是高值斑块区呈现零星分布的特点，并没有表现出条带状的分布格局。这从另一个角度确认了其变异函数所揭示的方向不太明显的特征。全磷的高值斑块区是水田分布的区域；低值斑块区恰恰是旱地、新开垦旱地和一些鱼塘的点缀分布区。最上部的低值区为一条河堤被新开垦成的旱地，这是其表现相对低的主要原因。这也体现出了土地利用方式的不同造成了全磷含量空间分布上的差异。 2.6 土壤全磷的标准差分布图
Kriging方差是与变异函数的结构和样点的分布有关系的，并且能提供Kriging插值结果的置信水平。Kriging的标准差（KSD）是Kriging方差的平方根，被认为是内插象素值的标准差，也是样本数据的标准差和大范围取样插值结果内在不确定性的综合反映。因此，其值大小可用于评价Kriging插值精度，值越小表明Kriging插值结果越可靠。
图7是土壤全磷Kriging插值的标准差分布图，其体现了模型误差的水平，其估值标准差变化于0.1052～0.1415 g·kg-1，平均值为0.1119 g·kg-1。在采样点的附近，KSD值较小；远离采样点的区域，KSD值一般较大。围绕采样点周围的近似椭圆形的KSD值是由变异函数的椭圆形结构造成的。KSD值的大小与采样点密度存在着明显的对应关系，这是由于KSD值仅与采样点密度、采样点的布设结构及实验函数模型有关，而与采样点实测结果的大小没有直接关系[8]。
2.7 土壤全磷含量的风险性评价
为了进一步了解洞庭湖垸典型区全磷含量的空间分布特点，本文基于ArcGIS Desktop 8.3软件，采用概率克立格方法并结合Box-Cox转化的数据，选取二阶趋势参数并考虑各向异性，对不同含量水平下全磷的概率分布做了风险性评价。值得注意的是，Box-Cox转化也应当应用于临界值的设定中。
1210 中 国 农 业 科 学 38卷
1.08 km2，占整个研究区面积的28.9%。其次是0.20～0.32，面积为0.74 km2，占整个研究区面积的19.8%。低的全磷含量水平主要以小概率区间分布，这说明研究区范围内全磷的整体水平是较高的。全磷含量≤1.0 g·kg-1的平均概率为0.5109，其空间分布以0.2～0.4、0.5～0.65和0.65～0.8的概率为主（图8-b），面积分别为0.85、0.82和0.81 km2，占整个研究区面积的22.8%、21.9%和21.8%。从图8-b也可以看出，0.0～0.2和0.8～1.0这两个概率区间所占的面积最小，其余
图7 土壤全磷的Kriging标准差分布图
Fig.7 Kriging standard deviation (KSD) of soil total P
4个概率区间占的面积相当。全磷含量>1.05 g·kg-1的平均概率为0.3263(图8-c)，其空间分布以0.25～0.40较低的概率区间为主，面积为1.05 km2，占整个研究区面积的28.3%。这部分田块区域成为农业面源污染的严重区。因此，开展全磷的风险性评价研究，能从区域尺度上发现农田土壤养分管理过程中可能存在的问题，这对整个洞庭湖平原区农业生态环境管理提供理论基础，将有助于农业新一轮的大发展。
土壤全磷含量的风险性评价如图8所示，从整体上看，全磷含量≤0.9 g·kg-1的平均概率为0.2480。其概率区间以0.0～0.1的分布最广（图8-a），面积为
图8 土壤全磷不同含量水平下的概率分布
Fig.8 Probability distribution of soil total P at different levels
3 讨论
正如传统的统计学研究一样，在地统计学中也需要变量服从正态分布。即使没有严格的要求，但是太高的偏度、峰度及异常值将影响变异函数的结构和Kriging插值的结果。数据集中的异常值往往能使变异函数展示奇怪的现象，然而数据转化能消弱极值间的差异[11]。对数转化被广泛应用于呈右偏态的数据集。然而，在环境科学领域中并不一直遵循对数正态分布，Box-Cox转化则经常用于这一领域。本文中全磷的数据经过Box-Cox转化后很好的符合了正态分布，为后
面获得较好的变异函数结构和准确的插值图打下了基础。目前，研究出了不考虑数据分布类型新的处理数据的方法，如非参数统计方法或称无分布方法。近年来提出的非参数地质统计学理论及其主要方法——指示Kriging方法是一种更为精确和有效的估计方法。这需要另文研究，以比较它们之间对插值结果精度的影响。
笔者应用地统计学和GIS技术相结合的方法，研究了全磷的空间变异特性。由于受成土因素的影响，区域土壤性质的空间分布常呈现明显的趋势特征和异向性分布[12]。但由于过去常用的地统计学软件都无趋
6期 路 鹏等：洞庭湖平原区土壤全磷含量地统计学和GIS分析 1211
势分析和异向性轴轴向自动搜寻功能，土壤性质数据在变异函数和Kriging插值建模的参数分析上均存在困难，其空间分布规律不能很好地得到体现。运用ArcGIS8.3软件的地统计学模块的趋势分析、半变异协方差函数云工具和异向性轴轴向自动搜寻功能，可以方便地获得全磷异向性分布特征参数，和趋势效应特征，并确定了全磷的空间方向效应，可以准确地选定各向异性的方向，进而揭示全磷的空间异质性。从而克服了以往的确定各向异性的主观性和片面性，以及插值结果异向性显示不明显的缺陷。
在进行土壤养分的风险性评价时，分布概率图中临界值的设定是一个值得商榷的问题。本研究是借鉴土壤第二次普查的养分分级标准，结合研究区内养分含量的分布设定的。但是，临界值到底是多大的时候，才有可能造成农业的面源污染，才被认为是“危险的”？这与作物的生产力、土壤的性质和施肥的种类和方式有关[13]，准确的确定其值的大小，还需要进一步研究。
假设Kriging估值遵循正态分布，则超过某一临界值的概率能被计算出来。一些研究者[14]争论由于普通克立格方差不依赖于数据值，并且不能用来作为局部估值精度的衡量标准。然而，有研究者[9]持不同的观点，认为这个方法应用于这个研究领域提供了一个合理的、公正的超过某一临界值的指示概率图，而不是用来计算强烈的非线性的地统计学和模拟两者之一。
本文的研究区域约为3.7 km2，采样密度较大（每公顷约2个），这将得到较准确的插值图。本文的分析方法可以作为确定具体实施治理时田块水平的辅助参考，而对于分析评价农业非点源污染问题的研究尺度还不够大，但是为非点源污染的重点区和控制区的有效识别以及进行相关的评价提供了一种思路和方法。然而，其它平原农区上的研究[15]，采样密度较小，主要集中在较大尺度上研究土壤速效磷、钾养分的时空分布上。地统计学方法本身就是一种先进的空间分析方法，通过变异函数可以确定变量因子的空间变异程度及空间变异尺度。因为空间异质性与观测的尺度密切相关，它是尺度的函数。所以，将来在亚热带地区研究地块（field）、流域（watershed）和区域（region）不同尺度上土壤养分的空间变异特征，分析不同尺度上空间信息的相似性和相异性，这些将为空间信息的尺度外推提供有益的启示。
4 结论
本文应用地统计学和GIS相结合的方法，以洞庭湖平原区为例，研究了其典型景观单元全磷含量的空间变异特征和风险性评价。主要结论如下：
全磷的原始数据存在较强的偏斜效应，并可能存在异常值。经过Box-Cox转化成功地使数据集服从正态分布并消弱了异常值的负面影响。全磷的好的变异函数结构揭示了好的空间格局将展示在其分布图中，也说明了当前的采样密度和设计方案能恰当的揭示这样的空间格局。
在步长间隔60 m下，土壤全磷具有中等程度的空间相关性，其相关距离为338 m。全磷在NE70˚、NE160˚两个方向上具有几何异向性结构特点。在NE25˚、NE115˚两个方向上具有带状异向性结构特点。但其方向性效应并不太显著，其空间分布总的趋势是高值斑块区呈现零星分布的特点，并没有表现出条带状的分布格局。这从另一个角度确认了其变异函数所揭示的方向不太明显的特征。在经过不同趋势阶数土壤全磷Kriging插值误差的综合比较的基础上，表明趋势效应参数宜选取二阶。基于Kriging插值的概率分布图和克立格标准差图为土壤养分的风险性评价和制定决定性的管理策略提供了十分有用的信息。
References
[1]
李元沅, 李迪华, 刘劲凡, 周伟民. 洞庭湖区稻麻土壤生态系统养分平衡状况的研究. 湖南农学院学报, 1991, 17(增刊): 209-217. Li Y Y, Li D H, Liu J F, Zhou W M. Nutrient balance of soil ecological system in rice and ramie fields in Dongting Lake region. Journal of Hunan Agricultural College, 1991, 17(suppl.): 209-217. (in Chinese) [2]
刘付程, 史学正, 于东升, 潘贤章. 太湖流域典型地区土壤磷素含量的空间变异特征. 地理科学, 2003, 23(1): 77-81.
Liu F C, Shi X Z, Yu D S, Pan X Z. Characteristics of spatial variability of total phosphorus in soil of the typical area of Taihu Lake watershed. Scientia Geographica Sinica, 2003, 23(1): 77-81. (in Chinese) [3]
朱益玲, 刘洪斌, 江希流. 江津市紫色土中N、P养分元素区域空间变异性研究. 环境科学, 2004, 25(1): 138-143.
Zhu Y L, Liu H B, Jiang X L.Investigation of the spatial variability of nitrogen and phosphorus in purple soils in Jiangjin city, Sichuan, China. Environmental Science, 2004, 25(1): 138-143. (in Chinese) [4]
Goovaerts P. Geostatistical modeling of uncertainty in soil science.
1212
中 国 农 业 科 学 38卷
Geoderma, 2001, 103: 3-26. [5]
白由路, 金继运, 杨俐苹, 梁呜早. 基于GIS的土壤养分分区管理模型研究. 中国农业科学, 200l, 34(1): 46-50.
Bai Y L, Jin J Y, Yang L P, Liang M Z. Research on the subarea management model of soil nutrients by GIS. Scientia Agricultura Sinica, 200l, 34(1): 46-50. (in Chinese)
[6] 高祥照, 胡克林, 郭 焱, 李保国, 马韫韬, 杜 森, 王运华. 土壤养分与作物产量的空间变异特征与精确施肥. 中国农业科学, 2003, 35(6): 660-666.
Gao X Z, Hu K L, Guo Y, Li B G, Ma Y T, Du S, Wang Y H. Spatial variability of soil nutrients and crop yield and site-specific fertilizer management. Scientia Agricultura Sinica, 2003, 35(6): 660-666. (in Chinese)
[7] 王政权. 地统计学及在生态学中的应用. 北京: 科学出版社, 1999.
Wang Z Q. Geo-statistics and Its Application in Ecology. Beijing: Science Press, 1999. (in Chinese)
[8] McGrath D, Zhang C S, Carton O T. Geostatistical analyses and hazard assessment on soil lead in Silvermines area, Ireland. Environmental Pollution, 2004, 127: 239-248.
[9] 胡克林, 张凤荣, 吕贻忠, 王 茹, 徐 艳. 北京市大兴区土壤重金属含量的空间分布特征. 环境科学学报, 2004, 24(3): 463-468. Hu K L, Zhang F R, LÜ Y Z, Wang R, Xu Y. Spatial distribution of concentrations of soil heavy metals in Daxing county, Beijing. Acta Scientiae Circumstantiae, 2004, 24(3): 463-468. (in Chinese)
[10] Goovaerts P. Geostatistics in soil science: state-of-the-art and
perspectives. Geoderm, 1999, 89: 1-45.
[11] Gringarten E, Deutsch C V. Teacher’s aide: variogram interpretation
and modeling. Mathematical Geology, 2001, 33(4): 507-534. [12] 黄元仿, 周志宇, 苑小勇, 张红艳. 干旱荒漠区土壤有机质空间
变异特征. 生态学报, 2004, 24(12): 2 776-2 781.
Huang Y F, Zhou Z Y, Yuan X Y, Zhang H Y. Spatial variability of soil organic matter content in an arid desert area. Acta Ecologica Sinica, 2004, 24(12): 2 776-2 781. (in Chinese)
[13] 郭胜利, 党廷辉, 郝明德. 施肥对半干旱地区小麦产量、NO-
3-N
累积和水分平衡的影响. 中国农业科学, 2005, 38(4): 754-760. Guo S L, Dang T H, Hao M D. Effects of fertilization on wheat yield, NO-3-N accumulation and soil water content in semi-arid area of China. Scientia Agricultura Sinica, 2005, 38(4): 754-760. (in Chinese)
[14] Deutsch C V, Journel A G. GSLIB: Geostatistical Software Library
and Users Guide( second ed). New York: Oxford University Press, 1997.
[15] 齐 伟, 徐 艳, 张凤荣. 黄淮海平原农区县城土壤养分平衡评价
方法及其应用. 中国农业科学, 2004, 37(2): 238-243.
Qi W, Xu Y, Zhang F R. Study on evaluation methods of soil nutrient balance and application at county level in Huang-huai-hai plain. Scientia Agricultura Sinica, 2004, 37(2): 238-243. (in Chinese)
（责任编辑 李云霞）