CN104881850A - 一种多时相遥感影像的厚云自动去除方法 - Google Patents

Abstract

本发明公开了一种多时相遥感影像的厚云自动去除方法，由以下步骤组成：步骤1：采集遥感卫星厚云影像和T幅与其同区域的多时相影像；步骤2：厚云区域检测，得到厚云区域指示模板 _；步骤3：自动选择参考影像：自动选择一幅多时相影像作为参考影像；步骤4：采用泊松方程修复方法去除厚云，得到初步去云结果；步骤5：将所述初步去云结果和参考影像纳入变分模型，再次去除厚云，得到最终的去云结果。本发明能通过厚云影像和多幅多时相影像在梯度值之间的均方根误差（）确定参考影像。它无需人工交互，自动地去除厚云及其阴影,它兼顾了原始影像的像素亮度与参考影像的梯度信息，对像素值有较好保真性。

Description

一种多时相遥感影像的厚云自动去除方法

技术领域

本发明涉及一种遥感影像的厚云自动去除方法，尤其是一种多时相遥感影像的厚云自动去除方法，属于遥感影像预处理技术领域。

背景技术

遥感影像中,云覆盖是造成遥感数据缺乏的重要因素之一。大量的遥感影像由于云覆盖的干扰，降低了感兴趣信息的清晰度，从而降低了利用率。有效地减少或去除云的影响，是提高遥感数据利用率的一个重要途径，也是遥感影像预处理中的一个重要问题。遥感影像厚云去除可以恢复不完整的影像，增加遥感影像数据来源、降低数据成本，为促进遥感影像军民两用提供技术支持，以获取较大的经济效益。

由于厚云区域缺乏可利用的信息，直接去除比较困难。目前国内外去除厚云的方法主要有单幅影像的厚云去除方法和多时相影像的厚云去除方法。单幅影像的厚云去除方法主要运用影像修复和合成技术。它虽然能够产生合理的视觉效果，但是无法保证影像像素值的真实性。多时相影像的厚云去除方法是利用同一区域不同时间的影像，根据影像间的时间和空间相关性来去除厚云。但是现有多时相影像厚云去除方法缺乏有效稳定的选择参考影像的标准。同时，这些方法需要很多人工干预，难以同时去除厚云和云的阴影。所以这些方法无法得到实际应用。

发明内容

本发明的目的在于针对已有技术存在的缺陷，提供一种能够有效的选择参考影像，自动去除厚云及其阴影的多时相遥感影像的厚云自动去除方法。

为实现上述目的，本发明采用下述技术方案：

一种多时相遥感影像的厚云自动去除方法，由以下步骤组成：

步骤1：采集遥感卫星厚云影像I_raw(i,j,k),1≤k≤K和T幅与其同区域的多时相影像I_raw,t(i,j,k),1≤i≤M,1≤j≤N,1≤t≤T,1≤k≤K,；其中i和j分别为图像中像素的行坐标和列坐标，k为图像的波段编号，M、N和K分别代表厚云影像的行数、列数和波段数。

步骤2：厚云区域检测，得到厚云区域指示模板mask(i,j)；

步骤3：自动选择参考影像：自动选择一幅多时相影像作为参考影像I_ref(i,j,k)；

步骤4：采用泊松方程修复方法去除厚云，得到初步去云结果I₀(i,j,k)；

步骤5：将所述初步去云结果I₀(i,j,k)和参考影像I_ref(i,j,k)纳入变分模型，再次去除厚云，得到最终的去云结果I'₀(i,j,k)。

上述步骤2由以下具体步骤组成：

步骤2-1：将厚云影像I_raw(i,j,k)的蓝光波段图像I(i,j)与阈值T_C比较，计算厚云区域初始指示模板mask0(i,j)：

$\mathrm{mask}0(i,j)=\left\{\begin{array}{c}1,I(i,j)\&GreaterEqual;T_C\\ 0,I(i,j)<T_C\end{array}\right.,1\&le;i\&le;M,1\&le;j\&le;N---\left(1\right)$

式中mask0(i,j)为1表示对应像素属于厚云区域，mask0(i,j)为0表示对应像素属于无云区域；

步骤2-2：对所述厚云区域初始指示模板mask0(i,j)进行形态学处理，去除孤立的斑点并填补细小的空洞，得到厚云区域指示模板mask(i,j)。

上述步骤2-1中阈值T_C的计算方法为：

T_C＝m+λσ    (2)

其中，m为所述厚云影像I_raw(i,j,k)的蓝光波段图像I(i,j)的像素均值：

$m=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}I(i,j)---\left(3\right)$

σ为厚云影像I_raw(i,j,k)的蓝光波段图像I(i,j)的像素标准方差：

$\mathrm{\&sigma;}=\sqrt{\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{(I(i,j)-m)}^2}---\left(4\right)$

λ为权重系数：

$\mathrm{\&lambda;}=\frac{(\max \left(I\right)-\min \left(I\right))}{2m}---\left(5\right)$

上述步骤3由以下具体步骤组成：

步骤3-1：计算厚云影像的各波段的图像I_raw(i,j,k),1≤k≤K在x方向的梯度值g_x,0(x,y,k)和y方向的梯度值g_y,0(x,y,k)：

$\begin{array}{c}{g}_{x,0}(x,y,k)=I(i+1,j,k)-I(i,j,k)\\ g_{y,0}(x,y,k)=I(i,j+1,k)-I(i,j,k)\end{array},1\&le;i\&le;M,1\&le;j\&le;N,1\&le;k\&le;K---\left(6\right)$

步骤3-2：计算所述各多时相影像I_raw,t(i,j,k),1≤t≤T在x方向的梯度值g_x,t(x,y,k)和y方向的梯度值g_y,t(x,y,k):

$\begin{array}{c}{g}_{x,t}(x,y,k)=I_t(i+1,j,k)-I_t(i,j,k)\\ g_{y,t}(x,y,k)=I_t(i,j+1,k)-I_t(i,j,k)\end{array},1\&le;i\&le;M,1\&le;j\&le;N,1\&le;t\&le;T---\left(7\right)$

步骤3-3：计算所述厚云影像I_raw(i,j,k),1≤k≤K与各多时相影像I_raw,t(i,j,k),1≤t≤T在无云区域对应位置梯度值之间的均方根误差RMSE_t：

${\mathrm{RMSE}}_t=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{(I_{\mathrm{raw}}(i,j,k)-I_t(i,j,k))}^2}{\mathrm{MNK}}},\mathrm{mask}(i,j)>0---\left(8\right)$

步骤3-4：选择RMSE_t最小的多时相影像I_t，作为所述厚云影像I_raw的去云参考影像I_ref(i,j,k)。

上述步骤4中所述初步去云结果I₀(i,j,k)中无云区域像素值与所述厚云影像I_raw(i,j,k)的无云区域像素值相等，其厚云区域图像的像素值由下式计算得到：

$\left\{\begin{array}{c}{I}_0(i,j,k)=\underset{I_1(i,j,k)}{\arg \min }\underset{\mathrm{\&Omega;}}{\&Integral;}{|\&dtri;I_1(i,j,k)-\&dtri;I_{\mathrm{ref}}(i,j,k)|}^2\\ I_0(i,j,k){|}_{\&PartialD;\mathrm{\&Omega;}}=I_{\mathrm{raw}}(i,j,k){|}_{\&PartialD;\mathrm{\&Omega;}}\end{array}\right.,\mathrm{mask}(i,j)>0---\left(9\right)$

其中，Ω代表厚云区域内部；代表厚云区域边界，▽表示梯度算子,I₁表示I₀在区域Ω的像素。

上述步骤4中厚云区域图像的像素值可以根据以下递推公式计算得到：

$\left\{\begin{array}{c}{X}^{(l+1)}=X^{\left(l\right)}+{\mathrm{\&mu;}}^{\left(l\right)}{p}^{\left(l\right)}\\ {\mathrm{\&mu;}}^{\left(l\right)}=\frac{<r^{\left(l\right)},{(C,C^T)}^{-1}\&CenterDot;r^{\left(l\right)}>}{<p^{\left(l\right)},A\&CenterDot;p^{\left(l\right)}>}\\ r^{(l+1)=r^{\left(l\right)}-A\&CenterDot;p^{\left(l\right)}}\\ {\mathrm{\&beta;}}^{\left(l\right)}=\frac{<r^{(l+1)},{(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{(l+1)}>}{<r^{\left(l\right)},{(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{\left(l\right)}>}\\ p^{(l+1)}={(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{(l+1)}+{\mathrm{\&beta;}}^{\left(l\right)}{p}^{\left(l\right)}\\ X^{\left(0\right)}={(C\&CenterDot;C^T)}^{-1}\&CenterDot;B\\ r^{\left(0\right)}=B-A\&CenterDot;X^{\left(0\right)}\\ p^{\left(0\right)}={(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{\left(0\right)}\end{array}\right.---\left(10\right)$

其中

式中p_i和p_j分别代表初步去云结果I₀(i,j,k)中第i个和第j个厚云区域像素；N(p_i)代表所述像素p_i的四邻域像素集；

$B=\left[b_i\right],b_i=\underset{q_i\&Element;N\left(p_i\right)\&cap;\&PartialD;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_0\left(q_i\right)+\underset{q_i\&Element;N\left(p_i\right)}{\mathrm{\&Sigma;}}(I_{\mathrm{ref}}\left(p_i\right)-I_{\mathrm{ref}}\left(q_i\right))---\left(12\right)$

式中I₀(q_i)表示初步去云结果I₀(i,j,k)中像素q_i的像素值，I_ref(p_i)I_ref(q_i)分别表示所述去云参考影像I_ref(i,j,k)中像素p_i和q_i的像素值，像素q_i为N(p_i)与Ω交集内的第i个像素；C为A的Choleskey分解，A＝CC^T，其中C为一个下三角矩阵，其下三角部分具有和A完全相同的结构。

上述步骤5中的变分模型代价函数为：

$J\left[{I^{\&prime;}}_0(i,j,k)\right]=\underset{T}{\&Integral;\&Integral;}\frac{\mathrm{\&alpha;}}{2}{\left(\right|\&dtri;{I^{\&prime;}}_0(i,j,k)-\&dtri;I_{\mathrm{ref}}(i,j,k)\left|\right)}^2+\frac{1}{2}{\left(\right|{I^{\&prime;}}_0(i,j,k)-I_0(i,j,k)\left|\right)}^2\mathrm{dxdy}---\left(13\right)$

其中α为用于权衡两个约束项贡献度参数，α≥0。

上述步骤5中的变分模型代价函数最小值的求解方法为采用下述迭代公式计算：

$\left\{\begin{array}{c}{\left({I^{\&prime;}}_0(i,j,k)\right)}^{n+1}={\left({I^{\&prime;}}_0(i,j,k)\right)}^n-\mathrm{\&Delta;t}[{\left({I^{\&prime;}}_0(i,j,k)\right)}^n-I_0(i,j,k)-{\mathrm{\&alpha;}}^n(\&dtri;\&CenterDot;{\left({I^{\&prime;}}_0(i,j,k)\right)}^n)\\ -\&dtri;\&CenterDot;\left(I_{\mathrm{ref}}(i,j,k)\right)\left)\right]\\ {I^{\&prime;}}_0(i,j,k)=I_0(i,j,k)\end{array}\right.---\left(14\right)$

其中Δt是迭代步长，α是正则化参数，αⁿ的计算公式为：

${\mathrm{\&alpha;}}^n=\frac{\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\&dtri;\&CenterDot;\left(\frac{\&dtri;{\left(I_0^{\&prime;}\right)}_{i,j}}{|\&dtri;{\left(I_0^{\&prime;}\right)}_{i,j}|}\right)(I_0^{\&prime;}(i,j)-I_0(i,j))}{\mathrm{var}\_n}---\left(15\right)$

式中var_n为设定的标准参数。

本发明的有益效果在于：

1、本发明能通过厚云影像和多幅多时相影像在梯度值之间的均方根误差(RMSE)确定参考影像。它无需人工交互，自动地去除厚云及其阴影

2、本发明兼顾了原始影像的像素亮度与参考影像的梯度信息，对像素值有较好保真性。

附图说明

图1为本发明的流程图；

图2为本发明步骤2的流程图；

图3为本发明步骤3的流程图。

图4为本发明的遥感卫星采集的厚云影像；

图5为本发明的遥感卫星采集的参考影像；

图6为本发明检测出的厚云区域的指示模板；

图7为本发明的初步去云结果；

图8为本发明的最终去云结果。

具体实施方式

实施例1：

如图1所示，一种多时相遥感影像的厚云自动去除方法，由以下步骤组成：

步骤1：采集遥感卫星厚云影像I_raw(i,j,k),1≤k≤K和T幅与其同区域的多时相影像I_raw,t(i,j,k),1≤i≤M,1≤j≤N,1≤t≤T,1≤k≤K,；其中i和j分别为图像中像素的行坐标和列坐标，k为图像的波段编号，M、N和K分别代表厚云影像的行数、列数和波段数；；在本实施例中M、N和K分别为400、400和7；

步骤2：厚云区域检测，得到厚云区域指示模板mask(i,j)；

步骤3：自动选择参考影像：自动选择一幅多时相影像作为参考影像I_ref(i,j,k)；

步骤4：采用泊松方程修复方法去除厚云，得到初步去云结果I₀(i,j,k)；

步骤5：将所述初步去云结果I₀(i,j,k)和参考影像I_ref(i,j,k)纳入变分模型，再次去除厚云，得到最终的去云结果I'₀(i,j,k)；

如图2所示，上述步骤2由以下具体步骤组成：

步骤2-1：将厚云影像I_raw(i,j,k)的蓝光波段图像I(i,j)与阈值T_C比较，计算厚云区域初始指示模板mask0(i,j)：

$\mathrm{mask}0(i,j)=\left\{\begin{array}{c}1,I(i,j)\&GreaterEqual;T_C\\ 0,I(i,j)<T_C\end{array}\right.,1\&le;i\&le;M,1\&le;j\&le;N---\left(1\right)$

式中mask0(i,j)为1表示对应像素属于厚云区域，mask0(i,j)为0表示对应像素属于无云区域；

步骤2-2：对所述厚云区域初始指示模板mask0(i,j)进行形态学处理，去除孤立的斑点并填补细小的空洞，得到厚云区域指示模板mask(i,j)。

上述步骤2-1中阈值T_C的计算方法为：

T_C＝m+λσ    (2)

其中，m为所述厚云影像I_raw(i,j,k)的蓝光波段图像I(i,j)的像素均值：

$m=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}I(i,j)---\left(3\right)$

σ为厚云影像I_raw(i,j,k)的蓝光波段图像I(i,j)的像素标准方差：

$\mathrm{\&sigma;}=\sqrt{\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{(I(i,j)-m)}^2}---\left(4\right)$

λ为权重系数：

$\mathrm{\&lambda;}=\frac{(\max \left(I\right)-\min \left(I\right))}{2m}---\left(5\right)$

如图3所示，上述步骤3由以下具体步骤组成：

步骤3-1：计算厚云影像的各波段的图像I_raw(i,j,k),1≤k≤K在x方向的梯度值g_x,0(x,y,k)和y方向的梯度值g_y,0(x,y,k)：

$\begin{array}{c}{g}_{x,0}(x,y,k)=I(i+1,j,k)-I(i,j,k)\\ g_{y,0}(x,y,k)=I(i,j+1,k)-I(i,j,k)\end{array},1\&le;i\&le;M,1\&le;j\&le;N,1\&le;k\&le;K---\left(6\right)$

步骤3-2：计算所述各多时相影像I_raw,t(i,j,k),1≤t≤T在x方向的梯度值g_x,t(x,y,k)和y方向的梯度值g_y,t(x,y,k):

$\begin{array}{c}{g}_{x,t}(x,y,k)=I_t(i+1,j,k)-I_t(i,j,k)\\ g_{y,t}(x,y,k)=I_t(i,j+1,k)-I_t(i,j,k)\end{array},1\&le;i\&le;M,1\&le;j\&le;N,1\&le;t\&le;T---\left(7\right)$

步骤3-3：计算所述各厚云影像I_raw(i,j,k),1≤k≤K与各多时相影像I_raw,t(i,j,k),1≤t≤T在无云区域对应位置梯度值之间的均方根误差RMSE_t：

${\mathrm{RMSE}}_t=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{(I_{\mathrm{raw}}(i,j,k)-I_t(i,j,k))}^2}{\mathrm{MNK}}},\mathrm{mask}(i,j)>0---\left(8\right)$

步骤3-4：选择RMSE_t最小的多时相影像I_t，作为所述厚云影像I_raw的去云参考影像I_ref(i,j,k)。

上述步骤4由以下具体步骤组成：

步骤4-1：所述初步去云结果I₀(i,j,k)中无云区域像素值与所述厚云影像I_raw(i,j,k)的无云区域像素值相等，其厚云区域图像的像素值由下式计算得到：

$\left\{\begin{array}{c}{I}_0(i,j,k)=\underset{I_1(i,j,k)}{\arg \min }\underset{\mathrm{\&Omega;}}{\&Integral;}{|\&dtri;I_1(i,j,k)-\&dtri;I_{\mathrm{ref}}(i,j,k)|}^2\\ I_0(i,j,k){|}_{\&PartialD;\mathrm{\&Omega;}}=I_{\mathrm{raw}}(i,j,k){|}_{\&PartialD;\mathrm{\&Omega;}}\end{array}\right.,\mathrm{mask}(i,j)>0---\left(9\right)$

其中，Ω代表厚云区域内部；代表厚云区域边界，▽表示梯度算子,I₁表示I₀在区域Ω的像素。

上述步骤4-1中厚云区域图像的像素值可以根据以下递推公式计算得到：

$\left\{\begin{array}{c}{X}^{(l+1)}=X^{\left(l\right)}+{\mathrm{\&mu;}}^{\left(l\right)}{p}^{\left(l\right)}\\ {\mathrm{\&mu;}}^{\left(l\right)}=\frac{<r^{\left(l\right)},{(C,C^T)}^{-1}\&CenterDot;r^{\left(l\right)}>}{<p^{\left(l\right)},A\&CenterDot;p^{\left(l\right)}>}\\ r^{(l+1)=r^{\left(l\right)}-A\&CenterDot;p^{\left(l\right)}}\\ {\mathrm{\&beta;}}^{\left(l\right)}=\frac{<r^{(l+1)},{(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{(l+1)}>}{<r^{\left(l\right)},{(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{\left(l\right)}>}\\ p^{(l+1)}={(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{(l+1)}+{\mathrm{\&beta;}}^{\left(l\right)}{p}^{\left(l\right)}\\ X^{\left(0\right)}={(C\&CenterDot;C^T)}^{-1}\&CenterDot;B\\ r^{\left(0\right)}=B-A\&CenterDot;X^{\left(0\right)}\\ p^{\left(0\right)}={(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{\left(0\right)}\end{array}\right.---\left(10\right)$

其中

式中p_i和p_j分别代表初步去云结果I₀(i,j,k)中第i个和第j个厚云区域像素；N(p_i)代表所述像素p_i的四邻域像素集；

$B=\left[b_i\right],b_i=\underset{q_i\&Element;N\left(p_i\right)\&cap;\&PartialD;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_0\left(q_i\right)+\underset{q_i\&Element;N\left(p_i\right)}{\mathrm{\&Sigma;}}(I_{\mathrm{ref}}\left(p_i\right)-I_{\mathrm{ref}}\left(q_i\right))---\left(12\right)$

式中I₀(q_i)表示初步去云结果I₀(i,j,k)中像素q_i的像素值，I_ref(p_i)I_ref(q_i)分别表示所述去云参考影像I_ref(i,j,k)中像素p_i和q_i的像素值，像素q_i为N(p_i)与Ω交集内的第i个像素；C为A的Choleskey分解，A＝CC^T，其中C为一个下三角矩阵，其下三角部分具有和A完全相同的结构。

上述步骤5中的变分模型代价函数为：

$J\left[{I^{\&prime;}}_0(i,j,k)\right]=\underset{T}{\&Integral;\&Integral;}\frac{\mathrm{\&alpha;}}{2}{\left(\right|\&dtri;{I^{\&prime;}}_0(i,j,k)-\&dtri;I_{\mathrm{ref}}(i,j,k)\left|\right)}^2+\frac{1}{2}{\left(\right|{I^{\&prime;}}_0(i,j,k)-I_0(i,j,k)\left|\right)}^2\mathrm{dxdy}---\left(13\right)$

其中α为用于权衡两个约束项贡献度参数，α≥0。

上述步骤5中的变分模型代价函数最小值的求解方法为采用下述迭代公式计算：

$\left\{\begin{array}{c}{\left({I^{\&prime;}}_0(i,j,k)\right)}^{n+1}={\left({I^{\&prime;}}_0(i,j,k)\right)}^n-\mathrm{\&Delta;t}[{\left({I^{\&prime;}}_0(i,j,k)\right)}^n-I_0(i,j,k)-{\mathrm{\&alpha;}}^n(\&dtri;\&CenterDot;{\left({I^{\&prime;}}_0(i,j,k)\right)}^n)\\ -\&dtri;\&CenterDot;\left(I_{\mathrm{ref}}(i,j,k)\right)\left)\right]\\ {I^{\&prime;}}_0(i,j,k)=I_0(i,j,k)\end{array}\right.---\left(14\right)$

其中Δt是迭代步长，α是正则化参数，αⁿ的计算公式为：

${\mathrm{\&alpha;}}^n=\frac{\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\&dtri;\&CenterDot;\left(\frac{\&dtri;{\left(I_0^{\&prime;}\right)}_{i,j}}{|\&dtri;{\left(I_0^{\&prime;}\right)}_{i,j}|}\right)(I_0^{\&prime;}(i,j)-I_0(i,j))}{\mathrm{var}\_n}---\left(15\right)$

式中var_n为设定的标准参数。在本实施例中，var_n为519.84。

Claims (8)

1.一种多时相遥感影像的厚云自动去除方法，其特征在于：由以下步骤组成：

步骤1：采集遥感卫星厚云影像I_raw(i,j,k),1≤k≤K和T幅与其同区域的多时相影像I_raw,t(i,j,k),1≤i≤M,1≤j≤N,1≤t≤T,1≤k≤K,；其中i和j分别为图像中像素的行坐标和列坐标，k为图像的波段编号，M、N和K分别代表厚云影像的行数、列数和波段数；

步骤2：厚云区域检测，得到厚云区域指示模板mask(i,j)；

步骤3：自动选择参考影像：自动选择一幅多时相影像作为参考影像I_ref(i,j,k)；

步骤4：采用泊松方程修复方法去除厚云，得到初步去云结果I₀(i,j,k)；

步骤5：将所述初步去云结果I₀(i,j,k)和参考影像I_ref(i,j,k)纳入变分模型，再次去除厚云，得到最终的去云结果I'₀(i,j,k)。

2.根据权利要求1所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步骤2由以下具体步骤组成：

步骤2-1：将厚云影像I_raw(i,j,k)的蓝光波段图像I(i,j)与阈值T_C比较，计算厚云区域初始指示模板mask0(i,j)：

$\mathrm{mask}0(i,j)=\left\{\begin{array}{c}1,I(i,j)\&GreaterEqual;T_C\\ 0,I(i,j)<T_C\end{array}\right.,1\&le;i\&le;M,1\&le;j\&le;N---\left(1\right)$

式中mask0(i,j)为1表示对应像素属于厚云区域，mask0(i,j)为0表示对应像素属于无云区域；

步骤2-2：对所述厚云区域初始指示模板mask0(i,j)进行形态学处理，去除孤立的斑点并填补细小的空洞，得到厚云区域指示模板mask(i,j)。

3.根据权利要求2所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步骤2-1中阈值T_C的计算方法为：

T_C＝m+λσ   (2)

其中，m为厚云影像I_raw的蓝光波段图像I的像素均值：

$m=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}I(i,j)---\left(3\right)$

σ为厚云影像I_raw的蓝光波段图像I的像素标准方差：

$\mathrm{\&sigma;}=\sqrt{\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{(I(i,j)-m)}^2}---\left(4\right)$

λ为权重系数：

$\mathrm{\&lambda;}=\frac{(\max \left(I\right)-\min \left(I\right))}{2m}---\left(5\right).$

4.根据权利要求1所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步骤3由以下具体步骤组成：

步骤3-1：计算厚云影像的各波段的图像I_raw(i,j,k),1≤k≤K在x方向的梯度值g_x,0(x,y,k)和y方向的梯度值g_y,0(x,y,k)：

g_x,0(x,y,k)＝I(i+1,j,k)-I(i,j,k)

                                ,1≤i≤M,1≤j≤N,1≤k≤K   (6)

g_y,0(x,y,k)＝I(i,j+1,k)-I(i,j,k)

步骤3-2：计算所述各多时相影像I_raw,t(i,j,k),1≤t≤T在x方向的梯度值g_x,t(x,y,k)和y方向的梯度值g_y,t(x,y,k):

g_x，t(x,y,k)＝I_t(i+1,j,k)-I_t(i,j,k)

                                  ,1≤i≤M,1≤j≤N,1≤t≤T   (7)

g_y,t(x,y,k)＝I_t(i,j+1,k)-I_t(i,j,k)

步骤3-3：计算所述各厚云影像I_raw(i,j,k),1≤k≤K与各多时相影像I_raw,t(i,j,k),1≤t≤T在

无云区域对应位置梯度值之间的均方根误差RMSE_t：

${\mathrm{RMSE}}_t=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{(I_{\mathrm{raw}}(i,j,k)-I_t(i,j,k))}^2}{\mathrm{MNK}}},\mathrm{mask}(i,j)>0---\left(8\right)$

步骤3-4：选择RMSE_t最小的多时相影像I_t，作为所述厚云影像I_raw的去云参考影像I_ref(i,j,k)。

5.根据权利要求1所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步骤4中所述初步去云结果I₀(i,j,k)中无云区域像素值与所述厚云影像I_raw(i,j,k)的无云区域像素值相等，其厚云区域图像的像素值由下式计算得到：

$\left\{\begin{array}{c}{I}_0(i,j,k)=\underset{I_1(i,j,k)}{\arg \min }\underset{\mathrm{\&Omega;}}{\&Integral;}{|\&dtri;I_1(i,j,k)-\&dtri;I_{\mathrm{ref}}(i,j,k)|}^2\\ I_0(i,j,k){|}_{\&PartialD;\mathrm{\&Omega;}}=I_{\mathrm{raw}}(i,j,k){|}_{\&PartialD;\mathrm{\&Omega;}}\end{array}\right.,\mathrm{mask}(i,j)>0$

其中，Ω代表厚云区域内部；代表厚云区域边界，▽表示梯度算子,I₁表示I₀在区域Ω的像素。

6.根据权利要求5所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步骤4中厚云区图像的像素值根据以下递推公式计算得到：

$\left\{\begin{array}{c}{X}^{(l+1)}=X^{\left(l\right)}+{\mathrm{\&mu;}}^{\left(l\right)}{p}^{\left(l\right)}\\ {\mathrm{\&mu;}}^{\left(l\right)}=\frac{<r^{\left(l\right)},{(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{\left(l\right)}>}{<p^{\left(l\right)},A\&CenterDot;p^{\left(l\right)}>}\\ r^{(l+1)}=r^{\left(l\right)}-A\&CenterDot;p^{\left(l\right)}\\ {\mathrm{\&beta;}}^{\left(l\right)}=\frac{<r^{(l+1)},{(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{(l+1)}>}{<r^{\left(l\right)},{(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{\left(l\right)}>}\\ p^{(l+1)}={(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{(l+1)}+{\mathrm{\&beta;}}^{\left(l\right)}{p}^{\left(l\right)}\\ X^{\left(0\right)}={(C\&CenterDot;C^T)}^{-1}\&CenterDot;B\\ r^{\left(0\right)}=B-A\&CenterDot;X^{\left(0\right)}\\ p^{\left(0\right)}={(C\&CenterDot;C^T)}^{-1}\&CenterDot;r^{\left(0\right)}\end{array}\right.---\left(10\right)$

其中

式中p_i和p_j分别代表初步去云结果I₀(i,j,k)中第i个和第j个厚云区域像素；N(p_i)代表所述像素p_i的四邻域像素集；

$B=\left[b_i\right],b_i=\underset{q_i\&Element;N\left(p_i\right)\&cap;\&PartialD;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_0\left(q_i\right)+\underset{q_i\&Element;N\left(p_i\right)}{\mathrm{\&Sigma;}}(I_{\mathrm{ref}}\left(p_i\right)-I_{\mathrm{ref}}\left(q_i\right))---\left(12\right)$

式中I₀(q_i)表示初步去云结果I₀(i,j,k)中像素q_i的像素值，I_ref(p_i)I_ref(q_i)分别表示所述去云参考影像I_ref(i,j,k)中像素p_i和q_i的像素值，像素q_i为N(p_i)与Ω交集内的第i个像素；C为A的Choleskey分解，A＝CC^T，其中C为一个下三角矩阵，其下三角部分具有和A完全相同的结构。

7.根据权利要求1所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步骤5中的变分模型代价函数为：

$J\left[{I^{\&prime;}}_0(i,j,k)\right]=\underset{\mathrm{\&Gamma;}}{\&Integral;\&Integral;}\frac{\mathrm{\&alpha;}}{2}{\left(\right|\&dtri;{I^{\&prime;}}_0(i,j,k)-\&dtri;I_{\mathrm{ref}}(i,j,k)\left|\right)}^2+\frac{1}{2}{\left(\right|{I^{\&prime;}}_0(i,j,k){-I}_0(i,j,k)\left|\right)}^2\mathrm{dxdy}---\left(13\right)$

其中α为用于权衡两个约束项贡献度参数，α≥0。

8.根据权利要求7所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步骤5中的变分模型代价函数最小值的求解方法为采用下述迭代公式计算：

$\left\{\begin{array}{c}{\left({I^{\&prime;}}_0(i,j,k)\right)}^{n+1}={\left({I^{\&prime;}}_0(i,j,k)\right)}^n-\mathrm{\&Delta;t}[{\left({I^{\&prime;}}_0(i,j,k)\right)}^n-I_0(i,j,k)-{\mathrm{\&alpha;}}^n(\&dtri;\&CenterDot;{\left(\right({I^{\&prime;}}_0(i,j,k))}^n)\\ -\&dtri;\&CenterDot;\left(I_{\mathrm{ref}}(i,j,k)\right)\left)\right]\\ {I^{\&prime;}}_0(i,j,k)=I_0(i,j,k)\end{array}\right.$

其中Δt是迭代步长，α是正则化参数，αⁿ的计算公式为：

${\mathrm{\&alpha;}}^n=\frac{\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\&dtri;\&CenterDot;\left(\frac{\&dtri;{\left(I_0^{\&prime;}\right)}_{i,j}}{|\&dtri;{\left(I_0^{\&prime;}\right)}_{i,j}|}\right)(I_0^{\&prime;}(i,j)-I_0(i,j))}{\mathrm{var}\_n}---\left(15\right)$

式中var_n为设定的标准参数。