CN103366092B - 基于状态转移的工程风险监测系统及方法 - Google Patents

Abstract

本发明揭示了一种基于状态转移的工程风险监测系统及方法，所述系统包括原始数据处理模块、状态转移可视化表达模块、图形特征提取模块、自动学习模块、图像自动辨识模块。状态转移可视化表达模块根据经过原始数据处理模块处理的数据生成工程中各状态的状态转移图；图形特征提取模块从上述状态转移图中提取状态转移图的设定特征；自动学习模块根据提取的图形特征计算每一特征值与历史平均值的比值并设定阈值；图像自动辨识模块根据图形特征提取模块提取状态转移图的设定特征及自动学习模块设定的阈值判断工程是否处于风险状态。本发明能够针对性地解决地下工程灾害的不确定性和突发性问题，且以可视化方式自动监测工程风险，具有直观性和智能性。

Description

基于状态转移的工程风险监测系统及方法

技术领域

本发明属于计算机及测量技术领域，涉及一种工程风险监测系统，尤其涉及一种基于状态转移的工程风险监测系统；同时，本发明还涉及一种基于状态转移的工程风险监测方法。

背景技术

由于地下工程结构和施工方法的特殊性，特别是施工开挖过程中周边环境（介质）的不确定性和复杂性，建立精确的数学模型对其进行灾害分析与预警非常困难。常见的地下工程安全分析方法包含以下三类：

（1）定性为主的风险评估

定性为主（或结合部分定量计算）的风险分析技术方法常用在风险评估领域，确定风险等级。包括专家调查法、过程危险分析、危险与可操作性、故障模式分析法、层次分析法、故障树/事件树分析法和模糊推理等。

（2）以统计和概率为基础的定量风险分析

以数理统计和概率分析为基础的定量分析方法，已成为当前风险定量分析最主要的方法，并在一些动态评价场合应用。主要方法有决策树、蒙特卡罗法、贝叶斯法及控制区间和记忆模型等。

（3）基于时序数据的预测技术。

随着监测技术发展，以时序数据为基础的安全性分析也成为一个重要的方向，包括马尔可夫链、人工神经网络、灰色理论、非线性自回归方法和多步预报模型等。

但上述三类地下工程风险分析方法，不能针对性地解决地下工程灾害的不确定性和突发性问题，无法捕捉未知灾害前兆。对于采样频率不同、尺度不同且高维海量的信息的利用上述方法亦远远不够。

有鉴于此，如今迫切需要设计一种新的工程安全分析方法，以克服现有方法的上述缺陷。

发明内容

本发明所要解决的技术问题是：提供一种基于状态转移的工程风险监测系统，能够针对性地解决地下工程灾害的不确定性和突发性问题，而且以可视化的方式自动监测工程风险，具有直观性和智能性。

此外，本发明还提供一种基于状态转移的工程风险监测方法，能够针对性地解决地下工程灾害的不确定性和突发性问题，而且以可视化的方式自动监测工程风险，具有直观性和智能性。

为解决上述技术问题，本发明采用如下技术方案：

一种基于状态转移的工程风险监测系统，所述系统包括：

原始数据处理模块，用以获取原始数据，并将原始数据经过无量纲化处理；

状态转移可视化表达模块，与所述原始数据处理模块连接，用以根据经过原始数据处理模块处理的数据，生成工程中各状态的状态转移图；

图形特征提取模块，与所述状态转移可视化表达模块连接，用以获取所述状态转移可视化表达模块生成的状态转移图，并从中提取状态转移图的设定特征；

自动学习模块，与所述图像特征提取模块连接，用以获取图像特征提取模块提取的图形特征值，计算每一特征值与历史平均值的比值，并设定阈值；

图像自动辨识模块，与所述图形特征提取模块、自动学习模块连接，用以根据所述图形特征提取模块提取状态转移图的设定特征以及自动学习模块设定的阈值判断工程是否处于风险状态。

作为本发明的一种优选方案，所述状态转移可视化表达模块采用同心半圆坐标系建立表达多维平面坐标，过程如下：

N维欧几里得空间R^N同心半圆坐标系由N个等距同心半圆坐标轴组成。每一个半圆为一条坐标轴，每条轴线都对应于一个多维数据集的一个分量，从圆心开始，从内到外依次编号：d₀，d₁，…，d_i，…，d_N,单条周轴线圆心角范围为0～180°，按逆时针方向递增；

假定一个N维时序数据的在k时刻的第i个分量为通过公式1和公式2函数转换映射为第i条轴线上的一个点从圆心开始，依次连接形成一条折线，用于表达一组独立的n维数据；

${\mathrm{\θ}}_i^k=\frac{s_i^k}{{\max }_i-{\min }_i}$ 公式(1)

$\left\{\begin{array}{c}{P_i}^k.x=i*u*\cos \left({\mathrm{\θ}}_i^k\right)\\ {P_i}^k.y=i*u*\sin \left({\mathrm{\θ}}_i^k\right)\end{array}\right.$ 公式(2)

其中，u表示两个相邻轴半径的差值。

作为本发明的一种优选方案，所述原始数据处理模块将原始数据经过无量纲化处理；

无量纲化处理通过数学变换来消除原始变量量纲影响的方法；将数据全部缩放到0～π之间来消除原始变量的量纲影响；具体过程如下：

假定在k时刻，第i个分量通过公式1将其缩放限制在0～π之间，并映射到坐标轴点；

${\mathrm{\θ}}_i^k=\frac{s_i^k}{{\max }_i-{\min }_i}$ 公式(1)

其中， $s_i^k\∈[{\min }_i,{\max }_i],{\min }_i=\min \{s_i^k,0\≤k\≤|S\left|\right\},{\max }_i=\max \{s_i^k,0\≤k\≤|S\left|\right\}.$

作为本发明的一种优选方案，所述状态转移可视化表达模块是将数据进行向空间转换、可视化表达的过程；具体包括：

同心半圆坐标系建立单元，用以建立同心半圆坐标系，将预处理的数据根据公式映射到该平面坐标系中；

状态折线绘制单元，用以绘制同一个同心半圆坐标轴下所有状态折线；

状态转移图生成单元，用直线将相邻两根状态折线的末端相连形成一个封闭区域并进行填色，生成一幅状态转移图；

系列状态转移图绘制单元，用以根据预设的周期和递进参数绘制系列状态转移图。

作为本发明的一种优选方案，所述自动学习模块根据图形特征提取模块，提取状态转移图的特征变量值，计算当前特征值与平均值比值的过程；平均值分别按历史数据中数值大小前5%、10%、20%、50%的数值计算；

所述图像自动辨识模块用以实现状态转移图风险识别的主要过程；具体包括：

（1）根据图形特征提取模块，获取系列状态转移图的特征变量值并绘制相应的指标曲线；

（2）选择是否进入自动学习模块；若进入自动学习模块，则获取特征变量当前值与历史平均值的比值；若该比值在设定阈值外即判断为异常；

若不进入自动学习模块，则计算特征变量当前值与前一值的比值；若该比值在设定阈值外即判断为异常；

（3）比值超过阈值越大，判断状态转移图越为异常，此时系统越有可能处于异常状态。

作为本发明的一种优选方案，所述系统通过一组时序数据为S表示工程中的各状态变量数据，那么S^k就表示在k时刻工程的状态，并且每一个状态都由n个代表工程状态的内部变量所组成；

S={S⁰，S¹，...，S^k-1，S^k，S^k+1，...}；

$S^k=\{s_1^k,s_2^k,...,s_i^k,..,s_n^k\};$

其中， $s_i^k\∈[{\min }_i,{\max }_i],{\min }_i=\min \{s_i^k,0\≤k\≤|S\left|\right\},{\max }_i=\max \{s_i^k,0\≤k\≤|S\left|\right\};$

假定某工程有n个可测分变量，那么在k时刻，其状态变量S^k表示为：

$S^k=\{s_1^k,s_2^k,.....,s_n^k\}$

定义第i个可测分变量值的量程范围分别为min_i和max_i，那么有：

$s_i^k\∈[{\min }_i,{\max }_i]$

所述状态转移可视化表达模块包括第一状态分量映射模块，用以将第k时刻第i个分量通过公式1将其缩放限制在0～π之间，并映射到坐标轴点；

${\mathrm{\θ}}_i^k=\frac{s_i^k}{{\max }_i-{\min }_i}$ 公式(1)

所述状态转移可视化表达模块包括第二状态分量映射模块，用以选择平面上一定点O，以O为原点坐标，建立一个同心半圆坐标系；对于任意一个状态变量的分量通过公式2将其映射为平面上的一点

$\left\{\begin{array}{c}{P_i}^k.x=i*u*\cos \left({\mathrm{\θ}}_i^k\right)\\ {P_i}^k.y=i*u*\sin \left({\mathrm{\θ}}_i^k\right)\end{array}\right.$ 公式(2)

其中，u表示两个相邻轴半径的差值；

对于k时刻的状态变量S^k，通过公式3，映射为平面上的点并将其依次连接，得到一条折线L_k来反映k时刻的状态变量；而系统的每一个状态则可以通过这些唯一的折线来表示，因此有：

$L^k=\{O,P_1^k,P_2^k,P_3^k,...,{P_i}^k,...P_n^k\}1\&le;i<n$ 公式(3)

在平面直角坐标系中，首先画出n个等距同心半圆，其为N维欧几里得空间R^N的同心半圆坐标轴；即使这些坐标轴的长度并不一样，但是每个半圆都有相同的角度范围，即0～π；每一个半圆代表一根坐标轴，且方向为逆时针方向；坐标轴d_i代表了第i维的数据，并通过公示将数据点缩放限制在0～π之间；

所述状态转移图用来描述一段时间内系统状态转移的过程；从k₁时刻到k₂时刻，系统的状态顺序集为{S^k1，S^k1+1，……S^k2-1，S^k2}，对应的折线集为{L^k1，L^k1+1，……L^k2-1，L^k2}；

所述状态转移可视化表达模块还包括：

状态折线绘制模块，用以在同一个同心半圆坐标轴下画出所有状态折线；

封闭区域形成模块，对于相邻两根状态折线，用直线将折线的末端相连形成一个封闭区域，A^k是由折线L^k和折线L^k+1以及他们的末端连线围成的封闭区域；

填色模块，用以对封闭区域A¹,A²,……,A^k-1进行一次填色，当区域相互重叠的时候，重叠区域的颜色就会加深，颜色的深度取决于这些区域重叠的次数。

作为本发明的一种优选方案，所述状态转移可视化表达模块用以生成工程中一系列状态的状态转移图；一组有序的时间序列图根据周期T以及递进D绘制产生系列轨迹图，0≤D≤T；其中第一张图所选取的原始数据系列为{X₁,X₂,X₃,...,X_T}，第二张图所选取的原始数据系列为{X_D+1,X_D+2,...,X_D+T}，那么第N张图所选取的原始数据系列为{X_D*(N-1)+1,X_D*(D-1)+2,...,X_D*(N-1)+T}不同的轨迹图代表了系不同时期的属性。

作为本发明的一种优选方案，所述图形特征提取模块提取状态转移图的设定特征包括：状态转移图的质心与质心间距、图形熵、图形面积、单位面积熵变率中的一个或多个；

质心为物质系统上被认为质量集中于此的一个假想点，质心的位置矢量是质点组中各个质点的位置矢量根据其对应质量加权平均之后的平均矢量；质心的变化反映图形特征的变化，对于每一幅状态转移图的质心用公式4表示：

$\stackrel{\&OverBar;}{X}=\frac{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N(P_{V(i.j)}*i)}{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N{P}_{V(i.j)}}$ 公式(4)

$\stackrel{\&OverBar;}{Y}=\frac{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N(P_{V(i.j)}*j)}{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N{P}_{V(i.j)}}$

其中，M*N表示图像在X方向M个像素点，Y方向N个像素点；Pv(i,j)表示一幅图像中像素点(i,j)的灰度值；

由公式4可知，图像质心与图像中像素点的灰度值明显相关；为了更易发现状态的变化，通过缩小灰度值范围求出同一副图的不同质心，而质心间距用于衡量图形质心的变化，定义如公式5所示；

$L_{i,i+1}=\sqrt{{(X_i-X_{i+1})}^2+{(Y_i-Y_{i+1})}^2}$ 公式(5)

其中，第i幅图的质心坐标为(X_i,Y_i)；第i+1幅图的质心坐标为(X_i+1,Y_i+1)；

图形熵是用来估计不同图像之间信息的变化量的指标，同时也反映了不同状态下系统的能量变化；系统能量越均匀，状态无序，熵值越大；系统能量不均匀，状态有序，熵值越小；公式6将图形熵定义为：

$\left\{\begin{array}{c}E={\mathrm{\&Sigma;}}_{i=0}^{G}d\left(i\right){\log }_{2}d\left(i\right)\\ d\left(i\right)=p\left(i\right)/N\end{array}\right.$ 公式(6)

其中，E表示图形熵；G表示8进制图像直方图中灰度值，取值范围为0-255；d(i)表示标准化后每一个灰度值出现的频率；N表示整幅图像的总像素点个数；p(i)表示灰度值为i的像素点个数；

图形面积变量包含一个值，统计的是图像像素点的个数，当图像的面积发生变化时，图像的像素点个数也会发生变化，像素点变化越大，代表着图像的面积特征变化越明显；定义如公式7：

$\left\{\begin{array}{c}A={\mathrm{\&Sigma;}}_{i=1}^X{\mathrm{\&Sigma;}}_{j=1}^{Y}F\left(\mathrm{Pv}(i,j)\right)\\ F\left(a\right)=\left\{\begin{array}{cc}1& a\&NotEqual;255\\ 0& a=255\end{array}\right.\end{array}\right.$ 公式(7)

其中，X表示一幅图像中水平方向的像素点个数；Y表示一幅图像中垂直方向的像素点个数；P_V(i,j)表示一幅图像中像素点(i,j)的灰度值；F()表示一幅图像中的某一个像素点是否为白色，若a=255为白色；a表示像素点的灰度值，取值范围为0到255；

对于两幅相邻的状态转移图，如果他们的图形熵和图形面积的值分别是H_i，S_i，H_i-1和S_i-1，定义单位面积熵变率ECRA如公式8所示：

${\mathrm{ECRA}}_i={\mathrm{\&Delta;HA}}_i=\frac{{\mathrm{\&Delta;H}}_i}{A_i}=\frac{\left|\frac{H_i-H_{i-1}}{H_i}\right|}{A_i}=\frac{|H_i-H_{i-1}|}{H_{i-1}{A}_i}$ 公式(8)

其中，ΔH_i表示第i幅状态转移图的熵变值，ΔHA_i表示第i幅状态转移图的单位面积熵变值；

所述图形特征提取模块根据上述定义分别计算每幅状态转移图的图形特征变量值。

一种基于状态转移的工程风险监测方法，所述方法包括如下步骤：

原始数据处理步骤；获取原始数据，并将原始数据经过无量纲化处理；

状态转移可视化表达步骤；根据经过原始数据处理模块处理的数据生成工程中各状态的状态转移图；

图形特征提取步骤；获取所述状态转移可视化表达模块生成的状态转移图，并从中提取状态转移图的设定特征；

自动学习步骤；获取图像特征提取模块提取的图形特征值，计算每一特征值与历史平均值的比值，并设定阈值；

图像自动辨识步骤；根据所述图形特征提取模块提取状态转移图的设定特征以及自动学习模块设定的阈值判断工程是否处于风险状态。

作为本发明的一种优选方案，所述原始数据处理步骤中，将原始数据经过无量纲化处理；

无量纲化处理通过数学变换来消除原始变量量纲影响的方法；将数据全部缩放到0～π之间来消除原始变量的量纲影响；具体过程如下：

假定在k时刻，第i个分量通过公式1将其缩放限制在0～π之间，并映射到坐标轴点；

${\mathrm{\&theta;}}_i^k=\frac{s_i^k}{{\max }_i-{\min }_i}$ 公式(1)

其中， $s_i^k\&Element;[{\min }_i,{\max }_i],{\min }_i=\min \{s_i^k,0\&le;k\&le;|S\left|\right\},{\max }_i=\max \{s_i^k,0\&le;k\&le;|S\left|\right\}.$

作为本发明的一种优选方案，所述状态转移可视化表达步骤中，将数据进行向空间转换、可视化表达的过程；具体包括如下步骤：

建立同心半圆坐标系，将预处理的数据根据公式映射到该平面坐标系中；

绘制同一个同心半圆坐标轴下所有状态折线；

用直线将相邻两根状态折线的末端相连形成一个封闭区域并进行填色，生成一幅状态转移图；

根据预设的周期和递进参数绘制系列状态转移图。

作为本发明的一种优选方案，所述方法自动学习步骤中，自动学习模块根据图形特征提取步骤提取状态转移图的特征变量值，计算当前值与平均值比值的过程；平均值分别按历史数据中数值大小前5%、10%、20%、50%的数值计算；

所述图像自动辨识步骤实现状态转移图风险识别的主要过程；具体包括如下步骤：

（1）根据图形特征提取模块获取系列状态转移图的特征变量值并绘制相应的指标曲线；

（2）选择是否进入自动学习模块；若进入自动学习模块，则获取特征变量当前值与历史平均值的比值；若该比值在设定阈值外即判断为异常；

若不进入自动学习模块，则计算特征变量当前值与前一值的比值；若该比值在设定阈值外即判断为异常；

（3）比值超过阈值越大，判断状态转移图越为异常，此时系统越有可能处于异常状态。

作为本发明的一种优选方案，所述方法通过一组时序数据为S表示工程中的各状态变量数据，那么S^k就表示在k时刻工程的状态，并且每一个状态都由n个代表工程状态的内部变量所组成；

S={S⁰,S¹,...,S^k-1,S^k,S^k+1,...}；

$S^k=\{s_1^k,s_2^k,...,s_i^k,..,s_n^k\};$

其中， $s_i^k\&Element;[{\min }_i,{\max }_i],{\min }_i=\min \{s_i^k,0\&le;k\&le;|S\left|\right\},{\max }_i=\max \{s_i^k,0\&le;k\&le;|S\left|\right\};$

假定某工程有n个可测分变量，那么在k时刻，其状态变量S^k表示为：

$S^k=\{s_1^k,s_2^k,.....,s_n^k\}$

定义第i个可测分变量值的量程范围分别为min_i和max_i，那么有：

$s_i^k\&Element;[{\min }_i,{\max }_i]$

所述状态转移可视化表达步骤包括第一状态分量映射步骤，将第k时刻第i个分量通过公式1将其缩放限制在0～π之间，并映射到坐标轴点；

${\mathrm{\&theta;}}_i^k=\frac{s_i^k}{{\max }_i-{\min }_i}$ 公式(1)

所述状态转移可视化表达步骤包括第二状态分量映射步骤，选择平面上一定点O，以O为原点坐标，建立一个同心半圆坐标系；对于任意一个状态变量的分量通过公式2将其映射为平面上的一点

$\left\{\begin{array}{c}{P_i}^k.x=i*u*\cos \left({\mathrm{\&theta;}}_i^k\right)\\ {P_i}^k.y=i*u*\sin \left({\mathrm{\&theta;}}_i^k\right)\end{array}\right.$ 公式(2)

其中，u表示两个相邻轴半径的差值；

对于k时刻的状态变量S^k，通过公式3，映射为平面上的点并将其依次连接，得到一条折线L_k来反映k时刻的状态变量；而系统的每一个状态则可以通过这些唯一的折线来表示，因此有：

$L^k=\{O,P_1^k,P_2^k,P_3^k,...,{P_i}^k,...P_n^k\}1\&le;i<n$ 公式(3)

在平面直角坐标系中，首先画出n个等距同心半圆，其为N维欧几里得空间R^N的同心半圆坐标轴；即使这些坐标轴的长度并不一样，但是每个半圆都有相同的角度范围，即0～π；每一个半圆代表一根坐标轴，且方向为逆时针方向；坐标轴d_i代表了第i维的数据，并通过公示将数据点缩放限制在0～π之间；

所述状态转移图用来描述一段时间内系统状态转移的过程；从k₁时刻到k₂时刻，系统的状态顺序集为{S^k1，S^k1+1，……S^k2-1，S^k2}，对应的折线集为{L^k1，L^k1+1，……L^k2-1，L^k2}；

所述状态转移可视化表达步骤还包括：

状态折线绘制步骤，在同一个同心半圆坐标轴下画出所有状态折线；

封闭区域形成步骤，对于相邻两根状态折线，用直线将折线的末端相连形成一个封闭区域，A^k是由折线L^k和折线L^k+1以及他们的末端连线围成的封闭区域；

填色步骤，对封闭区域A¹,A²,……,A^k-1进行一次填色，当区域相互重叠的时候，重叠区域的颜色就会加深，颜色的深度取决于这些区域重叠的次数。

作为本发明的一种优选方案，所述状态转移可视化表达步骤生成工程中一系列状态的状态转移图；一组有序的时间序列图根据周期T以及递进D绘制产生系列轨迹图，0≤D≤T；其中第一张图所选取的原始数据系列为{X₁,X₂,X₃,...X_T}，第二张图所选取的原始数据系列为{X_D+1,X_D+2,...,X_D+T}，那么第N张图所选取的原始数据系列为{X_D*(N-1)+1,X_D*(D-1)+2,...,X_D*(N-1)+T}，不同的轨迹图代表了系统不同时期的属性。

作为本发明的一种优选方案，所述图形特征提取步骤提取状态转移图的设定特征包括：状态转移图的质心与质心间距、图形熵、图形面积、单位面积熵变率中的一个或多个；

质心为物质系统上被认为质量集中于此的一个假想点，质心的位置矢量是质点组中各个质点的位置矢量根据其对应质量加权平均之后的平均矢量；质心的变化反映图形特征的变化，对于每一幅状态转移图的质心用公式4表示：

$\stackrel{\&OverBar;}{X}=\frac{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N(P_{V(i.j)}*i)}{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N{P}_{V(i.j)}}$ 公式(4)

$\stackrel{\&OverBar;}{Y}=\frac{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N(P_{V(i.j)}*j)}{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N{P}_{V(i.j)}}$

其中，M*N表示图像在X方向M个像素点，Y方向N个像素点；Pv(i,j)表示一幅图像中像素点(i,j)的灰度值；

由公式4可知，图像质心与图像中像素点的灰度值明显相关；为了更易发现状态的变化，通过缩小灰度值范围求出同一副图的不同质心，而质心间距用于衡量图形质心的变化，定义如公式5所示；

$L_{i,i+1}=\sqrt{{(X_i-X_{i+1})}^2+{(Y_i-Y_{i+1})}^2}$ 公式(5)

其中，第i幅图的质心坐标为(X_i,Y_i)；第i+1幅图的质心坐标为(X_i+1,Y_i+1)；

图形熵是用来估计不同图像之间信息的变化量的指标，同时也反映了不同状态下系统的能量变化；系统能量越均匀，状态无序，熵值越大；系统能量不均匀，状态有序，熵值越小；公式6将图形熵定义为：

$\left\{\begin{array}{c}E={\mathrm{\&Sigma;}}_{i=0}^{G}d\left(i\right){\log }_{2}d\left(i\right)\\ d\left(i\right)=p\left(i\right)/N\end{array}\right.$ 公式(6)

其中，E表示图形熵；G表示8进制图像直方图中灰度值，取值范围为0-255；d(i)表示标准化后每一个灰度值出现的频率；N表示整幅图像的总像素点个数；p(i)表示灰度值为i的像素点个数；

图形面积变量包含一个值，统计的是图像像素点的个数，当图像的面积发生变化时，图像的像素点个数也会发生变化，像素点变化越大，代表着图像的面积特征变化越明显；定义如公式7：

$\left\{\begin{array}{c}A={\mathrm{\&Sigma;}}_{i=1}^X{\mathrm{\&Sigma;}}_{j=1}^{Y}F\left(\mathrm{Pv}(i,j)\right)\\ F\left(a\right)=\left\{\begin{array}{cc}1& a\&NotEqual;255\\ 0& a=255\end{array}\right.\end{array}\right.$ 公式(7)

其中，X表示一幅图像中水平方向的像素点个数；Y表示一幅图像中垂直方向的像素点个数；P_V(i,j)表示一幅图像中像素点(i,j)的灰度值；F()表示一幅图像中的某一个像素点是否为白色，若a=255为白色；a表示像素点的灰度值，取值范围为0到255；

对于两幅相邻的状态转移图，如果他们的图形熵和图形面积的值分别是H_i，A_i，H_i-1和A_i-1，定义单位面积熵变率ECRA如公式8所示：

${\mathrm{ECRA}}_i={\mathrm{\&Delta;HA}}_i=\frac{{\mathrm{\&Delta;H}}_i}{A_i}=\frac{\left|\frac{H_i-H_{i-1}}{H_i}\right|}{A_i}=\frac{|H_i-H_{i-1}|}{H_{i-1}{A}_i}$ 公式(8)

其中，ΔH_i表示第i幅状态转移图的熵变值，ΔHA_i表示第i幅状态转移图的单位面积熵变值；

所述方法还包括：图形特征提取模块根据上述定义分别计算每幅状态转移图的图形特征变量值。

本发明的有益效果在于：本发明提出的基于状态转移的工程风险监测系统及方法，借鉴空间重构技术和可视化数据分析技术在表达复杂系统规律和特征的优势，从状态变化的角度分析系统的异常变化，从而发现工程的潜在风险。该方法与一般地下工程风险监测方法相比，不仅能够针对性地解决地下工程灾害的不确定性和突发性问题，而且以可视化的方式自动监测工程风险，具有直观性和智能性。并且本发明也隐含地给出了一种突发事件风险监测的一般方法。

附图说明

图1为本发明基于状态转移的工程风险监测方法的流程图。

图2为同心半圆坐标轴的示意图。

图3为画出不同状态下的折线的示意图。

图4为对状态转移区域进行填色的示意图。

图5为长江隧道状态转移图。

图6为长江隧道灾变前兆面积与图形熵信号示意图。

图7为长江隧道灾变前兆ECRA信号示意图。

具体实施方式

下面结合附图详细说明本发明的优选实施例。

实施例一

请参阅图1，本发明揭示了一种基于状态转移的工程风险监测系统，所述系统包括：原始数据处理模块、状态转移可视化表达模块、图形特征提取模块、自动学习模块、图像自动辨识模块。

原始数据处理模块用以获取原始数据，并将原始数据经过无量纲化处理。状态转移可视化表达模块与所述原始数据处理模块连接，用以根据经过原始数据处理模块处理的数据，生成工程中各状态的状态转移图。图形特征提取模块与所述状态转移可视化表达模块连接，用以获取所述状态转移可视化表达模块生成的状态转移图，并从中提取状态转移图的设定特征。自动学习模块与所述图像特征提取模块连接，用以获取图像特征提取模块提取的图形特征值，计算每一特征值与历史平均值的比值，并设定阈值。图像自动辨识模块与所述图形特征提取模块、自动学习模块连接，用以根据所述图形特征提取模块提取状态转移图的设定特征以及自动学习模块设定的阈值判断工程是否处于风险状态。

所述原始数据处理模块的无量纲化处理通过数学变换来消除原始变量量纲影响的方法；将数据全部缩放到0～π之间来消除原始变量的量纲影响。具体过程如下：

假定在k时刻，第i个分量通过公式1将其缩放限制在0～π之间，并映射到坐标轴点；

${\mathrm{\&theta;}}_i^k=\frac{s_i^k}{{\max }_i-{\min }_i}$ 公式(1)

其中， $s_i^k\&Element;[{\min }_i,{\max }_i],{\min }_i=\min \{s_i^k,0\&le;k\&le;|S\left|\right\},{\max }_i=\max \{s_i^k,0\&le;k\&le;|S\left|\right\}.$

所述状态转移可视化表达模块是将数据进行向空间转换、可视化表达的过程；具体包括：同心半圆坐标系建立单元、状态折线绘制单元、状态转移图生成单元、系列状态转移图绘制单元。同心半圆坐标系建立单元用以建立同心半圆坐标系，将预处理的数据根据公式映射到该平面坐标系中；状态折线绘制单元用以绘制同一个同心半圆坐标轴下所有状态折线；状态转移图生成单元，用直线将相邻两根状态折线的末端相连形成一个封闭区域并进行填色，生成一幅状态转移图；系列状态转移图绘制单元用以根据预设的周期和递进参数绘制系列状态转移图。

所述自动学习模块根据图形特征提取模块，提取状态转移图的特征变量值，计算当前值与平均值比值的过程；平均值可以为一个或多个，本实施例中，平均值分别按历史数据中数值大小前5%、10%、20%、50%的数值计算，相应地，阈值也有多个数值。优选地，历史数据为设定时间内正常状态下的历史数据。

所述图像自动辨识模块用以实现状态转移图风险识别的主要过程；具体包括：

（1）根据图形特征提取模块，获取系列状态转移图的特征变量值并绘制相应的指标曲线；

（2）选择是否进入自动学习模块；若进入自动学习模块，则获取特征变量当前值与历史平均值的比值；若该比值在设定阈值外即判断为异常；若不进入自动学习模块，则计算特征变量当前值与前一值的比值；若该比值在设定阈值外即判断为异常；

（3）比值超过阈值越大，判断状态转移图越为异常，此时系统越有可能处于异常状态。

具体地，所述系统首先通过一组时序数据为S表示工程中的各状态变量数据，那么S^k就表示在k时刻工程的状态，并且每一个状态都由n个代表工程状态的内部变量所组成；

S={S⁰,S¹,...,S^k-1,S^k,S^k+1,...}；

$S^k=\{s_1^k,s_2^k,...,s_i^k,..,s_n^k\};$

其中， $s_i^k\&Element;[{\min }_i,{\max }_i],{\min }_i=\min \{s_i^k,0\&le;k\&le;|S\left|\right\},{\max }_i=\max \{s_i^k,0\&le;k\&le;|S\left|\right\}.$

假定某工程有n个可测分变量，那么在k时刻，其状态变量S^k表示为：

$S^k=\{s_1^k,s_2^k,.....,s_n^k\}$

定义第i个可测分变量值的量程范围分别为min_i和max_i，那么有：

$s_i^k\&Element;[{\min }_i,{\max }_i]$

所述状态转移可视化表达模块包括第一状态分量映射模块，用以将第k时刻第i个分量通过公式1将其缩放限制在0～π之间，并映射到坐标轴点；

${\mathrm{\&theta;}}_i^k=\frac{s_i^k}{{\max }_i-{\min }_i}$ 公式(1)

所述状态转移可视化表达模块包括第二状态分量映射模块，用以选择平面上一定点O，以O为原点坐标，建立一个同心半圆坐标系；对于任意一个状态变量的分量通过公式2将其映射为平面上的一点

$\left\{\begin{array}{c}{P_i}^k.x=i*u*\cos \left({\mathrm{\&theta;}}_i^k\right)\\ {P_i}^k.y=i*u*\sin \left({\mathrm{\&theta;}}_i^k\right)\end{array}\right.$ 公式(2)

其中，u表示两个相邻轴半径的差值；

对于k时刻的状态变量S^k，通过公式3，映射为平面上的点并将其依次连接，得到一条折线L_k来反映k时刻的状态变量；而系统的每一个状态则可以通过这些唯一的折线来表示，因此有：

$L^k=\{O,P_1^k,P_2^k,P_3^k,...,{P_i}^k,...P_n^k\}1\&le;i<n$ 公式(3)

如图2所示，在平面直角坐标系中，首先画出n个等距同心半圆，其为N维欧几里得空间R^N的同心半圆坐标轴；即使这些坐标轴的长度并不一样，但是每个半圆都有相同的角度范围，即0～π；每一个半圆代表一根坐标轴，且方向为逆时针方向；坐标轴d_i代表了第i维的数据，并通过公示将数据点缩放限制在0～π之间。

所述状态转移图用来描述一段时间内系统状态转移的过程；从k₁时刻到k₂时刻，系统的状态顺序集为{S^k1，S^k1+1，……S^k2-1，S^k2}，对应的折线集为{L^k1，L^k1+1，……L^k2-1，L^k2}；

所述状态转移可视化表达模块还包括：状态折线绘制模块、封闭区域形成模块、填色模块。

状态折线绘制模块用以在同一个同心半圆坐标轴下画出所有状态折线；

对于相邻两根状态折线，封闭区域形成模块用直线将折线的末端相连形成一个封闭区域，A^k是由折线L^k和折线L^k+1以及他们的末端连线围成的封闭区域；

填色模块用以对封闭区域A¹,A²,……,A^k-1进行一次填色，当区域相互重叠的时候，重叠区域的颜色就会加深，颜色的深度取决于这些区域重叠的次数。

通常，工程中有一系列的状态；所述状态转移可视化表达模块对应地需要生成工程中一系列状态的状态转移图。一组有序的时间序列图根据周期T以及递进D绘制产生系列轨迹图，0≤D≤T；其中第一张图所选取的原始数据系列为{X₁,X₂,X₃,...,X_T}，第二张图所选取的原始数据系列为{X_D+1,X_D+2,...,X_D+T}，那么第N张图所选取的原始数据系列为{X_D*(N-1)+1,X_D*(D-1)+2,...,X_D*(N-1)+T}，不同的轨迹图代表了系统不同时期的属性。

所述图形特征提取模块提取状态转移图的设定特征包括：状态转移图的质心与质心间距、图形熵、图形面积、单位面积熵变率中的一个或多个。

质心为物质系统上被认为质量集中于此的一个假想点，质心的位置矢量是质点组中各个质点的位置矢量根据其对应质量加权平均之后的平均矢量；质心的变化反映图形特征的变化，对于每一幅状态转移图的质心用公式4表示：

$\stackrel{\&OverBar;}{X}=\frac{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N(P_{V(i.j)}*i)}{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N{P}_{V(i.j)}}$ 公式(4)

$\stackrel{\&OverBar;}{Y}=\frac{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N(P_{V(i.j)}*j)}{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N{P}_{V(i.j)}}$

其中，M*N表示图像在X方向M个像素点，Y方向N个像素点；Pv(i,j)表示一幅图像中像素点(i,j)的灰度值。

由公式4可知，图像质心与图像中像素点的灰度值明显相关；为了更易发现状态的变化，通过缩小灰度值范围求出同一副图的不同质心，而质心间距用于衡量图形质心的变化，定义如公式5所示；

$L_{i,i+1}=\sqrt{{(X_i-X_{i+1})}^2+{(Y_i-Y_{i+1})}^2}$ 公式(5)

其中，第i幅图的质心坐标为(X_i,Y_i)；第i+1幅图的质心坐标为(X_i+1,Y_i+1)。

图形熵是用来估计不同图像之间信息的变化量的指标，同时也反映了不同状态下系统的能量变化；系统能量越均匀，状态无序，熵值越大；系统能量不均匀，状态有序，熵值越小；公式6将图形熵定义为：

$\left\{\begin{array}{c}E={\mathrm{\&Sigma;}}_{i=0}^{G}d\left(i\right){\log }_{2}d\left(i\right)\\ d\left(i\right)=p\left(i\right)/N\end{array}\right.$ 公式(6)

其中，E表示图形熵；G表示8进制图像直方图中灰度值，取值范围为0-255；d(i)表示标准化后每一个灰度值出现的频率；N表示整幅图像的总像素点个数；p(i)表示灰度值为i的像素点个数。

图形面积变量包含一个值，统计的是图像像素点的个数，当图像的面积发生变化时，图像的像素点个数也会发生变化，像素点变化越大，代表着图像的面积特征变化越明显；定义如公式7：

$\left\{\begin{array}{c}A={\mathrm{\&Sigma;}}_{i=1}^X{\mathrm{\&Sigma;}}_{j=1}^{Y}F\left(\mathrm{Pv}(i,j)\right)\\ F\left(a\right)=\left\{\begin{array}{cc}1& a\&NotEqual;255\\ 0& a=255\end{array}\right.\end{array}\right.$ 公式(7)

其中，X表示一幅图像中水平方向的像素点个数；Y表示一幅图像中垂直方向的像素点个数；P_V(i,j)表示一幅图像中像素点(i,j)的灰度值；F()表示一幅图像中的某一个像素点是否为白色，若a=255为白色；a表示像素点的灰度值，取值范围为0到255；

对于两幅相邻的状态转移图，如果他们的图形熵和图形面积的值分别是H_i，A_i，H_i-1和A_i-1，定义单位面积熵变率ECRA如公式8所示：

${\mathrm{ECRA}}_i={\mathrm{\&Delta;HA}}_i=\frac{{\mathrm{\&Delta;H}}_i}{A_i}=\frac{\left|\frac{H_i-H_{i-1}}{H_i}\right|}{A_i}=\frac{|H_i-H_{i-1}|}{H_{i-1}{A}_i}$ 公式(8)

其中，ΔH_i表示第i幅状态转移图的熵变值，ΔHA_i表示第i幅状态转移图的单位面积熵变值。

以上介绍了本发明基于状态转移的工程风险监测系统的组成，本发明在揭示上述系统的同时，还揭示一种基于状态转移的工程风险监测方法；所述方法包括如下步骤：

【步骤S1】原始数据处理步骤；获取原始数据，并将原始数据经过无量纲化处理。

【步骤S2】状态转移可视化表达步骤；根据经过原始数据处理模块处理的数据生成工程中各状态的状态转移图。

首先，所述方法通过一组时序数据为S表示工程中的各状态变量数据，那么S^k就表示在k时刻工程的状态，并且每一个状态都由n个代表工程状态的内部变量所组成；

S={S⁰,S¹,...,S^k-1,S^k,S^k+1,...}；

$S^k=\{s_1^k,s_2^k,...,s_i^k,..,s_n^k\};$

其中， $s_i^k\&Element;[{\min }_i,{\max }_i],{\min }_i=\min \{s_i^k,0\&le;k\&le;|S\left|\right\},{\max }_i=\max \{s_i^k,0\&le;k\&le;|S\left|\right\};$

假定某工程有n个可测分变量，那么在k时刻，其状态变量S^k表示为：

$S^k=\{s_1^k,s_2^k,.....,s_n^k\}$

定义第i个可测分变量值的量程范围分别为min_i和max_i，那么有：

$s_i^k\&Element;[{\min }_i,{\max }_i]$

所述状态转移可视化表达步骤包括第一状态分量映射步骤，将第k时刻第i个分量通过公式1将其缩放限制在0～π之间，并映射到坐标轴点；

${\mathrm{\&theta;}}_i^k=\frac{s_i^k}{{\max }_i-{\min }_i}$ 公式(1)

所述状态转移可视化表达步骤包括第二状态分量映射步骤，选择平面上一定点O，以O为原点坐标，建立一个同心半圆坐标系；对于任意一个状态变量的分量通过公式2将其映射为平面上的一点

$\left\{\begin{array}{c}{P_i}^k.x=i*u*\cos \left({\mathrm{\&theta;}}_i^k\right)\\ {P_i}^k.y=i*u*\sin \left({\mathrm{\&theta;}}_i^k\right)\end{array}\right.$ 公式(2)

其中，u表示两个相邻轴半径的差值；

对于k时刻的状态变量S^k，通过公式3，映射为平面上的点并将其依次连接，得到一条折线L_k来反映k时刻的状态变量；而系统的每一个状态则可以通过这些唯一的折线来表示，因此有：

$L^k=\{O,P_1^k,P_2^k,P_3^k,...,{P_i}^k,...P_n^k\}1\&le;i<n$ 公式(3)

在平面直角坐标系中，首先画出n个等距同心半圆，其为N维欧几里得空间R^N的同心半圆坐标轴；即使这些坐标轴的长度并不一样，但是每个半圆都有相同的角度范围，即0～π；每一个半圆代表一根坐标轴，且方向为逆时针方向；坐标轴d_i代表了第i维的数据，并通过公示将数据点缩放限制在0～π之间；

所述状态转移图用来描述一段时间内系统状态转移的过程；从k₁时刻到k₂时刻，系统的状态顺序集为{S^k1，S^k1+1，……S^k2-1，S^k2}，对应的折线集为{L^k1，L^k1+1，……L^k2-1，L^k2}；

所述状态转移可视化表达步骤还包括：

状态折线绘制步骤，如图3所示，在同一个同心半圆坐标轴下画出所有状态折线；

封闭区域形成步骤，对于相邻两根状态折线，用直线将折线的末端相连形成一个封闭区域，A^k是由折线L^k和折线L^k+1以及他们的末端连线围成的封闭区域；

填色步骤，如图4所示，对封闭区域A¹,A²,……,A^k-1进行一次填色，当区域相互重叠的时候，重叠区域的颜色就会加深，颜色的深度取决于这些区域重叠的次数。

通常，工程中有一系列的状态；所述状态转移可视化表达步骤生成工程中一系列状态的状态转移图；一组有序的时间序列图根据周期T以及递进D绘制产生系列轨迹图，0≤D≤T；其中第一张图所选取的原始数据系列为{X₁,X₂,X₃,...,X_T}，第二张图所选取的原始数据系列为{X_D+1,X_D+2,...,X_D+T}，那么第N张图所选取的原始数据系列为{X_D*(N-1)+1,X_D*(D-1)+2,...,X_D*(N-1)+T}，不同的轨迹图代表了系统不同时期的属性。

【步骤S3】图形特征提取步骤；获取所述状态转移可视化表达模块生成的状态转移图，并从中提取状态转移图的设定特征。

所述图形特征提取步骤提取状态转移图的设定特征包括：状态转移图的质心与质心间距、图形熵、图形面积、单位面积熵变率中的一个或多个。

质心为物质系统上被认为质量集中于此的一个假想点，质心的位置矢量是质点组中各个质点的位置矢量根据其对应质量加权平均之后的平均矢量；质心的变化反映图形特征的变化，对于每一幅状态转移图的质心用公式4表示：

$\stackrel{\&OverBar;}{X}=\frac{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N(P_{V(i.j)}*i)}{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N{P}_{V(i.j)}}$ 公式(4)

$\stackrel{\&OverBar;}{Y}=\frac{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N(P_{V(i.j)}*j)}{{\mathrm{\&Sigma;}}_{i=1}^M{\mathrm{\&Sigma;}}_{j=1}^N{P}_{V(i.j)}}$

其中，M*N表示图像在X方向M个像素点，Y方向N个像素点；Pv(i,j)表示一幅图像中像素点(i,j)的灰度值；

由公式4可知，图像质心与图像中像素点的灰度值明显相关；为了更易发现状态的变化，通过缩小灰度值范围求出同一副图的不同质心，而质心间距用于衡量图形质心的变化，定义如公式5所示；

$L_{i,i+1}=\sqrt{{(X_i-X_{i+1})}^2+{(Y_i-Y_{i+1})}^2}$ 公式(5)

其中，第i幅图的质心坐标为(X_i,Y_i)；第i+1幅图的质心坐标为(X_i+1,Y_i+1)；

图形熵是用来估计不同图像之间信息的变化量的指标，同时也反映了不同状态下系统的能量变化；系统能量越均匀，状态无序，熵值越大；系统能量不均匀，状态有序，熵值越小；公式6将图形熵定义为：

$\left\{\begin{array}{c}E={\mathrm{\&Sigma;}}_{i=0}^{G}d\left(i\right){\log }_{2}d\left(i\right)\\ d\left(i\right)=p\left(i\right)/N\end{array}\right.$ 公式(6)

其中，E表示图形熵；G表示8进制图像直方图中灰度值，取值范围为0-255；d(i)表示标准化后每一个灰度值出现的频率；N表示整幅图像的总像素点个数；p(i)表示灰度值为i的像素点个数；

图形面积变量包含一个值，统计的是图像像素点的个数，当图像的面积发生变化时，图像的像素点个数也会发生变化，像素点变化越大，代表着图像的面积特征变化越明显；定义如公式7：

$\left\{\begin{array}{c}A={\mathrm{\&Sigma;}}_{i=1}^X{\mathrm{\&Sigma;}}_{j=1}^{Y}F\left(\mathrm{Pv}(i,j)\right)\\ F\left(a\right)=\left\{\begin{array}{cc}1& a\&NotEqual;255\\ 0& a=255\end{array}\right.\end{array}\right.$ 公式(7)

其中，X表示一幅图像中水平方向的像素点个数；Y表示一幅图像中垂直方向的像素点个数；P_v(i,j)表示一幅图像中像素点(i,j)的灰度值；F()表示一幅图像中的某一个像素点是否为白色，若a=255为白色；a表示像素点的灰度值，取值范围为0到255；

对于两幅相邻的状态转移图，如果他们的图形熵和图形面积的值分别是H_i，A_i，H_i-1和H_i-1，定义单位面积熵变率ECRA如公式8所示：

${\mathrm{ECRA}}_i={\mathrm{\&Delta;HA}}_i=\frac{{\mathrm{\&Delta;H}}_i}{A_i}=\frac{\left|\frac{H_i-H_{i-1}}{H_i}\right|}{A_i}=\frac{|H_i-H_{i-1}|}{H_{i-1}{A}_i}$ 公式(8)

其中，ΔH_i表示第i幅状态转移图的熵变值，ΔHA_i表示第i幅状态转移图的单位面积熵变值。

【步骤S4】自动学习步骤；根据所述图像特征提取模块提取的图形特征值，计算每一特征值与历史平均值的比值，并设定阈值。

【步骤S5】图像自动辨识步骤；根据所述图形特征提取模块提取状态转移图的设定特征以及自动学习模块设定的阈值判断工程是否处于风险状态。

【风险评判的方法背景及特点】

在自然界和人类社会活动中，除了渐变的和连续光滑的变化现象外，还存在着大量的突然变化和跃迁现象，如岩石的破裂、桥梁的崩塌、地震、海啸、细胞的分裂、生物的变异、人的休克、情绪的波动、战争、市场变化、企业倒闭、经济危机等。这种现象被称为“突变”,其通常发生在非线性动力系统中，是系统连续性状态突然中断导致质变的过程，是系统从某一个稳定态到另一个稳定态的转化。

地下工程的风险常常具有这种特征，由于地下工程结构和施工方法的特殊性，特别是施工开挖过程中周边环境(介质)的不确定性和复杂性，其灾害发生往往具有“突变”的特性。与常规的预测对象不同，突变现象不具有渐变、量变特点，而是出现了不连续和质变的显著特征，因此采用微积分等传统的数据模型或数学工具很难描述或预测这种飞跃和不连续过程。

本发明是一种基于灾变理论和状态演化的可视化风险监测方法，该方法借鉴空间重构技术和可视化数据分析技术在表达复杂系统规律和特征的优势，从状态变化的角度来反映整个系统中的变化，而不单纯以时间作为主轴来考虑问题。通过图形来反映一段时间内随着系统从一个状态到另一个状态转移的情况，不仅能直观的反应系统状态变化的过程，更能深层次探索系统状态变化的原因，从而针对性地解决地下工程灾害的不确定性和突发性问题。

【风险评判的方法原理与思想】

本发明风险判断的方法是在目前突变研究的相关理论上发展得到的。

突变理论：由法国数学家托姆(R.Thom)在1969年提出。突变理论认为：系统所处的状态，可用一组参数描述。当参数在某个范围内变化，该函数值有不止一个极值时，系统必然处于不稳定状态。系统从一种稳定状态进入不稳定状态，随参数的再变化，又使不稳定状态进入另一种稳定状态，那么，系统状态就在这一刹那间发生了突变。

耗散理论：由比利时物理化学家普利高津在1969年提出。耗散理论认为：当系统离开平衡态的参数达到一定阈值时，系统将会出现“行为临界点”，在越过这种临界点后系统将离开原来的热力学无序分支，发生突变而进入到一个全新的稳定有序状态，即“耗散结构”。耗散结构是在与外界交换物质和能量的过程中，通过能量耗散和内部非线性动力学机制的作用，经过突变而形成并持久稳定的宏观有序结构。因此，突变的出现与能量变化（熵变）密切关联。

突变早期预警征兆研究：突变理论和耗散理论均指出突变往往发生在非线性系统的分叉处的临界点，近年来Scheffer[1]和Dakos[2]等指出，在分叉出现之前，反映系统状态的时间序列数据将会出现数据波动幅度下降、频度放缓和数据分布不对称等结构性异常变化特征。

[1]Scheffer M,Bascompte J,Brock W A,et al.Early-warning signals forcritical transitions[J].Nature.2009,461(7260):53-59.

[2]Dakos V,Scheffer M,Van Nes E H,et al.Slowing down as an earlywarning signal for abrupt climate change[J].Proceedings of the NationalAcademy of Sciences.2008,105(38):14308-14312.

地下工程的主体作用对象是岩土介质。地下工程岩体经多次开挖，岩体本身所具有的平衡结构受到严重的干扰，为了维持或恢复原有自身的平衡，在开挖强度不是太大的情况下，岩体将通过应力和变形等方式自动地对其本身结构进行一系列地调整，最后达到新的平衡。。但是，当人工开挖强度过大，对岩体的干扰破坏作用，超过了岩体自我调节的限度时，再加上其它外界因素（如地下水、地震等）的影响，岩体将随着时间的推移由原来的稳定状态走向失稳，也就是岩体系统出现了失稳，即地下工程系统发生了突变（灾害性事故）。

在地下工程灾害性事故演化过程中，整个系统符合耗散结构理论和灾变理论的描述，其灾变的原因是当系统处于某一临界点时，伴随着一定能量的转换和熵变，岩体远离平衡态，引起围岩体失稳。

因此，本发明的思路是将一组反映地下工程特征的时间序列数据按照一定的规则绘制成图像，结合熵变、分叉、波动频度和不对称性等突变前兆特征，通过定量化的手段提取图像中反映系统结构变化的特征值，从而达到进行地下工程风险早期预警的目的。

【风险评判定量计算指标含义】

本发明给出了图像质心、质心间距、、图形面积、图形熵和单位面积熵变率五个定量计算指标，评判系统是否出现了突变前兆，保证了突变前兆判断可以通过计算机自动实现。这多个定量指标与系统突变特征的关系如下：

1、质心、质心间距：质心体现了整个系统的数据波动结构特征和不对称性，质心间距反映了数据结构变化，这与突变前兆特征中的不对称性存在明显关联，当系统中的状态数据值向某一领域漂移时，其图形质心必然发生移动。质心位置特征和质心间距增大是系统出现突变的前兆信号，系统可能从一个稳定状态向另一个状态发生跃变。

2、图形面积：该特征变量体现了系统的波动幅度和频度，与突变前兆特征中的波动幅度和分叉特性存在明显关联，图像面积突然减少或增大是系统出现突变的前兆信号，系统可能从一个稳定状态向另一个状态发生跃变。

3、图像熵和单位面积熵变率：熵是用来描述、表征系统不确定程度。图像熵体现系统能量特征，图像熵的变化体现了系统能量耗散的过程，单位面积熵变率反映了系统的熵变，高熵变是系统出现突变的前兆信号，系统可能从一个稳定状态向另一个状态发生跃变。

【状态转移图绘制原理】

本发明提出了同心半圆坐标系，将时序多维数据集映射到二维平面，忽略了时间序列数据的时间性，通过图像宏观结构和微观纹理变化使得突变征兆更加方便和明显的体现出来。

同心半圆坐标系：是平行坐标技术和星形图技术综合后的改进。它的基本思想是将n维数据归一化处理后通过转移函数，映射到等距离的半圆轴上，每条轴线都对应于一个属性维，从圆心开始，依次连接属性值在轴线上映射点，形成一条折线，用于表达一组独立的n维数据。

数据时序关系的表达：前后两组n维数据所映射的折线按逆时针方向封闭成一个多边形区域，并进行涂色处理，以表达数据在时间上的转移过程。当数据转移过程中，所经过的区域的颜色深度，与经过的次数成正比。这样通过颜色可以显示数据所经过的范围以及数据变化过程的规律性。

与其它多维平面坐标相比，同心半圆坐标系不同于以时间轴为中心的坐标系，它更加注重状态的变化过程，将原本难以直观描述的系统状态，通过状态轨迹图谱法将单维或多维的数据绘制在一张图片上，以图形的特征量来判断当前的系统状态。这极大地压缩数据，提高了高频采样监测数据的处理，同时由于绘制的图反映了数据变化过程，对于数据的缺失或异常数据的处理具有较强的鲁棒性。

实施例二

本发明以长江隧道工程为例，说明方法的应用效果。

整个长江隧道实时采集到的数据维度众多，本实验利用盾构密封舱五维压力值进行风险预测。表1说明了长江隧道工程实验的基本信息。实验选取了长江隧道上行线2008年5月1日凌晨0:02分至2008年5月22日下午5:42分共10563条舱内压力数据作为数据样本，采样间隔为3分钟。

表1.长江隧道实验基本信息

实验中，以每240条记录作为一个周期绘制一张状态转移图，后一张图与前一张图的步长间隔为5条，共绘制2068幅状态转移图。图5选取了从2008/5/153:18到2008/5/170:33的20幅连续的状态转移图。图6、图7分别表示三种指标的判别曲线，其中虚线是实际在现场发现土体出现明显下陷的已知危险点，小圆圈代表了灾变征兆信号。图6实线代表了图形面积判别曲线，点线代表了图形熵判别曲线，图7实线代表了ECRA判别信号。

根据图5可以看出，从第1311幅图到第1314幅图是正常状态的图，图形密度分布均衡，图形形状相似。但是从第1316幅图到第1320幅图，与之前的图形明显有差异，图形的轮廓明显变大。对比第1483到第1489幅这组图，第1491幅图（5/170:03）有明显的差别，这意味着这里可能是一个风险信号。

从图6可以看出，面积曲线在（5/175:33）发出了明显的信号，而图形熵值在（5/194:33）也发出了强烈的信号。从图7可以看出，在灾变发生之前，共出现了三次明显的ECRA信号，分别是5/155:03，5/170:03和5/1914:00。其中前两次信号均早于实际灾害发生。

表2说明了灾害发生前ECRA峰值与历史平均值的关系。根据表2，可以看出前两个信号的ECRA值分别是前50%历史峰值平均值的11倍和14倍。因此，这两个信号均是可信的，这两个时间段工程极有可能处于风险状态。

实施例三

本实施例与实施例一的区别在于，本实施例中，所述状态转移可视化表达模块采用同心半圆坐标系建立表达多维平面坐标，过程如下：

N维欧几里得空间R^N同心半圆坐标系由N个等距同心半圆坐标轴组成。每一个半圆为一条坐标轴，每条轴线都对应于一个多维数据集的一个分量，从圆心开始，从内到外依次编号：d₀，d₁，…，d_i，…，d_N,单条周轴线圆心角范围为0～180°，按逆时针方向递增；

假定一个N维时序数据的在k时刻的第i个分量为通过公式1和公式2函数转换映射为第i条轴线上的一个点从圆心开始，依次连接形成一条折线，用于表达一组独立的n维数据；

${\mathrm{\&theta;}}_i^k=\frac{s_i^k}{{\max }_i-{\min }_i}$ 公式(1)

$\left\{\begin{array}{c}{P_i}^k.x=i*u*\cos \left({\mathrm{\&theta;}}_i^k\right)\\ {P_i}^k.y=i*u*\sin \left({\mathrm{\&theta;}}_i^k\right)\end{array}\right.$ 公式(2)

其中，u表示两个相邻轴半径的差值。

综上所述，本发明提出的基于状态转移的工程风险监测系统及方法，借鉴空间重构技术和可视化数据分析技术在表达复杂系统规律和特征的优势，从状态变化的角度分析系统的异常变化，从而发现工程的潜在风险。该方法与一般地下工程风险监测方法相比，不仅能够针对性地解决地下工程灾害的不确定性和突发性问题，而且以可视化的方式自动监测工程风险，具有直观性和智能性。并且本发明也隐含地给出了一种突发事件风险监测的一般方法。

这里本发明的描述和应用是说明性的，并非想将本发明的范围限制在上述实施例中。这里所披露的实施例的变形和改变是可能的，对于那些本领域的普通技术人员来说实施例的替换和等效的各种部件是公知的。本领域技术人员应该清楚的是，在不脱离本发明的精神或本质特征的情况下，本发明可以以其它形式、结构、布置、比例，以及用其它组件、材料和部件来实现。在不脱离本发明范围和精神的情况下，可以对这里所披露的实施例进行其它变形和改变。

Claims (13)

1.一种基于状态转移的工程风险监测系统，其特征在于，所述系统包括：

原始数据处理模块，用以获取原始数据，并将原始数据经过无量纲化处理；

状态转移可视化表达模块，与所述原始数据处理模块连接，用以根据经过原始数据处理模块处理的数据，生成工程中各状态的状态转移图；

图形特征提取模块，与所述状态转移可视化表达模块连接，用以获取所述状态转移可视化表达模块生成的状态转移图，并从中提取状态转移图的设定特征；

自动学习模块，与所述图形特征提取模块连接，用以获取图形特征提取模块提取的图形特征，计算当前特征值与设定历史平均值的比值，并设定阈值；

图像自动辨识模块，与所述图形特征提取模块、自动学习模块连接，用以根据所述图形特征提取模块提取状态转移图的设定特征以及自动学习模块设定的阈值判断工程是否处于风险状态；

所述状态转移可视化表达模块是将数据进行向空间转换、可视化表达的过程；具体包括：

同心半圆坐标系建立单元，用以建立同心半圆坐标系，将预处理的数据根据公式映射到该同心半圆坐标系中；

状态折线绘制单元，用以绘制同一个同心半圆坐标轴下所有状态折线；

状态转移图生成单元，用直线将相邻两根状态折线的末端相连形成一个封闭区域并进行填色，生成一幅状态转移图；

系列状态转移图绘制单元，用以根据预设的周期和递进参数绘制系列状态转移图。

2.根据权利要求1所述的基于状态转移的工程风险监测系统，其特征在于：

所述状态转移可视化表达模块采用同心半圆坐标系建立表达多维平面坐标，过程如下：

N维欧几里得空间R^N同心半圆坐标系由N个等距同心半圆坐标轴组成；每一个半圆为一条坐标轴，每条轴线都对应于一个多维数据集的一个分量，从圆心开始，从内到外依次编号：d₀，d₁，…，d_i，…，d_N,单条周轴线圆心角范围为0～180°，按逆时针方向递增；

假定一个N维时序数据的在k时刻的第i个分量为通过公式1和公式2函数转换映射为第i条轴线上的一个点从圆心开始，依次连接形成一条折线，用于表达一组独立的n维数据；

其中，u表示两个相邻轴半径的差值，

3.根据权利要求1所述的基于状态转移的工程风险监测系统，其特征在于：

所述原始数据处理模块将原始数据经过无量纲化处理；

无量纲化处理通过数学变换来消除原始变量量纲影响的方法；将数据全部缩放到0～π之间来消除原始变量的量纲影响；具体包括：

假定在k时刻，第i个分量通过公式1将其缩放限制在0～π之间，并映射到坐标轴点；

其中，

4.根据权利要求1所述的基于状态转移的工程风险监测系统，其特征在于：

所述自动学习模块根据图形特征提取模块，提取状态转移图的特征变量值，计算当前特征值与平均值比值的过程；平均值分别按历史数据中数值大小前5％、10％、20％、50％的数值计算；

所述图像自动辨识模块用以实现状态转移图风险识别的主要过程；具体包括：

(1)根据图形特征提取模块，获取系列状态转移图的特征变量值并绘制相应的指标曲线；

(2)选择是否进入自动学习模块；若进入自动学习模块，则获取特征变量当前值与历史平均值的比值；若该比值在设定阈值外即判断为异常；

若不进入自动学习模块，则计算特征变量当前值与前一值的比值；若该比值在设定阈值外即判断为异常；

(3)比值超过阈值越大，判断状态转移图越为异常，此时系统越有可能处于异常状态。

5.根据权利要求1所述的基于状态转移的工程风险监测系统，其特征在于：

所述系统通过一组时序数据为S表示工程中的各状态变量数据，那么S^k就表示在k时刻工程的状态，并且每一个状态都由n个代表工程状态的内部变量所组成；

S＝{S⁰,S¹,…,S^k-1,S^k,S^k+1,…}；

$S^k=\{s_1^k,s_2^k,\mathrm{...},s_i^k,\mathrm{..},s_n^k\};$

其中，

假定某工程有n个可测分变量，那么在k时刻，其状态变量S^k表示为：

$S^k=\{s_1^k,s_2^k,\mathrm{.....},s_n^k\}$

定义第i个可测分变量值的量程范围分别为min_i和max_i，那么有：

$s_i^k\&Element;\&lsqb;{\min }_i,{\max }_i\&rsqb;$

所述状态转移可视化表达模块包括第一状态分量映射模块，用以将第k时刻第i个分量通过公式1将其缩放限制在0～π之间，并映射到坐标轴点；

所述状态转移可视化表达模块包括第二状态分量映射模块，用以选择平面上一定点O，以O为原点坐标，建立一个同心半圆坐标系；对于任意一个状态变量的分量通过公式2将其映射为平面上的一点

其中，u表示两个相邻轴半径的差值；

对于k时刻的状态变量S^k，通过公式3，映射为平面上的点并将其依次连接，得到一条折线L_k来反映k时刻的状态变量；而系统的每一个状态则可以通过这些唯一的折线来表示，因此有：

在平面直角坐标系中，首先画出n个等距同心半圆，其为N维欧几里得空间R^N的同心半圆坐标轴；即使这些坐标轴的长度并不一样，但是每个半圆都有相同的角度范围，即0～π；每一个半圆代表一根坐标轴，且方向为逆时针方向；坐标轴d_i代表了第i维的数据，并通过公式将数据点缩放限制在0～π之间；

所述状态转移图用来描述一段时间内系统状态转移的过程；从k₁时刻到k₂时刻，系统的状态顺序集为{S^k1，S^k1+1，……S^k2-1，S^k2}，对应的折线集为{L^k1，L^k1+1，……L^k2-1，L^k2}；

所述状态转移可视化表达模块还包括：

状态折线绘制模块，用以在同一个同心半圆坐标轴下画出所有状态折线；

封闭区域形成模块，对于相邻两根状态折线，用直线将折线的末端相连形成一个封闭区域，A^k是由折线L^k和折线L^k+1以及他们的末端连线围成的封闭区域；

填色模块，用以对封闭区域A¹,A²,……,A^k-1进行一次填色，当区域相互重叠的时候，重叠区域的颜色就会加深，颜色的深度取决于这些区域重叠的次数。

6.根据权利要求1所述的基于状态转移的工程风险监测系统，其特征在于：

所述状态转移可视化表达模块用以生成工程中一系列状态的状态转移图；

一组有序的时间序列图根据周期T以及递进D绘制产生系列轨迹图，0≤D≤T；其中第一张图所选取的原始数据系列为{X₁,X₂,X₃,…,X_T}，第二张图所选取的原始数据系列为{X_D+1,X_D+2,…,X_D+T}，那么第N张图所选取的原始数据系列为{X_D*(N-1)+1,X_D*(N-1)+2,…,X_D*(N-1)+T}，不同的轨迹图代表了系统不同时期的属性。

7.根据权利要求1所述的基于状态转移的工程风险监测系统，其特征在于：

所述图形特征提取模块提取状态转移图的设定特征包括：状态转移图的质心、质心间距、图形熵、图形面积、单位面积熵变率中的一个或多个；

质心为物质系统上被认为质量集中于此的一个假想点，质心的位置矢量是质点组中各个质点的位置矢量根据其对应质量加权平均之后的平均矢量；质心的变化反映图形特征的变化，对于每一幅状态转移图的质心用公式4表示：

其中，M*N表示图像在X方向M个像素点，Y方向N个像素点；Pv(i,j)表示一幅图像中像素点(i,j)的灰度值；

由公式4可知，图像质心与图像中像素点的灰度值明显相关；为了更易发现状态的变化，通过缩小灰度值范围求出同一副图的不同质心，而质心间距用于衡量图形质心的变化，定义如公式5所示；

其中，第i幅图的质心坐标为(X_i,Y_i)；第i+1幅图的质心坐标为(X_i+1,Y_i+1)；

图形熵是用来估计不同图像之间信息的变化量的指标，同时也反映了不同状态下系统的能量变化；系统能量越均匀，状态无序，熵值越大；系统能量不均匀，状态有序，熵值越小；公式6将图形熵定义为：

其中，E表示图形熵；G表示8进制图像直方图中灰度值，取值范围为0-255；d(i)表示标准化后每一个灰度值出现的频率；N表示整幅图像的总像素点个数；p(i)表示灰度值为i的像素点个数；

图形面积变量包含一个值，统计的是图像像素点的个数，当图像的面积发生变化时，图像的像素点个数也会发生变化，像素点变化越大，代表着图像的面积特征变化越明显；定义如公式7：

其中，X表示一幅图像中水平方向的像素点个数；Y表示一幅图像中垂直方向的像素点个数；P_v(i,j)表示一幅图像中像素点(i,j)的灰度值；F()表示一幅图像中的某一个像素点是否为白色，若a＝255为白色；a表示像素点的灰度值，取值范围为0到255；

对于两幅相邻的状态转移图，如果他们的图形熵和图形面积的值分别是H_i，A_i，H_i-1和A_i-1，定义单位面积熵变率ECRA如公式8所示：

其中，ΔH_i表示第i幅状态转移图的熵变值，ΔHA_i表示第i幅状态转移图的单位面积熵变值；

所述图形特征提取模块根据上述定义分别计算每幅状态转移图的图形特征变量值。

8.一种权利要求1至7之一所述基于状态转移的工程风险监测系统的监测方法，其特征在于，所述方法包括如下步骤：

原始数据处理步骤；获取原始数据，并将原始数据经过无量纲化处理；

状态转移可视化表达步骤；根据经过原始数据处理模块处理的数据生成工程中各状态的状态转移图；

图形特征提取步骤；获取所述状态转移可视化表达模块生成的状态转移图，并从中提取状态转移图的设定特征；

自动学习步骤；获取所述图形特征提取模块提取的图形特征，计算每一特征值与历史平均值的比值，并设定阈值；

图像自动辨识步骤；根据所述图形特征提取模块提取状态转移图的设定特征以及自动学习模块设定的阈值判断工程是否处于风险状态；

所述状态转移可视化表达步骤中，将数据进行向空间转换、可视化表达的过程；具体包括如下步骤：

建立同心半圆坐标系，将预处理的数据根据公式映射到同心半圆坐标系中；

绘制同一个同心半圆坐标轴下所有状态折线；

用直线将相邻两根状态折线的末端相连形成一个封闭区域并进行填色，生成一幅状态转移图；

根据预设的周期和递进参数绘制系列状态转移图。

9.根据权利要求8所述的监测方法，其特征在于：

所述原始数据处理步骤中，将原始数据经过无量纲化处理；

无量纲化处理通过数学变换来消除原始变量量纲影响的方法；将数据全部缩放到0～π之间来消除原始变量的量纲影响；具体过程如下：

假定在k时刻，第i个分量通过公式1将其缩放限制在0～π之间，并映射到坐标轴点；

其中，

10.根据权利要求8所述的监测方法，其特征在于：

所述自动学习步骤中，自动学习模块根据图形特征提取步骤提取状态转移图的特征变量值，计算当前特征值与平均值比值的过程；平均值分别按历史数据中数值大小前5％、10％、20％、50％的数值计算；

所述图像自动辨识步骤实现状态转移图风险识别的主要过程；具体包括如下步骤：

(1)根据图形特征提取模块获取系列状态转移图的特征变量值并绘制相应的指标曲线；

(2)选择是否进入自动学习模块；若进入自动学习模块，则获取特征变量当前值与历史平均值的比值；若该比值在设定阈值外即判断为异常；

若不进入自动学习模块，则计算特征变量当前值与前一值的比值；若该比值在设定阈值外即判断为异常；

(3)比值超过阈值越大，判断状态转移图越为异常，此时系统越有可能处于异常状态。

11.根据权利要求8所述的监测方法，其特征在于：

所述方法通过一组时序数据为S表示工程中的各状态变量数据，那么S^k就表示在k时刻工程的状态，并且每一个状态都由n个代表工程状态的内部变量所组成；

S＝{S⁰,S¹,…,S^k-1,S^k,S^k+1,…}；

$S^k=\{s_1^k,s_2^k,\mathrm{...},s_i^k,\mathrm{..},s_n^k\};$

其中，

假定某工程有n个可测分变量，那么在k时刻，其状态变量S^k表示为：

$S^k=\{s_1^k,s_2^k,\mathrm{.....},s_n^k\}$

定义第i个可测分变量值的量程范围分别为min_i和max_i，那么有：

$s_i^k\&Element;\&lsqb;{\min }_i,{\max }_i\&rsqb;$

所述状态转移可视化表达步骤包括第一状态分量映射步骤，将第k时刻第i个分量通过公式1将其缩放限制在0～π之间，并映射到坐标轴点；

所述状态转移可视化表达步骤包括第二状态分量映射步骤，选择平面上一定点O，以O为原点坐标，建立一个同心半圆坐标系；对于任意一个状态变量的分量通过公式2将其映射为平面上的一点

其中，u表示两个相邻轴半径的差值；

对于k时刻的状态变量S^k，通过公式3，映射为平面上的点并将其依次连接，得到一条折线L_k来反映k时刻的状态变量；而系统的每一个状态则可以通过这些唯一的折线来表示，因此有：

在平面直角坐标系中，首先画出n个等距同心半圆，其为N维欧几里得空间R^N的同心半圆坐标轴；即使这些坐标轴的长度并不一样，但是每个半圆都有相同的角度范围，即0～π；每一个半圆代表一根坐标轴，且方向为逆时针方向；坐标轴d_i代表了第i维的数据，并通过公式将数据点缩放限制在0～π之间；

所述状态转移图用来描述一段时间内系统状态转移的过程；从k₁时刻到k₂时刻，系统的状态顺序集为{S^k1，S^k1+1，……S^k2-1，S^k2}，对应的折线集为{L^k1，L^k1+1，……L^k2-1，L^k2}；

所述状态转移可视化表达步骤还包括：

状态折线绘制步骤，在同一个同心半圆坐标轴下画出所有状态折线；

封闭区域形成步骤，对于相邻两根状态折线，用直线将折线的末端相连形成一个封闭区域，A^k是由折线L^k和折线L^k+1以及他们的末端连线围成的封闭区域；

填色步骤，对封闭区域A¹,A²,……,A^k-1进行一次填色，当区域相互重叠的时候，重叠区域的颜色就会加深，颜色的深度取决于这些区域重叠的次数。

12.根据权利要求8所述的监测方法，其特征在于：

所述状态转移可视化表达步骤生成工程中一系列状态的状态转移图；

一组有序的时间序列图根据周期T以及递进D绘制产生系列轨迹图，0≤D≤T；其中第一张图所选取的原始数据系列为{X₁,X₂,X₃,…,X_T}，第二张图所选取的原始数据系列为{X_D+1,X_D+2,…,X_D+T}，那么第N张图所选取的原始数据系列为{X_D*(N-1)+1,X_D*(N-1)+2,…,X_D*(N-1)+T}，不同的轨迹图代表了系统不同时期的属性。

13.根据权利要求8所述的监测方法，其特征在于：

所述图形特征提取步骤提取状态转移图的设定特征包括：状态转移图的质心与质心间距、图形熵、图形面积、单位面积熵变率中的一个或多个；

质心为物质系统上被认为质量集中于此的一个假想点，质心的位置矢量是质点组中各个质点的位置矢量根据其对应质量加权平均之后的平均矢量；质心的变化反映图形特征的变化，对于每一幅状态转移图的质心用公式4表示：

其中，M*N表示图像在X方向M个像素点，Y方向N个像素点；Pv(i,j)表示一幅图像中像素点(i,j)的灰度值；

由公式4可知，图像质心与图像中像素点的灰度值明显相关；为了更易发现状态的变化，通过缩小灰度值范围求出同一副图的不同质心，而质心间距用于衡量图形质心的变化，定义如公式5所示；

其中，第i幅图的质心坐标为(X_i,Y_i)；第i+1幅图的质心坐标为(X_i+1,Y_i+1)；

图形熵是用来估计不同图像之间信息的变化量的指标，同时也反映了不同状态下系统的能量变化；系统能量越均匀，状态无序，熵值越大；系统能量不均匀，状态有序，熵值越小；公式6将图形熵定义为：

其中，E表示图形熵；G表示8进制图像直方图中灰度值，取值范围为0-255；d(i)表示标准化后每一个灰度值出现的频率；N表示整幅图像的总像素点个数；p(i)表示灰度值为i的像素点个数；

图形面积变量包含一个值，统计的是图像像素点的个数，当图像的面积发生变化时，图像的像素点个数也会发生变化，像素点变化越大，代表着图像的面积特征变化越明显；定义如公式7：

其中，X表示一幅图像中水平方向的像素点个数；Y表示一幅图像中垂直方向的像素点个数；P_v(i,j)表示一幅图像中像素点(i,j)的灰度值；F()表示一幅图像中的某一个像素点是否为白色，若a＝255为白色；a表示像素点的灰度值，取值范围为0到255；

对于两幅相邻的状态转移图，如果他们的图形熵和图形面积的值分别是H_i，A_i，H_i-1和A_i-1，定义单位面积熵变率ECRA如公式8所示：

其中，ΔH_i表示第i幅状态转移图的熵变值，ΔHA_i表示第i幅状态转移图的单位面积熵变值；

所述方法还包括：图形特征提取模块根据上述定义分别计算每幅状态转移图的图形特征变量值。