本发明涉及土壤学研究领域，具体涉及一种土壤属性插值的估算方法。
背景技术：
土壤连续属性(土壤养分含量、重金属含量)的分布特征和定量分布信息是进行土壤生产力评价及环境综合评估的基础。近年来，随着精准农业的战略的实施，对土壤属性信息的精度提出了更高的要求。因此如何提高土壤属性的预测精度，为农业生产、生态评价提供精确的土壤属性信息，已成为土壤学研究的热点问题。
以克里金方法为代表的经典地统计学方法在许多学科领域都得到了广泛的应用，但就其方法本身来说还存在一定的缺陷，制约了预测结果的精确度，如：
(1)经典的克里金方法仅考虑变量自身的采样数据，而忽略了与变量相关的其他信息(如环境因素、专家经验、历史资料等)的利用，造成了预测结果的降低，在采样点稀疏的情况下还会造成预测结果出现非常不合理的情况，产生了较大的不确定性；
(2)克里金方法具有一定的平滑效应，针对短距离内变化剧烈的区域，克里金的平滑效应会将降低或完全消除这些变化在空间上的反应，导致信息丢失。
(3)克里金方法要求区域变量满足高斯分布，而在实践中通过数据变换单点的高斯分布容易满足，多点的高斯分布却难以满足。另外数据中往往存在一些离群点与异常点，这些特殊的采样点的存在也将极大的影响模型的精度。
(4)克里金方法为局部估计方法，对估计值的整体空间相关性考虑不够，它仅保证了数据的估计局部最优而不能保证数据估计的全体最优。这就要求进行克里金方法预测的采样数据需要在空间上分布均匀且达到一定的采样密度。
新兴的bme方法虽然在空间属性预测上取得了较好的成果，但是其发展时间较短，缺少相应的可视化软件的支持，如主流的空间分析及制图软件arcgis和qgis都不支持该方法。bme法的另一大缺点是在软数据的获取上，要搜集组织出与待测变量相关软数据需要对研究区域和待测变量的特性有较深的了解，因而在同一区域的空间变量预测可能由研究者的对该区域的了解不同造成最终结果的不一致。总而言之，bme方法虽然可以综合利用各方面的数据，但是对这些数据的来源、类型、组织方式等都没有具体的要求，该方法的具体实施还需要对具体研究区域和预测变量进行分析，可以说其实施过程比经典的克里金法要复杂的多，从而限制了bme方法在土壤、环境等领域的利用。
技术实现要素：
本发明的目的在于提供一种土壤属性插值的估算方法，解决目前的土壤属性插值的估算采用单纯的克里金法，估算精度较低的问题。
为解决上述的技术问题，本发明采用以下技术方案：
一种土壤属性插值的估算方法，包括以下步骤：
s1：利用克里金法生成土壤属性待测变量的先验概率分布；
s2：利用环境因素与预测点的相关关系生成软数据；
s3：利用bme法更新土壤属性待测变量的先验概率分布；
s4：利用土壤属性待测变量的先验概率分布进行土壤空间属性制图。
更进一步的方案是，上述s1步骤中利用克里金法生产土壤属性待测变量的先验概率分布的具体方法是：
s101：按照s4步骤中的制图精度要求将土壤研究区域划分为row×col的网格，设土壤研究区域有n个采样点，则第k(k<n)个采样点的属性值可表示为χij，其中i和j为采样点在row×col的网格中的位置，i<row,j<col；对于待测点pij，使用实验变差函数确定待测点周围各点的权重值，计算出每个待测点的属性估值
s102：利用每个待测点的属性估值得到每个待测点的先验概率分布函数。
更进一步的方案是，上述s101步骤中计算每个待测点
更进一步的方案是，上述s101步骤中通过以下线性方程组求得观测点属性值x(pn)的权重λn：
其中p0表示待观测点的位置；pi和pj为待测点周围各观测点的空间位置，i、j取值为1～n；
更进一步的方案是，上述s101步骤中设定实验变差函数的取值只与空间中的相对位置有关，因此使用以下公式计算实验变差函数：
其中，n(h)为空间中距离为h的观测变量的对数；pi为观测点在空间中的位置；x(pi)为待测点pi周围各观测点的属性值；
计算出实验变差函数的值后使用变差函数理论模型对实验变差函数值进行拟合，得到一条最优的实验变差函数。
更进一步的方案是，上述s102步骤中利用每个待测点的属性估值得到每个待测点的先验概率分布函数的公式是：
其中，g表示由克里金法得到的概率分布；μ为克里金的预测结果；σ2为预测方差。
更进一步的方案是，上述s2步骤中利用环境因素与预测点的相关关系生成软数据的具体方法是：
s201：将实验区域划分为n个属性值的取值区间，计算第k个属性取值区间对该点上的单个环境因素的隶属度mk,si,j，计算公式为
s202：计算环境因子对待测点属性值的决定系数rj2，
s203：在研究区域随机生成分布均匀的一定数量的软数据点，对于某一个软数据点psi，由s201步骤中确定的隶属度函数确定某个属性值取值区间对该软数据点psi上的单个环境因子的隶属度mk,si,j，并结合环境因子对待测点属性值的决定系数，确定该软数据点psi的取值范围对某个属性值取值区间的相似度
更进一步的方案是，上述s3步骤中，利用bme法更新土壤属性待测变量的先验概率分布采用下述公式：
其中，fg为先验概率密度函数；fs为软数据概率密度函数。
更进一步的方案是，上述s4步骤中利用土壤属性待测变量的先验概率分布进行土壤空间属性制图的具体方法是：
根据s3步骤中获得的土壤属性待测变量的先验概率分布函数，确定各待测点的预测属性值，预测属性值由下述公式确定：
其中，
按照上述公式计算出各待测点的预测属性值后，得到属性分布图。
更进一步的方案是，上述s2步骤中，对于描述型环境因子数据，对每一个描述型环境因子赋予唯一的数字编号，根据编号的数量将采样点的属性值划分为与编号数量相等的取值区间，计算每个取值区间对各描述型环境因子的隶属度。
与现有技术相比，本发明的有益效果是：
1、本发明将操作简单、成熟的克里金法与精确、合理的贝叶斯最大熵法进行结合，使得插值过程更加简便，插值结果更加合理。
2、本发明的方法能够有效的利用各方面的数据，从而提高插值精度。
3、本发明的插值过程中考虑到了环境因素对土壤属性值的影响，将其作为软数据参与到插值的过程中来，解决了克里金法在采样点稀疏或分布不均时插值结果不稳定、不合理的问题，降低了对数据的需求。
4、本发明以环境因子与属性值的关系定量的产生软数据，降低了软数据搜集、组织的难度。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
具体实施例：
1、使用克里金法生成待测变量的先验概率分布
按照制图的精度要求将研究区域划分为row×col(这里以二维平面插值为例)的格网，设该区域有n个采样点，则第k(k<n)个采样点的属性值可表示为χij，其中i<row,j<col；对于待测点pij，其属性的估计值
式中：x(pn)为待测点pij周围各观测点的属性值；pn为待测点周围各观测点的空间位置；n为待测点周围观测点的数量；λn为观测点属性值x(pn)的权重。
上式中确定待测点估计值的关键在于确定λn，可以通过解以下线性方程组求得：
式中：p0表示待测点的位置；
式中：n(h)为空间中距离为h的观测变量的对数。pi为观测点在空间中的位置。
计算出实验变差函数的值以后，使用理论模型对实验变差函数值进行拟合，得到一条最优的实验变差函数。
使用变差函数确定待测点周围各点的权重值，从而确定待测点的估计值。使用这种方法得到格网中每个点的属性值及其先验概率密度函数：
式中：g表示由克里金法得到的概率分布；μ为克里金的预测结果；σ2为预测方差。
2、利用相关数据生成软数据
可以利用以下3种方法生成软数据：
(1)利用环境因素与预测点的相关关系
假设预测变量受环境因素的影响，那么在空间点pij的属性分布就可以在一定程度上用该点的环境因
子eij(eij,1,…,eij,m)来表示，按照这一推断我们就可以在区域中利用某一点的环境因子来近似的表达待测点区域化变量的分布。
首先根据制图精度的要求，将采样点的属性值划分为n个区间，对于第k个区间，其属性值含量为：
式中：χmin，χmax为采样点属性值的最小取值和最大取值。
划分出属性值的取值区间后，就可以定量的分析属性取值区间与环境因子的关系。方法是将区域中与待测变量取值有关的环境因素按照同样的方法划分为n个区间，那么对于属性值的第k个区间，它对某一环境因子(例如第j个环境因子)取值区间的隶属度可用以下数据序列表示：
式中：ck为属性值属于第k个区间的采样点数；c1为属性值属于第k个区间且该点的第j个环境因子属于第1个环境因子区间的采样点数。将上式的计算结果以隶属度为纵坐标，环境因子的区间序号为横坐标用折线表示出来，采用一条平滑的曲线拟合该折线，从而得到采样点属性值范围对第j个环境因子范围的隶属度函数。使用该方法最终得到n×m个隶属度函数。
另外，考虑到每个环境因子对待测点属性值的影响程度不同，这里还需要计算这m个环境因子对待测点属性值的决定系数：
首先由采样点属性值与对应的环境因子的取值建立线性函数关系，按照上式计算环境因子j的决定系数。式中：
按照上述方法计算出m个决定系数：
(r12,…,rm2)
随后对决定系数进行归一化：
这样以来，便可以在研究区域随机生成分布均匀的一定数量的软数据点，对于某一个软数据点psi，可由之前确定的隶属度函数确定某个属性值取值区间对该点上单个环境因子的隶属度mk,si,j，并考虑到环境因子对属性值的影响程度，使用下式确定该点取值范围对属性值某一取值区间的相似度：
对上式计算的sk,si进行归一化处理，即
fs(psi)＝((s1,si,n1),…,(sk,si,nk),…,(sn,si,nn))
式中:fs(psi)表示由环境数据确定属性值概率密度分布；nk表示属性值第k个取值区间。
上述关于环境因子的计算主要是数值型环境因子，对于描述型(范畴型)环境因子，需要对该方法进行一些改进，如对土壤碱解氮进行空间统计插值时，根据以往的经验我们可以知道碱解氮的含量与土类或是区域的施肥措施存在一定的联系，而土壤类型和施肥措施属于描述型的环境因子数据，无法直接套用第三项的计算方法，这里的做法是对于每一个土类或者施肥措施，赋予其唯一的数字编号，根据编号的数量将采样点的属性值划分为与编号数量相等的取值区间，这样就可以计算每个取值区间对土类编号的隶属度，从而进行后续的计算。在这一步需要注意的是，由于环境因子中描述性数据的存在，使得我们划分属性取值区间时就不能任意选择区间的数量，土壤属性的取值区间与环境因子的取值区间应统一为描述型环境因子数据的最小编号数量，这样就可以将数值型与描述型的环境因子以相同的方式进行计算，从而拓展了环境因子的选取范围。
(2)使用区域历史图件生成软数据
利用研究区域的历史图件也可生成软数据，例如进行土壤某种属性插值时，往往可以获取到土壤属性的历史分布图件，则可以据此生成软数据点psi，其取值范围为该点所在历史土壤属性分布图的分级范围。以此方法生成区间软数据集。
(3)使用其他与预测变量相关的数据生成软数据
如缺少上述两种方法所需要的数据，也可以使用其他与预测属性值相关的数据生成软数据，如对土壤有机质进行空间统计插值时，由以往的经验可以得知土壤的有机质含量与土地利用\覆被存在一定的联系，因而可以利用土地利用\覆被图或直接利用遥感影像，由土壤有机质含量与土地利用\覆被的关系生成均匀的区间软数据点。
3、使用bme法更新待测变量的先验概率分布
根据贝叶斯条件概率的原理，我们可以在考虑软硬数据的情况下更新由克里金法得到的先验概率密度函数。如预测空间变量χ0在p0位置处考虑软硬数据的完全概率分布为：
式中：fg为先验概率密度函数，由克里金法得到；fs为软数据概率密度函数，由第二步确定。
该公式给出了在空间任意位置pij属性值χij的完全概率分布，基于该公式可以对其进行预测。
4、空间属性制图
有了上一步得到的预测变量完全概率分布，就可以在特定的位置上对属性值进行预测，需要制图时，由于第1步已经将研究区域划分为row×col的格网，那么我们可以将每个格网的节点作为一个像素点，像素值即为预测的属性值。预测的属性值可由下式确定：
式中：
尽管这里参照本发明的解释性实施例对本发明进行了描述，但是，应该理解，本领域技术人员可以设计出很多其他的修改和实施方式，这些修改和实施方式将落在本申请公开的原则范围和精神之内。更具体地说，在本申请公开和权利要求的范围内，可以对主题组合布局的组成部件和/或布局进行多种变型和改进。除了对组成部件和/或布局进行的变形和改进外，对于本领域技术人员来说，其他的用途也将是明显的。