基于p范数回归模型的光学遥感图像时间序列云检测方法与流程

本发明涉及云检测领域，尤其是涉及基于p范数回归模型的光学遥感图像时间序列云检测方法。

背景技术：

传统的云检测算法是通过设置像素的阈值、使用像素周围的信息、使用空间信息等方法。但是这些方法过于简单且有局限性，难以广泛的使用。

很多传统的去云算法都是针对特定场景的，当脱离这个场景后，就不适用了，而且没办法检测到不同的云，比如说当云为薄云时，就有可能检测不出来。

技术实现要素：

本发明的目的在于：针对现有技术存在的问题，提供一种基于p范数回归模型的光学遥感图像时间序列云检测方法，解决传统云检测方法有局限性，只能针对特定场景检测的问题。

本发明的发明目的通过以下技术方案来实现：

基于p范数回归模型的光学遥感图像时间序列云检测方法，该方法包括步骤：

(1)预处理：选择感兴趣的区域；

(2)基于p范数的云检测回归建模；

(3)通过观测值与回归值的偏差检测云。

步骤(1)包括：获得固定区域一定时间内如一年的图像，本专利中采用22个时刻采集得到的光学卫星Landsat-8OLI图像，用D＝{S1，S2...，S22}表示一个区域的子图像数据集，用Dm，n＝{S1(m，n)，S2(m，n)...，S22(m，n)}来表示一个处于(m，n)处22个时刻的像素值序列。

步骤(2)包括：

1)假定：yi＝f(xi)+εi，i＝1，2...，n，其中f(*)是回归函数，yi是在采样时刻xi处给出的一个像素的反射测量，n是采样的数量；

2)假设εi遵循p范数分布来模拟云的影响，且服从参数为期望μ，方差σ²和正实数p，则εi的密度为：其中Γ(.)是Gamma函数，E(εi)＝μ；

3)假设期望μ等于零，在服从p范数分布的噪声εi的影响下，从测量值yi中找到一种的未知的像素点的变化f(*)。

步骤(2)包括：假设f(*)具有3阶局部平滑性；假设xi是在x附近的采样时刻，则表达式f(xi)的泰勒展开式写为：

其中，

β0＝f(x)，β1＝f′(x)，β2＝f″(x)，β3＝f″′(x)，β＝[β0，β1，β2，β3]^T，a＝[1，(xi-x)，(xi-x)²，(xi-x)³].因此，有：

又由于εi服从p范数分布，因此yi的密度表示成：

并且样本的似然函数：

则β的似然估计的最大值由下式计算得到：通过求解这个优化问题，得到回归系数

步骤(3)包括：像素的真实值通过如下的回归方程估计：

步骤(4)包括：给定回归模型，用s表示误差或内点的标度，通过下式来检测云：其中T是门限值，

T为2.5。

与现有技术相比，本发明方法简单，检测结果可靠，能将非常薄的云检测出来。

附图说明

图1为本发明的方法流程图；

图2在m，n像素点上的观测时间序列Dm，n和估计的回归函数f(*)；

图3为原始遥感图像；

图4为云检测结果图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

实施例

如图1～图4所示，本发明方法包括：

1.预处理。

2.基于p范数的云检测回归建模。

3.云检测。

第1个步骤预处理：

我们首先在Landsat-8OLI数据集中对所有图像进行配准，以获得固定区域一年的有云图像。在22个时刻中，用D＝{S1，S2...，S22}表示一个区域的子图像数据集,用Dm，n＝{S1(m，n)，S2(m，n)...，S22(m，n)}来表示处于(m，n)处22个时刻的像素值序列。

第2个步骤基于p范数的云检测回归建模：

基于云及其阴影总是在时间序列上引起像素反射率突然变化的现象，我们可以将它们视为具有p范数的噪声分布。p范数分布的噪声是单峰并且对称的，通过选择适当的p值，理论模型可以产生只有少量大的非零值和大量接近于零的较小值的噪声。因此对于我们的问题，这种分布的拟合程度要比正态分布更优。

对于给定有云图像的时间序列来寻找每个像素对应的变化模型，这实际上是一种在p范数分布噪声下的鲁棒回归问题。

假定：

yi＝f(xi)+εi，i＝1，2...，n (1)

其中f(*)是回归函数，yi是在采样时刻xi处反射率的观测值，n是时间序列长度。在一个确定的采样时刻，当它被云污染时，测量值yi可能偏离像素真实的反射率。假设εi服从期望为μ、方差为σ²的p范数分布来模拟云的影响，p为正实数。那么εi的密度函数可以表示为：

其中Γ(.)为Gamma函数，E(εi)＝μ。为了简单起见，我们可以假设期望μ等于零。从而问题变为：我们要在服从p范数分布的噪声εi的影响下，从测量值yi中找到一种未知的像素值的变化函数f(*)。

由于像素的反射率总是随着时间序列而缓慢而连续地变化，我们可以假设f(*)具有3阶局部平滑性。假设xi是在x附近的采样时刻，则表达式f(xi)的泰勒展开式可写为：

其中，

β0＝f(x)，β1＝f′(x)，β2＝f″(x)，β3＝f″′(x)，β=[β0，β1，β2，β3]^T，α＝[1，(xi-x)，(xi-x)²，(xi-x)³].

因此，我们有：

又由于εi服从p范数分布，因此yi的密度函数可以表示成：

并且样本的似然函数为

则β的似然估计的最大值，即可以由下式计算得到：

通过求解这个优化问题，我们可以得到回归系数于是，像素的真实值可以通过如下的回归方程估计得到

第3个步骤云检测：

通过观测值与回归值的偏差检测云。对每个像素点上的像素值时间序列与回归值相减得到一个误差序列，该误差可以有效地模拟为高斯噪声，而异常值则是与云对应的。给定回归模型，用s表示正常值的尺度，则云可以通过下式来检测：

其中T是门限值。通常情况下，选择T为2.5，可以将98％的高斯分布点识别为正常值，这也是在我们的实验中使用的。为了有效地区分异常值和高斯噪声，我们必须正确地估计正常值的尺度。我们使用鲁棒的中值尺度估计来获得s的值如下：

其中n是样本点的数量。

图2给出了m，n处的Dm，n序列和估计得到的回归函数曲线结果，其中鲁棒回归函数用水平方向的曲线表示，有云观测值用”*”表示。注意，时刻13、14处像素值分别为79、76，这两个点是异常值，它们对应于厚云；而在时刻17、21处像素值分别为39、38，这两个异常值对应于薄云。图3为有云图像，图4为用本专利方法检测得到的云图像。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，应当指出的是，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。