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

一种多时相遥感影像的厚云自动去除方法
【技术领域】
[0001] 本发明设及一种遥感影像的厚云自动去除方法，尤其是一种多时相遥感影像的厚 云自动去除方法，属于遥感影像预处理技术领域。
【背景技术】
[0002] 遥感影像中，云覆盖是造成遥感数据缺乏的重要因素之一。大量的遥感影像由于 云覆盖的干扰，降低了感兴趣信息的清晰度，从而降低了利用率。有效地减少或去除云的影 响，是提高遥感数据利用率的一个重要途径，也是遥感影像预处理中的一个重要问题。遥感 影像厚云去除可W恢复不完整的影像，增加遥感影像数据来源、降低数据成本，为促进遥感 影像军民两用提供技术支持，W获取较大的经济效益。
[000引由于厚云区域缺乏可利用的信息，直接去除比较困难。目前国内外去除厚云的方 法主要有单幅影像的厚云去除方法和多时相影像的厚云去除方法。单幅影像的厚云去除方 法主要运用影像修复和合成技术。它虽然能够产生合理的视觉效果，但是无法保证影像像 素值的真实性。多时相影像的厚云去除方法是利用同一区域不同时间的影像，根据影像间 的时间和空间相关性来去除厚云。但是现有多时相影像厚云去除方法缺乏有效稳定的选择 参考影像的标准。同时，该些方法需要很多人工干预，难W同时去除厚云和云的阴影。所W 该些方法无法得到实际应用。

【发明内容】

[0004] 本发明的目的在于针对已有技术存在的缺陷，提供一种能够有效的选择参考影 像，自动去除厚云及其阴影的多时相遥感影像的厚云自动去除方法。
[0005] 为实现上述目的，本发明采用下述技术方案：
[0006] -种多时相遥感影像的厚云自动去除方法，由W下步骤组成：
[0007] 步骤1;采集遥感卫星厚云影像Ir"(i,j,k), 1《k《K和T幅与其同区域的多时 相影像Iraw.t(1，j，k)，1《i《M，1《j《N，1《t《T，1《k《K，；其中i和j分别为图 像中像素的行坐标和列坐标，k为图像的波段编号，M、N和K分别代表厚云影像的行数、列 数和波段数。
[000引步骤2 ;厚云区域检测，得到厚云区域指示模板mask(i，如；
[0009] 步骤3;自动选择参考影像；自动选择一幅多时相影像作为参考影像Iuf(i，j，k);
[0010] 步骤4 ;采用泊松方程修复方法去除厚云，得到初步去云结果I。(i，j，k);
[0011] 步骤5 ;将所述初步去云结果1。(1，j，k)和参考影像luf(i，j，k)纳入变分模型，再 次去除厚云，得到最终的去云结果I'。(i，j，k)。
[0012] 上述步骤2由w下具体步骤组成：
[001引步骤2-1 ;将厚云影像j，k)的藍光波段图像I(iJ)与阔值Te比较，计算厚 云区域初始指示模板mask0(i，j);
[0014] (1)
[0015] 式中maskO(i,j)为1表示对应像素属于厚云区域，maskO(i,j)为0表示对应像 素属于无云区域；
[0016] 步骤2-2 ;对所述厚云区域初始指示模板maskO(i，j)进行形态学处理，去除孤立 的斑点并填补细小的空洞，得到厚云区域指示模板mask(i，j)。
[0017] 上述步骤2-1中阔值Tc的计算方法为：
[001引Tc=m+Ao似
[001引其中，m为所述厚云影像j，k)的藍光波段图像I(i，如的像素均值：
[0020]
(3)
[0021] 0为厚云影像j，k)的藍光波段图像I(i，j)的像素标准方差：
[0022]
(4)
[0023] A为权重系数；
[0024]
(5)
[00巧]上述步骤3由W下具体步骤组成：
[0026] 步骤3-1 ;计算厚云影像的各波段的图像Ir"(i,j,k)，1《k《K在X方向的梯度 值吊，〇 (X，y，k)和y方向的梯度值gy，〇 (X，y,k):
[0027]
(6)
[002引步骤3-2 ;计算所述各多时相影像j，k)，1《t《T在X方向的梯度值gx，t(X，y,k)和y方向的梯度值gy，t(X，y,k):
[002引
（7)
[0030] 步骤3-3 ;计算所述厚云影像Ir"(i，j，k)，1《k《K与各多时相影像 (i，j，k)，1《t《T在无云区域对应位置梯度值之间的均方根误差RMSEt;
[00引]
W
[0032] 步骤3-4 ;选择RMSEt最小的多时相影像It，作为所述厚云影像的去云参考影 像 (i,j,k)。
[0033] 上述步骤4中所述初步去云结果1。(1，j，k)中无云区域像素值与所述厚云影像 It"(i，j，k)的无云区域像素值相等，其厚云区域图像的像素值由下式计算得到：
[0034] (9)
[00对其中，Q代表厚云区域内部；洗代表厚云区域边界，▽表示梯度算子，Ii表示I。 在区域Q的像素。
[0036] 上述步骤4中厚云区域图像的像素值可W根据W下递推公式计算得到：
[0040]式中P郝Pj.分别代表初步去云结果I。(1，j，k)中第i个和第j个厚云区域像素；N(Pi)代表所述像素Pi的四邻域像素集；
[00 …
（12)
[004引式中1。站）表示初步去云结果I0(i，j，k)中像素di的像素值，Iref(Pi)Irefhi)分 别表示所述去云参考影像Iaf(i，j，k)中像素Pi和的像素值，像素为N(Pi)与Q交集 内的第i个像素；C为A的化olesk巧分解，A=CCT，其中C为一个下S角矩阵，其下S角 部分具有和A完全相同的结构。
[0043] 上述步骤5中的变分模型代价函数为：
[0044]
(13)
[0045] 其中a为用于权衡两个约束项贡献度参数，a> 0。
[0046] 上述步骤5中的变分模型代价函数最小值的求解方法为采用下述迭代公式计算：
[0047]
[0048] 其中At是迭代步长，a是正则化参数，an的计算公式为；
[004引
山）
[0050] 式中var_n为设定的标准参数。
[0051] 本发明的有益效果在于：
[0052] 1、本发明能通过厚云影像和多幅多时相影像在梯度值之间的均方根误差（MSE) 确定参考影像。它无需人工交互，自动地去除厚云及其阴影
[0053] 2、本发明兼顾了原始影像的像素亮度与参考影像的梯度信息，对像素值有较好保 真性。
【附图说明】
[0054] 图1为本发明的流程图；
[0055] 图2为本发明步骤2的流程图；
[0056] 图3为本发明步骤3的流程图。
[0057] 图4为本发明的遥感卫星采集的厚云影像；
[0058] 图5为本发明的遥感卫星采集的参考影像；
[0059] 图6为本发明检测出的厚云区域的指示模板；
[0060] 图7为本发明的初步去云结果；
[0061] 图8为本发明的最终去云结果。
【具体实施方式】
[0062] 实施例1 ;
[0063] 如图1所示，一种多时相遥感影像的厚云自动去除方法，由W下步骤组成：
[0064] 步骤1 ;采集遥感卫星厚云影像Ir"(i，j，k)，1《k《K和T幅与其同区域的多时 相影像Iraw.t(1，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 ;
[006引步骤2 ;厚云区域检测，得到厚云区域指示模板mask(i，如；
[0066] 步骤3 ;自动选择参考影像；自动选择一幅多时相影像作为参考影像Iuf(i，j，k);
[0067] 步骤4 ;采用泊松方程修复方法去除厚云，得到初步去云结果1。(1，j，k);
[006引步骤5 ;将所述初步去云结果1。(1，j，k)和参考影像luf(i，j，k)纳入变分模型，再 次去除厚云，得到最终的去云结果r0(i，j，k);
[006引如图2所示，上述步骤2由W下具体步骤组成：
[0070] 步骤2-1 ;将厚云影像j，k)的藍光波段图像I(iJ)与阔值Te比较，计算厚 云区域初始指示模板maskO(i，j);
[0071]
(1)
[0072] 式中maskO(i,j)为1表示对应像素属于厚云区域，maskO(i,j)为0表示对应像 素属于无云区域；
[0073] 步骤2-2 ;对所述厚云区域初始指示模板maskO(i，j)进行形态学处理，去除孤立 的斑点并填补细小的空洞，得到厚云 区域指示模板mask(i，j)。
[0074] 上述步骤2-1中阔值Tc的计算方法为：
[00巧]Tc= m+入0似
[0076] 其中，m为所述厚云影像j，k)的藍光波段图像I(i，j)的像素均值：
[0077]
(3)
[007引 0为厚云影像j，k)的藍光波段图像I(i，如的像素标准方差：
[008引如图3所示，上述步骤3由W下具体步骤组成：
[0083] 步骤3-1 ;计算厚云影像的各波段的图像j，k)，1《k《K在X方向的梯度 值吊，〇 (X，y，k)和y方向的梯度值gy，〇 (X，y,k):
[0084]
(6)
[00财步骤3-2 ;计算所述各多时相影像j，k)，1《t《T在X方向的梯度值gx，t(X，y,k)和y方向的梯度值gy，t(X，y,k):
[0086]
(7)
[0087] 步骤3-3 ;计算所述各厚云影像j，k)，1《k《K与各多时相影像 (i，j，k)，1《t《T在无云区域对应位置梯度值之间的均方根误差RMSEt;
[0088] (8)
[0089] 步骤3-4 ;选择RMSEt最小的多时相影像It，作为所述厚云影像的去云参考影 像 (i,j,k)。
[0090] 上述步骤4由w下具体步骤组成：
[0091] 步骤4-1 ;所述初步去云结果Iu(i，j，k)中无云区域像素值与所述厚云影像 It"(i，j，k)的无云区域像素值相等，其厚云区域图像的像素值由下式计算得到：
[0092]
巧）
[0093] 其中，Q代表厚云区域内部；过)_代表厚云区域边界，▽表示梯度算子，Ii表示I。 在区域Q的像素。
[0094] 上述步骤4-1中厚云区域图像的像素值可W根据W下递推公式计算得到：
[009引式中P郝Pj.分别代表初步去云结果I。(1，j，k)中第i个和第j个厚云区域像素；N(Pi)代表所述像素Pi的四邻域像素集；
[009引
。2)
[0100]式中1。姑）表示初步去云结果1。(1，j，k)中像素Qi的像素值，I 分 别表示所述去云参考影像Iaf(i，j，k)中像素Pi和的像素值，像素为N(Pi)与Q交集 内的第i个像素；C为A的化olesk巧分解，A=CCT，其中C为一个下S角矩阵，其下S角 部分具有和A完全相同的结构。
[0101] 上述步骤5中的变分模型代价函数为：
[0102]
(13)
[0103] 其中a为用于权衡两个约束项贡献度参数，a> 0。
[0104] 上述步骤5中的变分模型代价函数最小值的求解方法为采用下述迭代公式计算：
[0105]
[0106] 其中At是迭代步长，a是正则化参数，an的计算公式为；
[0108] 式中var_n为设定的标准参数。在本实施例中，var_n为519. 84。
[0107] / , (15)
【主权项】
1. 一种多时相遥感影像的厚云自动去除方法，其特征在于：由以下步骤组成： 步骤1 :采集遥感卫星厚云影像IMW(i，j，k)，1 < k < K和T幅与其同区域的多时相影 像Wi，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 :自动选择参考影像：自动选择一幅多时相影像作为参考影像IMf(i，j，k); 步骤4 :采用泊松方程修复方法去除厚云，得到初步去云结果Itl (i，j，k); 步骤5 :将所述初步去云结果IJi，j，k)和参考影像IMf (i，j，k)纳入变分模型，再次去 除厚云，得到最终的去云结果Γ C1 (i，j，k)。2. 根据权利要求1所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步 骤2由以下具体步骤组成： 步骤2-1 :将厚云影像IMW(i，j，k)的蓝光波段图像I(i，j)与阈值Tc比较，计算厚云区 域初始指示模板maskO(i，j):式中maskO(i, j)为1表示对应像素属于厚云区域，maskO(i, j)为O表示对应像素属 于无云区域； 步骤2-2 :对所述厚云区域初始指示模板maskO (i，j)进行形态学处理，去除孤立的斑 点并填补细小的空洞，得到厚云区域指示模板mask(i, j)。3. 根据权利要求2所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步 骤2-1中阈值T。的计算方法为： Tc= m+ λ σ ⑵ 其中，m为厚云影像1_的蓝光波段图像I的像素均值：σ为厚云影像1_的蓝光波段图像I的像素标准方差：λ为权重系数：4. 根据权利要求1所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步 骤3由以下具体步骤组成： 步骤3-1 :计算厚云影像的各波段的图像IMW(i，j，k)，1彡k彡K在X方向的梯度值 gx，Q(x, y, k)和 y 方向的梯度值 gy，Q(x, y, k): gx；0 (x, y, k) =I (i+1, j, k) -I (i, j, k) ,I ^ I ^ μ, I ^ j ^ N, I ^ k ^ K (6) gy,〇 (x, y, k) =I (i, j+1, k) -I (i, j, k) 步骤3-2 :计算所述各多时相影像Iraw> t (i, j, k), I < t < T在x方向的梯度值 gx，t (X, y, k)和 y 方向的梯度值 gy，t (X, y, k): gx, t (x, y, k) = It (i+1, j, k) -It (i, j, k) ,I ^ i ^ M, I ^ j ^ N, I ^ t ^ T (7) gy；t(x, y, k) = It (i, j+1, k)-It(i, j, k) 步骤3-3 :计算所述各厚云影像Iraw (i, j, k), I < k < K与各多时相影像 Iraw，t(i, j，k)，l 彡 t彡T在 无云区域对应位置梯度值之间的均方根误差RMSEt:步骤3-4 :选择RMSEt最小的多时相影像I t，作为所述厚云影像1_的去云参考影像 Iref (i，j，k) 〇5. 根据权利要求1所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步 骤4中所述初步去云结果Itl(Hk)中无云区域像素值与所述厚云影像I MW(i，j，k)的无云 区域像素值相等，其厚云区域图像的像素值由下式计算得到：其中，Ω代表厚云区域内部；?Ω代表厚云区域边界，▽表示梯度算子，I1表示I ^在区 域Ω的像素。6. 根据权利要求5所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步 骤4中厚云区图像的像素值根据以下递推公式计算得到：式中PdP P j分别代表初步去云结果10 (i，j，k)中第i个和第j个厚云区域像素；N(Pi) 代表所试俛麦η.的四部±或俛麦隹.式中1〇(1)表示初步去云结果I〇(i，j，k)中像素qi的像素值，IjpDUqi)分别表 示所述去云参考影像IMf(i，j，k)中像素pJPq ^勺像素值，像素CiiSN(Pi)与Ω交集内的 第i个像素；C为A的Choleskey分解，A = CCT，其中C为一个下三角矩阵，其下三角部分 具有和A完全相同的结构。7. 根据权利要求1所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步 骤5中的变分模型代价函数为：其中α为用于权衡两个约束项贡献度参数，α彡〇。8. 根据权利要求7所述的多时相遥感影像的厚云自动去除方法，其特征在于：所述步 骤5中的变分模型代价函数最小值的求解方法为采用下述迭代公式计算：其中At是迭代步长，α是正则化参数，αη的计算公式为：式中var_n为设定的标准参数。
【专利摘要】本发明公开了一种多时相遥感影像的厚云自动去除方法，由以下步骤组成：步骤1：采集遥感卫星厚云影像和T幅与其同区域的多时相影像；步骤2：厚云区域检测，得到厚云区域指示模板；步骤3：自动选择参考影像：自动选择一幅多时相影像作为参考影像；步骤4：采用泊松方程修复方法去除厚云，得到初步去云结果；步骤5：将所述初步去云结果和参考影像纳入变分模型，再次去除厚云，得到最终的去云结果。本发明能通过厚云影像和多幅多时相影像在梯度值之间的均方根误差（）确定参考影像。它无需人工交互，自动地去除厚云及其阴影,它兼顾了原始影像的像素亮度与参考影像的梯度信息，对像素值有较好保真性。
【IPC分类】G06T5/00
【公开号】CN104881850
【申请号】CN201510274174
【发明人】聂龙保, 黄微, 张婷婷, 孟新知, 叶分晓
【申请人】上海大学
【公开日】2015年9月2日
【申请日】2015年5月26日