内蒙古河套灌区土壤水盐空间变异性的ik法阈值选择

内蒙古河套灌区土壤水盐空间变异性的ik法阈值选择

指导思想的不规则性（缩写ik）法是条库尔提出的一种非参数估计方法。史密斯和霍夫最初用于土壤和生态科学。IK法与其他插值估计方法比较有明显优点:它不依赖空间现象的平稳性;也不要求区域化变量服从某种分布假设,由于对随机变量经过特有的指示变换后,使随机函数的贡献都一样,可自然消除特异值(Outlier)的影响;可对非取样点处的不确定性进行条件概率分布函数(ConditionalCumulativeDistributionFunction简写CCDF)估计,尤其是当人们感兴趣的不是某些点的具体估值,而是大于或小于某阈值(CutoffValue)或某范围内物质品位在空间的分布概率时,可简捷获得在一定风险条件下未知量分布概率的估值,符合地学科学家用人工智能思维处理空间问题的思路。IK法自创立以来,在国外广泛用于地质、土壤、水文(Neuman)等领域,成为一种直接估计局部区域品位分布的数学地质新方法。20世纪90年代开始在我国地学及土壤学科应用,也有一些报导:如地质矿产方面侯景儒、土壤水文方面李保国等,水库工程地质方面吴蓉等,农业水土工程方面刘全明等、徐英等均有一定成效,但多为局部估计的初步试探,而对IK法的运用规则及区域性IK法阈值及其函数关系的系统研究尚少见。本文以黄河大型灌区(河套)土壤水盐监测实验的水土资源空间变异为案例,重点分析IK方程三项主要结构因子:指示阈值、指示变异函数和条件概率分布函数的变化趋势、相互关系和取值规则,为IK法的应用提供参考。1ik平均估计的确定设在某实验区D上取样并测定其品位观测值,若该区域上的边界门限(即阈值)为Z,则D内每一个样品点x∈D上定义一个Z的阶梯函数,即指示函数:i(x;Ζ)={1当x点上的观测值Ζ(x)≤Ζ0当x点上的观测值Ζ(x)>Ζ(1)D内任一区域A∈D,低于阈值Z的观测值Z(x)所占A区的比例表示为:ϕ(A;Ζ)=1A∫Ai(x;Ζ)dx∈(2)式中ϕ(A;Z)是关于Z(x)和Z的二元函数,即小于阈值Z的全部i(x;Z),(x∈A)的平均值。指示函数i(x;Z)在给定阈值Z的条件下服从二项分布,期望值是:E{Ι(x;Ζ)}=1⋅Ρrob{Ζ(x)≤Ζ}+0⋅Ρrob{Ζ(x)>Ζ}=Ρrob{Ζ(x)≤Ζ}=F(Ζ)(3)式中F(Z)称观测值Z(x)的分布函数在阈值Z处的数值。当I(x+h;Z)和I(x;Z)是被矢量h分隔的两个随机变量,即两个指示变换的数据时,则指示变异函数(IndicatorVariogrom)定义为:γΙ(h;Ζ)=12E{[Ι(x+h;Ζ)-Ι(x;Ζ)]2}(4)用指示变异函数可描述IK的空间变异程度,并可获得条件概率分布函数(CCDF)的估计值:F*[z|(n)]=prob*[Ζ≤z|(n)]=1nn∑α=1i(Ζα;Ζ)‚∈[0‚1](5)这个估计能被接受的前提是:各样本观测值Zα=zα是独立的,它们每一个对Z的作用是相同的,否则应采用不等权估计:F*[z|(n)]=prob*[Ζ≤z|(n)]=n∑α=1λα(z)⋅i(za；z)α=1,⋯,n(6)可以看出指示值i(zα;z)是用Zα-zα单独估计Z时的ccdf:F*[z|(n)=zα]=prob*[Ζ≤z|Ζ=Ζα]=i(zα;z)α=1,⋯,n(7)故式(6)可以写成:F*[z|(n)]=prob*[Ζ≤z|(n)]=n∑α=1λα(z)⋅F*(z|Ζ=Ζα)(8)上式说明:用多个观测值对CCDF做出的估计等于用单个观测值进行估计时所得CCDF的加权平均。指示函数值的方差S2(Z),即指示函数的基台值C0,可从γ*I(h;Z)的图形求得,故分布函数F*(z|(n))可从解以下二次方程得到F2(z|(n))-F(z|(n))-C0=0(9)同普通Kriging(OK)法一样,再运用线性无偏最优插值和估计方差最小的OK要求,结合指示函数的有关条件建立类似普通克立格(OK)方程的IK方程组式(10),可获得待估域A的IK平均估计。{n∑β=1λβ¯γi(xα,xβ;Ζ)+μ=¯γi(xα,A;Ζ)n∑α=1λα=1(α=1,2,⋯,n)(10)用式(10)可求出这些权系数λα及未知点的i*(x;Z)估值,即Z(x)≤Z在该点出现的概率和整个CCDF的估值,可不必剔除特异值而获得具有一定风险条件下的未知Z(x)的估计量及空间分布。2案例分析2.1土样及测定方法研究区域设在黄河河套灌区内的沙壕渠试验区。试验区采样面积为4km2,用手持GPS接收机定位。以间隔h=200m(采样点数为121个,中等田间尺度)的网状进行采样,采样分两层取土,深度为0～20、20～40cm。土样的测定分析项目为含水率θ和电导率EC、pH值。其中EC的测定采用水土比5∶1浸提,电导法测定水溶性盐总量。采样点布置见图1。2.2数据的比较及预处理IK法关键是阈值的合理确定,目前还没有从理论上研究水盐的阈值合理选择问题。阈值的选择主要是根据研究目标的专业门限要求、估计误差和允许风险大小进行选择,在缺乏经验时,有人建议用中位数作为初步估计参考值,经作者初步实践认为:中位数并非最佳阈值,估计结果具有平滑性,且接近OK法,也可能会产生IK法与OK法无明显区别的误导,为此应进行不同方案的比较及探索研究。表1水分数据的偏态系数接近于零,表明数据分布接近于正态分布的左右对称,变异系数小反映水分变异性不大;表2盐分数据的偏态系数均大于零,表明数据分布呈右偏态,变异系数大反映盐分变异性较大;水盐数据的峰度系数均大于零,表明数据分布比正态分布高耸而狭窄。按照均值的3倍标准差邻域判断,发现盐分数据有明显特异值存在,0～20cm土层水分有特异值存在,适宜用IK法估计。对于大、中尺度样本应根据不同作物、不同土壤和气候及地理条件,统一选定合理阈值,如水分阈值可选用田间持水量、凋萎系数、适宜含水量等;盐分阈值可采用FAO推荐的作物全生育期平均土壤浓度临界值、作物能容忍的土壤最大盐分值等。本文综合考虑数据的分位数特点,水分初步采用20%、22%、25%、27%4个阈值;盐分初步采用0.2、0.4、0.6、0.8、1.0、1.2mS/cm6个阈值进行参数的分析与比较,旨在探索选择合理水盐阈值的原则。下文仅列出0～40cm土层土壤的水盐计算成果。3阈值与函数关系的研究3.1水盐阈值对盐胁迫下自适应的影响水盐指示条件概率的计算采用ARCGIS9.2软件的地质统计学模块完成。表3中水分预测平均条件概率随着阈值的增加由99.2%减少至18.2%,反映了超过大阈值的概率逐渐减少;RMSE结果表明随着阈值的增加,均方根误差由9.5%增加至22.4%、51.1%,后又减小至39.0%,有先增加、后减小的变化趋势;表3中盐分预测平均条件概率随着阈值的增加由82.0%减少至1.8%,反映了超过大阈值的概率逐渐减少;RMSE结果表明随着阈值的增加,均方根误差由38.4%增加至52.2%,后减小至9.4%,也有相似的先增加、后减小的趋势。从图2可以直观地看出水盐阈值与指示概率均值的关系:水分预测平均概率随着阈值的增加而减少,减小的速度不同,超过22%阈值对应的指示概率均值减小速度加大。盐分预测平均概率随着阈值的增加而减少,减小的速度在0.6mS/cm阈值处变慢。水分预测概率均方根误差随着阈值的增加,先增加后减少,在25%～26%(接近中位数25.5%)区间概率预测误差最大,接近50%,这将会降低风险评价与IK法估计的可靠性,在20%阈值处概率预测误差最小,可见单从提高概率预测精度考虑,水分指示阈值应该尽量选择低于中位数的偏小值。盐分预测概率均方根误差随着阈值的增加,先增加后减少,在0.4mS/cm(接近中位数0.38mS/cm)处最大超过了50%,也将会降低风险评价与IK法估计的精度,可见单从提高概率预测精度考虑,盐分指示阈值应该尽量选择超过中位数的偏大值。图3—6为不同阈值下土壤水盐指示概率预测图。分析后发现高阈值的大概率区均包含在低阈值的大概率区内,其分布位置随阈值的增加而逐渐减少和集中,可用于不同阈值目标下的概率与风险评价。分析水分阈值25%、27%的指示概率预测图,含水率超过阈值25%、27%的概率在试验区中部及左上角较大,含水率超过27%的概率最大值为50%,有4个明显大概率区均包含在超过25%的大概率区域内。分析盐分阈值0.4、0.6mS/cm相应的指示概率预测图,可以了解盐分积聚的可能性,如电导率超过阈值0.6mS/cm的概率最大为36%,分布在试验区北部两侧及西南角,且均包含在超过阈值0.4mS/cm的大概率区域内。3.2中性关系变异函数性变异性函数的检测由表3可见:当水分阈值由小逐渐增大时,指示变异函数值有一种由小到大、再到小的趋势,转折点在25%处,可见接近中位数阈值的指示变异函数值最大,其指示变异函数结构性最差,所反映的指示概率的空间变化性最大,由此得到的指示概率预测误差最大,这与上文RMSE在阈值接近中位数处最大的结果相一致。盐分指示阈值与指示变异函数的关系同水分,指示变异函数值转折点在0.4mS/cm处。3.3阈值对结构性分析结果的影响从表5分析水分指示变异函数的结构性,阈值22%的指示变异函数结构性表现为纯随机,因为其C0/(C0+C)=1;结构性较强对应的阈值为20%、25%、27%,随着阈值增大,结构性有由强到弱,再由弱到强的变化趋势,为非降函数,其他两层指示变异函数的结构性分析具有相同的结论。从表6分析盐分指示变异函数的结构性,无纯随机现象,随着阈值增大,结构性有由强到弱,再由弱到强的变化趋势,0.4～0.6mS/cm区间内阈值对应的指示变异函数结构性最差,其指标值接近于1。4与其他方法插值对比本文选用含水率25%、电导率0.6mS/cm分别作为水盐的指示阈值进行IK法区域水盐估计,首先借助于IK方程组求解ϕ(A;Z)的线性估计量ϕ*(A;Z);然后计算区域A平均值的估计值:Z**(A;Ζe)=n∑e=1[ϕ*(A;Ze)×Z*(A;Ze)],最后与OK法、传统反距离加权平均法(简写IDW法)的插值结果做了对比。图7—8为土壤水盐IK法插值图,可以反映试验区的水盐分布条件。水分分布比较连续,有3处明显湿润区;盐分在西南、东北及试区中部(对角线方向)有明显条带状积聚。从表7可以分析3种方法的水分预测均值与实测值基本相同,在121个采样点上的预测平均误差为1.6%,说明了IK法同样具有精确的插值功能。IK法的估计方差、变异系数最大,但却小于采样点对应参数,说明IK法对水分普通克立格插值平滑效应有一定程度的消弱。从表8可以分析3种方法的盐分预测均值与实测值基本相同,在121个采样点上的预测平均误差为0.230mS/cm,说明了IK法同样具有精确的插值功能。IK法的估计方差、变异系数大于OK法,但却小于采样点对应参数,说明IK法对盐分普通克立格插值平滑效应有一定程度的消弱。5加指示概率估计的优缺点(1)IK法的难度在于阈值的合理选择,不仅涉及研究目标和专业门限要求,也与研究允许风险概率及变异程度有关,目前尚无理论分析方法,多数情况仍然需要试算探求。(2)经过模型与参数分析,发现水盐指示Kriging平均概率与指示阈值、指示变异函数的关系为:随着阈值的增加指示Kriging平均概率减少;当阈值由小逐渐增大时,指示变异函数值有一种由小到大,再到小的趋势。(3)作者认为选择中位数作阈值并非最佳估计,接近中位数阈值的指示变异函数值最大,其指示变异函数结构性最差,所反映的指示概率的空间变化性最大,由此得到的指示概率预测误差最大。本