多元指示克立格法在黄土丘陵沟壑区的应用

多元指示克立格法在黄土丘陵沟壑区的应用

干旱和半干旱地区耕地的土壤水分和养分是限制作物生长的主要因素之一。监测和预测土壤水分和养分的时空分布对指导这些地区的节水灌溉、碱土壤改良和农业可持续发展具有重要意义。由于“盐随水来,盐随水走”,土壤水分和盐分的时空分布及其动态并不是相互独立的过程,而是具有一定的协同相关性,即为协同区域化变量,因此,国内外的学者们也应用协同克立格法对土壤水分和盐分的空间分布进行了估计。然而,无论是常用的普通克立格,还是单元指示克立格,或协同克立格,只能给出某一时期土壤水分或盐分的空间分布等值线图,却不能同时从土壤的水分、盐分含量及其它特性几方面综合评价土壤的质量,如果再加上时域特征,分析起来就更为困难了,这就需要一种可以将这些时空变量整合在一起形成一个综合变量,并考虑综合变量空间结构性的方法,多元指示克立格正是这样一种方法。多元指示克立格法是利用多元指示变换(Multiple-VariableIndicatorTransform)整合各种指示变量(例如,水土资源的物理性质、化学性质、生物性质及其它描述其质量的变量都可以作为水土资源质量指示变量),利用非参数地质统计学——指示克立格方法(IndicatorKriging)估计满足某一阈值集合(由各个指示变量的阈值组成)条件下综合指示变量(CombinedIndicator)条件概率的一种新方法。该方法1993年由Smith等提出,此后,Halvorson等应用该方法分析了灌木-草原生态系统中土壤特性的空间变异模型,但将该方法应用于较大面积不同时期土壤水盐特性的分析少有报道。本文应用多元指示克立格法,结合单元指示克立格对河套灌区约55hm2农田的两个时期土壤水分和盐分的空间变异性进行了综合分析和评价,这是该法在中尺度农田的土壤质量评价中应用的一个很好的例子,对应用多元地质统计学理论综合各种影响因素及其空间分布,从而制定科学的农业管理措施,具有一定的借鉴意义和参考价值。1基本原则和方法1.1传统克立格估计算法土壤具有高度的空间异质性,其特性的观测数据中经常会出现一些特异值,它们比所研究的全部观测值的算术平均值或中位数值要高(或低)得多,但并非误差所致,特异值在全部观测值中虽只有少部分,但它对全部观测值的统计结果影响很大,往往使观测值偏离正态分布,造成变异函数的不稳健,从而影响克立格估计结果。通常的解决办法是通过非线性变换(如,平方根变换,对数变换,反正弦变换等)“平滑”这些特异值,对变换后的数据进行(地质)统计分析,最后再进行逆变换。但经非线性变换后的数据不能直接应用普通克立格等线性估计技术,除非有办法将克立格估计的数据转换回原尺度,否则就会产生有偏估计。但目前只有对数变换可以应用对数正态克立格加以转换,其他变换没有相应的转换方法。另一个途径是剔除特异值或用估计邻域和影响系数等方法处理特异值,但这些方法都在一定程度上破坏了变量空间变异的真实信息。第三种途径就是应用指示克立格法进行分析。指示克立格是条件普通克立格的非参数形式,因为它能抑制特异值对变异函数的影响,因此是处理有偏数据的有力工具。假设在区域D上采样测定了某变量,若该变量在区域D上的阈值或边界值为z,则在x∈D处可定义如下阶梯函数(即指示函数):I(x‚z)={?1Z(x)≥z?0Z(x)＜z(1)在给定的阈值Z的条件下,记Z(x)≥z(即I(x,z)=1)的概率为F(z)=Prob[I(x‚z)=1](2)则指示函数I(x,z)的期望值为E{I(x‚z)}=1×F(z)+0×[1−F(z)]=F(z)(3)式(3)表明x处Z(x)≥z(即I(x,z)=1)的概率等于指示变换值的平均值。当I(x+h,z)和I(x,z)为被矢量h分隔的两个指示变换数据点时,则指示变异函数可定义为γI(h‚z)=12E{[I(x+h‚z)−I(x‚z)]2}(4)可见用式(4)求得的指示变异函数,可通过普通克立格估计得到的未知点的I(x,z)值,即为Z(x)≥z在该点出现的概率,即F∗(Z)=∑α=1nλαI(xα)(5)式中F*为Z(x)≥z的概率F的克立格估计;n是在待估点x的邻域中的有效点数;xα是第α个采样点的位置;λα是克立格权系数。用式(5)可以估计Z(x)≥z的条件概率图。式(4)表明,指示变异函数是用指示函数计算得到的,计算中所有的指示变换值1(或0)对变异函数的贡献都是一样的,因此,不受特异值的影响。须注意的是,指示函数除式(1)外,还有其它形式,一般是出现预期结果变换为1,否则为0,选用何种形式,按具体问题而定。1.2指标阈值集的确定在水土资源的分析评价中,影响决策的指标往往不只一个,如何将众多的指标整合为一个综合指标,是非常重要的。多元指示变换是实现这一目标的有用的方法,其步骤为:①根据问题的要求(如,农业生产管理或环境要求等),为参与分析的所有指标(即变量)独立地确定一个阈值,每个阈值可是一个临界值,也可以是一个范围,所有指标的阈值构成阈值集,阈值原则上是可任选的;②为每个指标确定指示函数,并根据指示函数对相应的指标进行二态指示变换;③进行指标整合,确定综合指标,即当某一点上有m(m≤n,n为变量数)个指标的指示变换值为1,则该点的综合指示变换值为1,否则为0。经多元指示变换后得到的由综合指示变换值(即1或0)组成的数据集即为综合指标,可以用作评价相应采样点上水土资源质量优劣的指标。然而,这一方法更大的价值在于可以利用综合指标的指示变异函数(求法与单变量指示变异函数相同),通过普通克立格法估计未采样点的值,进而绘制整个研究区域上所有或部分指标(即变量)同时满足相应标准(或要求)的概率图。由于这种方法可以整合多个指标(即变量),得到综合克立格概率分布图,故称之为多元指示克立格(MVIK)。多元指示克立格虽然只是单元指示克立格的简单外延,却具有更广泛的用途。该方法本身并不定义水土资源的质量,而是整合那些代表水土资源质量的标准(即阈值),因而根据要解决的问题选择合适的标准很关键。2研究区土壤水分、土壤样品采集与分析试验区设在黄河河套灌区长胜试验站,该地区气候干旱,地下水埋深浅,水平排水不畅,是我国盐渍土发育的典型地区。试验区是正在改良的盐渍化土地,试验区总面积为1200hm2,重点试验区为规则的长方形,面积为55hm2,1998年秋浇前(9月中旬)与1999年夏灌前(4月下旬)沿东西、南北方向每隔100m用剖面法采集1m深度(分5层)的土样,样品分析成果包括土壤全盐量(占干土重%)、土壤水分含量、EC值、pH值以及土壤8个主要离子含量。本文以重点试验区秋浇前与夏灌前0～100cm内土壤全盐量及土壤水分含量的平均值为例,应用指示克立格法分析土壤水盐特性的空间分布,讨论多元指示克立格法的适用性。3结果分析3.1夏灌前水分变异情况表1为土壤水分和盐分的统计特征值,水分的变异性明显小于盐分。若按一般对变异系数值Cv的评估,当Cv≤0.1时,为弱变异性,0.1≤Cv≤1.0为中等变异性,则夏灌前水分属弱变异性,秋浇前水分属中等偏弱变异,盐分都属于中等变异性。水分的偏度都较接近于0,峰度也较接近于3,所以水分接近正态分布,但盐分的峰度和偏度都很大程度偏离了3和0,即严重偏离正态分布。实际上,盐分的频率分布图(图1只给出了夏灌前盐分的频率分布图)明显向左偏倚,不符合正态分布,经对数变换后,大致接近正态分布,但右侧由于特异值存在仍有长尾现象。3.2指挥克立格分析(1)夏灌前土壤盐分的行为统计特征(表1和图1)表明,土壤盐分不仅呈非正态分布,而且存在特异值,两者都会影响变异函数的稳健性,进而影响克立格估计结果;此外,尽管各种克立格方法(包括指示克立格和普通克立格等)能给出单个变量,即某一时期土壤盐分的空间分布图,但有时人们更感兴趣的不是某一点处盐分的具体含量,而是大于(或小于)的某一阈值,或处于某一范围内的盐分含量在空间上的分布概率。指示克立格法不仅能削弱有偏分布和特异值对变异函数及估计结果的影响,而且能有效地估计区域化变量的分布概率,成为解决这类问题的有力工具。以夏灌前盐分为例,进行单元指示克立格分析。其计算的第一步(也是比较关键的一步)是选择阈值,对于夏灌前盐分本文选0.3%,且置大于0.3%的含盐量为1,否则为0。选0.3%为阈值的原因是:河套灌区盐渍化土地夏灌前未播种,这时如知道哪些是盐渍化发生的高概率区,哪些是低概率区,就可为种植作物种类的选择提供依据(即盐渍化发生的低概率区种植耐盐性较差的玉米等作物,反之种植耐盐性较好的小麦、向日葵等);而按一般盐土的分级,1m深度内土层全盐量小于0.2%或0.3%为非盐渍化土壤,否则为盐渍化土壤,考虑本试验区为正在改良的盐渍化土壤,因而选较大的0.3%为阈值。图2和图3分别为夏灌前土壤盐分数据的指示变异函数图和指示克立格估计的盐分含量大于0.3%的概率分布图。图4是夏灌前土壤盐分的指示克立格估值分布图。指示克立格估计值分布图(图4)实际是根据条件概率计算得到的盐分含量期望值的分布图,由于该图的估计不受特异值和数据有偏分布影响,无需处理特异值,也无需对数据进行正态变换,能避免普通克立格估计中由于这些处理而引起的估计误差。实际上,交叉验证的结果(略)表明,用指示克立格法对夏灌前各实测点盐分进行估值,所得误差均方比普通(对数正态)克立格估计(即先进行对数正态变换,再用影响系数法处理特异值,然后计算变异函数,最后用对数正态克立格法估值)的误差均方缩小27%。另一方面,由于普通克立格强烈的“平滑效应”,使得所得盐分分布图(未给出)中左侧的重度盐渍化(全盐量大于0.7%)未完全反映出来,而实际试验区的西南角(图4中的左上角)与盐荒地相邻,盐渍化较为严重。说明指示克立格所得盐分分布图较普通克立格更能反映实际情况。但有一点须承认:指示克立格也有一定程度的平滑效应,只是不如普通克立格明显。从图3能清楚地看出夏灌前盐分含量大于0.3%的高概率区域占很大一部分,对这些区域可选择较耐盐的作物种植,且在夏灌时(播种前)对盐分含量大的高概率区(图3中的深色部分)要适当增加灌水量以冲洗盐分。对含盐量超过0.3%的低概率区(图3中的浅色部分)为争取农时可酌情考虑夏灌前提早播种。(2)夏灌前盐渍化情况的调查与防治上面是只考察一个时期土壤盐分的情形,如果要同时考察两个或两个以上时期土壤的盐分,或同时考察土壤水盐的分布时,就需要用到多元指示克立格。这里所谓的多元指示克立格是(单元)指示克立格的简单外延,为的是将多个诊断指标(即变量)整合为一个指标,以便于分析和估值。下面以秋浇前和夏灌前两个时期的盐分为指标诊断试验区盐分的分布情况,并说明多元指示克立格具体计算步骤:①分别为每个指标选择各自的阈值;本例两个时期的土壤盐分阈值都取0.4%,原因是:对于河套灌区夏灌前和秋浇前是两个较特殊的时期,秋浇前反映的是经整个作物生长期人工灌溉干预下土壤水分和盐分的分布状况,夏灌前非饱和带中土壤水分和盐分已经完成了冻融季节的再分配过程,基本代表冬春自然过程影响下的分布状态,如果某区域两个时期含盐量都较大,则基本上可以认为该区域盐渍化程度较严重;而按一般盐土的分级,1m深度内土层全盐量大于0.4,则表明盐渍化已达中度以上,中度盐渍化以上的土壤中作物正常生长就会受到抑制;②为每个指标建立指示函数,按指示函数实施指示变换;本例两指标选相同的指示函数形式,即式(1);③将多个指标合并为一个综合指标(或多个指示变换值合并为一个综合指示变换值),即当某采样点上m(m≤n,n表示变量数,m依问题不同而定)个指示变换值如果为1,则综合指标(或综合指示变换值)为1,否则综合指标为0;本例m=n=2,即采样点上两个时期的盐分含量都大于0.4%时,综合指标为1,否则为0;④计算综合指标的指示变异函数(图5(a));⑤进行普通克立格插值,得到综合指标的概率图(图5(b))。图5(b)中颜色较深的部分即为秋浇前和夏灌前两个时期盐分含量超过0.4%的高概率区,即这些区域盐渍化达中度以上的概率较大,需要采取适当的措施进行盐碱化防治。图6为秋浇前水分(大于中值21.3%变换为1,否则为0)和盐分含量(小于0.4%变为1,否则为0)的综合指示变异函数图与概率图。图6(b)的深色部分表示土壤水分含量较大(大于中值21.3%)、盐分含量较小(小于0.4%)的高概率区,所以,在这些区域秋季的储水灌溉(目的是造墒压盐)灌水定额可以低于其它区域,从而节约灌水量,减小地下水位抬高的程度。制约作物生长的土壤特性不仅是水分和盐分,还有各种养分含量、pH值等,多元指示克立格允许同时考虑这些因素,绘制综合概率图,可以为土壤或其它资源环境区域信息的合理化管理决策提供有力依据,它为综合、全面地评价土壤及其它资源质量提供了新的方法。4多变量量化方法在水土资源特性的数据分析中,往往由于特异值的存在使数据偏离正态分布,而不能得到稳健的变异函数,影响克立格估计(产生有偏估计)或条件模拟结果的精度。指示克立格法是解决这类问题的有效工具。本例中,夏灌前土壤盐分近似符合对数正态分布,但存在特异值,经指示克立格估值,其结果的误差均方比普通(对数)克立格缩小27%,且指示克立格估计的分布图比普通克立格更能反映真实情况。特别是当变量存在临界值或某一合理取值范围时,指示克立格估计的条件概率图能给出各点取值的风险概率(如本文的夏灌前