环渤海低平原水土盐分与水位埋深的空间变异及协同克立格估值

周在明,张光辉,王金哲,严明疆

中国地质科学院水文地质环境地质研究所,河北石家庄 050803

环渤海低平原水土盐分与水位埋深的空间变异及协同克立格估值

周在明,张光辉,王金哲,严明疆

中国地质科学院水文地质环境地质研究所,河北石家庄 050803

针对环渤海低平原水土盐分和水位埋深,通过选取0～20cm深度内的127个代表性土样进行土壤全盐量测定,并对130个水井的水位埋深及128个地下水样的矿化度进行测定。综合运用普通克立格(Ordinary Kriging)与协同克立格(CoKriging)方法与GIS技术研究水土盐分及水位埋深的空间分布状况。结果表明,土壤盐分和地下水矿化度的变异系数分别为0.91和 0.87,属于中等空间变异强度,而水位埋深的变异系数为1.06,属于强变异强度。单变量空间相关程度均属于中等,交互变量的空间相关性较强,其空间自相距离(135.6～223.2km)高于单变量(40.6～144.8km)。Kriging与CoKriging分析表明,环渤海低平原土壤盐分与地下水矿化度分布规律基本一致,而与水位埋深相反。即土壤盐分和地下水矿化度自内陆平原向东部滨海平原逐渐增加,不同地面高程上存在分布差异,并随着高程的降低呈增加的趋势。盐分含量以唐山—天津—沧州—东营—滨州一线较高,矿化度以天津、沧州部分地区居高,水位在黄河三角洲较浅。通过CoKriging插值的均方根误差与Kriging相比减少了0.29%～4.78%,而预测值与实测值的相关系数提高了2.03%～20.58%。

环渤海低平原;土壤盐分;地下水矿化度;水位埋深;协同克立格

水盐和水位埋深对地下水的利用程度及土壤盐分有一定的影响(Kahlown et al.,2005),土壤盐分、盐渍化对农业可持续生产造成重大威胁(Cetin et al.,2003),土地盐化程度的增加不但减少了作物产量而且限制了作物种类的选择(Singh et al.,2009),干旱半干旱地下水浅埋区、滨海低平原区尤为显著(Il’ichev et al.,2008;Guswa,2002)。对地下水和土壤盐分含量的准确估计和分布研究在区域和全球水土环境及相关研究中是一项重要的课题(Liu et al.,2005)。对于水土而言其可溶盐的分布在一定程度上受到地形的影响(Salama et al.,1999),地形特征对地下水及土壤形成与利用具有重要作用(Lin et al.,2007;Wu et al.,2008)。因此在地下水和土壤盐分含量的插值估算中区域地形因素的影响不容忽视。

协同克立格(CoKriging)插值方法充分利用了待求变量和易测变量的相关性(Odeh et al.,1995),借助辅助变量提高主变量的预测精度,在地下水和土壤盐分空间插值中得到广泛的应用,是公认的有效插值方法(Triantafilis et al.,2001;Darwish et al.,2007)。Kim等应用水化学、Kriging与CoKriging方法对首尔以南的五松地区地下水硝酸盐污染进行了评价,阐明了硝酸盐的空间分布及与其它物理化学参数的关系,并解释了 CoKriging在复杂地下水化学数据插值中的优越性(Kim et al.,2009)。Eldeiry等综合应用地面验证和遥感数据对美国阿肯色流域土壤盐分含量的Kriging和CoKriging插值方法进行了比较,认为CoKriging方法对玉米地的土壤盐分估计精度最高,其次是小麦和苜蓿(Eldeiry et al.,2010)。José等应用CoKriging方法对西班牙东南部塞奎拉河域土地,以土壤水电导率为辅助变量进行了 50cm深度的盐分预测研究,较好地反应了作物根系层的盐渍化状况(José et al.,2010)。在国内相关的研究以西北干旱内陆和黄河三角洲为主,结果表明与普通克立格相比 CoKriging既能够减少样点数量又能提高估值精度(赵成义等,2003;姚荣江等,2006)。

针对整个环渤海低平原水土盐分问题,进行大尺度区域布点,结合地形因素对水土盐分的空间分布格局进行研究的尚不多见。为此,本文以环渤海低平原为研究对象,以研究区高程为辅助变量应用协同克立格插值方法对水土盐分及水位埋深的空间分布进行研究,为提高水土盐分含量的插值估算精度,为地下水和土地资源管理利用和农业生产的合理布局提供科学依据。

1 材料与方法

环渤海低平原位于华北东部,包括河北、山东内陆低平原和滨海低平原两大类型,含天津、沧州等80多个县市。由黄河、海河、滦河冲积而成,地势低平,大部分海拔50 m以下,滨海区10 m左右。总面积9.6×104km2;总耕地面积6000多万亩,是我国重要的粮棉、果蔬产区。

环渤海低平原属于欧亚大陆东岸暖温带半干旱季风气候区,冬春寒冷干燥,夏季炎热多雨。多年平均气温12.2℃,平均年降水量500～600mm,全区年内降水分配不均,主要集中在 6—9月,占全年降水量的60%～80%,年均蒸发量900～1400mm,干燥度达1.5左右。

1.1 数据来源

根据该区的土地利用类型和地下水状况,全区共测量水井130个,布设水样点128个、土样点127个(图1),取样点间距约 20km。取样范围为(N36°03′～N39°35′;E114°36′～E119°28′),取 样 面 积8.97×104km2。取样时间为2010年4—5月,每个点采用GPS定位,土样深度为0～20cm。

水位埋深用皮尺与测绳测定,水样矿化度采用重量法测定。土样全盐量测定过程如下,称取过2mm筛的风干土试样50～100 g,按土水比1:5配制浸出液。采用蒸干法,添加15%双氧水溶液与2%碳酸钠溶液进行测定,详细操作方法参见《土壤农业化学分析方法》(鲁如坤等,1999)。

以1:5万地形图(等高距20 m)与SRTM高程数据为基础,用实测高程点进行高程验证(图2),最终确定高程点1120个,平均间距约6.7km,覆盖研究区主要地貌单元。利用均方根误差描述高程精度,Zi为高程数据真值,zi为计算值,n为误差个数,计算得RMSE=6.52,说明所得到的高程空间插值结果满足数据分析要求。

图1 环渤海低平原取样点分布图Fig.1 Distribution of soil and groundwater samples in the low plain around the Bohai Sea

图2 环渤海低平原高程分布图Fig.2 Digital Elevation Model of the low plain around the Bohai Sea

1.2 普通克立格与协同克立格方法

在地统计学中,半方差函数的一些重要参数如块金值、基台值和变程等可以用来表示区域化变量在一定尺度上的空间变异和相关程度,它是研究土壤特性空间变异的关键,也是克立格插值的精度要素(Yang et al.,2008;姜勇等,2005)。在本征平稳假设下,半方差计算公式为(式1)

式中,r(h)为步长h的半方差函数,N(h)是间距为h的计算对数,Z(xi)和Z(xi+h)分别是区域化变量Z(xi)和Z(xi+h)在空间位置xi和xi+h处的实测值。

两个随机变量的协同区域化可以用交互半方差函数来表示(式2)

式中,rij(h)是两个变量的交互半方差值,N(h)是具有相同间距h的变量Zi(x)和Zj(x)的离散点的数目。如果两个变量是正相关的,那么变量Zi从xa到xa+h的增加或减少会引起Zj的增加或减少,交互半方差就是正值。

如果变异函数和相关分析的结果表明某一属性的空间相关性存在,则可以利用普通克立格进行插值(式 3),

式中,z*(x0)是待估点处的估计值,z(xi)是实测值,λi是分配给每个实测值的权重且∑λi=1。n是参与点估值的实测值的数目。

协同克立格是普通克立格的扩展形式,它要用到两个或两个以上的变量,其中一个是主变量,其它的作为辅助变量,将主变量的自相关性和主辅变量的交互相关性结合起来用于无偏最优估值中。其公式为(式4),

式中,z*(x0)是待估点x0处的估计值,z1(xi)和z2(xj)分别是主变量z1和辅助变量z2的实测值,λi和λj分别是分配给主变量z1和辅助变量z2的实测值的权重,且。n和p是参与x0点估值的主变量z1和辅助变量z2的实测值数目。

1.3 数据处理

采用SPSS 16.0软件进行水土全盐量及水位埋深的统计分析,变异函数及协同变异函数模型的拟合采用地学统计软件 GS+7.0,克立格与协同克立格空间插值应用ArcGIS 9.3软件。

2 结果分析

2.1 描述性统计特征

根据测定结果对地下水矿化度、水位埋深和土壤盐分含量进行统计(表1)。环渤海低平原区耕层土壤盐分含量均值为1.08 g/kg属于轻度盐渍化土(王遵亲等,1993),地下水矿化度平均值为2.06 g/L属于微咸水,水位埋深平均为16.61 m。由表1可见,环渤海低平原区的地下水埋深属于强变异强度(变异系数为1.06),地下水矿化度与土壤盐分的空间变异强度属于中等,但也较大(变异系数分别为0.87和0.91)。这主要是由于研究区独特的水文、地形结构以及复杂的人类活动所致。由于地下水过量开采出现区域性地下水位变化不均,内陆低平原浅层地下水位下降幅度较大,东部滨海平原区浅层地下水下降幅度较小,地下水超采是地下水位埋深空间变异加大的主要原因。从内陆到滨海地下水矿化度逐渐增大,水质变差,再加上人为排灌的影响,其共同作用使得平原区地下水矿化度空间变异相对较大。而研究区范围较大且存在地势差异、土地耕种方式差异、田块灌溉制度差异等因素的共同作用导致耕层土壤盐分变异性较大。

表1 环渤海低平原土壤盐分、地下水矿化度和水位埋深统计特征值Table 1 Descriptive statistics of soil salinity,TDS and GWL in the low plain around the Bohai Sea

2.2 高程与观测项目的相关性

根据相应取样点的高程分别与土壤盐分、地下水矿化度和水位埋深进行 person相关分析发现,其相关性分别满足 0.05显著性水平下的显著负相关,相关系数−0.253;0.01显著性水平下的极显著负相关,相关系数−0.394;和 0.01显著性水平下极显著的正相关,相关系数为0.322。这说明土壤盐分含量、地下水矿化度、水位埋深在不同高程上存在分布差异,并随着高程的增大,土壤盐分含量和地下水矿化度呈减少的趋势,水位埋深则呈现增加的趋势。在研究区内高程与土壤盐分、地下水矿化度、水位埋深受区域化现象或空间过程的影响,属于协同区域化变量,且应用偏度、峰度联合计算法进行检验发现,通过对数转换后 P＞0.05,因此区域高程属于对数正态分布类型。

2.3 普通克立格与协同克立格插值

根据半方差函数理论,计算得到半方差函数拟合模型,其决定系数在0.58～0.96之间均达到显著性水平(表2)。可见,土壤盐分、矿化度和水位埋深半方差和交互半方差函数可分别应用球状模型和指数模型进行拟合。由块金值与基台值的比值C0/Sill可知,土壤盐分、矿化度和水位埋深单变量空间相关程度均属于中等,其值分别为41.78%、25.68%和40.00%,对交互变量而言,除土壤盐分的空间相关性属于中等外,矿化度和水位埋深的空间相关性均较强。此外在空间自相关距离上交互变量(135.6km、200.2km 和 223.2km)也高于单变量(40.6km、144.8km和90.0km)。总体上看,土壤盐分、地下水矿化度和水位埋深的交互变量半方差函数理论模型较之单变量都有一定程度的改变。

分别根据半方差和交互半方差函数模型对研究区土壤盐分、矿化度和水位埋深数据进行插值,获得其空间分布图(图3)。本研究中 Kriging 和CoKriging是基于环渤海低平原区80多个县市的实测数据做出的,在这种范畴广、数据点较少的条件下,我们制作分布图时,尽量通过参数的选择使均方根标准误(RMSSE)接近于 1,以保证预期误差的变异尽可能地小,从而使得到的分布图是我们现有数据资料的最佳和最优效的预期结果。

表2 土壤盐分、地下水矿化度和水位埋深的半方差、交互半方差函数模型Table 2 Models of semivariogram and cross-semivariogram for soil salinity,TDS and GWL

由图3可见,两种插值方法所获得的分布图在整体趋势和具体斑块形状上基本相似,表明协同克立格和普通克立格都能够较好地反应研究区土壤盐分、矿化度和水位埋深的空间分布情况。环渤海低平原土壤盐分含量,自内陆平原向东部滨海平原逐渐增加,盐分含量较高的地区出现在唐山—天津—沧州—东营—滨州一线。盐分含量小于1 g/kg的非盐化土以内陆平原为主,分布在保定—衡水—邢台—邯郸一线。这与地下水矿化度的分布规律基本吻合,地下水矿化度大于5 g/L咸水主要分布在东部滨海天津、沧州地区,1～5 g/L的微咸水分布于广大内陆平原区,其盐分含量也是逐渐向滨海地区增大,而水位埋深与之相反,在黄河三角洲较浅。

从地形因子看,土壤盐分和地下水矿化度高的地区地势平坦,海拔较低,水位较浅;与之相反的是海拔与地势起伏相对较大的内陆区,土壤盐分和地下水矿化度不大。从成土母质上看,滨海平原的土体一部分是由河流入海冲积成的三角洲如黄河三角洲,一部分是海积平原。质地多为粉砂、细砂组合,粗粉砂含量多达 60%以上。这种砂壤土体构型使得水盐迁移量高于内陆平原的粘壤土和亚砂土的组合。从农业活动看,河北平原区的灌溉水源以抽取深层地下淡水为主,山东平原区以引黄灌溉为主,并建有完善的水利工程,有合理的耕作制度,因此农田长期处于脱盐状态,盐分累积小。滨海平原区农田主要依靠自然降雨,长期的地面蒸发使得盐化加重。

图3 环渤海低平原土壤盐分(a、a’)、矿化度(b、b’)和水位埋深(c、c’)空间分布(a、b、c 普通克立格插值;a’、b’、c’协同克立格插值)Fig.3 Spatial distribution of soil salinity (a,a’),TDS (b,b’)and GWL (c,c’)in the low plain around the Baohai Sea(a,b ,c-Ordinary Kriging;a’,b’,c’-Cokriging)

2.4 两种插值方法预测精度比较

用评价方法(协同克立格CK)的均方根误差相对于参考方法(普通克立格 OK)的均方根误差减少的百分数(RRMSE)表示预测精度的提高程度(式6)。

用RR表示评价方法相对于参考方法相关系数的提高程度(式7)。

式中,RMSEOK和ROK分别是参考方法的预测均方根误差及预测值与实测值间的相关系数,RMSEOK和RCK分别表示评价方法的预测均方根误差及预测值与实测值间的相关系数。

由表3可见,土壤盐分、矿化度和水位埋深通过协同克立格插值的均方根误差与普通克立格插值相比分别减少了0.29%、2.18%和4.78%;而预测值与实测值的相关系数分别提高了 20.58%、2.03%和11.31%。表明在相同的取样条件下,协同克立格插值由于融合了更丰富的空间信息,预测的精度要高于普通克立格插值,可优化样点的插值精度。

表3 土壤盐分、地下水矿化度和水位埋深的普通克立格与协同克立格精度比较Table 3 Comparison of prediction accuracy of soil salinity,TDS and GWL between Ordinary Kriging (OK)and CoKriging(CK)

3 结论

环渤海低平原区土壤盐分和地下水矿化度属于中等空间变异强度(Cv值分别为0.91和0.87)而水位埋深属于强变异强度(Cv值为1.06)。半方差和交互半方差函数可分别应用球状和指数模型进行拟合。单变量空间相关程度均属于中等,交互变量的空间相关性较强,空间自相距离(135.6～223.2km)明显高于单变量(40.6～144.8km)。

克立格与协同克立格空间插值表明,环渤海低平原土壤盐分与地下水矿化度分布规律基本一致,而与水位埋深相反。即土壤盐分和地下水矿化度自内陆平原向东部滨海平原逐渐增加,不同地面高程上存在分布差异,并随着高程的降低呈增加的趋势。盐分含量以唐山—天津—沧州—东营—滨州一线较高,矿化度以天津、沧州部分地区居高,水位在黄河三角洲较浅。由于高程与上述观测项目存在很好的相关性,通过协同克立格插值的均方根误差与普通克立格插值相比减少了 0.29%～4.78%,而预测值与实测值的相关系数提高了2.03%～20.58%。

姜勇,庄秋丽,梁文举,施春健,欧伟.2005.空间变异在土壤性质长期定位观测及取样中的应用[J].土壤通报,36(4):531-535.

鲁如坤.1999.土壤农业化学分析方法[M].北京:中国农业科学技术出版社.

王遵亲,祝寿泉,俞仁培.1993.中国盐渍土[M].北京:科学出版社.

姚荣江,杨劲松,刘广明.2006.土壤盐分和含水量的空间变异性及其 CoKriging估值——以黄河三角洲地区典型地块为例[J].水土保持学报,20(5):133-138.

赵成义,王玉潮,李子良,李国振.2003.田块尺度下土壤水分和盐分的空间变异性[J].干旱区研究,20(4):252-256.

CETIN M,KIRDA C.2003.Spatial and temporal changes of soil salinity in a cotton field irrigated with low-quality water[J].Journal of Hydrology,272:238-249.

DARWISH K H M,KOTB M M,ALI R.2007.Mapping soil salinity using collocated cokriging in Bayariya,Oasis,Egypt[C]//Proceedings of the 5th International Symposium on Spatial Data Quality,ITC Enschede,The Netherlands.

ELDEIRY A A,GARCIA L A.2010.Comparison of ordinary Kriging,regression Kriging,and Cokriging techniques to estimate soil salinity using Landsat images[J].Journal of Irrigation and Drainage Engineering,136(6):355-364.

GUSWA A J.2002.Models of soil moisture dynamics inecohydrology:a comparative study[J].Water Resources Research,38(9):1-15.

IL’ICHEV A T,TSYPKIN G G,PRITCHARD D,RICHARDSON C N.2008.Instability of the salinity profile during the evaporation of saline groundwater[J].The Journal of Fluid Mechanics,614:87-104.

JIANG Yong,ZHUANG Qiu-li,LIANG Wen-ju,SHI Chun-jian,OU Wei.2005.Application of sptial variability in long-term site-specific observatory study of soil properties and sampling strategy[J].Chinese Journal of Soil Science,36(4):531-535(in Chinese with English abstract).

JOSÉ M P,FERNANDO V,JOSÉ L R.2010.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,doi:10.1002/jpln.200900221.

KAHLOWN M A,ASHRAF M,ZIAUL H.2005.Effect of shallow groundwater table on crop water requirements and crop yields[J].Agricultural Water Management,76:24-35.

KIM K H,YUN S T,CHOI B Y,CHAE G T,JOO Y S,KIM K,KIM H S.2009.Hydrochemical and multivariate statistical interpretations of spatial controls of nitrate concentrations in a shallow alluvial aquifer around oxbow lakes (Osong area,central Korea)[J].Journal of Contaminant Hydrology,107:114-127.

LIN Y S,LIN Y W,WANG Y,CHEN Y G,HSU M L,CHIANG S H,CHEN Z S.2007.Relationship between topography and spatial variations in groundwater and soil morphology within the Taoyuan-Hukou Tableland,Northwestern Taiwan[J].Geomorphology,90(1-2):36-54.

LIU X Y,PETERSON J,ZHANG Z Y,CHANDRA S.2005.Improving soil salinity prediction with high resolution DEM derived from LIDAR data[J].International Archives of the Photogrammetry,Remote Sensing and Spatial Information Sciences,36(7):41-43.

LU Ru-kun.1999.Analytical methods of soil agricultural chemistry[M].Beijing:China Agricultural Science Press(in Chinese).

ODEH I O A,MCBRATNEY A B,CHITTLEBOROUGH D J.1995.Further results on prediction of soil properties from terrain attributes:Heterotopic cokriging and regression-kriging[J].Geoderma,67(3-4):215-226.

SALAMA R B,OTTO C J,FITZPATRICK R W.1999.Contributions of groundwater conditions to soil and water salinization[J].Hydrogeology Journal,7:46-64.

SINGH R B,CHAUHAN CPS,MINHAS P S.2009.Water production functions of wheat (Triticum aestivum L.)irrigated with saline and alkali waters using double-line source sprinkler system[J].Agricultural Water Management,96:736-744.

TRIANTAFILIS J,ODEH I O A,MCBRATNEY A B.2001.Five geostatistical models to predict soil salinity from electromagnetic induction data across irrigated cotton[J].Soil Science Society of America Journal,65:869-878.

WANG Zun-qin,ZHU Shou-quan,YU Ren-pei.1993.Saline Soil in China[M].Beijing:Science Press(in Chinese).

WU W,FAN Y,WANG Z Y,LIU H B.2008.Assessing effects of digital elevation model resolutions on soil-landscape correlations in a hilly area[J].Agriculture,Ecosystems and Environment,126:209-216.

YANG J,HUANG Z C,CHEN T B,LEI M,ZHENG Y M,ZHENG G D,SONG B,LIU Y Q,ZHANG C S.2008.Predicting the probability distribution of Pb-increased lands in sewage-irrigated region:A case study in Beijing,China[J].Geoderma,147:192-196.

YAO Rong-jiang,YANG Jin-song,LIU Guang-ming.2006.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,20(5):133-138(in Chinese with English abstract).

ZHAO Cheng-yi,WANG Yu-chao,LI Zi-liang,LI Guo-zhen.2003.Study on the spatial variability of soil moisture content and salt content in the field scale[J].Arid Zone Research,20(4):252-256(in Chinese with English abstract).

Spatial Variability of Soil Salinity,Total Dissolved Solid and Groundwater Depth Based on Cokriging in the Low Plain around the Bohai Sea

ZHOU Zai-ming,ZHANG Guang-hui,WANG Jin-zhe,YAN Ming-jiang
Institute of Hydrogeology and Environmental Geology,Chinese Academy of Geological Sciences,Shijiazhuang,Hebei050803

Taking into account the problem of soil salinity,total dissolved solid (TDS)of groundwater and groundwater level (GWL)existient in the low plain around the Bohai Sea,the authors collected soil samples at the depths of 0～20cm from 127 sites in the plain,determined the soil salinity,and measured shallow groundwater depth in 130 sites andTDS ofwater samples in 128 sites.Classical statistical and geostatistical methods combined with GIS technique were used to analyze spatial variability of soil salinity,TDS and GWL.The results show that spatial variability of soil salinity and TDS belongs to the moderate degree while GWL belongs to the strong degree with Cv being 0.91,Cv being 0.87,and Cv being 1.06.The semivariogram and cross-semivariogram of soil salinity,TDS and GWL were fitted by the spherical and the exponential model respectively.Spatial correlation of single variables belongs to the moderate degree,while cross variable belongs to the strong degree,with spatial correlation distance 135.6km to 223.2km longer than 40.6km to 144.8km.Ordinary Kriging and CoKriging maps show that soil salinity and TDS spatial distribution are almost the same,increasing from the inland plain to the east coastal plain,spatial distribution difference exists in different altitudes,and the soil salinity and TDS rise with the decrease of the altitude;nevertheless,things are just the opposite for GWL.Soil salinity is high in Tangshan-Tianjin-Cangzhou-Dongying-Binzhou,TDS is high in some areas of Tianjin and Cangzhou,and GWL is low in the Yellow River Delta.In comparison with the ordinary Kriging with the same sampling numbers,the root-mean-square error produced by CoKriging decreases by 0.29%～4.78%,while the correlation coefficient between the predicted value and the measured value increases by 2.03%～20.58%.

low plain around the Bohai Sea;soil salinity;total dissolved solid;groundwater level;CoKriging

P641.51;P322.3

A

10.3975/cagsb.2011.04.14

本文由国家科技支撑计划项目(编号:2007BAD69B02;2009BADA3B05)和河北省科技厅重点基础研究项目(编号:08966711D)联合资助。

2011-03-16;改回日期:2011-04-15。责任编辑:魏乐军。

周在明,男,1980年生。博士研究生。主要从事土壤水盐管理的研究工作。E-mail:tougaozhou@163.com。