一种多光谱光学遥感图像数据去除薄云的方法

暂无

一种多光谱光学遥感图像数据去除薄云的方法\n技术领域\n[0001] 本发明属于光学遥感图像去云技术领域，具体涉及一种多光谱光学遥感图像数据去除薄云的方法。\n背景技术\n[0002] 大气对太阳能的散射和吸收作用对光学卫星图像的获取存在着不同程度的影响。\n由于大气校正技术的发展，广泛用于大气校正的软件和算法包括ACORN–Atmospheric CORrection Now(InSpec,2002)，ATCOR–the ATmospheric CORrection program(Thiemann and Hermann 2002)，ATREM–the ATmospheric REMoval program[Center for the Study of Earth from Space(CSES),University of Colorado)],FLAASH–Fast Line–of–sight Atmospheric Analysis of Spectral Hypercubes(Research Systems,Inc.,2003)。然而，很多遥感卫星图像被云覆盖，云对光学卫星图像的影响限制了卫星图像的使用。由于云在空间上的不确定性，因波段的个数有限，对多波段遥感图像的去云仍面临着很大的挑战。尤其对于常年有云覆盖的地区，干净可用的遥感卫星数据很有限的情况下，发明一种有效的去薄云方法是非常有必要的。\n[0003] 云大致可以分为卷云、层云和积云。卷云属于高云，由细小稀疏的冰晶组成，卷云比较薄，对可见光和近红外太阳辐射几乎是透明的。层云属于中层云，由小的水和冰珠组成，具有一定的水平结构，当其厚度较薄时，短波红外的太阳辐射是有可能被穿透的。积云分布高度更低，且轮廓分明，由水汽凝结而成，光学波段的太阳辐射是几乎不能穿透它，因此在光学遥感图像中呈高亮且覆盖了地物信息。\n[0004] 去薄云的方法可分为单幅图像法和多时相或多传感器法。单幅图像的去云方法如同态滤波法和小波变换等，具体实施方法参见[1]冯春，马建文，戴芹，陈雪.一种改进的遥感图像薄云快速去除方法[J].国土资源遥感，2004,04:1–3+18–77；[2]张波，季民河，沈琪.基于小波变换的高分辨率快鸟遥感图像薄云去除[J].遥感信息，2011,03:38–43。现有的单幅图像去云方法缺点主要在于将整幅图像在频率域上去除低频成分，不可避免的去除了一些背景信息；选取截止频域时候是凭经验选取，需要不断尝试；而且要把图像转换为频率域不适合处理数量较大的图像。基于多时相或多传感器的方法，如多时相图像叠加去云法，具体处理过程参见：[3]Chao–Hung Lin and Po–Hung Tsai,“Cloud removal from multitemporal  satellite images  using information cloning,”IEEE \nTran.Geoscience and remote sensing,vol.51,pp.232–241,January 2013。多时相或多传感器的去云方法是基于多时相的遥感卫星数据，而这在常年有云覆盖的地区几乎不可行的。虽然不同传感器之间的图像是可以构成多时相的图像集，但如何解决不同传感器之间的图像配准和辐射差异仍是有待解决的问题。因此，多时相或多传感器的方法对数据要求苛刻，实用的可行性差。\n发明内容\n[0005] 本发明的目的是为了解决上述已有去薄云技术中存在的问题，提出了一种多光谱光学遥感图像数据去除薄云的方法，该方法充分利用了多光谱光学遥感数据单幅图像自身数据的优势，摆脱了多数据去云方法数据要求苛刻等诸多束缚，克服了现有单幅影像去薄云方法滤波作用整幅图像等缺点。从而恢复在多光谱光学遥感数据中薄云覆盖区域处地物光谱信息，提高多光谱光学遥感数据的使用质量和图像的应用能力。\n[0006] 为了方便描述本发明的内容，首先作以下定义：\n[0007] 定义1、辐亮度\n[0008] 辐亮度是单位立体角和单位面积上的能量。单位是W·cm-2·sr-1·μm-1，其中sr(球面度)是立体角的单位，μm(微米)是波长单位。\n[0009] 定义2、地表反射率\n[0010] 地表反射率是指地表物体向各个方向上反射的太阳总辐射通量与到达该物体表面上的总辐射通量之比。\n[0011] 定义3、大气校正\n[0012] 光学遥感的大气校正是去除遥感数据中的大气效应，获取地表反射率的过程。大气校正主要包括两部分：大气参数估计和地表反射率反演。对于水平均匀的大气和朗伯体地面，地表反射率rλ通过以下公式得到。\n[0013]\n[0014] 其中，Lλ为表观光谱辐亮度，Lp是大气路径辐射，S是大气的半球反照率，F0是在大气顶部垂直于太阳光束入射的太阳能通量密度的 τ(μs)和τ(μv)是太阳到地面和地表到传感器的总透过率。μs和μv是太阳角和观测角的余弦值。详见文献“定量遥感”，梁顺林等编著，科学出版社。\n[0015] 定义4、大气校正软件模块\n[0016] 传统的基于辐射传输模型的大气校正软件模块可以通过用户提供基础大气特征信息或者特点的大气吸收波段计算特定时间的大气散射和吸收特性。通过特定的大气特征参数就可以获取地表反射率。广泛用于大气校正的软件和算法包括：ACORN–Atmospheric CORrection Now(InSpec,2002)，ATCOR–the ATmospheric CORrection program(Thiemann and Hermann 2002)，ATREM–the ATmospheric REMoval program[Center for the Study of Earth from Space(CSES),University of Colorado)],FLAASH–Fast Line–of–sight Atmospheric Analysis of Spectral Hypercubes(Research Systems,Inc.,2003)等。这些大气校正软件模块通常需要用户提供：\n[0017] ·遥感图像经纬度信息，\n[0018] ·遥感图像获取日期和时刻，\n[0019] ·遥感图像的海报高度,\n[0020] ·卫星传感器的高度，\n[0021] ·大气模型(如：中纬度-夏季，中纬度-冬季，热带)，\n[0022] ·辐射校正的辐亮度数据(如：数据单位必须是W·cm-2·sr-1·μm-1)，\n[0023] ·传感器特定的波段信息(如：波段的半高全宽)，\n[0024] ·遥感图像获取时的大气能见度等。\n[0025] 定义5、云数据\n[0026] 由卫星自上而下观测到的地球上的云层覆盖和地表面特征的图像。利用云数据可以识别不同的云特征，确定它们的范围，区域，例如：Landsat8的QA波段。\n[0027] 定义6、二值化处理\n[0028] 任何利用云数据提取云二值图的方法在本发明专利中都适用，例如：Landsat8的QA波段，QA波段的14,15数据位用于描述图像云信息，当有云的可信度高于33％时，设置该区域为有云区域，将其设置成1，小于33％时，设置该区域为无云区域设置成0，从而得到二值图。\n[0029] 定义7、高斯平滑\n[0030] 高斯滤波器是一类根据高斯函数的形状来选择权值的线性平滑滤波器。高斯平滑滤波器对于抑制服从正态分布的噪声非常有效。\n[0031] 二维零均值高斯函数计算公式为：\n[0032] 其中σ表示高斯分布参数。详见参考文献“图像降噪的自适\n应高斯平滑滤波”，谢勤岚等，计算机工程与应用。\n[0033] 本发明提供了一种多光谱光学遥感图像数据去除薄云的方法，该方法的步骤如下：\n[0034] 步骤1：数据准备\n[0035] 本发明提供的数据包括：多光谱卫星原始图像数据Qcal，它是一个有m列，n行的矩阵，包括j个可见光和近红外波段数据V和k个短波红外波段数据S，其中m,n,j,k均为正整数；特定波段的倍率修正因子ML，特定波段的增量修正因子AL；云数据，记为Dc。\n[0036] 步骤1中所述参数或数据均为已知。\n[0037] 步骤2：对多光谱光学遥感数据进行大气校正预处理\n[0038] 对步骤1中的多光谱卫星原始图像数据Qcal利用公式L＝MLQcal+AL计算得到表观辐亮度数据L。其中ML是特定波段的倍率修正因子，Qcal是步骤1提供的多光谱卫星原始图像数据，AL是步骤1提供的特定波段的增量修正因子。\n[0039] 将上述计算得到的表观辐亮度数据L通过经典的大气校正软件模块进行传统的大气校正处理，获得大气校正结果，即得到j个可见光和近红波段的反射率矩阵RV和k个短波红外波段的反射率矩阵RS，i,j均为正整数。\n[0040] 步骤3：提取二值图矩阵Dc'(0–无云，1–有云)\n[0041] 在遥感图像成像区域中，将步骤1中的云数据Dc进行传统的二值化方法处理得到矩阵Dc'，矩阵Dc'由“0”和“1”两类元素组成，矩阵Dc'中0代表无云区域像素点的像素值，1代表有云区域像素点的像素值。另外在矩阵Dc'中，定义无云区域像素点的总数为P，有云区域像素点的总数为Q，P和Q均为正整数。\n[0042] 步骤4：有云区域像素点和无云区域像素点配对及替换\n[0043] 操作a：计算短波红外波段的反射率矩阵\n[0044] 利用步骤2得到的短波红外波段的反射率矩阵RS和步骤3得到的二值图矩阵Dc'，通过公式RSc＝RS*Dc'计算得到有云区域像素点短波红外波段的反射率矩阵RSc，*表示相乘。矩阵RSc是三维矩阵。在矩阵RSc中，第i层，第q个有云区域像素点(q＝1,2,3,…,Q)的反射率值表示为RSc(aq,bq,i)，其中aq表示矩阵RSc中第q个有云区域像素点的行号，bq表示矩阵RSc中第q个有云区域像素点的列号，i表示矩阵RSc的层数，Q为有云区域像素点总数，a，b，i，h和q均为正整数。\n[0045] 利用步骤2得到的短波红外波段的反射率矩阵RS和步骤3得到的二值图矩阵Dc'，通过公式RSu＝RS*～Dc'计算得到无云区域所有像素点短波红外波段的反射率矩阵RSu，*表示相乘，～表示取反。在矩阵RSu中，第i层无云区域所有像素点的可见光和近红外波段的反射率值表示为RSu(x,y,i)，其中x表示矩阵RSu的行号，y表示矩阵RSu的列号，i表示矩阵RSu的层数，x，y和i均为正整数，1≤x≤m,1≤y≤n。\n[0046] 操作b：寻找与有云区域像素点对应的无云区域像素点位置\n[0047] 首先定义与第q个有云区域像素点对应的无云区域像素点位置为\n[0048] 对矩阵RSc中每一层的第q个有云区域像素点RSc(aq,bq,i)，通过公式d(x,y,i)＝|RSc(aq,bq,i)-RSu(x,y,i)|，计算得到矩阵d(x,y,i)，其中i＝(1,2,…,k)，“||”表示取绝对值；\n然后，通过公式 得到矩阵d(x,y)；最后，通过公式\n计算得到无云区域像素点位置 其中符号“{A|B}”表示满足条件B的A的所有取值，“min()”表示计算最小值。与无云区域像素点位置 对应的有云区域像素点位置为aq,bq。有云区域像素点位置与对应的无云区域像素点位置的地物在短波红外波段有相似的光谱特征。\n[0049] 对所有Q个有云区域像素点位置确定对应的无云区域云像素点位置，共得到Q个aq,bq与 像素点位置对，Q为有云区域像素点总数。\n[0050] 操作c：计算可见光和近红外波段的反射率矩阵\n[0051] 利用步骤2得到的可见光和近红外波段的反射率矩阵RV和步骤3得到的二值图矩阵Dc'，通过公式RVc＝RV*Dc'计算得到有云区域可见光和近红外波段的反射率矩阵RVc；利用可见光和近红波段的反射率矩阵RV和二值图矩阵Dc'，通过公式RVu＝RV*(～Dc')计算得到无云区域可见光和近红外波段的反射率矩阵RVu，*表示相乘，～表示取反。\n[0052] 操作d：替换有云区域像素值\n[0053] 定义：有云区域像素点在可见光和近红外波段去云后的反射率矩阵为RVc'；\n[0054] 在无云区域可见光和近红外波段的反射率矩阵RVu中，取出无云区域像素点位置的反射率值，记为 在有云区域可见光和近红外波段的反射率矩阵RVc中，\n取出有云区域像素点位置(aq,bq)的反射率值，记为RVc(aq,bq)。对每一个有云区域像素点位置的反射率值RVc(aq,bq)用对应的无云区域像素点位置的反射率值 替换，最后得到所有有云区域像素点在可见光和近红外波段去云后的反射率矩阵RVc'。\n[0055] 步骤5：平滑滤波和确定去云后结果\n[0056] 对步骤4中得到的矩阵RVc'进行传统的高斯平滑滤波处理，得到有云区域可见光和近红外波段矩阵去云后结果RVc”。通过公式R＝RVc”∪RVu计算，最终得到整幅图像可见光和近红外波段数据去云后的结果R，其中符号“∪”表示“并”操作。\n[0057] 本发明提供的一种多光谱光学遥感图像数据去除薄云的方法，该方法充分利用了多光谱光学遥感数据单幅图像自身数据的优势，摆脱了多数据去云方法数据要求苛刻等诸多束缚，和现有单幅影像去薄云方法滤波作用整幅图像等缺点。本发明与现有的去薄云的技术相比，在处理难度上，相对于单幅数据方法要进行时域变换、频域变换和频率域滤波等一系列复杂处理，本发明操作简单，更有优势，也更有效率。\n附图说明\n[0058] 图1为本发明方法流程图；\n[0059] 图2为多光谱Landsat8卫星原始图像数据ML和AL参数表；\n[0060] 图3为本发明具体实施案例中用多光谱Landsat8卫星原始图像数据Qcal的波段4(红)、3(绿)、2(蓝)合成的真彩色图像；\n[0061] 图4为Landsat8卫星云数据Dc，即为QA波段；\n[0062] 图5为多光谱Landsat8卫星原始图像数据进行大气校正后结果用波段4(红)、3(绿)、2(蓝)合成的真彩色图像；\n[0063] 图6为有云和无云区域的二值图，图中白色区域为有云区域，黑色区域为无云区域；\n[0064] 图7为利用本发明的具体实施方法处理获得的Landsat8多光谱光学遥感数据去除薄云处理后结果用波段4(红)、3(绿)、2(蓝)合成的真彩色图像。\n具体实施方式\n[0065] 为了对本发明的技术特征、目的和效果有更加清楚的理解，现说明本发明的具体实施方式。\n[0066] 参考图1处理流程图，本发明以Landsat8多光谱光学遥感数据为例，具体阐述这种多光谱光学遥感数据的薄云去除方法。\n[0067] 该实施案例的步骤主要分为5步：\n[0068] 步骤1：参数准备：\n[0069] 本发明提供的数据包括：Landsat-8卫星多光谱原始图像数据Qcal，有118列，118行，其中包括5个可见光和近红外波段数据V和2个短波红外波段数据S，其波段4(红)、3(绿)、2(蓝)合成的真彩色图像见附图3，；特定波段的倍率修正因子ML和特定波段的增量修正因子AL系数见附图2；云数据为Landsat8的QA波段，见附图4，记为Dc。\n[0070] 步骤2：对多光谱光学遥感数据进行大气校正预处理\n[0071] 对步骤1中的Landsat-8卫星多光谱原始图像数据Qcal利用公式L＝MLQcal+AL计算得到表观辐亮度数据L。其中ML是特定波段的倍率修正因子，Qcal是步骤1提供的Landsat-8卫星多光谱原始图像数据，AL是步骤1提供的特定波段的增量修正因子。\n[0072] 将上述计算得到的表观辐亮度数据L通过经典的大气校正软件模块进行传统的大气校正处理，获得大气校正结果，即得到j个可见光和近红波段的反射率矩阵RV和k个短波红外波段的反射率矩阵RS。\n[0073] 大气校正结果用波段4(红)、3(绿)、2(蓝)合成的真彩色图像如附图5所示。\n[0074] 步骤3：提取二值图矩阵Dc'(0–无云，1–有云)\n[0075] 在遥感图像成像区域中，将步骤1中的云数据Dc进行传统的二值化方法处理得到矩阵Dc'，矩阵Dc'由“0”和“1”两类元素组成，矩阵Dc'中0代表无云区域像素点的像素值，1代表有云区域像素点的像素值。另外在矩阵Dc'中，定义无云区域像素点的总数为P，有云区域像素点的总数为Q，P和Q均为正整数。\n[0076] 如图6所示，图中白色区域为有云区域，黑色区域为无云区域。\n[0077] 步骤4：有云区域像素点和无云区域像素点配对及替换\n[0078] 操作a：计算短波红外波段的反射率矩阵\n[0079] 利用步骤2得到的短波红外波段的反射率矩阵RS和步骤3得到的二值图矩阵Dc'，通过公式RSc＝RS*Dc'计算得到有云区域像素点短波红外波段的反射率矩阵RSc，*表示相乘。矩阵RSc是三维矩阵。在矩阵RSc中，第i层，第q个有云区域像素点(q＝1,2,3,…,Q)的反射率值表示为RSc(aq,bq,i)，其中aq表示矩阵RSc中第q个有云区域像素点的行号，bq表示矩阵中第q个有云区域像素点RSc的列号，i表示矩阵RSc的层数，Q为有云区域像素点总数，a，b，i，h和q均为正整数。\n[0080] 利用步骤2得到的短波红外波段的反射率矩阵RS和步骤3得到的二值图矩阵Dc'，通过公式RSu＝RS*～Dc'计算得到无云区域所有像素点短波红外波段的反射率矩阵RSu，*表示相乘，～表示取反。在矩阵RSu中，第i层无云区域所有像素点的可见光和近红外波段的反射率值表示为RSu(x,y,i)，其中x表示矩阵RSu的行号，y表示矩阵RSu的列号，i表示矩阵RSu的层数，x，y和i均为正整数，1≤x≤m,1≤y≤n。\n[0081] 操作b：寻找与有云区域像素点对应的无云区域像素点位置\n[0082] 首先定义与第q个有云区域像素点对应的无云区域像素点位置为\n[0083] 对矩阵RSc中每一层的第q个有云区域像素点RSc(aq,bq,i)，通过公式d(x,y,i)＝|RSc(aq,bq,i)-RSu(x,y,i)|，计算得到矩阵d(x,y,i)，其中i＝(1,2,…,k)，“||”表示取绝对值；然后，通过公式 得到矩阵d(x,y)；最后，通过公式\n计算得到无云区域像素点位置 其中符号“{A|B}”表\n示满足条件B的A的所有取值，“min()”表示计算最小值。与无云区域像素点位置 对应的有云区域像素点位置为aq,bq。有云区域像素点位置与对应的无云区域像素点位置的地物在短波红外波段有相似的光谱特征。\n[0084] 对所有Q个有云区域像素点位置确定对应的无云区域云像素点位置，共得到Q个aq,bq与 像素点位置对，Q为有云区域像素点总数。\n[0085] 操作c：计算可见光和近红外波段的反射率矩阵\n[0086] 利用步骤2得到的可见光和近红外波段的反射率矩阵RV和步骤3得到的二值图矩阵Dc'，通过公式RVc＝RV*Dc'计算得到有云区域可见光和近红外波段的反射率矩阵RVc；利用可见光和近红波段的反射率矩阵RV和二值图矩阵Dc'，通过公式RVu＝RV*(～Dc')计算得到无云区域可见光和近红外波段的反射率矩阵RVu，*表示相乘，～表示取反。\n[0087] 操作d：替换有云区域像素值\n[0088] 定义：有云区域像素点在可见光和近红外波段去云后的反射率矩阵为RVc'；\n[0089] 在无云区域可见光和近红外波段的反射率矩阵RVu中，取出无云区域像素点位置的反射率值，记为 在有云区域可见光和近红外波段的反射率矩阵RVc中，\n取出有云区域像素点位置(aq,bq)的反射率值，记为RVc(aq,bq)。对每一个有云区域像素点位置的反射率值RVc(aq,bq)用对应的无云区域像素点位置的反射率值 替换，最后得到所有有云区域像素点在可见光和近红外波段去云后的反射率矩阵RVc'。\n[0090] 步骤5：平滑滤波和最终去云结果\n[0091] 对步骤4中得到的矩阵RVc'进行传统的高斯平滑滤波处理，得到有云区域可见光和近红外波段矩阵去云后结果RVc”。通过公式R＝RVc”∪RVu计算，最终得到整幅图像可见光和近红外波段数据去云后的结果R，其中“∪”表示“并”操作。\n[0092] 如图7所示为本发明具体实施案例获取结果的用波段4,3,2合成的真彩色图像。