CN115618720A - 一种基于海拔高度的土壤盐渍化预测方法及*** - Google Patents

Abstract

本公开涉及土壤盐渍化预测技术领域，提供了一种基于海拔高度的土壤盐渍化预测方法及***，包括：获取待预测区域的海拔高度；基于海拔高度，采用logistic回归模型，对待预测区域的土壤盐渍化二元状况进行预测；其中，logistic回归模型的构建方法为：基于采集的若干个采样点的海拔高度和土壤盐渍化二元状况，通过逐分量的Metropolis‑Hastings抽样方法从后验分布中产生随机数，并结合提议分布按分量进行逐个更新logistic回归参数，通过似然函数更新先验分布。实现土壤盐渍化的随机预测预报。

Description

一种基于海拔高度的土壤盐渍化预测方法及***

技术领域

本公开属于土壤盐渍化预测技术领域，尤其涉及一种基于海拔高度的土壤盐渍化预测方法及***。

背景技术

本部分的陈述仅仅是提供了与本公开相关的背景技术信息，不必然构成在先技术。

盐渍化的发生和盐渍土的形成，主要是受自然、人为等多种因素的影响，尤其受水分和盐分在土体中迁移与再分配的水盐运动过程所控制。土壤水盐运动过程及其调控机理是目前盐渍土研究的核心问题之一，盐渍化土地的产生是气候、土壤、地形地貌和人为条件等因素综合作用的结果。土壤盐渍化预测预报，从广义来说，是根据一个地区土壤、地下水水盐动态监测资料、地形地貌、河流灌溉综合状况对已经发生盐化土壤的发展方向及非盐渍化土壤产生盐渍化(俗称次生盐渍化)的可能性及程度进行预测预报；从狭义上来讲，指次生盐渍化的预测预报。总体来看，盐渍化预测预报研究可分为三个层次：

(1)通过对自然环境条件及土壤盐渍化发生、发展规律的研究进行定性预报。通过对预报区与己发生盐渍化区自然条件的比较分析，根据专家知识和经验进行土壤盐渍化可能性预报。常用的方法有地理相似法和专家预报法两种，地理相似法主要是根据已掌握的有关资料，作出定性预报，提出产生次生盐渍化的可能性，该方法难以作出较为详细的定量预报，因此是一种比较粗略的预报方法；专家预报法主要根据专家经验和当地的自然条件特点，对盐渍化可能性进行预估和初步预测。

(2)结合水盐平衡、概率统计等方法进行半定量预报。水盐平衡、概率统计等方法是区域水盐动态研究由定性向定量研究的过渡研究的重要方法，由于这些方法要求输入的数据比较简单，容易获得，在区域土壤盐渍化预测预报中常被采用。水盐平衡法是以质量守衡为理论依据，区域水盐平衡可以定量地表征区域土壤盐渍化发展的方向，区域水盐平衡研究不足之处是在揭示区域水分和盐分运动的作用机理方面较薄弱，在阐述区域土壤水分和盐分之间的内在关系方面也存在不足，不能很好地对区域水分和盐分的分布状况做出精确预报。由于土壤属性本身的变异性和影响水盐运动的各种因素如降水、蒸发、地下水位和地下水埋深等状况的随机性，区域水盐运动具有明显的随机性，因此概率统计法和成因分析法等方法也常用来研究水盐运动的随机特征，该类方法考虑了水盐运动的随机性特点，具有一定的灵活性，但可移植性差，多用于研究方面，在实际的盐渍化土地改良应用方面，仍需加强。

(3)在区域水盐动态研究的基础上，建立数学模型，借助电子计算机，对土壤盐渍化进行定量预报。“盐随水去，盐随水来”，土壤盐渍化和土壤潜在盐渍化受到复杂因素的影响，地下水位和埋深对于土壤盐分的预测预报固然重要，但地下水位和埋深的实时动态监测的困难在一定程度上对土壤盐渍化和潜在盐渍化的预测预报还未普及，地形地貌对于盐分的预测预报更具有现实意义，海拔高度获取比较容易，对于土壤盐渍化预测预报具有较强的可行性，现有技术中基于海拔高度识别和动态监测随机性特点明显的土壤盐渍化和潜在盐渍化，尚缺乏有效的解决方案。

发明内容

为了解决上述背景技术中存在的技术问题，本公开提供一种基于海拔高度的土壤盐渍化预测方法及***，用逐分量的Metropolis-Hastings算法从后验分布中产生随机数，提高马尔科夫链(Markov Chain)的混合效率，在Metropolis-Hastings算法中按分量进行逐个更新，不需要考虑调节参数，实现土壤盐渍化的随机预测预报。

为了实现上述目的，本公开采用如下技术方案：

本公开的第一个方面提供一种基于海拔高度的土壤盐渍化预测方法，其包括：

获取待预测区域的海拔高度；

基于海拔高度，采用logistic回归模型，对待预测区域的土壤盐渍化二元状况进行预测；

其中，logistic回归模型的构建方法为：基于采集的若干个采样点的海拔高度和土壤盐渍化二元状况，通过逐分量的Metropolis-Hastings抽样方法从后验分布中产生随机数，并结合提议分布按分量进行逐个更新logistic回归参数，通过似然函数更新先验分布。

进一步地，所述从后验分布中产生随机数的具体方法为：

(301)令β＝(β₀，β₁)^T，β₀和β₁为第t-1次迭代得到的logistic回归参数；

(302)从第一提议分布产生候选点β′₀；

(303)令β′＝(β′₀，β₁)^T，计算第一接受概率；

(304)以第一接受概率接受β＝β′；否则β保持不变；

(305)从第二提议分布产生候选点β′₁；

(306)令β′＝(β₀，β′₁)^T，计算第二接受概率；

(307)以第二接受概率接受β＝β′；否则β保持不变；

(308)令t＝t+1，返回(301)，直到t达到设定值，输出β。

进一步地，所述第一接受概率为：

α₀(β，β′₀)＝min{1，A}

其中，

β＝(β₀，β₁)^T，β₀和β₁为logistic回归参数，β′₀为第一提议分布产生候选点，y为所有样本点的土壤盐渍化二元状况向量，π()表示先验分布，f()表示目标分布。

进一步地，所述第二接受概率为：

α₁(β，β′₁)＝min{1，B}

其中，

β＝(β₀，β₁)^T，β₀和β₁为logistic回归参数，β′₁为第二提议分布产生候选点，y为所有样本点的土壤盐渍化二元状况向量，π()表示先验分布，f()表示目标分布。

进一步地，所述提议分布使得产生的马氏链满足不可约、正常返、非周期。

进一步地，所述先验分布为独立的正态分布。

本公开的第二个方面提供一种基于海拔高度的土壤盐渍化预测***，其包括：

模型构建模块，其被配置为：基于采集的若干个采样点的海拔高度和土壤盐渍化二元状况，通过逐分量的Metropolis-Hastings抽样方法从后验分布中产生随机数，并结合提议分布按分量进行逐个更新logistic回归参数，通过似然函数更新先验分布；

数据获取模块，其被配置为：获取待预测区域的海拔高度；

预测模块，其被配置为：基于海拔高度，采用logistic回归模型，对待预测区域的土壤盐渍化二元状况进行预测。

进一步地，所述先验分布为独立的正态分布。

本公开的第三个方面提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现如上述所述的一种基于海拔高度的土壤盐渍化预测方法中的步骤。

本公开的第四个方面提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现如上述所述的一种基于海拔高度的土壤盐渍化预测方法中的步骤。

与现有技术相比，本公开的有益效果是：

本公开提供了一种基于海拔高度的土壤盐渍化预测方法，其基于海拔高度通过马尔科夫链蒙特卡洛(Markov Chain Monte Carlo，MCMC)逐分量方法实现土壤盐渍化的随机预测预报。

本公开提供了一种基于海拔高度的土壤盐渍化预测方法，该方法仅需要相对较少的海拔高度采样就可以得到是否盐渍化的较准确的近似，根据先验分布、似然函数和提议分布而获得的目标后验分布把对后验求积分的过程转化成了对采样链所构成的向量进行求和的过程，土壤盐渍化风险预测后验高密度图用一种简洁、直观的方式提供了土壤盐渍化预测的可能性。

本公开提供了一种基于海拔高度的土壤盐渍化预测方法，其用逐分量的Metropolis-Hastings算法从后验分布中产生随机数，提高了链的混合效率；在Metropolis-Hastings算法中，按分量进行逐个更新，自动执行推断、不需要考虑调节参数。

附图说明

构成本公开的一部分的说明书附图用来提供对本公开的进一步理解，本公开的示意性实施例及其说明用于解释本公开，并不构成对本公开的不当限定。

图1是本公开实施例一的一种基于海拔高度的土壤盐渍化预测方法的流程图；

图2(a)是本公开实施例一的逐分量M-H算法logistic回归参数β₀的样本路径图；

图2(b)是本公开实施例一的逐分量M-H算法logistic回归参数β₁的样本路径图；

图3(a)是本公开实施例一的逐分量M-H算法logistic回归参数β₀的遍历均值图；

图3(b)是本公开实施例一的逐分量M-H算法logistic回归参数β₁的遍历均值图；

图4(a)是本公开实施例一的逐分量M-H算法logistic回归参数β₀的自相关图；

图4(b)是本公开实施例一的逐分量M-H算法logistic回归参数β₁的自相关图；

图5是本公开实施例一的土壤盐渍化预报β₀分量参数的后验密度和99％，95％，50％最高密度区域示意图；

图6是本公开实施例一的土壤盐渍化预报β₁分量参数的后验密度和99％，95％，50％最高密度区域示意图。

具体实施方式

下面结合附图与实施例对本公开作进一步说明。

应该指出，以下详细说明都是例示性的，旨在对本公开提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本公开所属技术领域的普通技术人员通常理解的相同含义。

术语解释：

地理相似法：通过与预报地区环境条件相似的己盐渍化地区的情况相比较进行预报。采用这种方法，必需对预报地区与已产生盐渍化地区的条件，进行全面的调查分析和详细的比较研究，从而进行预报。需要考虑地理区域的主要环境因素、自然条件和人类活动条件。自然条件特征包括：(1)气象、水文条件，(2)地形、地貌条件，(3)地质及水文地质状况，(4)土壤状况(包括土壤类型、分布，土壤物理性质、土壤物理化学和生物特性及其空间上的分布)，(5)地球化学特征。人类活动条件主要指：(1)土地利用情况，(2)农业技术措施，(3)土壤改良措施，(4)水利技术措施(包括灌溉、排水、蓄水等措施)。

水盐平衡方法：通过对预报地区水盐平衡计算，预报土壤盐渍化发生的可能性，一般常用的盐分平衡计算公式是：

Δs＝[P+I+R+G+W+F]-[L_p+L_i+r+g+u]

式中，Δs为某一地区盐分总贮量的变化；P为由大气带来的盐量(大气降水和风力搬运等)；I为由灌溉水带入的盐量；R为由地表水带入的盐量(径流、洪水、沥涝)；G为地表以下输入的盐量(地下水、深层土壤水)；W为由当地风化过程增加的盐量；F为施入肥料和化学改良剂带来的盐量；L_p为大气降水所淋失的盐量；L_i为灌溉或冲洗水所淋失的盐量；r为由地表径流带走的盐量；g为由土壤水水平出流带走的盐量；u为随农业收获物带走的盐量。

基于先验和后验信息的贝叶斯推理方法：是基于贝叶斯定理发展起来的用于***地阐述和解决随机问题的概率方法，贝叶斯统计最困难的是在寻求后验分布过程中所遇到的积分困难。与后验分布有关的积分是很难用数值方法去计算的，尤其在高维的情形，MCMC方法提供了一个有效的途径去解决这类问题，即蒙特卡罗(Monte Carlo)方法用平均值近似积分；马尔科夫链(Markov Chain)解决样本的抽样问题。

Metropolis-Hastings算法(M-H算法)，是最常用的马尔科夫链蒙特卡洛方法之一，该方法从一般的后验分布中进行抽样，抽样策略是建立一个不可约的、非周期的马氏链，并且其平稳分布即为感兴趣的目标后验分布，核心问题是确定从当前值转移到下一个值的规则，其主要任务是生成满足上述正则条件的一个马氏链{x；i＝1,...,J}。

实施例一

基于水盐运动机理的数学模拟研究主要根据流体动力学原理，对包气带和非包气带水盐运动进行数值模拟，并借助于计算机技术对方程进行求解。该方法的优点在于它可作出研究区任何一点、任何深度、任何时段的地下水、土壤盐分和水分分布预报。

作为制约土壤表层盐分含量的关键因素，地下水位的降低使表层土壤脱盐成为可能，而土壤表层盐分含量降低使盐碱地成为可垦地，开垦以后的土地投入增加和不同用地类型投入的差异促使养分平衡向养分积累方向发展，进而土壤养分大幅度上升和不同土地类型养分含量产生差异。虽然土壤肥力质量明显提升，土壤质量有了较大提高，但目前导致土壤盐渍化的生物气候条件并没有改变，地下水位的变化仍然制约着土壤表层的盐分含量，土壤潜在盐渍化仍然存在，盐渍化、旱涝灾害的反复性和农业生产的不稳定性仍然存在，关注土壤潜在盐渍化等障碍因素的动态变化，对土壤盐渍化实现预测预报，不仅有利于认识肥力特征、土壤障碍及健康状况，而且有利于重塑中低产田土壤生态功能，提升中低产田土壤质量和可持续生产能力，满足对中低产田土壤资源开发利用，对农业的可持续发展具有重要意义。

基于此，结合地下水位动态调控和水盐运动的机理，制定切实可行的防治土壤盐渍化的措施是有必要的。地下水埋深和地下水矿化度、土壤的电导率(利用电导率和土壤盐分含量的回归关系，可以得出采样点的土壤中的盐分含量)是土壤水盐运动定量模拟的重要指标，这些指标的确定对于土壤水盐运动的时间和空间的动态模拟和预测预报有重要意义。事实上，研究区域土壤盐分的空间格局与该区的地形地貌特征、土壤类型和河流的分布状况有着密切的联系，一般地，在地势较高的区域，地下水埋深相对较深，盐渍化威胁的概率相对较小。由于土壤属性本身的变异性和影响水盐运动的各种因素如耕作、管理措施、种植制度、降水、蒸发、地下水状况等的随机性，使得区域水盐运动具有随机性和复杂性的特点，因此概率统计法被常用来研究水盐运动和土壤盐渍化的预测预报。

近些年来，虽然概率编程/贝叶斯编程计算使贝叶斯统计有了强劲的发展势头，未见运用贝叶斯推理基于海拔高度预测土壤盐渍化的报道。

本实施例提供了一种基于海拔高度的土壤盐渍化预测方法，如图1所示，包括以下步骤：

步骤1、获取第i个待预测区域的海拔高度x_i；

步骤2、基于第i个待预测区域的海拔高度x_i，采用logistic回归模型，预测得到第i个待预测区域的土壤盐渍化二元状况y_i。

logistic回归模型的构建过程为：

(1)采集若干个采样点的海拔高度x_i和土壤盐渍化二元状况y_i，i＝1，....，n，n表示采样点的个数，可以取值为100。

(2)选取如下的logistic回归模型：

y_i～B(1，π_i)，

其中，B表示二项式分布，π_i表示二项式分布的参数，y_i表示第i个采样点的土壤盐渍化二元状况(取值为1表示具有土壤盐渍化，取值为0表示不具有盐渍化)，x_i表示第i个采样点的海拔高度，β₀和β₁为回归参数。

logistic回归模型的似然函数为：

其中，

表示后验分布。

考虑β₀和β₁的先验分布π为独立的正态分布，即π(β₀，β₁)＝π₁(β₀)π₂(β₁)，而

j＝0，1；当

且

很大时，先验分布接近无信息先验，其中，N为正态分布，

为均值，

为方差，表征正态分布的参数。因此后验分布为：

其中，σ₀和σ₁为后验分布的参数。

(3)基于采集的若干个采样点的海拔高度x_i和土壤盐渍化二元状况y_i，通过逐分量的Metropolis-Hastings抽样方法从后验分布π(β₀，β₁|y)中产生随机数，并结合提议分布按分量进行逐个更新logistic回归参数，选取似然函数f(y|β₀，β₁)更新先验分布。

提议分布为：

此处β′＝(β′₀，β′₁)^T，而β＝(β₀，β₁)^T。

其中，

为提议分布的参数。

M-H算法中，按分量进行逐个更新，其优势在于应用方便，不需要考虑调节参数。

通过逐分量的Metropolis-Hastings抽样方法从后验分布π(β₀，β₁|y)中产生随机数的具体算法为：对t＝1，...，T：

(301)令

即令第t-1次迭代得到的β中的两个分量(logistic回归参数)β₀和β₁作为第t-1次迭代的β中的两个分量β₀和β₁；

(302)从第一提议分布

产生候选点β′₀；

(303)令β′＝(β′₀，β₁)^T，计算第一接受概率α₀(β，β′₀)＝min{1，A}；

其中，

y为所有样本点的土壤盐渍化二元状况向量，π()表示先验分布，f()表示目标分布；

(304)以第一接受概率α₀(β，β′₀)接受β＝β′，即将β的第一个分量β₀更新为β′₀；否则β保持不变；

(305)从第二提议分布

产生候选点β′₁；

(306)令β′＝(β₀，β′₁)^T，第二接受概率α₁(β，β′₁)＝min{1，B}；

其中，

(307)以第二接受概率α₁(β，β′₁)接受β＝β′，即将β的第二个分量β₁更新为β′₁；否则β保持不变；

(308)令β^(t)＝β，并t＝t+1，返回(301)，直到t达到设定值T，输出β。

用逐分量的M-H算法从后验分布中产生随机数，提高了链的混合效率。在M-H算法中，按分量进行逐个更新，其优势在于应用方便，不需要考虑调节参数。

本实施例中，用f()表示目标分布，以g()表示提议分布。设，希望从目标分布f()中抽样，M-H算法从初始值X₀出发，指定一个从当前值X_t转移到下一个值X_t+1的规则，从而产生马氏链{X_t＝0，1，...}。

具体来说，在给定当前值X_t，从一个提议分布g(·|x)产生一个候选点x′，若候选点x′被接受，则链从状态x′转移到链的下一时刻t+1，令X_t+1＝x′；否则，链停留在状态X_t，令X_t+1＝X_t，候选点x′是否被接受为链的下一个值，根据接受概率的表达式α(X_t，x′)＝min{1，A}来确定，其中：

需要注意的一点是上式分子和分母中的密度函数可分别用“密度函数的核”代替，故密度函数中的正则化常数因子可省略掉，以便简化计算。

M-H算法生成满足正则条件的马氏链迭代过程如下：

(1)选择一个提议分布g(·|X_t)；

(2)从提议分布中生成初始值x₀；

(3)对t＝1，2，...重复步骤(a)-(d)：

(a)从提议分布g(·|X_t)中产生一个候选值x′；

(b)从均匀分布U(0，1)中生成随机数U；

(c)计算接受概率，若U≤A，则接受x′且令X_t+1＝x′，否则令X_t+1＝X_t；

(d)增加t，返回到(a)。

上述步骤中提到的接受概率并非越大越好，因为这可能导致较慢的收敛性。当参数维数是1时，接受概率应略小于0.5是最优的，当维数大于5时，接受概率应降至0.25左右。为说明M-H抽样方法产生的马氏链具有平稳分布f，则可以通过说明此链的转移核(或转移概率)和f均满足平衡方程。此外，提议分布g的选择除了使得产生的马氏链满足不可约、正常返、非周期且具有平稳分布f等正则化条件外，还应满足：

(1)提议分布的支撑集包含目标分布的支撑集；

(2)容易从中抽样，常取为已知的分布，如正态或t分布等；

(3)提议分布应使接受概率容易计算；

(4)提议分布的尾部要比目标分布的尾部厚；

(5)新的候选点被拒绝的频率不高。

并且，当目标分布为后验分布时，前面的M-H算法可直接在贝叶斯框架下实现，只要将x用感兴趣的参数θ代替，将目标分布f(x)用后验分布π(θ|x)代替方法概括如下：

(1)选择一个提议分布g(·|θ_t)；

(2)从提议分布中生成初始值θ₀；

(3)对t＝1，2...重复步骤(a)-(d)：

(a)从提议分布g(·|θ_t)中产生一个候选值θ’；

(b)从均匀分布U(0，1)中生成随机数U；

(c)计算接受概率，若

则接受θ’且令θ_t+1＝θ′，否则令θ_t+1＝θ_t；

(d)增加t，返回到(a)。

在Metropolis抽样方法中，提议分布是对称的，即g(·|X_n)满足g(X|Y)＝g(Y|X)，因此接受概率为：

根据提议分布g的不同选择，逐分量的Metropolis-Hastings抽样方法是Metropolis-Hastings抽样方法一个变种之一，具有广泛的用途。

当状态空间为多维时，不整体更新Xn，而是对其分量进行逐个更新，即称为逐分量的M-H抽样方法(componentwise Metropolis-Hastings sampler)，这样做更方便和更有效率。

记为：Xn＝(X_n，1，...，X_n，k)，Xn，-i＝(X_n，1，....，X_n，i-1，....，X_n，k)，分别表示在第n步链的状态，及在第n步除第i个分量外其他分量的状态，f(x)＝f(x₁，...，x_k)为目标分布，f(x_i|x_-i)＝f(x)/∫f(x₁，...，x_k)dx_i表示X_i对其他分量的条件密度。

逐分量的M-H抽样方法由k步构成：令X_n，i表示在第n次迭代后Xn第i个分量的状态，则在第n+1步迭代的第i步中，使用M-H算法更新Xn，i做法如下：

首先，对i＝1，....，k，从第i个提议分布q_i(·|X_n，i，X^* _n，-i)中产生Y_i，这里X^* _n，-i＝(X_n+1，1，.....，X_n+1，i-1，...，X_n，k)；

然后以概率

若Y_i被接受，则令X_n+1，i＝Y_i；否则令X_n+1，i＝X_n，i。

图2(a)和图2(b)两个关于β₀和β₁的时间轨迹图显示链的混合较好。画出的迭代次数和样本值的曲线图，轨迹图中的数值出现在一定区域内，并没有呈现趋势性和周期性变化，样本路径图稳定下来，并且无法区别彼此(混合在一起)，可以认为达到收敛状态，链的混合较好。其中，Value of beta0表示参数β₀的值；Value of beta1表示参数β₁的值；Iterations表示迭代次数。

收敛特征的诊断对于MCMC模拟是非常重要的，样本路径图和遍历均值图都可以诊断MCMC的收敛性。样本路径图(trace plot)是将生成的MCMC按迭代次数采样得到。马氏链的样本路径图为了避免链写入目标分布的某个局部区域，通常从不同的初始点出发同时产生多个马氏链，经过一段时间运行后，如果他们的样本路径图都稳定下来了，具有很好的混合度(mixing)，看不到任何可识别的模式，包括向上或向下的曲线，曲线在某个值附近震荡，则认为抽样已经收敛。将所生成的马氏链的累计均值对迭代次数作图得到此链的遍历均值图(ergodic mean plot)，用于判定遍历均值是否达到收敛。达到平稳后的马氏链的遍历均值会趋于一条水平的直线。同样，为了避免链写入目标分布的某个局部区域，也可以考察从不同初始点出发的多条马氏链的遍历均值是否收敛。如果仅使用一条马氏链，则要求迭代次数足够多，使链能够到达支撑的每一部分。中间两个呈现的是β₀和β₁的遍历均值图，去掉前1000次迭代(即预烧期)后趋于平稳，也表明链的收敛性较好。

另外，自相关系数也是检验马氏链是否收敛的比较有效的方法，较小的相关系数表示较快的收敛性。土壤盐渍化预测链的自相关函数图中显示马氏链的抽样步长L＝25次迭代，也就是每间隔25个样本抽取一个，获得的样本自相关性较低，显示了较好的收敛性。

图2(a)和图2(b)、图3(a)和图3(b)、图4(a)和图4(b)三组图诊断结果都显示链的收敛性较好。去掉预烧期后由链获得的β的两个分量β₀和β₁的链的样本均值分别是1.8715，-0.0611，样本标准差分别是0.7862，0.0418。其中，Lag表示滞后，Thin-25 iterations表示抽样步长25次迭代。收敛性也是指算法生成的样本是否达到均衡分布，即样本是否由目标分布生成。土壤盐渍化预测链的自相关函数图中显示马氏链的抽样步长L＝25次迭代，也就是每间隔25个样本抽取一个，获得的样本自相关性较低，显示了较好的收敛性。

参数β₀和β₁的后验分布包含了全部信息，也是较好的、完全的估计，参考图5和图6土壤盐渍化预报β₀和β₁参数的后验密度和99％，95％，50％最高密度区域。其中，Densitydefault表示缺省状况下的密度分布，N＝50000表示迭代次数为50000次，Bandwidth＝0.1417表示带宽为0.1417，Bandwidth＝0.006143表示带宽为0.006143。对于参数后验分布式单峰的情况下，可能会有很多组满足下面条件的区间(H_α/2L_α/2)：

满足以上条件的最小区间(H_α/2L_α/2)称为1-α的最高密度区域(highest densityregion)或者bayesian incredible interval，而对于多峰后验分布，最高密度区间可有多于一个连续区间组成。

参数的收敛状况和参数后验分布的高密度区域均表明基于海拔高度可以实现对于土壤盐渍化的预测。

实施例二

本实施例提供了一种基于海拔高度的土壤盐渍化预测***，其具体包括如下模块：

模型构建模块，其被配置为：基于采集的若干个采样点的海拔高度和土壤盐渍化二元状况，通过逐分量的Metropolis-Hastings抽样方法从后验分布中产生随机数，并结合提议分布按分量进行逐个更新logistic回归参数，通过似然函数更新先验分布；

数据获取模块，其被配置为：获取待预测区域的海拔高度；

预测模块，其被配置为：基于海拔高度，采用logistic回归模型，对待预测区域的土壤盐渍化二元状况进行预测。

其中，先验分布为独立的正态分布。

此处需要说明的是，本实施例中的各个模块与实施例一中的各个步骤一一对应，其具体实施过程相同，此处不再累述。

实施例三

本实施例提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现如上述实施例一所述的一种基于海拔高度的土壤盐渍化预测方法中的步骤。

实施例四

本实施例提供了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现如上述实施例一所述的一种基于海拔高度的土壤盐渍化预测方法中的步骤。

本领域内的技术人员应明白，本公开的实施例可提供为方法、***、或计算机程序产品。因此，本公开可采用硬件实施例、软件实施例、或结合软件和硬件方面的实施例的形式。而且，本公开可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器和光学存储器等)上实施的计算机程序产品的形式。

本公开是参照根据本公开实施例的方法、设备(***)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(RandomAccessMemory，RAM)等。

以上所述仅为本公开的优选实施例而已，并不用于限制本公开，对于本领域的技术人员来说，本公开可以有各种更改和变化。凡在本公开的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。

Claims (10)

1.一种基于海拔高度的土壤盐渍化预测方法，其特征在于，包括：

获取待预测区域的海拔高度；

基于海拔高度，采用logistic回归模型，对待预测区域的土壤盐渍化二元状况进行预测；

其中，logistic回归模型的构建方法为：基于采集的若干个采样点的海拔高度和土壤盐渍化二元状况，通过逐分量的Metropolis-Hastings抽样方法从后验分布中产生随机数，并结合提议分布按分量进行逐个更新logistic回归参数，通过似然函数更新先验分布。

2.如权利要求1所述的一种基于海拔高度的土壤盐渍化预测方法，其特征在于，所述从后验分布中产生随机数的具体方法为：

(301)令β＝(β₀,β₁)^T，β₀和β₁为第t-1次迭代得到的logistic回归参数；

(302)从第一提议分布产生候选点β′₀；

(303)令β′＝(β′₀,β₁)^T，计算第一接受概率；

(304)以第一接受概率接受β＝β′；否则β保持不变；

(305)从第二提议分布产生候选点β′₁；

(306)令β′＝(β₀,β′₁)^T，计算第二接受概率；

(307)以第二接受概率接受β＝β′；否则β保持不变；

(308)令t＝t+1，返回(301)，直到t达到设定值，输出β。

3.如权利要求2所述的一种基于海拔高度的土壤盐渍化预测方法，其特征在于，所述第一接受概率为：

α₀(β,β′₀)＝min{1,A}

其中，

β＝(β₀,β₁)^T，β₀和β₁为logistic回归参数，β′₀为第一提议分布产生候选点，y为所有样本点的土壤盐渍化二元状况向量，π()表示先验分布，f()表示目标分布。

4.如权利要求2所述的一种基于海拔高度的土壤盐渍化预测方法，其特征在于，所述第二接受概率为：

α₁(β,β′₁)＝min{1,B}

其中，

β＝(β₀,β₁)^T，β₀和β₁为logistic回归参数，β′₁为第二提议分布产生候选点，y为所有样本点的土壤盐渍化二元状况向量，π()表示先验分布，f()表示目标分布。

5.如权利要求2所述的一种基于海拔高度的土壤盐渍化预测方法，其特征在于，所述提议分布使得产生的马氏链满足不可约、正常返、非周期。

6.如权利要求1所述的一种基于海拔高度的土壤盐渍化预测方法，其特征在于，所述先验分布为独立的正态分布。

7.一种基于海拔高度的土壤盐渍化预测***，其特征在于，包括：

模型构建模块，其被配置为：基于采集的若干个采样点的海拔高度和土壤盐渍化二元状况，通过逐分量的Metropolis-Hastings抽样方法从后验分布中产生随机数，并结合提议分布按分量进行逐个更新logistic回归参数，通过似然函数更新先验分布；

数据获取模块，其被配置为：获取待预测区域的海拔高度；

预测模块，其被配置为：基于海拔高度，采用logistic回归模型，对待预测区域的土壤盐渍化二元状况进行预测。

8.如权利要求7所述的一种基于海拔高度的土壤盐渍化预测***，其特征在于，所述先验分布为独立的正态分布。

9.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现如权利要求1-6中任一项所述的一种基于海拔高度的土壤盐渍化预测方法中的步骤。

10.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如权利要求1-6中任一项所述的一种基于海拔高度的土壤盐渍化预测方法中的步骤。

Images