利用卡尔曼滤波精确校正土壤侵蚀理论值的方法及系统

专利名称：利用卡尔曼滤波精确校正土壤侵蚀理论值的方法及系统

专利类型：发明专利

专利申请号：CN202111482323.5

专利申请（专利权）人：华中师范大学
权利人地址：湖北省武汉市洪山区珞瑜路152号

专利发明（设计）人：李畅,颜真珍

专利摘要：本发明公开了一种利用卡尔曼滤波精确校正土壤侵蚀理论值的方法及系统，收集得到月降雨量等数据，通过对遥感影像进行解译、对DEM数据进行提取获得各土壤侵蚀因子具体值，利用中国土壤侵蚀方程(CSLE)计算得到土壤侵蚀理论值，并收集获取野外站点实际土壤侵蚀量，采用贝叶斯抽样一致性算法抽取样本统计流域总体理论侵蚀值，在假设历史时刻土壤侵蚀校正值已知的基础上，建立起历史时刻校正值和理论值之间的关系，通过时间迭代更新，得到高准确度的增益值，从而推求出当前时刻的校正值，并最终能根据当前模型预测未来时刻土壤侵蚀理论值。本发明提高了CSLE计算的精度并拓宽了水土流失方程的应用范围，提高了实用性。

主权利要求：
1.一种利用卡尔曼滤波精确校正土壤侵蚀理论值的方法，其特征在于，包括以下步骤：步骤1：获取研究区水土流失方程各因子的高时空分辨率数据，计算研究区土壤侵蚀理论值；
步骤2：利用改进无迹卡尔曼滤波校正模型建立实际土壤侵蚀实测值和土壤侵蚀理论值相关性方程，并通过时间迭代进行更新拟合；
具体包括以下子步骤：
步骤2.1：将流域所有观测点土壤侵蚀理论值作为一个集合，计算其先验概率，抽取其中概率高于预设阈值的观测点理论值构建最小二乘法统计模型，采用贝叶斯抽样一致性算法进行抽样剔除粗差，抽取可靠的代表性样点，计算流域总体土壤侵蚀理论值；
步骤2.2：构建卡尔曼滤波时间更新模型，将前一时刻的土壤侵蚀理论值加权平均值作为输入量，当前时刻土壤侵蚀理论值的加权平均值作为输出量，进行时间递推；计算土壤侵蚀理论值的时间不确定性，使其达到阈值，满足精度要求；
步骤2.3：计算增益值K，将其作为参数，输入前一过程中更新后的土壤侵蚀理论值和实测值，输出理论土壤侵蚀理论值，并计算理论值和实测值之间的状态不确定性，尽可能提高状态更新精度；
步骤3：利用时间更新中经过调整参数最后获得的卡尔曼滤波时间更新方程，将当前时刻的理论土壤侵蚀量代入，获得未来时刻的土壤侵蚀理论值。
2.根据权利要求1所述的利用卡尔曼滤波精确校正土壤侵蚀理论值的方法，其特征在于：步骤1中，所述水土流失方程各因子值，包括与A＝R×K×L×S×B×E×T各因子的相关数据；对研究区高时空分辨率数据进行解译，获得植被覆盖度数据、土地类型数据和相关工程措施数据E；通过对研究区DEM影像的提取分析，得到坡度数据S、坡长数据L；
‑2 ‑1 ‑2 ‑1 ‑1
其中：A为年均土壤侵蚀量，单位thm a ；R为降雨侵蚀力，单位MJmmhm h a ；K为土
2 ‑2 ‑1 ‑1
壤可蚀性因子，单位thmhhm MJ mm ；L坡长因子、S坡度因子；B为生物措施因子；E为工程措施因子；T为耕作措施因子；L、S、B、E、T是无量纲。
3.根据权利要求2所述的利用卡尔曼滤波精确校正土壤侵蚀理论值的方法，其特征在于：针对中国水土流失方程CSLE，首先，x为关于理论土壤侵蚀量和时间的状态向量，对土壤侵蚀理论值各因子分别进行校正；
x1＝R
x2＝K
x3＝L
x4＝S
x5＝B
x6＝E
x7＝T；
引入对应的实测土壤侵蚀量，分别对7个因子进行校正，得到精确的土壤侵蚀理论值；
然后对土壤侵蚀理论值整体进行校正：
根据CSLE基本方程A＝R×K×L×S×B×E×T得到：
CSLE的状态转移矩阵：
4.根据权利要求2所述的利用卡尔曼滤波精确校正土壤侵蚀理论值的方法，其特征在于：步骤2.1中，从利用n个观测点的相关数据计算出的流域内各观测点某一时刻土壤侵蚀理论值中选取能够代表该时刻流域理论土壤侵蚀值的s个观测点理论值，并计算各点的权重，构建最小二乘法统计模型来拟合流域整体理论土壤侵蚀情况；
所述统计模型依据最小二乘法构建：
式中，D为流域总体土壤侵蚀理论值；ci为s个点中某点土壤侵蚀理论值，i∈1～s；wi为这s个观测点理论值拟合流域平均土壤侵蚀量的均值；为这s个观测点理论值的平均值；
由于同一时刻各观测点理论土壤侵蚀值基本符合高斯分布规律，因此将其看作一组符合高斯分布的数集；假设其先验概率是P1、P2......Pn，第一次从n个点中随机抽取s个点作为初始样本，根据这s个点的土壤侵蚀理论值构建统计模型，将全部n个点代入统计模型进行迭代计算，根据以下公式更新各点的先验概率值：其中，s为某次假设检验符合土壤侵蚀模型的点的个数，N为所有的观测点理论值，St、Pt分别为第t次检验的假设点集和先验概率，Pt‑1表示第t‑1次检验的先验概率，I表示局内点个数；
每次迭代后，会得到前一次迭代所构建统计模型的契合度和新的先验概率值，其中各点新先验概率根据如下公式计算：每次迭代后，可根据局内点个数是否增加判断统计模型契合度是否上升，并获得新的先验概率值，其中各点新先验概率根据如下公式计算：其中，Pi为第i点的先验概率，|ΔAi|为第i点理论土壤侵蚀值与流域理论土壤侵蚀拟合值的差的绝对值，m为设定的阈值；
每次迭代结束后，算得的各点新的先验概率，选出其中最高的s个点，重新构建统计模型，并根据本次迭代的次数依据下式调整迭代循环次数Q重复操作，直到统计模型拟合度达到阈值l以下，且所有s个点的先验概率都小于阈值m，取最后一次参与迭代的全部s个局内点作为流域代表性土壤侵蚀理论值，根据各点的最后一次参与迭代得到的权重计算流域的总体土壤侵蚀理论值；
其中，I为局内点个数，N为全部测量点个数，p为置信概率，limp＝1；设最后一次迭代的统计模型算得的流域土壤侵蚀理论值为C，设：其中，t为时间，C表示流域土壤侵蚀理论值关于时间的向量，反映土壤侵蚀理论值随时间的变化；
当s趋于n时，次数最多的点满足：lims→∞C0＝lims→nC0＝E；其中，C0表示初始时刻流域土壤侵蚀值，E为期望值；故当s→n时，用 来代替E[a0]表示初始时刻的期望值，有：对统计模型计算得到的流域土壤侵蚀理论值进行优化，将其作为一个因子，和过程噪声、观测噪声一起组成优化后的状态向量；
其中， 为k时刻优化后的流域理论土壤侵蚀量, 为 的期望值，Ck为k时间流域土壤侵蚀理论值，Vk为k时刻根据贝叶斯抽样一致性算法得到的模型的噪声，Nk为k时刻野外站点v n a观测过程中的噪声，C为不同时刻系统过程噪声，C为不同时刻系统观测噪声，C为不同时刻初步优化后的流域理论土壤侵蚀量的矩阵向量； 分别为Ck、Vk、Nk的转置，分别为初始时刻流域土壤侵蚀理论值的转置、初始时刻根据贝叶斯抽样一致性算法得到的模型噪声的转置、初始时刻野外站点观测过程中噪声的转置；根据优化后的土壤侵蚀理论值计算其协方差矩阵，获得不确定性：其中， 为初始时刻模型的总不确定性，即充分考虑了中国土壤侵蚀方程本身、抽样过程、观测过程的不确定性协方差矩阵；P0为初始时刻；Pv为初始时刻根据贝叶斯抽样一致性算法得到的模型的噪声协方差矩阵，反映抽样的不确定性；Pn为初始时刻野外站点观测过程的噪声协方差矩阵，反映野外观测过程的不确定性。
5.根据权利要求4所述的利用卡尔曼滤波精确校正土壤侵蚀理论值的方法，其特征在于：步骤2.2中，对从初始时刻到当前时刻每一个时间点的土壤侵蚀理论值C和野外站点土壤侵蚀实测值分别利用改进无迹卡尔曼滤波变化，对其进行时间更新，修正系统参数其中， 为通过k‑1时刻更新到k时刻的状态量， 为k‑1时刻的状态量， 表示k‑1时刻的过程噪声；最优模型G是一种非线性变换，通过代入代表性土壤侵蚀量点进行计算来确定该函数；
从初始时刻t＝0开始，进行迭代计算，反复调整参数，获得最优模型G，取最后一次迭代得到的G，作为时间最优状态量；
基于土壤侵蚀理论值的n个观测点，取这n个点的野外实际观测值Yi,i∈1～n,构建如下时间更新模型：其中，Yk‑1为k‑1时刻的野外观测值，Yk|k‑1为从k‑1时刻更新到k时刻的优化实测野外观测值， 为过去时刻t＝k‑1的观测噪声；H函数是通过实际土壤侵蚀量和理论土壤侵蚀量来代入确定；
根据步骤2.1中统计模型获得的各观测点权重w，将其作为Y的权重，代入时间更新模型计算，根据时间更新模型时间迭代更新，调整w的值；将新修正的w值代入下式计算最优野外站点实测土壤侵蚀量其中， 为推演得到的当前时刻最优实测值，即野外站点实测值经过滤波时间更新得到的；Yi,k|k‑1为i观测点经过时间更新得到的k时刻实测土壤侵蚀值， 为新修正的适用于Y的各观测点权重，i∈(1,……,s)。
6.根据权利要求5所述的利用卡尔曼滤波精确校正土壤侵蚀理论值的方法，其特征在于：步骤2.3中，构建经过时间更新后的土壤侵蚀理论值、和野外站点土壤侵蚀实测值之间的函数关系，利用协方差矩阵计算二者之间的偏差，再调整函数关系参数，反复重复这一过程，降低不确定性，直到不确定性低于阈值，则停止循环；其不确定性根据下式求得：其中， 表示当前时刻多观测点平均实测值和多观测点平均最优实测值之间的偏差；
其中，K为最优增益值，是后验协方差矩阵的函数 和 的函数；
构建土壤侵蚀理论值和土壤侵蚀实测值之间的函数关系，得到以下二式：
其中， 是最终的状态更新结果，即最优土壤侵蚀理论值，它由经过时间更新后的状态向量通过野外实测值的增益得到，是最终的输出值；Pk是 和 之间的偏差，通过迭代缩小其值，使之降低到阈值以下，使系统修正精度达到要求。
7.根据权利要求6所述的利用卡尔曼滤波精确校正土壤侵蚀理论值的方法，其特征在于：步骤3中，若未来时刻为K，取过去某时刻K‑t为参照标准，0≤t≤K，根据最后一次迭代更新得到的状态更新方程，如下：首先利用步骤2中构建的时间更新模型，将当前时刻的最优理论土壤侵蚀量 系统过程噪声 作为自变量输入G变化，输出即为未来时刻k+1的理论土壤侵蚀量 得到(1)式；根据2)式和(3)式利用时间迭代对Y进行时间更新，其中，Yi，k+1|k为从历史时刻更新到当前时刻的观测点实测土壤侵蚀量校正值， 为当前时刻系统的观测噪声， 为不同观测点实测值代表性权重，根据贝叶斯抽样模型获得，获得k+1时刻的实测土壤侵蚀值预测值再根据(3)式计算最优实测土壤侵蚀值，利用(4)式进行调参，使不确定性 降低到阈值以下；根据(6)式计算系统整体不确定性 其中 为各观测点理论值代表性权重，Ci，k+1|k为从历史时刻更新到当前时刻的观测点理论土壤侵蚀量校正值， 为k+1时刻的理论土壤侵蚀值预测值，再根据(7)计算未来时刻的最优增益值K，代入(8)计算最优k+1时刻理论土壤侵蚀量 并计算其相对于 的不确定性Pk+1，若大于阈值，则重复(1)式至(9)式的计算，直到不确定性降低到阈值以下为止；若不确定性小于阈值，则输出最终结果
8.根据权利要求1‑7任意一项所述的利用卡尔曼滤波精确校正土壤侵蚀理论值的方法，其特征在于：利用滤波修正过去时刻的土壤侵蚀理论值，和过去同一时刻野外站点土壤侵蚀实测值进行比较，检验卡尔曼滤波状态更新得到的优化值优化程度；
取过去某相关资料丰富时段通过土壤侵蚀方程计算得到的土壤侵蚀理论值的加权平均值 同一时刻野外站点实测土壤侵蚀量Yk+1|k及其加权平均值 代入经过步骤2得到的最终卡尔曼滤波模型如下式，获得的土壤侵蚀理论值优化结果将结果 和同一时段野外站点真实测量值 进行比较，判断改进无迹卡尔曼滤波校正模型优化结果是否达到理想状态，并根据检验结果决定是否重复步骤2，直到符合理想要求为止。
9.一种利用卡尔曼滤波精确校正土壤侵蚀理论值的系统，其特征在于，包括以下模块：模块1，用于获取研究区水土流失方程各因子的高时空分辨率数据，计算研究区土壤侵蚀理论值；
模块2，用于利用改进无迹卡尔曼滤波校正模型建立土壤侵蚀实测值和土壤侵蚀理论值相关性方程，并通过时间迭代进行更新拟合；
具体包括以下子模块：
模块2.1，用于将流域所有观测点土壤侵蚀理论值作为一个集合，计算其先验概率，抽取其中概率高于预设阈值的观测点理论值构建最小二乘法统计模型，采用贝叶斯抽样一致性算法进行抽样剔除粗差，抽取可靠的代表性样点，计算流域总体土壤侵蚀理论值；
模块2.2，用于构建卡尔曼滤波时间更新模型，将前一时刻的土壤侵蚀理论值加权平均值作为输入量，当前时刻土壤侵蚀理论值的加权平均值作为输出量，进行时间递推；计算土壤侵蚀理论值的时间不确定性，使其达到阈值，满足精度要求；
模块2.3，用于计算增益值K，将其作为参数，输入前一过程中更新后的土壤侵蚀理论值和实测值，输出理论土壤侵蚀理论值，并计算理论值和实测值之间的状态不确定性，尽可能提高状态更新精度；
模块3，用于利用时间更新中经过调整参数最后获得的卡尔曼滤波时间更新方程，将当前时刻的理论土壤侵蚀量代入，获得未来时刻的土壤侵蚀理论值。 说明书 : 利用卡尔曼滤波精确校正土壤侵蚀理论值的方法及系统技术领域[0001] 本发明属于地理信息科学与水土保持技术领域，涉及一种校正土壤侵蚀理论值的方法及系统，具体涉及一种利用改进卡尔曼滤波精确校正土壤侵蚀理论值的方法及系统。背景技术[0002] 水土流失问题是当前我国乃至全世界都面临的一个重大环境问题，已知的计算土壤侵蚀量的方程包括通用水土流失方程USLE、由刘宝元改进的适用于中国地区的中国水土流失方程CSLE等，在估算长江流域水土流失量上尚可，但涉及精确校正则表现欠佳，亟需准确性更高的计算方法。[0003] 土壤侵蚀理论值的精确校正属于时空环境下的非线性问题，而卡尔曼滤波能够通过时空双重更新来提高线性及非线性问题的精度。尤其是在解决非线性问题上，拓展无迹卡尔曼和无迹卡尔曼滤波具有显著优势；其中，无迹卡尔曼滤波对非线性函数的概率密度分布进行近似，用一系列确定样本来逼近状态的后验概率密度，而不是对非线性函数进行近似，无需求导计算雅可比矩阵，有效克服了拓展卡尔曼滤波精度低，计算繁琐等问题，具有准确性高，计算量小等优点。[0004] 而在初次样本的选取上，无迹卡尔曼滤波中所采用的sigma点抽样法主要根据高斯分布抽样，但如果样本存在较大粗差，会导致理论值精度产生较大误差，基于RANSAC随机抽样一致性的贝叶斯抽样一致性算法具有抗噪性强等优点，故采用贝叶斯抽样一致性算法代替sigma点抽样法，改进抽样样本，加强抗噪性，避免异常影响。发明内容[0005] 为了解决上述技术问题，本发明提供了一种效率较高，并可用于大流域的利用改进卡尔曼滤波基于中国水土流失方程获得较精确的土壤侵蚀理论值的方法及系统。[0006] 本发明的方法所采用的技术方案是：一种利用卡尔曼滤波精确校正土壤侵蚀理论值的方法，包括以下步骤：[0007] 步骤1：获取研究区水土流失方程各因子的高时空分辨率数据，计算研究区土壤侵蚀理论值；[0008] 步骤2：利用改进无迹卡尔曼滤波校正模型建立实际土壤侵蚀实测值和土壤侵蚀理论值相关性方程，并通过时间迭代进行更新拟合；[0009] 具体包括以下子步骤：[0010] 步骤2.1：将流域所有观测点土壤侵蚀理论值作为一个集合，计算其先验概率，抽取其中概率高于预设阈值的观测点理论值构建最小二乘法统计模型，采用贝叶斯抽样一致性算法进行抽样剔除粗差，抽取可靠的代表性样点，计算流域总体土壤侵蚀理论值；[0011] 步骤2.2：构建卡尔曼滤波时间更新模型，将前一时刻的土壤侵蚀理论值加权平均值作为输入量，当前时刻土壤侵蚀理论值的加权平均值作为输出量，进行时间递推；计算土壤侵蚀理论值的时间不确定性，使其达到阈值，满足精度要求；[0012] 步骤2.3：计算增益值K，将其作为参数，输入前一过程中更新后的土壤侵蚀理论值和实测值，输出理论土壤侵蚀理论值，并计算理论值和实测值之间的状态不确定性，尽可能提高状态更新精度；[0013] 步骤3：利用时间更新中经过调整参数最后获得的卡尔曼滤波时间更新方程，将当前时刻的理论土壤侵蚀量代入，获得未来时刻的土壤侵蚀理论值。[0014] 本发明的系统所采用的技术方案是：一种利用卡尔曼滤波精确校正土壤侵蚀理论值的系统，包括以下模块：[0015] 模块1，用于获取研究区水土流失方程各因子的高时空分辨率数据，计算研究区土壤侵蚀理论值；[0016] 模块2，用于利用改进无迹卡尔曼滤波校正模型建立土壤侵蚀实测值和土壤侵蚀理论值相关性方程，并通过时间迭代进行更新拟合；[0017] 具体包括以下子模块：[0018] 模块2.1，用于将流域所有观测点土壤侵蚀理论值作为一个集合，计算其先验概率，抽取其中概率高于预设阈值的观测点理论值构建最小二乘法统计模型，采用贝叶斯抽样一致性算法进行抽样剔除粗差，抽取可靠的代表性样点，计算流域总体土壤侵蚀理论值；[0019] 模块2.2，用于构建卡尔曼滤波时间更新模型，将前一时刻的土壤侵蚀理论值加权平均值作为输入量，当前时刻土壤侵蚀理论值的加权平均值作为输出量，进行时间递推；计算土壤侵蚀理论值的时间不确定性，使其达到阈值，满足精度要求；[0020] 模块2.3，用于计算增益值K，将其作为参数，输入前一过程中更新后的土壤侵蚀理论值和实测值，输出理论土壤侵蚀理论值，并计算理论值和实测值之间的状态不确定性，尽可能提高状态更新精度；[0021] 模块3，用于利用时间更新中经过调整参数最后获得的卡尔曼滤波时间更新方程，将当前时刻的理论土壤侵蚀量代入，获得未来时刻的土壤侵蚀理论值。[0022] 本发明根据已有的CSLE计算得到的土壤侵蚀理论值，和野外站点的实测数据，利用卡尔曼滤波实现对该理论值的自动校正。提高CSLE计算的精度并拓宽了水土流失方程的应用范围，提高了实用性。附图说明[0023] 图1是本发明实施例的流程图。[0024] 图2本发明实施例的改进卡尔曼滤波修正和预测过程、检验过程示意图具体实施方式[0025] 为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。[0026] 本发明的技术方案为利用改进卡尔曼滤波校正水土流失方程计算得到的理论值的方法，以中国水土流失方程(ChineseSoilLossEquation，CSLE)和利用贝叶斯抽样一致性算法(Bayessampleconsensus，BaySAC)改进的无迹卡尔曼滤波(theUnscentedKalmanfilter，UKF)为例，请见图1，本发明提供的一种利用卡尔曼滤波精确校正土壤侵蚀理论值的方法，包括以下步骤：[0027] 步骤1：获取研究区水土流失方程各因子的高时空分辨率数据，计算研究区土壤侵蚀理论值；[0028] 本实施例中，水土流失方程各因子的高时空分辨率数据，包括与A＝R×K×L×S×B×E×T各因子的相关数据；对研究区高时空分辨率数据进行解译，获得植被覆盖度数据、土地类型数据和相关工程措施数据E；通过对研究区DEM影像的提取分析，得到坡度数据S、坡长数据L；[0029] 其中：A为年均土壤侵蚀量，单位thm‑2a‑1；R为降雨侵蚀力，反映降雨强度对土壤的侵蚀作用，是土壤侵蚀的主导因子，基于日降雨量，以半月为时长进行计算，单位MJmm‑2 ‑1 ‑1hm h a ；K为土壤可蚀性因子，反映土壤的属性，采用Williams等提出的侵蚀力影响计算2 ‑2 ‑1 ‑1器来计算，单位thmhhm MJ mm ；L坡长因子、S坡度因子，反映地形地貌特征对土壤侵蚀的影响，采用刘宝元坡度公式计算S，并根据USLE通用手册中给出的公式计算L；B为生物措施因子，取值(0‑1)；E为工程措施因子，取值(0‑1)；T为耕作措施因子，取值(0‑1)；L、S、B、E、T是无量纲。[0030] 对遥感影像数据进行波段融合，计算植被指数(NDVI)，进而根据像元二分模型计算植被覆盖度(FVC)，参考经验值，计算B；对目标区的遥感影像进行室内解译并将其结果矢量化；根据解译结果并参考E、T经验值，提取目标区E、T因子。[0031] 本实施例针对中国水土流失方程CSLE，首先，设x为关于理论土壤侵蚀量和时间的状态向量，对土壤侵蚀理论值各因子分别进行校正；[0032] x1＝R[0033] x2＝K[0034] x3＝L[0035] x4＝S[0036] x5＝B[0037] x6＝E[0038] x7＝T；[0039] 引入对应的实测土壤侵蚀量，分别对7个因子进行校正，得到精确的土壤侵蚀理论值；[0040] 然后对土壤侵蚀理论值整体进行校正：[0041] 根据CSLE基本方程A＝R×K×L×S×B×E×T得到：[0042][0043][0044] CSLE的状态转移矩阵：[0045][0046] 步骤2：利用改进无迹卡尔曼滤波校正模型建立实际土壤侵蚀实测值和土壤侵蚀理论值相关性方程，并通过时间迭代进行更新拟合；[0047] 请见图2，本实施例中步骤2的具体包括以下子步骤：[0048] 步骤2.1：将流域所有观测点土壤侵蚀理论值作为一个集合，计算其先验概率，抽取其中概率高于预设阈值的观测点理论值构建统计模型，采用贝叶斯抽样一致性算法进行抽样剔除粗差，抽取可靠的代表性样点，计算流域总体土壤侵蚀理论值；[0049] 本实施例中，从利用n个观测点的相关数据计算出的流域内各观测点某一时刻土壤侵蚀理论值中选取能够代表该时刻流域理论土壤侵蚀值的s个观测点理论值，并计算各点的权重，构建统计模型来拟合流域整体理论土壤侵蚀情况；由于同一时刻各观测点理论土壤侵蚀值基本符合高斯分布规律，因此将其看作一组符合高斯分布的数集；假设其先验概率是P1、P2......Pn，第一次从n个点中随机抽取s个点作为初始样本，根据这s个点的土壤侵蚀理论值构建统计模型。[0050] 统计模型依据最小二乘法构建:[0051][0052][0053] 式中，C为流域总体土壤侵蚀理论值，ci为s个点中某点土壤侵蚀理论值(i∈1～s)，wi为这s个观测点理论值拟合流域平均土壤侵蚀量的均值。为这s个观测点理论值的平均值。[0054] 将全部n个点代入统计模型进行迭代计算，根据以下公式更新各点的先验概率值：[0055][0056] 其中，s为某次假设检验符合土壤侵蚀模型的点的个数，N为所有的观测点理论值，St、Pt分别为第t次检验的假设点集和先验概率，Pt‑1表示第t‑1次检验的先验概率，I表示局内点个数；[0057] 每次迭代后，会得到前一次迭代所构建统计模型的契合度和新的先验概率值，其中各点新先验概率根据如下公式计算：[0058] 每次迭代后，可根据局内点个数是否增加判断统计模型契合度是否上升，并获得新的先验概率值，其中各点新先验概率根据如下公式计算：[0059][0060] 其中，Pi为第i点的先验概率，|ΔAi|为第i点理论土壤侵蚀值与流域理论土壤侵蚀拟合值的差的绝对值，m为设定的阈值；[0061] 每次迭代结束后，算得的各点新的先验概率，选出其中最高的s个点，重新构建统计模型，并根据本次迭代的次数依据下式调整迭代循环次数Q重复操作，直到统计模型拟合度达到阈值l以下，且所有s个点的先验概率都小于阈值m，取最后一次参与迭代的全部s个局内点作为流域代表性土壤侵蚀理论值，根据各点的最后一次参与迭代得到的权重计算流域的总体土壤侵蚀理论值；[0062][0063] 其中，I为局内点个数(N个测量点中基于此式落在阈值范围内的点)，N为全部测量点个数，p为置信概率，limp＝1；设最后一次迭代的统计模型算得的流域土壤侵蚀理论值为C，设：[0064][0065] 其中，t为时间，C表示流域土壤侵蚀理论值关于时间的向量，反映土壤侵蚀理论值随时间的变化；[0066] 根据高斯分布的特点，当s趋于无穷大(在本说明实例中趋于n)时，一组数集中出现频率最高，次数最多的点满足：lims→∞C0＝lims→nC0＝E；其中，C0表示初始时刻流域土壤侵蚀值，E为期望值；故当s→n时，用 来代替E[a0]表示初始时刻的期望值，有：[0067][0068][0069] 由于观测时系统本身存在一定偏差，当同时考虑贝叶斯抽样一致性算法选取代表性点过程中和野外站点观测获取用于计算理论侵蚀值过程中的误差时，对统计模型计算得到的流域土壤侵蚀理论值进行优化，将其作为一个因子，和其他影响因子(过程噪声、观测噪声)一起组成优化后的状态向量；[0070][0071][0072][0073] 其中， 为k时刻优化后的流域理论土壤侵蚀量, 为 的期望值，Ck为k时间流域土壤侵蚀理论值，Vk为k时刻根据贝叶斯抽样一致性算法得到的模型的噪声，Nk为k时刻野外v n a站点观测过程中的噪声，C为不同时刻系统过程噪声，C为不同时刻系统观测噪声，C为不同时刻初步优化后的流域理论土壤侵蚀量的矩阵向量； 分别为Ck、Vk、Nk的转置，分别为初始时刻流域土壤侵蚀理论值的转置、初始时刻根据贝叶斯抽样一致性算法得到的模型噪声的转置、初始时刻野外站点观测过程中噪声的转置；[0074] 根据优化后的土壤侵蚀理论值计算其协方差矩阵，获得不确定性：[0075][0076] 其中， 为初始时刻模型的总不确定性，即充分考虑了中国土壤侵蚀方程本身、抽样过程、观测过程的不确定性协方差矩阵；P0为初始时刻；Pv为初始时刻根据贝叶斯抽样一致性算法得到的模型的噪声协方差矩阵，反映抽样的不确定性；Pn为初始时刻野外站点观测过程的噪声协方差矩阵，反映野外观测过程的不确定性。[0077] 步骤2.2：构建卡尔曼滤波时间更新模型，将前一时刻的土壤侵蚀理论值加权平均值作为输入量，当前时刻土壤侵蚀理论值的加权平均值作为输出量，进行时间递推；计算土壤侵蚀理论值的时间不确定性，使其达到阈值，满足精度要求；[0078] 本实施例中，对从初始时刻到当前时刻每一个时间点的土壤侵蚀理论值C和野外站点土壤侵蚀实测值分别利用改进无迹卡尔曼滤波变化，对其进行时间更新，修正系统参数[0079][0080] 其中， 为通过k‑1时刻更新到k时刻的状态量， 为k‑1时刻的状态量，表示k‑1时刻的过程噪声；最优模型G是一种非线性变换，通过代入代表性土壤侵蚀量点进行计算来确定该函数(此处假设已知函数，若未知，可以给出线性和非线性函数然后通过最小二乘法拟合一个函数)；[0081] 从初始时刻t＝0开始,进行迭代计算，反复调整参数，计算时间不确定性，使其小于阈值，获得最优模型G，取最后一次迭代得到的G，作为时间最优状态量；时间不确定性根据下式计算：[0082][0083] 式中， 表示时间更新状态量 的时间不确定性，Ci,k|k‑1为全部s个代表性观测点中的第i个观测点k时刻理论土壤侵蚀量。[0084] 基于土壤侵蚀理论值的n个观测点，取这n个点的野外实际观测值Yi,i∈1～n,构建如下时间更新模型：[0085][0086] 其中，Yk‑1为k‑1时刻的野外观测值，Yk|k‑1为从k‑1时刻更新到k时刻的优化实测野外观测值， 为过去时刻t＝k‑1的观测噪声；H函数是通过实际土壤侵蚀量和理论土壤侵蚀量来代入确定(此处假设已知函数，若未知，可以给出线性和非线性函数然后通过最小二乘法拟合一个函数)；[0087] 根据步骤2.1中统计模型获得的各观测点权重w，将其作为Y的权重，代入时间更新模型计算，根据模型时间迭代更新，调整w的值；将新修正的w值代入下式计算最优野外站点实测土壤侵蚀量[0088][0089] 其中， 为推演得到的当前时刻最优实测值，即野外站点实测值经过滤波时间更新得到的；Yi,k|k‑1为i观测点经过时间更新得到的k时刻实测土壤侵蚀值， 为新修正的适用于Y的各观测点权重，i∈(1,……,s)；[0090] 步骤2.3：计算增益值K(作为中间量代入计算)，将其作为参数，输入前一过程中更新后的土壤侵蚀理论值和实测值，输出理论土壤侵蚀理论值，并计算理论值和实测值之间的状态不确定性，尽可能提高状态更新精度；[0091] 本实施例中，构建经过时间更新后的土壤侵蚀理论值和野外站点土壤侵蚀实测值之间的函数关系，利用协方差矩阵计算二者之间的偏差，再调整函数关系参数，反复重复这一过程，降低不确定性，直到不确定性低于阈值，则停止循环；其不确定性根据下式求得：[0092][0093] 其中， 表示当前时刻多观测点平均实测值和多观测点平均最优实测值之间的偏差；[0094][0095][0096] 其中，K为最优增益值，是后验协方差矩阵的函数 和 的函数；[0097] 构建土壤侵蚀理论值和土壤侵蚀实测值之间的函数关系，得到以下二式：[0098][0099][0100] 其中， 是最终的状态更新结果，即最优土壤侵蚀理论值，它由经过时间更新后的状态向量通过野外实测值的增益得到，是最终的输出值；Pk是 和 之间的偏差，通过迭代缩小其值，使之降低到阈值以下，使系统修正精度达到要求。[0101] 步骤3：利用时间更新中经过调整参数最后获得的卡尔曼滤波时间更新方程，将当前时刻的理论土壤侵蚀量代入，获得未来时刻的土壤侵蚀理论值。[0102] 若对精度要求不高，到此即可输出结果，若要求更高精度，则需要再对实测土壤侵蚀量进行时间更新，通过状态更新，输出优化后的外来时刻理论土壤侵蚀量。[0103] 本实施例中，若未来时刻为k，取过去某时刻k‑t为参照标准，0≤t≤K，根据最后一次迭代更新得到的状态更新方程，如下：[0104][0105][0106][0107][0108][0109][0110][0111][0112][0113] 首先利用步骤2中构建的时间更新模型，将当前时刻的最优理论土壤侵蚀量 系统过程噪声 作为自变量输入G变化，输出即为未来时刻(k+1)的理论土壤侵蚀量 得到(1)式；根据2)式和(3)式利用时间迭代对Y进行时间更新，其中，Yi,k+1|k为从历史时刻更新到当前时刻的观测点实测土壤侵蚀量校正值， 为当前时刻系统的观测噪声， 为不同观测点实测值代表性权重，根据贝叶斯抽样模型获得，获得k+1时刻的实测土壤侵蚀值预测值 再根据(3)式计算最优实测土壤侵蚀值，利用(4)式进行调参，使不确定性降低到阈值以下；根据(6)式计算系统整体不确定性 其中 为各观测点理论值代表性权重，Ci,k+1|k为从历史时刻更新到当前时刻的观测点理论土壤侵蚀量校正值，为k+1时刻的理论土壤侵蚀值预测值，再根据(7)计算未来时刻的最优增益值K，代入(8)计算最优k+1时刻理论土壤侵蚀量 并计算其相对于 的不确定性Pk+1，若大于阈值，则重复(1)式至(9)式的计算，直到不确定性降低到阈值以下为止；若不确定性小于阈值，则输出最终结果[0114] 步骤4：利用滤波修正过去时刻的土壤侵蚀理论值，和过去同一时刻野外站点土壤侵蚀实测值进行比较，检验卡尔曼滤波状态更新得到的优化值优化程度。[0115] 本实施例中，取过去某相关资料丰富时段通过土壤侵蚀方程计算得到的土壤侵蚀理论值的加权平均值 同一时刻野外站点实测土壤侵蚀量Yk+1|k及其加权平均值代入经过步骤2得到的最终卡尔曼滤波模型如下式，获得的土壤侵蚀理论值优化结果[0116][0117] 将结果 和同一时段野外站点真实测量值 进行比较，判断改进无迹卡尔曼滤波校正模型优化结果是否达到理想状态，并根据检验结果决定是否重复步骤2，直到符合理想要求为止。[0118] 本发明根据已有的CSLE计算得到的土壤侵蚀理论值，和野外站点的实测数据，利用卡尔曼滤波实现对该理论值的自动校正。提高CSLE计算的精度并拓宽了水土流失方程的应用范围，提高了实用性。[0119] 应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

专利地区：湖北

专利申请日期：2021-12-07

专利公开日期：2024-06-18

专利公告号：CN114357719B