一种基于贝叶斯推理的土地利用数据同化方法

摘要

本发明涉及一种基于贝叶斯推理的土地利用数据同化方法，该方法的核心是将离散的类型分布变量对应到连续的类型分布密度，并将连续的类型分布密度作为中间变量，进而计算类型变量的后验分布概率。本发明将基于共轭先验的贝叶斯推理思想引入到土地利用动态模型中，通过融合观测数据，对离散的、多维的模型模拟的土地利用类型进行同化，提高土地利用动态模型的模拟预测精度。本发明在模型方面，首次利用多维土地利用/覆被变化动态模型对多种土地利用类型进行同化研究。本发明在方法上，首次引入一种同化策略，对离散的类型分布状态变量进行更新。

权利要求

1.一种基于贝叶斯推理的土地利用数据同化方法，其特征在于：所述同化方法包括以下步骤：第1步：利用历史变化数据集对土地利用动态模型进行校准，获得模型的参数值；

第2步：将研究区划分为若干个正方形网格，每个网格包含ｋ×ｋ个像元；选择其中的若干个正方形网格作为观测点；

第3步：运行土地利用动态模型，当模拟到观测年份时，得到土地利用类型转换概率，以及模拟值α，同时提取该同化年份的观测值ｃ；

第4步：运行基于共轭先验的贝叶斯推理公式，得到同化值；其中，基于共轭先验的贝叶斯推理公式如下：式中， 是同化后的土地利用类型 的分布密度值；是土地利用类型的观测数目；

代表模型中模拟的土地利用类型数目；N是所有土地利用类型的观测数目； 是模拟的所有土地利用类型数目；

第5步：根据同化值，对当前的模拟结果进行修正；

具体修正步骤如下：

S1：在观测点 ，基于同化值 和正方形方格的面积 ，计算得到各土地利用类型的变化数量 ，其计算公式如下： S2：依次寻找变化数量 的土地利用类型；

 S3：将 的土地利用类型 的转换概率从小到大逐像元减少个元胞，即像元；同时将这些像元的土地利用类型赋值为0；

S4：将所有 的土地利用类型按照步骤S3进行重新分配，最终得到观测点中需重新分配的所有元胞；

S5：在待分配的元胞中，将剩余的土地利用类型，即 的土地利用类型，依次增加 个元胞；

S6：重复步骤S1-S5，直至完成所有观测点的修正；

第6步：把修正后的结果再输入到模型中，作为初始值，进行一下年份的模拟同化；

当模拟到下一观测时刻，重复第3步～第5步，直到模拟同化年份结束。

2.如权利要求1所述的一种基于贝叶斯推理的土地利用数据同化方法，其特征在于：在第3步中的模拟值和观测值均指各土地利用类型对应的数目。

说明书

一种基于贝叶斯推理的土地利用数据同化方法

[0001] 摘要

[0002] 本发明发展了一种基于贝叶斯推理的土地利用数据同化方法，该方法的核心是将离散的类型分布变量对应到连续的类型分布密度，并将连续的类型分布密度作为中间变量，进而计算类型变量的后验分布概率。本发明将基于共轭先验的贝叶斯推理思想引入到土地利用动态模型中，通过融合观测数据，对离散的、多维的模型模拟的土地利用类型进行同化，提高土地利用动态模型的模拟预测精度。

技术领域

[0003] 本发明涉及地理信息科学技术领域，具体来说是一种基于贝叶斯推理的土地利用数据同化方法。

背景技术

[0004] 土地利用变化是人类活动对地球表层系统影响最直接的表现方式，通过改变地球表层系统的下垫面，影响了地球表层系统各大圈层的循环过程。因此，准确模拟和预测过去、现在、未来地球表层真实的土地利用变化格局和过程，可为土地利用规划、政策决策提供科学依据，也可服务于区域可持续发展。

[0005] 模型和观测是土地利用变化研究的两种基本手段。模型有助于深入了解土地利用变化成因，重建过去、预测未来发展趋势、评估环境影响以及支持土地利用规划决策分析。目前，元胞自动机（Cellular Automata，CA）已成为土地利用变化模拟最常用的工具之一。
但元胞自动机模型结构、静态参数以及边界条件、初始场引起的误差和不确定性会随着模型的运行不断的传递和积累；而观测能够快速准确地获取各种时空尺度上观测时刻的地表真实土地利用状况。但无法了解无观测时刻的状况，不能提供时间上连续的土地利用变化情况。如何将观测信息和模型有机融合，来量化和减少土地利用模拟和预测中的各种不确定性，获取时空连续的土地利用信息，是目前土地利用变化研究中面临的一重大挑战。

[0006] 数据同化（Data Assimilation），也称“模型-数据融合”（mode-data fusion），是一种基于模型的动力学框架，通过融合不同来源、不同分辨率的直接和间接观测，在考虑模型和各种观测误差的基础上不断地依靠观测而自动调整模型轨迹、优化模型状态，并可同步优化参数、定量评价和减少不确定性的一种方法。数据同化为融合观测和土地利用模型提供了理论框架，为该领域的研究方法注入了新鲜的血液。

[0007] 然而，无论国际还是国内，土地利用数据同化研究刚刚起步，理论和方法还存在不完善的地方。不同于水文、气象、和陆地过程模型，土地利用动态模型的输出为离散的类型变量（土地利用类型），缺少能够直接用于同化的状态变量。纵观国内外土地利用数据同化相关研究可知，（1）当前研究中使用的均为二维的土地利用动态模型，没有开展对多维的土地利用动态模型的同化。（2）绝大多数研究是基于粒子滤波和集合卡尔曼滤波算法对连续的模型参数进行优化，缺少对离散的、多维的类型分布状态变量（土地利用类型）的更新。

[0008] 鉴于此，有必要瞄准国际上土地利用研究中模型-观测融合方法。

发明内容

[0009] 针对上述的两个主要问题，本发明从方法论的角度提出了一种基于贝叶斯推理的土地利用数据同化策略，为提高土地利用动态模型的模拟预测精度，以及离散型变量的数据同化方法研究提供了较好的解决方案。

[0010] 土地利用动态模型的输出是一个离散的土地利用类型。基于贝叶斯推理方法对土地利用进行同化的核心思想是将离散的类型分布变量对应到连续的类型分布密度，并将连续的类型分布密度作为中间变量，来计算类型变量的后验分布概率。

[0011] 具体土地利用数据同化实施步骤如下：

[0012] S1：运行土地利用动态模型，得到同化年份的模拟值α，同时提取该同化年份的观测值ｃ。此处的模拟值和观测值均指各土地利用类型对应的数目。

[0013] S2：运用基于共轭先验的贝叶斯推理公式，得到同化值，即各观测点中各土地利用类型对应的分布密度值。其中，基于共轭先验的贝叶斯推理公式如下：

[0014]

[0015] 式中， 是同化后的土地利用类型 的分布密度值；是土地利用类型的观测数目； 代表模型中模拟的土地利用类型 数目； 是所有土地利用类型的观测数目； 是模拟的所有土地利用类型数目。

[0016] S3：根据同化值，对当前的模拟结果，即模拟的土地利用类型进行修正，获得该年份同化后的土地利用类型。

[0017] 本发明的有益效果是：

[0018] 本发明利用研究中模型-观测融合方法这一崭新的发展方向，从方法论的角度探索土地利用研究的新途径，在多维土地利用动态模型的基础上，融合观测数据，优化多维离散的类型分布状态变量，以其更准确地给出土地利用的空间和时间分布特征。

[0019] 本发明在模型方面，首次利用多维土地利用/覆被变化动态模型对多种土地利用类型进行同化研究。

[0020] 本发明在方法上，首次引入一种同化策略，对离散的类型分布状态变量进行更新。

附图说明

[0021] 图1为本发明实施土地利用数据同化方法的流程图。

[0022] 图2为对土地利用动态模型的模拟结果进行修正的流程图。

[0023] 图3为本发明的同化结果与模拟结果以及2001年和2009年真实土地利用空间分布对比示意图

具体实施方式

[0024] 下面结合附图和具体实施例对本发明做进一步阐述说明。

[0025] 实施例1

[0026] 本发明中的研究对象为位于黑河流域中游地区的张掖绿洲甘州区。甘州区是张掖绿洲的行政中心，总面积约3665平方公里。甘州区是2008年张掖绿洲启动的湿地保护政策的重点实施区域，并且该区域的城镇用地扩张及扩耕在张掖绿洲较为典型。本研究中所采用土地利用动态模型为多维逻辑回归马尔科夫元胞自动机模型（MLRMCA），采用的测试数据为：1992-1999年的甘州区土地利用历史变化数据集作为MLRMCA模型的校准数据集，以确定模型参数；1999年-2009年的甘州区土地利用数据集作为同化的观测及结果验证数据，每2年进行一次同化。以上土地利用数据集由TM/ETM+影像解译而成，采用相同的分类系统。所有区域的数据均采用同坐标系、同尺度的栅格数据文件。

[0027] 图1为本发明的方法流程图，主要包括以下几个步骤：

[0028] 第1步：利用历史变化数据集对土地利用动态模型进行校准，获得模型的参数值。

[0029] 第2步：将研究区划分为若干个正方形网格，每个网格包含ｋ×ｋ个像元（元胞）。选择其中的若干个正方形网格作为观测点。

[0030] 第3步：运行土地利用动态模型，当模拟到观测年份时，得到土地利用类型转换概率，以及模拟值α，同时提取该同化年份的观测值ｃ。此处的模拟值和观测值均指各土地利用类型对应的数目。

[0031] 第4步：运行基于共轭先验的贝叶斯推理公式，得到同化值。其中，基于共轭先验的贝叶斯推理公式如下：

[0032]

[0033] 式中， 是同化后的土地利用类型 的分布密度值；是土地利用类型的观测数目； 代表模型中模拟的土地利用类型 数目； 是所有土地利用类型的观测数目； 是模拟的所有土地利用类型数目。

[0034] 第5步：根据同化值 ，对当前的模拟结果，即模拟的土地利用类型进行修正，获得该观测年份同化后的土地利用类型。

[0035] 修正模拟结果的流程见图2，具体修正步骤如下：

[0036] S1：在观测点j，基于同化值 和正方形方格的面积 ，计算得到各土地利用类型的变化数量，其计算公式如下：

[0037] S2：依次寻找变化数量的土地利用类型。

[0038] S3：将的土地利用类型的转换概率从小到大逐像元减少

[0039] 个元胞（像元），同时将这些像元的土地利用类型赋值为0。

[0040] S4：将所有的土地利用类型按照步骤S3进行重新分配，最终得到观测点j中需重新分配的所有元胞。

[0041] S5：在待分配的元胞中，将剩余的土地利用类型，即的土地利用类型，依次增加个元胞。

[0042] S6：重复步骤S1-S5，直至完成所有观测点的修正。

[0043] 第6步：把修正后的结果再输入到模型中，作为初始值，进行一下年份的模拟同化。当模拟到下一观测时刻，重复第3步～第5步，直到模拟同化年份结束。

[0044] 图3为同化结果与模型模拟的结果对比。该图表明基于本发明的方法更新后的土地利用类型空间分布更加接近于真实情况，能够提高模型的模拟预测精度。

[0045] 具体的同化与模拟精度结果如下： 

[0046] （OA为总体精度，Kappa为Kappa系数，FM为性能指数）。