摘要：针对环渤海低平原土壤盐化特征，通过选取0~20 cm深度内的127个代表性土样进行土壤全盐量测定，结合dem，综合运用地质统计学方法与gis技术研究土壤盐分的空间分布状况。结果表明，土壤盐分空间变异强度属于中等，土壤盐分单变量和交互变量的空间相关程度均属于中等，空间自相关距离分别为40.60 km 和135.60 km。环渤海低平原土壤盐分含量在不同高程上存在分布差异，并随着高程的降低土壤盐分含量呈增加的趋势。自内陆平原向东部滨海平原土壤盐分含量逐渐增加，盐分含量较高的地区出现在唐山-天津-沧州-东营-滨州一线。非盐化土以内陆平原为主，分布在保定-衡水-邢台-邯郸一线。通过协同克立格插值的均方根误差与普通克立格插值相比减少了0.29%，而预测值与实测值的相关系数提高了20.58%。
关键词：环渤海低平原；土壤盐分；dem； 协同克立格
spatial variability of soil salinity and its estimation by co-kriging based on dem in the low plain around the bohai sea
zhou zai-ming, zhang guang-hui, wang jin-zhe, yan ming-jiang
( the institute of hydrogeology and environmental geology, chinese academy of geological sciences, shijiazhuang 0500061, china)
abstract: taking into account salinization existing in the low plain around the bohai sea at present, soil samples from 0 to 20 cm were collected from 127 sites in the plain, and soil salinity was analyzed. in additional, the digital elevation model (dem) of the study area was also considered. geostatistical method combined with gis was used to analyze spatial variability of soil salinity. results showed that: spatial variability belonged to moderate degree with cv=0.91. spatial correlation of single variable and cross variable belonged to moderate degree, spatial correlation distance was 40.60 km and 135.60 km respectively. soil salinity spatial distribution difference existed in different altitude as the altitude falling soil salinity rising. soil salinity content showed an increasing trend from inland plain (baoding-hengshui-xingtai-handan) to the east coastal plain (tangshan-tianjin-cangzhou-dongying-binzhou). compared with the ordinary kriging under the same sampling numbers, the root-mean-square error produced by co-kriging decrease by 0.29%, while the correlation coefficient between the predicted value and the measured value increased by 20.85%.
key words: low plain around the bohai sea; soil salinity; digital elevation model (dem); co-kriging
土壤盐分、盐渍化对农业可持续生产构成重大威胁［1］，土地盐化程度的加剧不但减少了作物产量而且限制了对作物种类的选择［2］，在干旱半干旱地下水浅埋区、滨海低平原区尤为显著［3-4］。wwW.11665.cOM对土壤盐分含量的准确估计和分布研究在区域和全球土壤环境及相关研究中是一项重要的课题［5］。对于土壤而言，其理化性质、 养分盐分变化与地形位置关系密切，地形特征对土壤形成及利用具有重要作用［6-7］。因此在土壤盐分含量的插值估算中区域地形因素的影响不容忽视。
在土壤盐分空间插值中，协同克立格借助辅助变量提高主变量的预测精度得到广泛应用［8］。eldeiry 等，综合应用地面验证和遥感数据对美国阿肯色流域土壤盐分含量的kriging和cokriging插值方法进行了比较，认为cokriging方法对玉米地的土壤盐分估计精度最高，其次是小麦和苜蓿［9］。josé 等应用cokriging 方法，对西班牙东南部塞奎拉河域土地，以土壤水电导率为辅助变量进行了50 cm深度的盐分预测研究，较好的反应了作物根系层的盐渍化状况［10］。在国内相关的研究以西北干旱内陆和黄河三角洲为主，应用表层盐分含量对底层进行插值［11-13］。
针对整个环渤海低平原土壤盐化问题，本文以研究区dem为辅助变量应用协同克立格插值方法对土壤盐分的空间分布进行研究，为提高土壤盐分含量的插值估算精度，为土地资源管理利用和农业生产的合理布局提供科学依据。
1 研究区概况
环渤海低平原位于华北东部，包括河北、山东内陆低平原和滨海低平原两大类型，含天津、沧州等80多个县市（图1）。由黄河、海河、滦河冲积而成，地势低平，大部分海拔50 m以下，滨海区10 m左右。总面积9.6×104 km2；总耕地面积400多万hm2，是我国重要的粮棉、果蔬产区。
图1 环渤海低平原行政区划与取样点分布
fig.1 administrative divisions and distribution of soil samples in the low plain around the bohai sea
2 数据与方法
2.1 数据来源
根据该区的土地利用类型和地下水状况，全区共布设127个土样点（图1），取样点间距约20 km。取样范围为（n 36°03′- n 39°35′；e 114°36′- e 119°28′），取样面积8.97×104 km2 。取样时间为2010年4月-5月，每个点采用gps定位，土层深度为0~20 cm。同时调查当地的土地耕作、灌溉制度等信息。
称取过2 mm筛的风干土试样50~100 g，按土水比1∶5配制浸出液。采用蒸干法，添加15%双氧水溶液与2%碳酸钠溶液进行了土壤易容盐总量的测定，详细操作方法参见《土壤农业化学分析方法》［14］。
2.2 应用软件
采用spss16.0软件进行土壤全盐量的统计分析，变异函数及协同变异函数模型的拟合采用地学统计软件gs+7.0，克立格与协同克立格空间插值应用arcgis9.3软件。
3 半方差函数
半方差函数的一些重要参数如块金值、基台值和变程等可以用来表示区域化变量在一定尺度上的空间变异和相关程度，它是研究土壤特性空间变异的关键，也是精确克立格插值的基础［15］。在本征平稳假设下，半方差计算公式为：
r（h）=12n（h）∑n（h）i=1［z（x.i）-z（x.i+h）］2（1）
式中：r（h）—间距h的半方差，也是以h为间距所有观测点的成对数目；n（h）—间距为h的计算对数；z（x.i）、z（x.i+h）—区域化变量z（x.i）和z（x.i+h）在空间位置x.i和x.i+h处的实测值。
两个随机变量的协同区域化可以用交互半方差函数来表示：
r.ij（h）=12n（h）∑n（h）a=1［z.i（x.a）-z.i（x.a+h）］*［z.j（x.a）-z.j（x.a+h）］（2）
式中:r.ij（h）-两个变量的交互半方差值;n（h）-具有相同间距h的变量z.i（x）和z.j（x）的离散点的数目。如果两个变量是正相关的，那么变量z.i从x.a到x.a+h的增加或减少会引起z.j的增加或减少，交互半方差就是正值。
4 研究区dem
数字高程模型dem，是以高程来描述地面特征，是进行地形分析和二维地形空间模拟的基础［16-17］。根据1∶5 万地形图（等高距20 m）生成tin模型，用实地采样高程点进行高程验证，利用均方根误差rmse描述dem精度rmse=∑ni=1（z.i-z.i）2/n。式中z.i为高程数据真值，z.i为计算值，n为误差个数。得rmse=6.52< 7，满足dem分析要求［18］。图2 环渤海低平原高程分布
fig.2 digital elevation model in the low plain around the bohai sea
5 结果分析
5.1 土壤盐分的统计特征
根据测试结果对土壤盐分含量进行统计（表1），环渤海低平原区耕层土壤盐分含量均值为1.08 g/kg属于轻度盐渍化土［19］。变异系数cv值反映样点的离散程度，cv<0.1为弱变异性；0.1≤cv≤1为中等变异性；cv>1为强变异性［20］。由表1可见，环渤海低平原区土壤盐分的空间变异强度属于中等（cv=0.91）。这主要是由于研究区范围较大且存在地势差异、土地耕种方式差异、田块灌溉制度差异等因素共同作用的结果。
项 目样点数偏度系数峰度系数均值标准差最小值最大值变异系数(cv)
土壤盐分/(g·kg-1)1273.029.481.080.980.435.900.91
根据cochran［21］的随机采样公式计算发现，本研究的取样点数满足90%置信水平和15%相对误差的精度要求，在此条件下土壤盐分的样点数为100。
由频率分布图（图3）可以看出，土壤盐分的分布呈现明显向左偏倚，大部分样点的土壤盐分在0.98~1.25 g/kg之间。图3 环渤海低平原土壤盐分含量频率分布
fig.3 frequency distribution of soil salinity in the low plain around the bohai sea
由于土壤盐分局部地区存在高值故在频率图右侧出现长尾现象。单样本k-s正态检验结果表明，不服从正态分布，也不服从对数正态分布。
5.2 土壤盐分与高程关系
根据相应取样点的高程信息进行person相关分析发现，在0.05显著性水平下通过双边检验（two-tailed）土壤盐分与高程之间存在显著的负相关，相关系数-0.253。这说明，土壤盐分含量在不同高程上存在分布差异，并随着高程的降低，土壤盐分含量呈增加的趋势。在研究区内高程与土壤盐分受区域化现象或空间过程的影响，属于协同区域化变量。
5.3 普通克立格和协同克立格插值
表2表明，土壤盐分含量半方差和交互半方差函数可分别应用球状模型和高斯模型进行拟合。由块金值与基台值的比值c.0/(c.0+c)可知，土壤盐分单变量和交互变量的空间相关程度均属于中等，其值分别为41.78%和27.95%。单变量与交互变量的空间自相关距离分别为40.60 km 和135.60 km，表明交互变量由于融合了辅助变量的空间信息，扩展了变量具有空间相关性的范围。因此在相同样点数量下，协同克立格插值比普通克立格具有更大的适用范围。
table2. models of semivariogram and cross-semivariogram for soil salinity
函数类型模型块金值c.0基台值c.0+cc.0/(c.0+c)(%)变程/kmr2
单变量球状0.4271.02241.7840.600.577
交互变量高斯-0.282-1.00927.95135.600.735环渤海低平原土壤盐分含量，自内陆平原向东部滨海平原逐渐增加，盐分含量较高的地区出现在唐山-天津-沧州-东营-滨州一线。盐分含量小于1 g/kg的非盐化土以内陆平原为主，分布在保定-衡水-邢台-邯郸一线。从dem及地形因子分析，土壤盐分含量较高的地区地势平坦，海拔较低；而土壤盐分含量低的区域则以海拔与地势起伏相对较大的内陆区为主。从成土母质上看，滨海平原的土体一部分是由河流入海冲积成的三角洲(如黄河三角洲)，一部分是海积平原。质地多为粉砂、细砂组合，粗粉砂含量多达60%以上。
图4 环渤海低平原土壤盐分空间分布
fig4. spatial distribution of soil salinity in the low plain around the baohai sea
(a. ordinary kriging; b. co-kriging)
这种砂壤土体构型使得水盐迁移量高于内陆平原的黏壤土和亚砂土的组合。从农业活动看，河北平原区的灌溉水源以抽取深层地下淡水为主，山东平原区以引黄灌溉为主，并建有完善的水利工程，有合理的耕作制度，因此农田长期处于脱盐状态，盐分累积小。滨海平原区农田主要依靠自然降雨，长期的地面蒸发使得盐化加重。
5.4 精度检验
用均方根误差减少的百分数（rrmse）表示预测精度的提高程度（式3）。用r.r表示相关系数的提高程度（式4）。式中，rmse.ok和r.ok分别是参考方法的预测均方根误差及预测值与实测值间的相关系数，rmse.ok和r.ck分别表示评价方法的预测均方根误差及预测值与实测值间的相关系数。
rrmse=（rmse.ok-rmse.ck）×100％/rmse.ok（3）
r.r=（r.ck-r.ok）/r.ok×100％（4）
由表3可见，通过协同克立格插值的均方根误差与普通克立格插值相比减少了0.29%，而预测值与实测值的相关系数提高了20.58%。表明在相同的取样条件下，协同克立格插值由于融合了更丰富的空间信息，预测的精度要高于普通克立格插值，可优化样点的 样点数
普通克立格ordinary kriging协同克立格co-krigingrmserrmser
rrmse(%)r.r(%) 1271.0340.5781.0310.6970.2920.58
6 结论
通过对环渤海低平原耕层0~20 cm土壤盐分统计分析表明，空间分布较为复杂属于偏态，大部分样点的土壤盐分在0.98~1.25 g/kg之间，空间变异强度属于中等。土壤盐分含量半方差和交互半方差函数可分别应用球状模型和高斯模型进行拟合。土壤盐分单变量和交互变量的空间相关程度均属于中等，单变量与交互变量的空间自相关距离分别为40.60 km 和135.60 km。
克立格与协同克立格插值表明，环渤海低平原土壤盐分含量，自内陆平原向东部滨海平原逐渐增加，土壤盐分含量在不同高程上存在分布差异，并随着高程的降低，土壤盐分含量呈增加的趋势。盐分含量较高的地区出现在唐山-天津-沧州-东营-滨州一线。盐分含量小于1 g/kg的非盐化土以内陆平原为主，分布在保定-衡水-邢台-邯郸一线。精度分析发现，通过协同克立格插值的均方根误差与普通克立格插值相比减少了0.29%，而预测值与实测值的相关系数提高了20.58%。
参考文献:
［1］ cetin m, kirda c. spatial and temporal changes of soil salinity in a cotton field irrigated with low-quality water［j］.journal of hydrology, 2003, 272: 238-249.
［2］ singh r b, chauhan cps, minhas p s. water production functions of wheat (triticum aestivum l.) irrigated with saline and alkali waters using double-line source sprinkler system［j］. agricultural water management, 2009, 96: 736-744.
［3］ il’ichev a t, tsypkin g g, pritchard d, et al. instability of the salinity profile during the evaporation of saline groundwater［j］. the journal of fluid mechanics, 2008, 614:87-104.
［4］ guswa a j. models of soil moisture dynamics inecohydrology: a comparative study［j］. water resources research, 2002, 38(9): 1-15.
［5］ liu x y, peterson j, zhang z y, et al. improving soil salinity prediction with high resolution dem derived from lidar data［j］. international archives of the photogrammetry, remote sensing and spatial information sciences, 2005, 36(7): 41-43.
［6］ lin y s, lin y w, wang y, et al. relationship between topography and spatial variations in groundwater and soil morphology within the taoyuan-hukou tableland, northwestern taiwan［j］. geomorphology, 2007, 90(1-2): 36-54.
［7］ wu w, fan y, wang z g, et al. assessing effects of digital elevation model resolutions on soil-landscape correlations in a hilly area［j］. agriculture, ecosystems and environment, 2008, （126）: 209-216.
［8］ odeh i o a, mcbratney a b, chittleborough d j. further results on prediction of soil properties from terrain attributes: heterotopic cokriging and regression-kriging［j］. geoderma, 1995, 67(3-4): 215-226.
［9］ eldeiry a a, garcia l a. comparison of ordinary kriging, regression kriging, and cokriging techniques to estimate soil salinity using landsat images［j］. journal of irrigation and drainage engineering, 2010, 136(6): 355-364.
［10］ josé m p, fernando v, josé l r. spatial evaluation of soil salininty using the wet sensor in the irrigated area of seuura river lowland［j］. journal of plant nutrition and soil science, 2010, doi: 10. 1002/jpln.200900221.
［11］ 赵成义, 王玉潮, 李子良, 等. 田块尺度下土壤水分和盐分的空间变异性［j］. 干旱区研究, 2003, 20(4):252-256.（zhao cheng-yi, wang yu-chao, li zi-liang, et al. study on the spatial variability of soil moisture content and salt content in the field scale［j］. arid zone research, 2003, 20(4):252-256.(in chinese)）
［12］ 杨劲松, 姚荣江, 刘广明, 等. 黄河三角洲地区土壤盐分的空间变异性及其cokriging估值［j］. 干旱区研究, 2006, 23(3): 439-445.(yang jin-song, yao rong-jiang, liu guang-ming, et al. spatial variability of soil salinity in the yellow river delta and its estimation by cokriging mthod［j］. arid zone research, 2006, 23(3): 439-445. (in chinese)）
［13］ 姚荣江, 杨劲松, 刘广明. 土壤盐分和含水量的空间变异性及其cokriging估值—以黄河三角洲的去典型地块为例［j］. 水土保持学报, 2006, 20(5): 133-138.(yao rong-jiang, yang jin-song, liu guang-ming. spatial variability of soil salinity and moisture and their estimations by cokriging method-a case study in characteristic field of yellow river delta［j］. journal of soil and water conservation, 2006,20(5):133-138.(in chinese)）
差值精度。
