CN103984938B - 一种遥感时间序列异常检测方法 - Google Patents

Abstract

本发明公开一种技术方案，可对遥感时间序列中存在的多个时间的异常同时进行检测，在显著性水平为α＝0.05时异常检测结果的整体置信度为95％、单个异常置信度大于99％，同时给出异常的时间、水平、程度和置信度信息。本发明的基本思路为：首先将原始遥感时间序列数据进行一阶季节差分，得到的差分序列作为季节误差序列；对季节误差序列的均值和标准差进行稳健化估计，并将季节误差序列进行z标准化，得到的标准分数序列作为时间序列异常检验统计量序列；在给定的显著性水平和多重检验总体误差率控制下，计算检验统计量临界值；将检验统计量序列的绝对值与临界值相比较，满足预定条件时将其标记为异常点；由各个异常点的检验统计量，计算其p值，并转换为置信水平。通过上述流程检测出遥感时间序列中的异常，并得到异常的时间、水平、程度和置信度。

Description

一种遥感时间序列异常检测方法

技术领域

本发明属于遥感图像处理和信息提取领域，涉及遥感图像变化检测和地表覆盖异常信息提取，具体涉及一种基于遥感时间序列数据进行地表覆盖异常变化的检测以及异常信息提取的方法。

背景技术

地表覆盖的状态随着时间推移而不断变化，这些变化既包括正常变化，即周期性季节变化(如植被的发芽到茂盛到落叶)，也包括非正常变化，如地表覆盖类型的变化(如植被变成建成区、林地变为农田、旱地变为水田)，和地表覆盖状态的异常变化(如森林砍伐、过度放牧、洪水淹没、过火、植被病虫害等)。

利用遥感影像进行地表覆盖变化的研究已经非常广泛。传统的多时相方法(如分类后比较法、差值法等)一般选择不同年份的相同或相近时相影像进行分析，若影像时相不相近，区分地物的正常变化和非正常变化比较困难。多时相方法一般只能检测所选定时相的影像间是否发生变化，并未对变化类型(地表覆盖类型变化和地表覆盖异常变化)进行区分，且变化发生的时间段和变化强度是无法得知的。

高时间分辨率的遥感时间序列影像为地表覆盖的高频率长期观测和实时监测提供了可能。利用遥感时间序列影像，不仅可以避免传统的多时相方法中因为云覆盖而导致的数据缺失或时相不匹配问题，而且可以在不同时间尺度上(年内季节性、年际变化和长期趋势)表征地物状态变化的模式。

遥感时间序列是指：对相同空间位置的对象进行固定频率的遥感观测而获得的表征对象状态某种特征的时间序列数据。遥感时间序列异常是指：由于被观测对象状态的异常变化，遥感时间序列表现出非规则性的变化。遥感时间序列异常检测是指：通过某种方法，检测出遥感时间序列中的上述非规则性变化。

现有的遥感时间序列异常检测方法主要分为两类：基于时间序列分解的方法和基于预测的方法。基于时间序列分解的方法的基本思想是用STL分解法将时间序列分解为周期项、趋势项和残余项，进而检测周期项和趋势项中的结构性变化。基于预测的方法的基本思想是用一种模型(如SARIMA、高斯过程、卡尔曼滤波器等)对部分时间序列进行拟合，然后用模型进行预测，将预测值与实际观测值进行对比以判断是否发生异常。但是这两类方法都有不同缺陷。一方面，基于时间序列分解的方法假设遥感时间序列具有相同的周期性变化和线性的趋势变化，而实际上由于气温、降水、光照等的年际变化，地物状态每年的季节性变化并不是相同的，而且年际变化也并非线性的。这就导致时间序列分解误差较大，最终导致周期项和趋势项的结构性变化检测出现较多错误。另一方面，基于预测的方法，一是在建模时假定时间序列是没有异常的，而实际上时间序列是否存在异常是未知的和待检测的，二是用于建模的时间序列长度的选择是主观的，这些都会导致模型的鲁棒性降低。三是基于预测的方法只能按时间顺序逐步检测，且只能检测到一次异常，其后模型将失效。因此，上述两大类方法都有各自明显的缺陷，鲁棒性低，适用性不强。

因此，需要一种新的遥感时间序列异常检测方法，在避免上述问题的同时提高异常检测的精度，同时提取异常的相关信息。

发明内容

针对现有遥感时间序列异常检测方法存在的复杂度和错检率较高、多异常检测率低、异常信息不可靠等缺陷，本发明公开一种新的技术方案，可对遥感时间序列中存在的多个时间的异常同时进行检测，在显著性水平为α＝0.05时异常检测结果的整体置信度为95％、单个异常置信度大于99％，同时给出异常的时间、水平、程度和置信度信息。

本发明的基本思路为：首先将原始遥感时间序列数据进行一阶季节差分，得到的差分序列作为季节误差序列；对季节误差序列的均值和标准差进行稳健化估计，并将季节误差序列进行z标准化，得到的标准分数序列作为时间序列异常检验统计量序列；在给定的显著性水平和多重检验总体误差率控制下，计算检验统计量临界值；将检验统计量序列的绝对值与临界值相比较，满足预定条件时将其标记为异常点；由各个异常点的检验统计量，计算其p值，并转换为置信水平。通过上述流程检测出遥感时间序列中的异常，并得到异常的时间、水平、程度和置信度。

技术方案：实现本发明思路的技术方案流程如图1所示，具体流程描述如下：

A.从遥感时间序列影像中逐像素地提取时间序列数据，将时间序列表示为{Y_t：t＝1，2，...，n}，时间序列的长度为n，季节周期为s；

B.由时间序列{Y_t}计算季节误差序列{a_t：t＝s+1，s+2，...，n}；

C.对季节误差序列{a_t}的均值μ和标准差σ进行估计；

D.对季节误差序列{a_t}进行变换，得到时间序列的异常检验统计量{λ_t：t＝s+1，s+2，...，n}；

E.获得异常检验统计量临界值λ⁽¹⁾、λ^(m)；

F.将异常检验统计量与步骤E中的异常检验统计量临界值相比，满足特定条件时判定为异常，其中异常判定的特定条件为：

当s+1≤t≤n-s时，

当t＞n-s时，|λ_t|＞λ^(m)，

满足以上条件的时间t记为异常时间T；

G.由{Y_t}、{a_t}、{λ_t}获得异常的相关信息：计算如下：

异常水平L_t＝a_t，t∈T，

异常程度D_t＝λ_t，t∈T，

异常置信度C_t＝1-P(|λ|＞|λ_t|)，t∈T；

H.输出遥感时间序列异常检测结果{Y_t，L_t，D_t，C_t：t∈T}。

进一步，在步骤B中，可以通过一阶季节差分计算{Y_t}的季节误差序列{a_t：t＝s+1，s+2，...，n}：

a_t＝Y_t-Y_t-s，t＞s。

进一步，在步骤C中，季节误差序列{a_t}的均值μ和标准差σ根据定义计算，或者进行稳健化估计，稳健化估计方式如下：

μ＝median(a_t)，

进一步，在步骤D中，对季节误差序列{a_t}进行变换的方式为对{a_t}进行z标准化：得到标准分数序列{z_t：t＝s+1，s+2，...，n}，以此标准分数序列{z_t}的取值作为时间序列的异常检验统计量{λ_t：t＝s+1，s+2，...，n}，即λ_t＝z_t。

进一步，在步骤E中，异常检验统计量临界值λ⁽¹⁾、λ^(m)为：

a)单一异常检验的异常检验统计量临界值λ⁽¹⁾用于检验已知异常时刻t₀上的检验统计量是否为异常值，λ⁽¹⁾＞0；

b)多个异常检验的异常检验统计量临界值λ^(m)用于检验未知时刻t上检验统计量λ_t是否为异常值，λ^(m)＞0，其中m＝n-s为序列λ_t的长度。

进一步，异常检验统计量临界值λ⁽¹⁾、λ^(m)为用户预设值，或者通过给定异常检验的显著性水平来进行显著性计算：

a)对于λ⁽¹⁾，在给定的显著性水平α下，异常检验统计量临界值的取值满足以下条件：

或

b)多个异常检验的异常检验统计量临界值λ^(m)用于检验未知时刻t上检验统计量λ_t是否为异常值，λ^(m)＞0，λ^(m)取值为预设值，或者也在给定的显著性水平α下进行显著性计算，对全部λ_t进行显著性检验，进行多重检验总体误差率控制，此时异常检验统计量临界值的取值满足以下条件：

进一步，多重检验总体误差率控制方法为Bonferroni、holm、hochberg、hommel、FDR、BY控制法以及它们改进的方法。

本发明的遥感时间序列异常检测方法，不需要对时间序列数据进行分解或者模型拟合，能在相对较低的时间复杂度和计算复杂度下实现较好的检测效果；可同时对多个时间的异常进行检测，并给出异常的时间、水平、程度和置信度信息；在显著性水平为α＝0.05时异常检测结果的整体置信度为95％、单个异常置信度大于99％；方法不需要观测对象类型和异常的先验知识，适用于检测遥感时间序列变化不规则或不明显的地表覆盖(如水、荒漠)的异常变化，尤其适用于具有明显周期性变化的遥感时间序列的异常检测。本发明的遥感时间序列异常检测方法，复杂度较低，可靠性较高，适应性较强，具有明显的应用价值。

附图说明

附图用来提供对本发明的进一步理解，并且构成本说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1是本发明的技术方案流程图；

图2是应用于本发明实施例的遥感NDVI时间序列数据及其异常检测过程；

图3a是应用于本发明实施例的试验区部分时相的NDVI影像；

图3b是根据本发明的遥感时间序列异常检测方法对试验区NDVI时间序列异常检测的部分结果图。

具体实施方式

现在结合附图和具体实施例，进一步阐明本发明，应理解具体实施例仅用于说明本发明而不用于限定本发明，对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

图1是本发明的遥感时间序列异常检测方法技术方案流程图，结合“发明内容”中的具体描述，本发明的具体实施方式包括以下步骤：

A.从遥感时间序列影像中提取某像素对应的时间序列，用{Y_t：t＝1，2，...，n}表示，时间序列的长度为n，季节周期为s；

B.计算{Y_t}的一阶季节差分序列(季节误差序列){a_t：t＝s+1，s+2，...，n}：

a_t＝Y_t-Y_t-s，t＞s；

C.对季节误差序列{a_t}的均值μ和标准差σ进行稳健化估计：

μ＝median(a_t)，

D.对{a_t}进行z标准化，得到标准分数序列{z_t：t＝s+1，s+2，...，n}：

用{z_t}的值作为时间序列异常检验统计量{λ_t：t＝s+1，s+2，...，n}，即λ_t＝z_t；

E.计算检验统计量临界值：

a)若检验已知异常时刻t₀上的检验统计量是否为异常值，在给定的显著性水平α下，检验统计量临界值的取值为当其满足以下条件时：

或

b)当检验未知时刻t上检验统计量λ_t是否为异常值时，需对全部λ_t进行检验，属于多重假设检验情况，需进行多重检验总体误差率控制。例如一种简单而保守的控制方法为Bonferroni校正，设序列λ_t的长度为m＝n-s，则检验统计量临界值的取值为当其满足以下条件时：

F.将异常检验统计量与上步骤中的临界值相比，满足特定条件时判定为异常。异常的判定条件为：

当s+1≤t≤n-s时，

当t＞n-s时，

满足以上条件的时间t记为异常时间T；

G.获得异常的相关信息，计算如下：

异常水平L_t＝a_t，t∈T，

异常程度D_t＝λ_t，t∈T，

异常置信度C_t＝1-P(|λ|＞|λ_t|)，t∈T。

H.输出遥感时间序列异常检测结果{Y_t，L_t，D_t，C_t：t∈T}

其中，步骤C中，季节误差序列{a_t}的均值μ和标准差σ可以由定义直接计算，也可以进行其它不同方式的稳健化估计。

步骤E和G中，时间序列异常检验统计量{λ_t}的分布假设为独立同分布(i.i.d.)且服从标准正态分布λ～N(0，1)，也可假设为服从t分布。

步骤E(b)中，多重检验总体误差率控制也可以使用其它方法，如holm、hochberg、hommel、FDR、BY控制法以及它们改进的方法。在给定显著性水平时，不同的控制方法在对检验错误率进行限制的情况下，对检测率进行不同程度的提高。

图2是应用于本发明实施例的遥感NDVI时间序列数据及其异常检测过程，下面根据实施例和图2对本发明的具体实施方式进行详细说明：

A.数据来源为MODIS卫星16天合成的250米分辨率NDVI时间序列影像，MODIS行列号为h26v04，覆盖中国东北和华北部分区域。影像时相从2000年2月至2014年2月，每年23景，共323景。提取一例NDVI时间序列，用{Y_t：t＝1，2，...，n}表示，时间序列的长度为n＝323，季节周期为s＝23；

B.计算{Y_t}的一阶季节差分序列(季节误差序列){a_t：t＝24，25，...，323}：

a_t＝Y_t-Y_t-23，t＞23；

C.对季节误差序列{a_t}的均值μ和标准差σ进行一种稳健化估计：

μ＝median(a_t)＝-0.004，

D.对{a_t}进行z标准化，得到标准分数序列{z_t：t＝24，25，...，323}：

用{z_t}的值作为时间序列异常检验统计量{λ_t：t＝24，25，...，323}，即λ_t＝z_t；

E.计算检验统计量临界值：

a)若检验已知异常时刻t₀上的检验统计量是否为异常值，在给定的显著性水平α＝0.05下，检验统计量临界值的取值为当其满足以下条件时：

或

计算得

b)当检验未知时刻t上检验统计量λ_t是否为异常值时，需对全部λ_t进行检验，属于多重假设检验情况，需进行多重检验总体误差率控制。使用一种简单而保守的控制方法Bonferroni校正，序列λ_t的长度为m＝n-s＝300，则检验统计量临界值的取值为当其满足以下条件时：

计算得

F.将异常检验统计量与上步骤中的临界值相比，满足特定条件时判定为异常。异常的判定条件为：

当24≤t≤300时，

当t＞300时，|λ_t|＞3.765，

满足以上条件的时间t记为异常时间T，T＝{217，218，219，310，311，312，313，314}；

G.获得异常的相关信息，计算如下：

异常水平L_t＝a_t，t∈T，L_t＝{-0.39，-0.47，-0.43，-0.38-0.53，-0.63，-0.59，-0.40}，

异常程度D_t＝λ_t，t∈T，D_t＝{-4.46，-5.40，-4.93，-4.34，-6.16，-7.21，-6.79，-4.63}，

异常置信度C_t＝1-P(|λ|＞|λ_t|)，t∈T，

C_t＝{99.99％，99.99％，99.99％，99.99％，100％，100％，100％，99.99％}；

H.输出遥感时间序列异常检测结果{Y_t，L_t，D_t，C_t：t∈T}。

图2显示了本实施例的主要步骤和数据，以及检测出的异常点。可以清楚的看出，利用本发明的遥感时间序列异常检测方法可以有效、准确地检测出遥感时间序列数据中的异常。

下面对本发明实施例的具体应用进行说明。

应用本发明实施例将黑龙江省黑龙江同江段区域作为试验区。由于2014年夏季黑龙江出现特大洪水并且2014年8月23日同江段出现决堤，造成试验区部分区域被洪水淹没。使用本实施例提供的方法对试验区的MODIS卫星16天合成的250米分辨率NDVI时间序列影像(影像时相从2000年2月至2014年2月，每年23景)进行异常检测，以检测试验区的异常淹没区域。

图3a显示了应用于本发明实施例的试验区部分时相的NDVI影像，图3b是根据本发明的遥感时间序列异常检测方法对试验区NDVI时间序列异常检测的部分结果图。

以上仅为本发明的优选实施例而已，并不用于限制本发明。对于本领域的技术人员来说，本发明可以有各种更改和变化。但凡在本发明的基本思路和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种遥感时间序列异常检测方法，该方法可对遥感时间序列中存在的多个时间的异常同时进行检测，同时给出异常的时间、水平、程度和置信度信息，其特征在于，包括以下步骤：

A.从遥感时间序列影像中逐像素地提取时间序列数据，将时间序列表示为{Y_t：t＝1，2，…，n}，时间序列的长度为n，季节周期为s；

B.由时间序列{Y_t}计算季节误差序列{a_t：t＝s+1，s+2，…，n}；

C.对季节误差序列{a_t}的均值μ和标准差σ进行估计；

D.对季节误差序列{a_t}进行变换，得到时间序列的异常检验统计量{λ_t：t＝s+1，s+2，…，n}；

E.获得单一异常检验的异常检验统计量临界值λ⁽¹⁾和多个异常检验的异常检验统计量临界值λ^(m)；

F.将异常检验统计量与步骤E中的异常检验统计量临界值相比，满足特定条件时判定为异常，其中异常判定的特定条件为：

当s+1≤t≤n-s时，

当t＞n-s时，|λ_t|＞λ^(m)，

满足以上条件的时间t记为异常时间T；

G.由{Y_t}、{a_t}、{λ_t}获得异常的相关信息，计算如下：

异常水平L_t＝a_t，t∈T，

异常程度D_t＝λ_t，t∈T，

异常置信度C_t＝1-P(|λ|＞|λ_t|)，t∈T；

H.输出遥感时间序列异常检测结果{Y_t，L_t，D_t，C_t：t∈T}。

2.根据权利要求1所述的遥感时间序列异常检测方法，其特征在于，在步骤B中，通过一阶季节差分计算{Y_t}的季节误差序列{a_t：t＝s+1，s+2，…，n}：

a_t＝Y_t-Y_t-s，t＞s。

3.根据权利要求1所述的遥感时间序列异常检测方法，其特征在于，在步骤C中，季节误差序列{a_t}的均值μ和标准差σ根据定义计算，或者进行稳健化估计，估计方式如下：

μ＝median(a_t)，

$\mathrm{\σ}=\sqrt{\frac{2}{\mathrm{\π}}}\·\frac{\mathrm{\Σ}|a_t-\mathrm{\μ}|}{n-s}.$

4.根据权利要求1所述的遥感时间序列异常检测方法，其特征在于，在步骤D中，对季节误差序列{a_t}进行变换的方式为对{a_t}进行z标准化：

$z_t=\frac{a_t-\mathrm{\μ}}{\mathrm{\σ}},$

得到标准分数序列{z_t：t＝s+1，s+2，…，n}，以此标准分数序列{z_t}的取值作为时间序列的异常检验统计量{λ_t：t＝s+1，s+2，…，n}，即λ_t＝z_t。

5.根据权利要求1所述的遥感时间序列异常检测方法，其特征在于，在步骤E中，异常检验统计量临界值λ⁽¹⁾、λ^(m)为：

a)单一异常检验的异常检验统计量临界值λ⁽¹⁾用于检验已知异常时刻t₀上的检验统计量是否为异常值，λ⁽¹⁾＞0；

b)多个异常检验的异常检验统计量临界值λ^(m)用于检验未知时刻t上检验统计量λ_t是否为异常值，λ^(m)＞0，其中m＝n-s为序列λ_t的长度。

6.根据权利要求5所述的遥感时间序列异常检测方法，其特征在于，异常检验统计量临界值λ⁽¹⁾、λ^(m)为用户预设值，或者通过给定异常检验的显著性水平来进行显著性计算得到：

a)对于λ⁽¹⁾，在给定的显著性水平α下，异常检验统计量临界值的取值满足以下条件：

${\∫}_{-{\mathrm{\λ}}_{\[\mathrm{\α}/2\]}^{\left(1\right)}}^{{\mathrm{\λ}}_{\[a/2\]}^{\left(1\right)}}\frac{1}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{t^2}{2}}dt=1-\frac{\mathrm{\α}}{2}$

或

b)多个异常检验的异常检验统计量临界值λ^(m)用于检验未知时刻t上检验统计量λ_t是否为异常值，λ^(m)＞0，λ^(m)取值为预设值，或者也在给定的显著性水平α下进行显著性计算，对全部λ_t进行显著性检验，进行多重检验总体误差率控制，此时异常检验统计量临界值的取值满足以下条件：

$P\left(\right|\mathrm{\λ}|>{\mathrm{\λ}}_{\[\mathrm{\α}/2m\]}^{\left(m\right)})=\frac{\mathrm{\α}}{2m}.$

7.根据权利要求5所述的遥感时间序列异常检测方法，其特征在于，多重检验总体误差率控制方法为Bonferroni、holm、hochberg、hommel、FDR、BY控制法。