专利名称：一种精细尺度下的动态风险及易损性预测方法
技术领域：
本发明一种精细尺度下的动态风险及易损性预测方法，用于对突发性公共灾害的预测，属于地球信息科学技术领域。
背景技术：
突发性公共灾害事件的监测预警在国家防灾减灾中占据着重要地位，准确及时的预警将会大幅减少生命财产的损失，提高防灾减灾的效率。建立及时、准确的监测预警系统一直是国家防灾减灾的重点；同时，由于突发公共灾害事件的突发性、随机性、影响因素的多样性及不确定性，准确及时的预报是比较困难的，而风险分析作为监测预警中的关键技术，其研究及相关技术的实现日趋重要。
当前风险分析技术主要是根据灾害的来源、发生的环境及承受者，将影响灾害风险的因素划分为致灾因子、孕灾环境及承灾体三方面，对这些方面的因子进行聚类或者模糊综合评判从而计算灾害风险的水平，是一种脆性分析方法。如对洪水灾害，通过分析直接诱因既大暴雨、孕灾环境既自然地理条件及承灾体既人口分布及房屋布局等，综合评判形成不同风险级别分区(低、中、高)。而易损性主要考虑了承灾体本身的损失特性，是风险分析的扩展到承灾体上，技术基本上同风险分析差异不大。
传统的风险分析及易损性分析方法存在着以下局限 (1)有的方法基于脆性(crisp)数值的分析方法，体现在风险水平的低、中、高或者类似级别的绝对划分，不确定性信息考虑不足(黄崇福，2001，自然灾害风险分析，北京北京师范大学出版社；魏一鸣等，2002，洪水灾害风险管理理论，北京科学出版社)； (2)还有方法是基于中尺度或大尺度下的分析结果，缺乏精细尺度下的风险概率及易损性估算，更难进行精细尺度下的空间定位，预测结果定位偏差大(王艳艳，2002，不同尺度的洪涝灾害损失评估模式述评，水利发展研究，第2卷，第12期)； (3)另外还有方法忽视灾害事件或者影响因子在空间上的相关性，但空间相关性常常对结果有着重要影响，且灾害事件作为地球表层事件，灾害事件或者诱发因子之间存在着一定的空间相关性(Li，L.，et al.2005，Typhoon insurance pricing with spatialdecision support tools，International Journal of Geographical Information Science，19(3)363-384；Tobler，W.R.，1979，Cellular Geography，Philosophy in Geography，DordrechtReidel)； (4)再有的方法缺乏根据实时资料灵活地对风险进行动态概率推理的能力，实时预测功能比较缺欠(William，J.P.，and Arthur，A.，1982，Natural Hazard Risk Assessmentand Public Policy，New YorkSpringer-Verlag New York Inc.)。
随着信息技术的发展，近年来日益发展的人工智能、数据挖掘及空间分析技术提供了更好的手段及新思路来探测日益增多的时空数据，从而可采用比原来基于少量历史数据的统计方法更有效的方法。人工智能方法是通过机器模拟人类的思考与推理，其分析结果更体现灵活性与客观性，如基于概率的不确定性推理的贝叶斯网络；数据挖掘则是从海量数据中找到最优解(模型、参数)，可为构造预测性能良好的贝叶斯网络提供学习算法；空间数据挖掘则是针对空间分布的对象的数据挖掘，由于空间对象本身具有空间属性(位置、空间相关性、拓扑关系及形状等)，空间数据挖掘则通过分析探索地物属性在空间上的相关性或变异性，为学习、预测及推理提供空间信息。
发明内容
本发明的技术解决问题克服现有针对现有的风险综合评判技术的不足，提供一种精细尺度下的动态风险及易损性预测方法，该方法扩展了训练数据的来源，提高分析的细节及精细的空间尺度定位功能，提高了预测的精度及效果，可进行实时的动态的灾害风险及易损性分析。
本发明通过以下的步骤得以实现一种精细尺度下的动态风险及易损性预测方法，步骤如下 (1)相关数据的采集及预处理按致灾因子、孕灾环境及承灾体三个采集预测因子X及灾害损失变量y，并进行预处理产生精细尺度下的训练样本数据集；所述的致灾因子是指直接导致灾害的可测量因子变量，所述的孕灾环境指在致灾因子存在的前提下导致灾害损失发生的自然地理环境，所述的承灾体特性指各种主要以人为主的承受灾害的对象遭受灾害所带来损失的性质； (2)采用融入空间相关性的核密度函数方法对步骤(1)的训练样本数据集衍生出新的训练样本数据，以弥补调查数据不足，并补充完整精细尺度下的训练样本数据； (3)预测因子的特征选择根据训练样本数据，采用三种方法衡量影响因子同目标因子的相关性，即线性相关系数、基于熵的信息获取率及粗糙集的属性重要性，结合主成分分析选择预测因子X； (4)对步骤(3)所选择的预测因子集中的连续变量进行离散化，采用模拟退火优化算法搜索贝叶斯网络模型，估计网络的条件概率参数； (5)利用生成的贝叶斯网络模型进行完全或不完全证据支持下的风险及易损性概率预测； (6)在地图上实时动态地定位所预测的风险概率及易损性，识别区域上的风险及潜在损失分布状况，为救灾减灾的实施提供辅助决策支持信息。
所述步骤(1)中的预处理方法为对致灾因子、孕灾环境及承灾体特性数据经预处理转换成统一分辨率的栅格数据集，分不同的情况采用以下步骤 a.对矢量数据预处理方法 对于点数据，采用k-最近邻方法或者Kriging方法插值算法，将点数据转换成栅格数据；对于线数据采用数字微分分析法DDA法或Bresenham方法转换成栅格数据；对于面数据采用扫描法；转换成统一分辨率的栅格数据； b.对于栅格的预处理方法 对于栅格数据，如果坐标或目标栅格数据不一致时，采用最近邻法进行重采样，形成转换成统一分辨率的栅格数据； 将矢量数据和栅格再栅格化后即形成统一分辨率的栅格数据集。
所述步骤(2)采用了融入空间相关性的核函数方法衍生样本点数据，弥补调查数据及历史资料的不足，补充完整精细尺度下的训练样本数据的步骤为 a.调查确定灾害事件及其损失情况； b.采用正态核密度函数计算估计目标变量的核函数 其中的影响宽度λ可随灾害强度的不同设置不同值，d(z，Zi)为空间上任意一点z与调查数据点Zi间的欧氏距离，z可为灾害事件y或者因子变量x； c.由各观察值点的核密度函数可推出空间任意一点的影响强度 其中，n为观察值数据点数，z为空间任意点，Zi为相应的观察值点数值，Kλ(z，Zi)为核函数。
所述步骤(3)采用三种方法衡量影响因子同目标因子的相关性具体如下 a.所述的线性相关性是采用统计学Pearson相关系数来衡量连续变量之间的线形相关性，但该方法忽略了变量之间可能存在的非线形关系及非连续变量相关关系； b.所述的基于信息熵的分类贡献函数，即Quinlan的Information Gain IG及Information Gain Ratio GR衡量辅助变量同目标变量之间的关系，选择对分类贡献较大的变量，IG及GR的定义如下 IG(y，x)＝H(y)-H(y|x) GR(y，x)＝(H(x)-H(y|x))/H(x) H(y)是熵的定义，而H(y|x)是针对类别x的条件熵的定义； c.所述的采用粗糙集属性重要性为通过计算辅助变量对目标变量的分类能力来衡量二个变量的相关性，选择重要性更大的特征变量 其中，X为因子变量集，xj是要计算属性重要性的因子变量，γ(X，y)是粗糙集中目标变量y对因子变量集X的依赖性； d.所述的采用主成分分析为选择独立性强的代表因子，删除信息冗余因子。
所述步骤(4)对选择的预测因子集中的连续变量进行离散化，采用模拟退火优化算法搜索贝叶斯网络模型，估计网络的条件概率参数的步骤如下 a.建立贝叶斯网络模型框架对贝叶斯网络B＝(BS，BP)，令目标变量y为灾害事件或相关灾害损失事件，则y状态空间定义为“灾害发生”与“灾害不发生”两种状况，或者“损失严重”与“损失不严重”两种状况，或者“损失最严重”、“损失严重”“损失一般”与“轻微损失或无损失”四种状况；若已知一些证据，即解释变量因子集X中的一些因子，通过计算边际概率而求得y的相应状态发生的后验概率，亦称信任函数 其中，为随机变量集V上的联合概率，xi既为影响因子集X中的随机变量，yj表示目标变量y取第j种状态情况； b.采用优化的离散化方法转换连续类型数据，作为因子变量的状态空间值输入贝叶斯网络模型进行计算，算法采用基于熵的衡量标准 其中，impurity(k，j，i)为当从j到i的训练样本被划分成k个区间时的最小的测量值impurity，最优的划分区间k则是impurity(k，j，i)，N为连续变量因子值的基，测量值impurity是每次划分的平均类熵，其最小使得用于编码离散化区间而进行预测的字节数最小； c.搜索最优的网络结构 采用贝叶斯质量记分函数作为标准，搜索策略是通过模拟金属的退火即硬化过程而找到网络结构的最优解，退火过程从开始加热金属致其融化与金属离子自由游动开始，温度缓慢有序地降低直到金属粒子在一个高结构化的格网中成型；一个潜在的网络相当一个解，对应一个同网络质量测量值有关的价值函数，通过退火即温度T的降低，搜索最优的网络结构； d.估计网络参数 得到网络拓扑结构后，采用基于Dirichlet分布的经典算法，叠代地估计条件概率参数值； e.风险概率及易损性预测 将所建立的贝叶斯网络用于一定证据支持下风险或易损性概率预测，采用适用于环路的信任传递算法即经典的Pearl信息传递算法实现，如所得到网络是单连接树，该算法则采用PolyTree的精确算法实现预测；而对包括环路的网络，该算法则采用为精度较高的近似算法实现预测。
所述(5)利用生成的贝叶斯网络模型进行完全或不完全证据支持下的风险及易损性概率预测方法为估算中数据输入为允许缺失数值的同分辨率及坐标的多维栅格数据，每栅格单元为一个预测对象，将每单元对应属性值输入模型，通过估计精细尺度下既一定分辨率的各个栅格单元的风险概率或易损性，从而估算区域上风险概率。
本发明与现有技术相比的优点在于本发明方法可融合多源、异构的数据进行风险(易损性)概率预测估计，突破了传统方法的局限，体现在以下几方面 (1)本发明是建立在精细尺度下既一定分辨率的栅格单元上的风险及易损性预测模型，特别是提出了在调查样本数据稀少时，结合机理可采用融合了空间信息的核密度函数方法来衍生因子及目标变量栅格数据，弥补调查样本缺乏的不足。核密度方法，可克服在进行精细尺度分析时数据稀少而不能进行的缺点，有利于精细尺度下的风险及易损性估算及预测，也克服了中或粗尺度下预测效果低，不能很好地进行空间定位的缺陷。
(2)通过采用优化的离散化方法离散化连续变量，使得本发明能融合不同数据类型(连续、离散及类别)采用贝叶斯网络模型进行风险及易损性概率预测，提高了预测因子的信息含量，这是传统风险分析方法所不具备的。
(3)采用结构灵活的贝叶斯网络模型来建立预测因子同目标变量之间的概率依赖关系，灵活多变的网络拓扑有助于融入不同来源及类型的数据进行风险概率的预测，为风险分析提供是一种不确定性信息，而不是“脆”的风险分级信息，体现在地图上则是风险概率的栅格图像，提供了比传统分析方法更客观的软信息。
(4)在网络的预测因子及训练过程中融合了空间信息，而贝叶斯网络灵活多变的网络拓扑也有助于通过概率依赖关系融入因子变量之间复杂的空间相关关系，算法考虑了空间相关性及空间变异性的影响，这也是传统方法所缺欠的。
(5)通过优化的模拟退火搜索算法可学到跟数据背景适应的预测模型，即使在不完备信息既“证据”情况下，也可通过环路Pearl的信息传播算法估计风险概率水平，比传统方法有更强的概率估算及预测功能。
图1为时间与灾害风险之间的关系图； 图2为本发明方法流程图； 图3为本发明具体实施方式
中结合的典型研究案例地图黑河流域地区； 图4为本发明方法从数据中搜索得到的优化的贝叶斯网络结构图； 图5为三种方法，既本发明方法、
Bayes及逻辑斯特方法的交叉验证结果的受试者特征(ROC)曲线的对比，其中本发明方法更靠近最优分类点，说明预测性能较好； 图6为三种方法，既本发明方法、
Bayes及逻辑斯特方法的由2006年的数据预测2007年的灾害风险结果的ROC曲线对比； 图7为本发明提出的贝叶斯网络模型预测的各栅格单元的风险概率图。
具体实施例方式 如图1所示，为时间与灾害风险之间的经验关系，说明了时间对减灾救灾的重要性，而本发明提出的方法通过预测与评估可适当提高减灾效率。
一、如图2所示，本发明方法的具体实现步骤为 1、数据集的采集及预处理，本发明通过融合精细尺度下多种类型的(包括三种连续、离散及类别)及多个来源的数据的空间数据集X进行建模预测。这些数据经过适当的预处理及转换，得到训练样本及测试数据集。
所搜集的数据可根据灾害的原理划分成三方面 (1)致灾因子，直接导致灾害的可测量因子变量，不同灾种有不同致灾因子，如洪水的致灾因子可为降雨，而台风的致灾因子则为风速。
(2)孕灾环境，在致灾因子存在的前提下导致灾害损失发生的自然地理环境，如高程、坡度、地貌、植被、土地利用状况等。
(3)承灾体属性，指各种承受灾害的对象(主要以人为主)遭受灾害所带来损失的性质，包括社会经济(如GDP)、收入、人口分布等。
这些不同来源的数据经过处理转换成统一栅格单元大小的栅格数据集。预处理涉及到以下几方面 (1)数据的插补及缺失数据分析 (a)点数据需要转换成栅格数据，转化考虑了空间相关性，采用k-最近邻方法或Kriging方法插值算法； (b)而如多维数据有某些缺失项，采用融合其他域变量信息的期望最大化，插补后数据更符合实际数据分布。
(2)矢量栅格化 对于矢量栅格化，不同类型的矢量数据有不同处理方法 (a)点数据采用前面提到的插值算法，既k-最近邻或者Kriging；线数据采用可采用DDA法(数字微分分析法)和Bresenham方法； (b)而对于面数据采用扫描法。
(3)栅格的再栅格化 对栅格的再栅格化，目的是形成同样栅格单元大小的数据集，也称重采样，既用原格网中的单元值填充新格网中的单元值，采用k-最近邻方法。
2、采用融入空间相关性的核密度函数方法衍生训练样本，补充完成精细尺度下的栅格数据集。根据空间自相关原理，既“越近越相关”，结合遥感影像、历史记录及调查资料，在样本数据少，获取昂贵情况下，可采用核函数法衍生样本数据，以弥补样本数据的不足；具有以下步骤 (1)调查确定灾害事件及其损失情况，或者相关属性的取值情况。
(2)采用核函数进行估算，该方法依据空间相关性机理，认为影响点或线，如泛滥的河流，对周边的影响随距离逐步递减，采用正态核密度函数 其中的影响宽度λ可随灾害强度的不同设置不同值，d(z，Zi)为空间上任意一点z与调查数据点Zi间的欧氏距离。z可为灾害事件y或者因子变量x当z为灾害事件y时，d(y，Yi)表示为空间上任意一点的灾害事件的强度或损失y同空间上历史或调查事件点的观察值Yi间的欧氏距离，可以将灾害事件强度或损失划分成不同的风险影响级别(如低、中、高)而获得不同级别风险水平的样本；当z为影响因子x时，d(x，Xi)表示空间上任意一点该因子的强度x同空间上历史或调查事件点的观察值Xi间的欧氏距离，如离开河流的远近，靠近河流的更容易遭受损失，远离的遭受损失的概率就低。
(3)由各观察值点的核密度函数可推出空间任意一点的影响强度 其中，n为观察值数据点数，z为空间任意点，Zi为相应的观察值点数值，Kλ(z，Zi)为式(1)所确定的核函数。根据对影响强度Density(z)的分级，由此将栅格单元划分成不同的风险影响级别或水平的训练样本。
核函数方法符合空间数据的基本性质既空间相关性，可弥补采用遥感数据、历史数据或实地调查资料的不足，产生更全面的训练样本数据。主要体现在两方面就风险影响因子而言，可在矢量栅格化中采用核密度方法衍生影响因子栅格图(如河流、城镇居民地)；而就目标变量(灾害事件而言)，可采用核密度方法衍生出目标变量的核心(如高、低二级，高为核心地带)或者多级影响(如高、中、底三级)的空间栅格单元作为训练样本。通过核函数方法，实现了采用少量样本数据衍生精细尺度下的栅格训练样本数据。
3、选择影响因子集，估计影响因子与目标因子之间的相关关系，选择有强相关性的因子，既特征选择。本发明采用了三种不同的方法分析影响因子同目标因子的相关关系，采用了主成分分析提取独立的变量，移除冗余信息。特征选择的步骤为 (1)采用统计学中的相关系数来衡量连续变量之间的线形相关性，但该方法忽略了变量之间可能存在的非线形关系及非连续变量相关关系； (2)采用了基于信息熵的分类贡献函数，既Quinlan的Information Gain(IG)及Information Gain Ratio(GR)，衡量辅助变量同目标变量间关系，选择对分类贡献较大的变量(即IG或GR较大)，IG及GR的定义如下 IG(y，x)＝H(x)-H(y|x) (3) GR(y，x)＝(H(x)-H(y|x))/H(x)(4) 其中，H(x)是熵的定义，而H(y|x)是针对类别y的条件熵的定义。
(3)采用粗糙集的属性重要性概念，通过计算辅助变量对目标变量的分类能力来衡量二个变量的相关性，选择重要性更大的特征变量当辅助变量对目标(决策变量)属性重要性强时，说明二者具有强相关性；而辅助变量对决策变量属性重要性弱时，说明二者相关性弱。
其中，X为因子变量集，xi是要计算属性重要性的因子变量，γ(X，y)是粗糙集中目标变量y对因子变量集X的依赖性。
一般情况下，如果因子变量在三种方法中都较小或者有一种方法计算值为0，就应该舍弃该因子变量。
(4)主成分分析提取独立的信息量丰富的因子变量，有以下几个步骤 (a)进行主成分分析，采用方差极大(varimax)旋转双来突出每个主成分； (b)从第一主成分开始，选择特征值大于1的主成分； (c)这些选出的特征值大于1的主成分中，其因子相应的载荷(loading)，如果loading最大的主成分其loading接近0.8(0.75-0.8)，则选择该因子变量代表该主成分；否则如果最大的loading小于0.7，则舍弃该主成分； (d)如果有1个因子的loading大于0.8，则选择该因子代表该成分； (e)如果有多个因子的loading大于0.8，则选择那些GR或者SA较大的且解译性强的因子变量。
4、建立基于贝叶斯网络的风险概率的不确定性推理模型 建立灾害风险及易损性评估的贝叶斯网络模型，要经过建立贝叶斯网络模型框架及确定目标函数状态空间、连续变量离散化、网络结构的学习、参数的学习等步骤，按照以下顺序进行 (1)建立贝叶斯网络模型框架，确定目标函数的状态空间。对一贝叶斯网络B＝(BS，BP)，令随机变量y为灾害事件或者相关的灾害损失事件(既目标变量)，则y的状态空间可以定义为“灾害发生”与“灾害不发生”两种状态，或者“损失严重”与“损失不严重”两种状态，或者“损失最严重”、“损失严重”、“损失一般”与“轻微损失或无损失”四种状态，若已知了一些“证据”，既部分因子变量X的数值，以通过计算边际概率而求得y的相应状态的发生(后验)概率(也称为信任函数) 其中，为随机变量集V上的联合概率，i既为影响因子变量集X中的随机变量，yj表示目标变量y取第j种状态情况。
(2)采用了优化的离散化方法来转换连续数据，获取得到概率推理模型既贝叶斯网络模型的连续变量的状态空间，便于在提高网络性能的同时融合连续数据与离散、类别数据一起进行搜索、估算及预测。
优化的离散化是在有训练样本下的一种监督离散化算法，它是以基于熵的测度条件下的求最优离散化区间方法。离散化后的连续变量数据作为贝叶斯网络随机变量的状态空间输入，可对网络的训练及预测性能有很大提高。该算法首先设置一个很大的可离散化区间数目，然后算法适应性地训练得到离散化区间。算法的衡量标准(目标函数) 其中，impurity(k，j，i)为当从j到i的训练样本被划分成k个区间时的最小的测量值(impurity)。最优的划分区间k则是impurity(k，j，i)，N为连续变量因子值的基。测量值(impurity)是每次划分的平均类熵，其最小使得用于编码离散化区间而进行预测的字节数最小，满足了信息论中的Oecam的最小描述长度(MDL)原则。本算法通过自回归叠代而求得最优解。
令S为栅格数据中获取的样本单元集，k为最大的离散化区间数，x为要离散化的连续变量，T为离散化x的区间划分值集，IT为初始的区间划分值集，则本算法的主要步骤为 (a)根据属性x的值排序S中的实例； (b)通过计算排序后的S中连续两个不同的x值的乘积的平方根得到初始的区间划分值，保存到IT中。
(c)令i＝0，n＝|IT|； (d)若n＜＝k，则T＝IT，计算结束；否则继续； (e)令j＝n到k，以-1作为步长，递减循环进行以下操作 (e.1)T0＝T，minx为x中最小值，maxx为x中最大值； (e.2)newimpurity＝GetImpurity(1，ITj，maxA，T0)+GetImpurity(k-1，minA，ITj，T0)； (e.3)若j＝n，则minImpurity＝newimpurity；i＝j；T＝T0； 否则，进入以下步骤 (e.3.1)若newimpurity＜minImpurity，则 minImpurity＝new_impurity；i＝j；T＝T0； (e.4)T←ITi； 算法中Getlmpurity(...)为根据式(7)所定义的函数，该函数可适应性地加入新的分离值到T中或修改T中的离散分离值。
(3)搜索最优的网络结构 在将所有的连续型随机变量转换成离散变量后，既可进行最优的贝叶斯网络拓扑结构的搜索。本方法采用了贝叶斯质量记分函数及通用的模拟退火算法搜索拓扑空间而得到最优的网络拓扑结构。
贝叶斯质量记分函数用于衡量所得网络是否满足质量最优，采用公式 其中，D训练样本P(BS，D)为后验概率，vi为随机变量节点(X中因子变量或者目标变量y)，πi为节点vi的父节点集，ri为
的基，qi为
为有序的，xik为节点vi的第k个值，
为节点vi的父节点集的第j种组合，Nijk为vi＝xik与的案例数，而Nij为的案例数。
模拟退火的搜索策略是通过模拟金属的退火即硬化过程而找到最优解。退火过程从开始加热金属致其融化与金属离子自由游动开始，温度缓慢有序地降低直到金属粒子在一个高结构化的格网中成型。搜索一个优化的BN，一个潜在的网络相当一个解，s∈S，S为温度状态集，s对应了一个同网络质量测量值(式(6))有关的价值函数f(s)，该函数即称为“能量水平”。退火的过程即随着温度T的降低，算法叠代地从当前解决方案s的邻近的解(一般通过加或减去一条幅而得到)中选择新解，r，当f(r)＜f(s)时，新解被接受，r取代s进入下一轮循环；否则，r将以概率exp((f(r)-f(s))/T)被接受(取代s)。
该算法的主要步骤 输入T0-开始温度；Tend-结束温度；calc_temp为冷却函数初始化s T←T0；k←0；sbest←s 重复以下步骤 从当前解s的邻近的解空间N(s)提取备选解，r； 如果exp((f(s)-f(r))/T)＞random(0-1)则 s←r 如果f(s)＞f(sbest)，则sbest←s k←k+1 T←calc_temp(T，k) 直致T＜Tend。
返回sbest为最优解。
注f(s)为同质量评估函数相关的当前解s价值函数；冷却函数calc_temp可根据一定准则设定。
(4)估计网络参数 在训练得到网络拓扑结构之后，需要估计网络的参数。网络参数的估计假定随机变量共有τ种状态，每一种状态用mi来表示，变量的状态函数符合Dirichlet分布 第i种状态的均值 变差 Dirichlet分布是beta分布的扩展，其均值及变差与样本数相关的性质使得该分布很适合BN参数的估计。下面给出了简要算法流程。
贝叶网络参数的估计算法 每个随机变量节点，ui ui的每种父节点联合状态，πi，分配一个针对ui的τ种状态取值的Dirichlet分布D(α1，...，αi，...ατ) 对每个ui节点进行以下操作 对所有随机变量的联合观察值ui∈U进行以下操作 (a)识别ui的第j种状态的取值； (b)更新D(α1，...，αj，...ατ)至D(α1，...，αj+1，...ατ)，其中父节点与训练实例的父节点相对应。
(5)概率预测方法的选择 当贝叶斯网络模型建立的时候，可以在所建立模型基础上在一定证据支持下进行风险或易损性的概率推理。在推理方面本发明考虑了所得到网络的多变性，采用了适用于环路的信任传递算法，既Pearl的信息传递算法实现。如所得到网络是单连接树，该算法采用PolyTree的精确算法实现预测；对包括环路的网络，该算法采用为精度较高的近似算法实现预测。
5、在完全或不完全证据支持条件下的风险或易损性概率估算 当建立了风险或者易损性的贝叶斯推理模型之后，就可以在精细尺度下的完全或不完全证据支持之下进行风险或者易损性的概率估算。数据的输入为允许缺失数值的具有相同分辨率及坐标的多维栅格数据，每个栅格单元为一个实施对象，将每个栅格单元对应的属性值输入到模型之中，最后得到在整个区域上的概率风险分布情况。根据贝叶斯满足的马尔可夫性质，在证据支持下的概率估算有以下四种类型，分别视不同情况予以采用 (1)因果类型，既从原因到结果的概率估算，根据当前情况推断灾害事件产生及其损失情况的可能性； (2)诊断类型，由结果到原因的概率估算，根据损失情况推断导致损失的原因； (3)因果间类型，一个结果多个原因，可根据结果及部分原因值推断其他原因发生概率。
(4)混合类型，即以上三种类型推理的混合，一般所搜索得到的网络结构都有较为复杂的拓扑结构，涉及到混合推理方式，通过采用Pearl的信息传递算法实现推理。
6.在地图上实时动态地定位所预测的风险概率及易损性，识别区域上的风险及潜在损失分布状况，为救灾减灾的实施提供辅助决策支持信息。
具体按照以下步骤进行 (1)收集尽可能多的当前解释变量既预测因子在区域上的数值，如预测洪水灾害风险，致灾因子可采集日降雨、日平均风速及日最大风速三个变量，孕灾环境可采集高程、坡度、坡向、地质类型、植被指数(ndvi)、河流、道路7个变量，而承灾体特性可采集城镇聚集地1个自变量，总共11个自变量。当然，对一些变化缓慢的变量如高程可采用历史数据；而且数据不一定采集全面，因为贝叶斯网络具有在缺值情况下的推理功能，即使只有一个自变量如日平均降雨，亦可预测灾害损失风险水平的空间分布，这是贝叶斯网络同其他模型相比的主要特点之一； (2)对所采集的数据，参照前面第1节与第2节对数据进行必要的预处理，如栅格矢量化，栅格的再栅格化等，生成同训练网络相同的分辨率的预测因子数据(之所以要求分辨率相同是为了消除尺度不一致导致的误差)； (3)将每个栅格单元对应的预测因子变量输入网络，通过网络的推理运算得出不同风险级别对应的发生概率水平；、 (4)选定某一风险水平(如“高风险”或者“高损失”)，得到各栅格单元的产生高风险的概率的空间分布影像图，即可对此影像矢量化，形成风险水平矢量分布图。
(5)根据以上产生的风险概率分布图，可准确及时地了解当前灾害损失在空间的分布情况，对于高风险区域，可采取相应的减灾措施及时施救，减少因施救不及时产生的损失后果，最大限度的减少损失。
(6)如某点预测不很准确，可将相关的数据重新输入网络，继续对网络进行学习训练，提高网络的计算功能及结果精度。
二、下面以黑河流域的洪水灾害作为典型案例，再详细阐述本发明的实施方式简述如下 1、时空数据的采集与预处理 根据域知识及目标变量采集相关的数据，在此主要是采用黑河流域地区的洪水灾害作为实验案例，因此本研究初步采集同洪水灾害相关的影响因子变量，从致灾因子、孕灾环境及承灾体属性三方面进行采集。
主要是搜集2006年7月份遭受洪水灾害前后的数据来建立模型的，在所采集的因子变量中，致灾因子包括日降雨、日平均风速及日最大风速，与孕灾环境相关的因素包括高程、坡度、坡向、地质类型、植被指数(ndvi)、河流及道路，而承灾体相关的因素主要是城镇聚集地数据。总共采集了11个影响因子。
对这些数据进行插值、缺值分析以补充数据的完整性。将不同来源、格式、分辨率及类型(类别、离散、连续)数据综合到同样分辨率及坐标的多维栅格数据集中。具体而言，对因子变量日降雨、日平均风速及日最大风速变量由测站数据采用Kriging插值算法生成栅格，而地质类型多边形矢量数据采用扫描法生成栅格数据，河流、道路、城镇聚集地数据则采用核密度函数方法衍生栅格，其他的栅格数据则通过最近邻“重采样”转化成统一分辨率及坐标的数据集中。
2、采用核密度函数方法衍生训练样本数据 采用了核密度函数方法来衍生训练样本数据，该算法考虑了空间距离的影响，选择正态(高斯)密度函数，既式(1)；再由各观察值点的核密度函数可推出空间任意一点的影响强度，既式(2)。
当分析的变量因子为灾害或其损失强度y时，我们将灾害损失强度划分为二类，既“强”与“弱”或者严重灾害“发生”与“不发生”；而在以下的贝叶斯网络模型计算中通过计算“强”或者“发生”对应的状态概率值而求得相应的风险概率值，其中“强”或“发生”对应了较强的核密度值，由此得到灾害事件或其损失强度的核密度衍生栅格数据，作为“结果已知”的训练样本。
当分析的变量因子为影响因子变量时，采取类似的方法获得衍生的属性数据。此类因子包括河流、道路、城镇聚集地数据三类数据，其中河流及道路为线状要素，采取核密度函数法按照离开道路的远近依次划分为“难”、“中”、“易”三类，对应了灾害发生时采取避免措施的难易；而河流则按照距离河流的远近将受洪水灾害影响级别划分为“低”、“中”、“高”，对应了受洪水灾害影响的高低，越靠近河流，影响越显著；而城镇聚集地为点状数据，按照距离远近划分为“低”、“中”、“高”，可解译为“越靠近城镇，灾害的影响越高”。
3、特征选择 采用了三种方法来进行特征选择既Pearson的统计相关性系数、基于熵的GR、粗糙集中的属性重要性数值，通过对比变量因子的这三项指标，保留值相对较大者，而删除值较小或为0的因子变量；同时按照主成分分析提取独立因子的方法，从备选因子中选择较为独立的影响因子。在本案例中，坡度的GR为0，且属性重要性较小，在最后的预测变量集中去除了该变量而保留了其他10个较独立的变量作为预测因子变量。
4、贝叶斯网络建模，包括以下步骤 (1)建立贝叶斯网络模型框架，确定目标变量状态空间。目标变量y可表示灾害或其损失强度，可将损失强度划分为二类，既“强”与“弱”或者严重灾害“发生”与“不发生”；而在贝叶斯网络模型中通过计算“强”或者“发生”对应的状态概率值而求得相应水平的风险或易损性概率。
(2)采用优化的离散化算法离散化连续变量 采用了基于熵的监督的离散化算法来离散化连续的因子变量，为贝叶斯网络B＝(BS，BP)提供变量的状态空间。离散化以公式(7)作为衡量标准，设置20作为最大的离散区间数，最后通过自回归叠代找出各连续变量的离散状态空间，这些变量包括日平均降雨、日平均风速、日最大风速、高程、坡度、植被指数。在离散化的变量中，日平均降雨被离散化成13个区间，日平均风速被离散化成3个区间，日最大风速被离散化成7个区间，高程被离散化成5个区间，坡度被离散化成5个区间，植被指数被离散化成7个区间。在离散化的变量中，每个变量所产生的离散区间数小于20，说明最大的离散区间数设置合理，最后将这些离散化的区间采用整数进行标记，便于下一步的处理。
(3)采用模拟退火的搜索策略搜索网络拓扑结构 采用了模拟退火策略从数据中按照预测误差最低原则搜索网络结构，使用了贝叶斯质量计分法则。在使用该算法搜索优化的网络结构中的相关参数的设置为初始温度为10，每次叠代中温度减少因子为0.99，最大叠代次数为1000，随机数种子为1，网络满足马尔可夫的条件独立性。我们采用致灾因子、孕灾环境及承灾体特性三方面共10个离散因子变量采用本算法进行运算，得到了网络结构BS如图4。
(4)对搜索得到的网络结构BS，在假定网络节点的条件概率分布符合Dirichlet分布前提下进行参数估计，得到条件概率参数估计值集BP。
(5)对训练得到的网络进行10次10倍交叉验证，经检验得到的估计精度对高风险单元的真正率为0.867，假正率为0.19，命中率为0.689，ROC面积指数为0.917，平衡指数为0.836；而对低风险单元的真正率为0.81，假正率为0.133，命中率为0.926，ROC面积指数为0.917，平衡指数为0.836。总的预测精度为0.83，而Kappa统计量为0.6349。结果表明预测达到了一定的精度。
5、采用其他的风险概率预测模型，既
Bayes及逻辑斯特模型，对连续型数据构成的记录数据表，建立相应的
Bayes及逻辑斯特预测模型。同时对同样的数据集，进行10次10倍的交叉验证，结果表明贝叶斯网络模型行为较优，图5显示了验证结果中的受试者特征曲线(英文为Receive Operating Characteristic curve，ROC)，是一种综合性地比较预测模型性能的方法，曲线越靠近左上角的完美分类点(0，1)(该点的假正率为0，真正率为1)，则说明性能越优。我们看到BN模型的ROC曲线相对于逻辑斯特高，与
Bayes有交叉，但大半在其上，说明其预测性能较优。
6、针对同一地区，采用三个不同的模型，采集2007年7月21号的因子变量数据，其中包括了相应的10个因子变量，用以上所建立的逻辑斯特回归模型、
Bayes及所搜索得到的贝叶斯网络模型进行预测该年度7月份的洪水灾害风险概率分布。
结果表明所得到的贝叶斯网络模型具有较好的预测精度，所得结果比较如下对高风险预测而言，逻辑斯特模型的真正率为0.007，假正率为0.002，命中率为0.625，ROC面积指数为0.813，平衡指数为0.294；
Bayes模型的真正率为0.329，假正率为0.043，命中率为0.766，ROC面积指数0.813，平衡指数为0.284；前面通过模拟退火算法所得到的贝叶斯网络模型的真正率为0.421，假正率为0.102，命中率为0.638，ROC面积指数0.854，平衡指数为0.30。比较表明贝叶斯网络模型的真正率、平衡指数及ROC面积都较高，说明其能识别更多的高风险区域，体现更好的预测性能。
三个方法预测各栅格单元的ROC曲线如图6，图中可见本发明方法的ROC曲线靠近最优分类点(0，1)更近，说明其预测性能更好。
图7展现了本发明提出的贝叶斯网络模型用于预测研究区域2007年7月的洪水灾害风险概率估计结果，图中栅格单元的分辨率为0.0118精度×0.0118纬度。风险的级别从无风险(图中空白区域)，到损失发生高概率值(图中黑色区域)，颜色越黑说明损失发生的概率越高。虽然即使是较好的贝叶斯网络模型预测精度不高，但基本上正确识别了栅格单元中灾害发生或损失严重的地方(图中颜色较深对应的部分)。
总之，本发明采用精细尺度下既一定分辨率的栅格单元作为训练样本，扩展了训练数据的来源，提高了分析的细节及精细的空间尺度定位功能；在精细尺度下既基于一定分辨率的栅格的多维时空数据中量化影响因子同灾害损失风险及易损性之间的关系，设计优化的离散化方法离散化连续变量，采用模拟退火算法搜索优化的网络模型，并采用此模型进行风险概率及易损性预测，提高了风险概率及易损性估算效果，并根据结果对风险及易损性进行精细的空间尺度定位，在应用中结果可为防灾及减灾提供重要的辅助决策支持信息，对提高减灾救灾效果及国家建立智能化的突发公共事件应急预警系统具有重要的理论意义及实用价值。
权利要求
1、一种精细尺度下的动态风险及易损性预测方法，其特征在于步骤如下
(1)相关数据的采集及预处理按致灾因子、孕灾环境及承灾体三个采集预测因子X及灾害损失变量y，并进行预处理产生精细尺度下的训练样本数据集；所述的致灾因子是指直接导致灾害的可测量因子变量，所述的孕灾环境指在致灾因子存在的前提下导致灾害损失发生的自然地理环境，所述的承灾体特性指各种主要以人为主的承受灾害的对象遭受灾害所带来损失的性质；
(2)采用融入空间相关性的核密度函数方法对步骤(1)的训练样本数据集衍生出新的训练样本数据，以弥补调查数据不足，并补充完整精细尺度下的训练样本数据；
(3)预测因子的特征选择根据训练样本数据，采用三种方法衡量影响因子同目标因子的相关性，即线性相关系数、基于熵的信息获取率及粗糙集的属性重要性，结合主成分分析选择预测因子X；
(4)对步骤(3)所选择的预测因子集中的连续变量进行离散化，采用模拟退火优化算法搜索贝叶斯网络模型，估计网络的条件概率参数；
(5)利用生成的贝叶斯网络模型进行完全或不完全证据支持下的风险及易损性概率预测；
(6)在地图上实时动态地定位所预测的风险概率及易损性，识别区域上的风险及潜在损失分布状况，为救灾减灾的实施提供辅助决策支持信息。
2、根据权利要求1所述的精细尺度下的动态风险及易损性预测方法，其特征在于所述步骤(1)中的预处理方法为对致灾因子、孕灾环境及承灾体特性数据经预处理转换成统一分辨率的栅格数据集，分不同的情况采用以下步骤
a.对矢量数据预处理方法
对于点数据，采用k-最近邻方法或者Kriging方法插值算法，将点数据转换成栅格数据；对于线数据采用数字微分分析法DDA法或Bresenham方法转换成栅格数据；对于面数据采用扫描法；转换成统一分辨率的栅格数据；
b.对于栅格的预处理方法
对于栅格数据，如果坐标或目标栅格数据不一致时，采用最近邻法进行重采样，形成转换成统一分辨率的栅格数据；
将矢量数据和栅格再栅格化后即形成统一分辨率的栅格数据集。
3、根据权利要求1所述的精细尺度下的动态风险及易损性预测方法，其特征在于所述步骤(2)采用了融入空间相关性的核函数方法衍生样本点数据，弥补调查数据及历史资料的不足，补充完整精细尺度下的训练样本数据的步骤为
a.调查确定灾害事件及其损失情况；
b.采用正态核密度函数计算估计目标变量的核函数
其中的影响宽度λ可随灾害强度的不同设置不同值，d(z，Zi)为空间上任意一点z与调查数据点Zi间的欧氏距离，z可为灾害事件y或者因子变量x；
c.由各观察值点的核密度函数可推出空间任意一点的影响强度
其中，n为观察值数据点数，z为空间任意点，Zi为相应的观察值点数值，Kλ(z，Zi)为核函数。
4、根据权利要求1所述精细尺度下的动态风险及易损性预测方法，其特征在于所述步骤(3)采用三种方法衡量影响因子同目标因子的相关性具体如下
a.所述的线性相关性是采用统计学Pearson相关系数来衡量连续变量之间的线形相关性，但该方法忽略了变量之间可能存在的非线形关系及非连续变量相关关系；
b.所述的基于信息熵的分类贡献函数，即Quinlan的Information Gain IG及Information Gain Ratio GR衡量辅助变量同目标变量之间的关系，选择对分类贡献较大的变量，IG及GR的定义如下
IG(y，x)＝H(y)-H(y|x)
GR(y，x)＝(H(x)-H(y|x))/H(x)
H(y)是熵的定义，而H(y|x)是针对类别x的条件熵的定义；
c.所述的采用粗糙集属性重要性为通过计算辅助变量对目标变量的分类能力来衡量二个变量的相关性，选择重要性更大的特征变量
其中，X为因子变量集，xj是要计算属性重要性的因子变量，γ(X，y)是粗糙集中目标变量y对因子变量集X的依赖性；
d.所述的采用主成分分析为选择独立性强的代表因子，删除信息冗余因子。
5、根据权利要求1所述的精细尺度下的动态风险及易损性预测方法，其特征在于所述步骤(4)对选择的预测因子集中的连续变量进行离散化，采用模拟退火优化算法搜索贝叶斯网络模型，估计网络的条件概率参数的步骤如下
a.建立贝叶斯网络模型框架对贝叶斯网络B＝(BS，BP)，令目标变量y为灾害事件或相关灾害损失事件，则y状态空间定义为“灾害发生”与“灾害不发生”两种状况，或者“损失严重”与“损失不严重”两种状况，或者“损失最严重”、“损失严重”“损失一般”与“轻微损失或无损失”四种状况；若已知一些证据，即解释变量因子集X中的一些因子，通过计算边际概率而求得y的相应状态发生的后验概率，亦称信任函数
其中，为随机变量集V上的联合概率，xi既为影响因子集X中的随机变量，yj表示目标变量y取第j种状态情况；
b.采用优化的离散化方法转换连续类型数据，作为因子变量的状态空间值输入贝叶斯网络模型进行计算，算法采用基于熵的衡量标准
其中，impurity(k，j，i)为当从j到i的训练样本被划分成k个区间时的最小的测量值impurity，最优的划分区间k则是impurity(k，j，i)，N为连续变量因子值的基，测量值impurity是每次划分的平均类熵，其最小使得用于编码离散化区间而进行预测的字节数最小；
c.搜索最优的网络结构
采用贝叶斯质量记分函数作为标准，搜索策略是通过模拟金属的退火即硬化过程而找到网络结构的最优解，退火过程从开始加热金属致其融化与金属离子自由游动开始，温度缓慢有序地降低直到金属粒子在一个高结构化的格网中成型；一个潜在的网络相当一个解，对应一个同网络质量测量值有关的价值函数，通过退火即温度T的降低，搜索最优的网络结构；
d.估计网络参数
得到网络拓扑结构后，采用基于Dirichlet分布的经典算法，叠代地估计条件概率参数值；
e.风险概率及易损性预测
将所建立的贝叶斯网络用于一定证据支持下风险或易损性概率预测，采用适用于环路的信任传递算法即经典的Pearl信息传递算法实现，如所得到网络是单连接树，该算法则采用PolyTree的精确算法实现预测；而对包括环路的网络，该算法则采用为精度较高的近似算法实现预测。
6、根据权利要求1所述的精细尺度下的动态风险及易损性计算方法，其特征在于所述(5)利用生成的贝叶斯网络模型进行完全或不完全证据支持下的风险及易损性概率预测方法为估算中数据输入为允许缺失数值的同分辨率及坐标的多维栅格数据，每栅格单元为一个预测对象，将每单元对应属性值输入模型，通过估计精细尺度下既一定分辨率的各个栅格单元的风险概率或易损性，从而估算区域上风险概率。
全文摘要
一种精细尺度下的动态风险及易损性预测方法，属于地球信息科学领域。主要特点在精细尺度下即一定分辨率的栅格格式基础上从多源异构的时空数据搜索优化的贝叶斯网络，结合域知识改进网络，由此进行灾害风险及易损性的不确定性估算即概率估算。本发明提出用核密度方法来根据样本衍生栅格训练样本；提出优化离散化方法离散化连续变量以便为网络提供离散的状态空间输入；采用模拟退火优化算法搜索优化的网络结构；采用精确推理结合近似推理方法预测风险及易损性概率。本发明提出的方法可在精细的空间尺度下实时定位灾害风险及易损性位置，估计其风险概率的空间分布，对提高减灾救灾效果及国家建立智能化的突发公共事件应急预警系统具有重要的理论意义及实用价值。
文档编号G06N5/04GK101354757SQ20081022205
公开日2009年1月28日 申请日期2008年9月8日 优先权日2008年9月8日
发明者李连发, 梁金龙 申请人:中国科学院地理科学与资源研究所