CN109308688B - 一种可见光和近红外波段厚云及阴影去除方法 - Google Patents
一种可见光和近红外波段厚云及阴影去除方法 Download PDFInfo
- Publication number
- CN109308688B CN109308688B CN201811119642.8A CN201811119642A CN109308688B CN 109308688 B CN109308688 B CN 109308688B CN 201811119642 A CN201811119642 A CN 201811119642A CN 109308688 B CN109308688 B CN 109308688B
- Authority
- CN
- China
- Prior art keywords
- cloud
- pixel
- image
- target image
- pixels
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 42
- 238000001514 detection method Methods 0.000 claims abstract description 7
- 238000012417 linear regression Methods 0.000 claims abstract description 7
- 238000012937 correction Methods 0.000 claims abstract description 5
- 238000001228 spectrum Methods 0.000 claims description 14
- 230000003595 spectral effect Effects 0.000 claims description 13
- 238000011160 research Methods 0.000 claims description 8
- 230000001419 dependent effect Effects 0.000 claims description 3
- 238000000611 regression analysis Methods 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 abstract description 12
- 238000013178 mathematical model Methods 0.000 abstract description 2
- 238000002310 reflectometry Methods 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 4
- 238000002474 experimental method Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000012952 Resampling Methods 0.000 description 1
- 230000003139 buffering effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G06T5/70—
-
- G06T5/80—
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10048—Infrared image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
- G06T2207/30192—Weather; Meteorology
Abstract
本发明公开了一种可见光和近红外波段厚云及阴影去除方法，目标晴空图像的构建主要是应用多时相图像的时空信息对云和云阴影覆盖像元值的进行估算，包括以下步骤：1)目标图像与参考图像的几何校正；2)对目标图像和参考图像进行云和阴影检测；3)云和阴影覆盖像元值(DN值、反射率等)估算。本发明根据云分布的动态变化性、遥感平台传感器的成像周期性以及下垫面在短时间内无明显变化等特点，基于多时相图像的时空信息，利用线性回归数学模型的逐波段逐像元进行定量估算云下像元值，建立了一套实现逐像元逐波段自动化厚云及其阴影去除的流程方法，提高了云和云阴影去除精度，实现全自动厚云及阴影去除。
Description
技术领域
本发明属于遥感应用技术领域，尤其涉及一种可见光和近红外波段厚云及阴影去除方法。
背景技术
在诸多遥感应用中，如土地覆盖分类、农业干旱监测和土地覆盖变化检测等，都需要用到晴空遥感图像。但是在光学遥感研究中，遥感图像普遍存在云覆盖的现象。而目前可见光和近红外(VNIR)波段的去云技术仍存在目标图像估算后的云覆盖像元与周围晴空像元光谱不一致的问题。
受到大气动力的影响，云层随着大气中的温度和压强的变化而移动，而不会在空中同一位置停留太长时间。通过周期性重访成像，遥感系统能够获取同一地点不同时间的图像。基于此，利用先前时间遥感图像重建云下地面信息成为可能。这是在可见光和近红外(VNIR)波段厚云及阴影去除的基本理论。云阴影常见于云层厚度、高度和光照条件的变化而生成不同形状的阴影，在云去除过程中，云阴影也可一同去除。基于以上条件，可知多时相遥感图像可用于去除遥感图像上厚云及阴影去除是可行的。
在厚云及阴影去除的过程中，我们将定义一个有云和阴影覆盖的图像为目标图像。云和阴影的去除，即我们将重建云覆盖区域地面丢失的信息，如此就必须借助于参考图像作为辅助信息，将用于信息重建的辅助图像定义为参考图像。参考图像与目标图像在相同的位置具有相同或极近似的下垫面。基于卫星轨道运行的周期性，我们可以获取相同位置一定周期的卫星图像。一般来说，从逻辑上可以假设，在其自然条件下或没有人为干扰的下垫面在短时间内(如传感器地面观测的几天间隔)不会发生较大变化。基于以上假设，使用先前成像日期获得的参考图像来重建目标图像云和云阴影覆盖下丢失信息是可行的。
厚云及阴影去除目前存在以下几方面问题：应用其它多种数据源，预处理工作较多，过程复杂；机器学习等方法可较好实现，但是局限于样本质量和数量。应用多时相图像进行云覆盖区地面信息重建，其基本要求是目标图像应保持待处理研究区与晴空区之间的光谱一致性。传统的去除云像元的方法通常是通过直接使用以前的参考图像来取代云覆盖区的地面信息，通常会导致重建像元和周围晴空像元值在可比性和一致性方面出现问题。
发明内容
本发明目的是解决上述问题，设计一种可见光和近红外波段厚云及阴影去除方法。根据云分布的动态变化性、遥感平台传感器的成像周期性以及下垫面在短时间内无明显变化等特点，利用线性回归等简单数学模型的逐波段逐像元应用基于多时相图像的时空信息进行定量估算云(及其阴影)覆盖像元值(DN值、反射率等)，建立了一种可见光和近红外波段厚云及阴影去除方法，提高了云和云阴影覆盖像元值估算的精度，实现全自动晴空图像构建。
为了实现上述目的，本发明的技术方案是：
一种可见光和近红外波段厚云及阴影去除方法，该方法的步骤是：
S1、采集同一研究区不同时相的卫星图像数据，定义一个有云和阴影覆盖的图像为目标图像T(0)，针对目标图像区域，选择若干不同时相的同种传感器影像作为参考图像，记T(-1)，T(-2)，…，T(j)；对目标图像和参考图像进行几何校正，确保它们相对应的像元在空间上处于相同的几何位置；
S2、读取图像T(0)，T(-1)，T(-2)，…，T(j)进行云检测，对应的云掩膜图像设为C(0)，C(-1)，C(-2)，…，C(j)，所有像元都被划分为晴空像元和云像元，晴空像元设为1，云像元设为0；
S3、确定目标图像T(0)中的云覆盖像元；
读取云掩膜图像C(0,i)来确定目标图像T(0)的第i个像元是晴空还是云覆盖，若值为1，那么该像元是晴空像元p；然后继续读取云掩膜图像下一个像元(i+1)的值C(0,i+1)，确定是晴空还是云覆盖；继续推进该过程直到找到一个云覆盖像元C(0,i)；
S4、确定参考图像相应的晴空像元；
对j＝-1，读C(-1,i)值来确定该像元是否为晴空；如果C(-1,i)＝0，则为云覆盖像元，就进入下一个参考图像j＝j-1，读取相应像元C(j,i)的云覆盖情况，以此类推，直到找到一个晴空像元C(j,i)＝1；
S5、在参考图像T(j)波段k围绕C(j,i)找光谱最相似的10个相邻像元；
定义SDi,k是参考图像波段k中周围晴空像元i与像元p的光谱差，公式如下：
SDi,k＝|Di,k-Dp,k| (1)
SDi,k和Dp,k是像元i和像元p在波段k的值；
在目标图像T(0)的云覆盖像元C(0,i)周围找相邻晴空像元，在参考图像T(j)中围绕C(j,i)找到与之相应的晴空像元；分别围绕C(0,i)和C(j,i)设置2个足够大的窗口，窗口边长像元数为奇数，如5*5，7*7，9*9，...，301*301，...，最大窗口设置为m,并检查两个窗口中每个像元对的值,在目标图像T(0)的云覆盖像元C(0,i)周围尽可能最小的窗口下找光谱差值SDi,k小于阈值n的晴空点，并将其光谱差值进行排序，取前10个光谱最相似的相邻晴空像元；调整窗口大小，直到在两个尽可能小的窗口中得到10个相应的符合阈值条件的晴空像元；如果超出最大窗口设置范围m仍然不能找到相应的10个像元，将阈值n逐渐从5，10，15，20，30，40，50…100…扩大，并继续逐窗口大小寻找,直到找到符合条件的点为止；
S6、建立波段k目标图像数据序列TIk和参考图像数据序列RIk；使用步骤S5的10个像元，分别获取它们在目标图像T(0)和参考图像T(j)的像元值来建立数据序列RIk和TIk；
TIk＝(TI1,k,TI2,k,TI3,k,TI4,k,TI5,k,TI6,k,TI7,k,TI8,k,TI9,k,TI10,k) (2)
TIk是目标图像种子像元在波段k的数据序列，TI1,k,TI2,k,TI3,k,TI4,k,TI5,k,TI6,k,TI7,k,TI8,k,TI9,k,TI10,k分别是种子像元1，2，3，4…9，10的像元值；
RIk＝(RI1,k,RI2,k,RI3,k,RI4,k,RI5,k,RI6,k,RI7,k,RI8,k,RI9,k,RI10,k)
(3)
RIk为参考图像种子像元在波段k的数据序列，RI1,k,RI2,k,RI3,k,RI4,k,RI5,k,RI6,k,RI7,k,RI8,k,RI9,k,RI10,k分别为种子像元1，2，3，4…9，10的像元值；
S7、通过最小二乘回归分析方法，建立数据序列TIk和RIk之间的光谱回归关系；
YTI,k＝a+bXRI,k (4)
其中，YTI,k是因变量，对应于目标图像波段k的像元值，XRI,k是自变量，对应于参考图像波段k的DN值，a和b分别是方程的回归系数；
S8、估算目标图像中云覆盖像元的值；用参考图像T(j)中C(i,j)＝1的晴空像元对应的像元值，实现目标图像中波段k云覆盖像元值的估算；
Dq,k＝＝a+bDp,k (5)
Dq,k是在目标图像云覆盖像元q在波段k的估算像元值，Dp,k是参考图像中波段k相应的晴空像元的真实像元值；
S9、循环步骤S5到步骤S8，估算目标图像其他波段云覆盖像元值；
S10、重复上述步骤S3到步骤S9，对所有云覆盖像元值进行估算，直到去除目标图像中的所有的云像元。
进一步的，在步骤S1中，参考图像的时间越接近目标图像越好，且保证所有参考图像的晴空像元的合集能够全部覆盖目标图像中的云、阴影和雪的分布范围。
进一步的，在步骤S2中，对目标图像与参考图像的进行几何校正时误差最好小于1个像元。
进一步的，所述目标图像与参考图像为采集Landsat 8 OLI的L1级卫星图像数据。
进一步的，所述利用目标云覆盖像元附近晴空像元的建立云像元逐波段线性回归方程。
进一步的，将目标图像和参考图像的QA波段像素值转换成16位二进制形式，直接读取QA波段文件，生成相应的云和云阴影掩膜文件。
进一步的，QA波段文件中，6-7字节和14-18字节分别对应了云阴影和云分布信息
本发明的可见光和近红外波段厚云及阴影去除方法，考虑目标图像与参考图像之间一定区域范围的光谱相关性，用参考图像的相应位置晴空像元值来替换目标图像云像元值。用数学公式表示，我们有
DTI,k＝DRI,k+SD (1)
DTI，k表示目标图像TI云像元波段k重建的像元值，DRI,k是参考图像RI波段k相应晴空像元像元值，SD是目标图像和参考图像在成像间隔内的光谱差别。因此，该方法的重要性在于如何估算该方法的SD。
如方程(1)所示，我们不能用参考图像相同位置晴空像元值直接替换目标图像云覆盖像元值q(Dq,k表示像元q在波段k像元值)，是由于原始像元与替换像元的光谱、辐射不可比性和不一致。
为了确保这种光谱和辐射的可比性和一致性，我们用目标图像和参考图像之间的相似性来估算云覆盖像元的差值，即公式(1)中SD值。可以通过考虑若干与云覆盖像元地理位置邻近的晴空像元在目标图像和参考图像之间的光谱差异，建立目标图像和参考图像之间的光谱相关性来估算云覆盖像元值。
首先，找到与目标图像云像元相对应的参考图像中的晴空像元。读取目标图像各像元并找到带有云像元时，继而转向参考图像1(成像日期距目标图像有最近)相同位置的像元，确定该像元是晴空像元还是云像元的。如果是云覆盖像元，继续读取另一个参考图像即参考图2(成像日期第二目标最近)相应位置的像元等，直到在参考图像上找到与目标图像云像元在同一地理位置的晴空像元。
下一步是确定目标图像与参考图像之间的光谱相关性。
在目标图像云覆盖像元周围找相邻晴空像元，在参考图像中围绕找到与之相应的晴空像元；分别围绕C(0,i)和C(j,i)设置2个足够大的窗口，窗口边长像元数为奇数，如5*5，7*7，...，101*101，...，并检查两个窗口中每个像元对的值,然后计算出参考图像中这些像元在波段k与像元p之间的光谱差：
SDi,k＝|Di,k-Dp,k| (2)
SDi,k是参考图像波段k中周围晴空像元i与像元p的光谱差，SDi,k和Dp,k是像元i和像元p波段k的值。
在目标图像的云覆盖像元周围尽可能最小的窗口下找光谱差值SDi,k小于阈值n的晴空点，并将其光谱差值进行排序，取前10个光谱最相似的相邻晴空像元。调整窗口大小，直到在两个尽可能小的窗口中得到10个相应的符合阈值条件的晴空像元；如果超出最大窗口设置范围m仍然不能找到相应的10个像元，将阈值n逐渐从5，10，15，20，30，40，50…100…扩大，并继续逐窗口大小寻找,直到找到符合条件的点为止。
从这些周围的像元中，找出与像元p的光谱差异最小的10个像元，也就是说具有最小DN差值SDi(参考上述公式2)，定义为种子像元：
RIk＝(RI1,k,RI2,k,RI3,k,RI4,k,RI5,k,RI6,k,RI7,k,RI8,k,RI9,k,RI10,k)
(3)
RIk为参考图像种子像元在波段k的数据序列，RI1,k,RI2,k,RI3,k,RI4,k,RI5,k,RI6,k,RI7,k,RI8,k,RI9,k,RI10,k分别为种子像元1，2，3，4…9，10的像元值；
同样的，在目标图像相应的种子像元定义为相关像元，如下:
TIk＝(TI1,k,TI2,k,TI3,k,TI4,k,TI5,k,TI6,k,TI7,k,TI8,k,TI9,k,TI10,k) (4)
TIk是目标图像种子像元在波段k的数据序列，TI1,k,TI2,k,TI3,k,TI4,k,TI5,k,TI6,k,TI7,k,TI8,k,TI9,k,TI10,k分别是种子像元1，2，3，4…9，10的像元值；
数据序列TIk为因变量，RIk为自变量，运用这两个变量我们可以建立如下线性回归方程：
YTI,k＝a+b XRI,k (5)
YTI,k是因变量，对应于目标图像波段k的像元值，XRI,k是自变量，对应于参考图像波段k的DN值，a和b分别是方程的回归系数。基于该方程，我们可以应用参考图像波段k的晴空像元p来估算相对应的目标图像云像元值，如下:
Dq,k＝＝a+b Dp,k (6)
Dq,k是在目标图像云覆盖像元q在波段k的估算像元值，Dp,k是参考图像中波段k相应的晴空像元的真实像元值。由于上述方法是通过线性回归建模(LRM)来估算云覆盖像元值，我们将其定义为用于可见光和近红外波段去云的LRM方法。
与现有技术相比，本发明的优势在于：1、不需要大量样本数据进行训练学习，对先前成像的依赖性不高，对成像间隔要求不是太高。2、所用的参考图像，不需要像应用其他传感器图像的复杂的处理过程，不会受到到下垫面类型影响，提高了云和云阴影像元估算效率和精度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为基于时空信息可见光近红外波段厚云及阴影去除的流程图；
图2为基于时空信息可见光近红外波段厚云及阴影去除的示意图；
图3为Landsat8 OLI图像南京附近实验区位置示意图(p120r38_20160429)；
图4为Landsat8 OLI参考图像1((p120r38_20160328)和参考图像2((p120r38_20160312)；
图5为基于时空信息的OLI图像可见光近红外波段厚云及阴影去除方法的示例图1(波段组合432，对应图4的a-b区)；
图6为基于时空信息的OLI图像可见光近红外波段厚云及阴影去除方法的示例图2(波段组合432，对应图4的c-d区)；
图7为基于时空信息的OLI图像可见光近红外波段厚云及阴影去除方法的示例图3(e波段组合432，f波段组合543，对应图4的e区)。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
如图1至图7所示，以2016年4月29日南京区域的有云覆盖的Landsat 8OLI图像可见光和近红外波段为例，结合附图对本发明作进一步详细的说明。
基于图1所示的流程，南京附近区域云区重建通过以下步骤完成：
步骤一：根据实现选定的南京研究区范围，采集Landsat 8OLI的L1级卫星图像数据。确定目标图像，并筛选若干可获取的成像时间距离目标图像最近的同传感器参考图像若干，保证筛选的参考图像的晴空区域集合能够在空间上覆盖所有目标图像所有的云、云阴影像元范围。最终确定目标图像和参考图像条带号为p120r38，其成像时间为20160429，其他两景参考图像成像时间分别是20160328和20160312。如图5所示，所选的两景参考图像的云覆盖区域较少，大部分为晴空图像，可满足目标图像厚云及阴影去除的要求。根据云区大小和下垫面的情况分为5个1000*1000像元的区块大小的研究区a-e，分布图见图4。
步骤二：对目标图像和参考图像进行几何校正；
对南京区域三个时相的OLI多光谱图像(30米)进行几何校正，以确保它们相对应的像元在空间上处于相同的几何位置。
采集1:10万地形图、全球Landsat7 ETM+正射图像图，建立1980西安平面坐标系、高斯－克吕格投影(6度分带)的图像数学基础。在此基础上，基于参考图与待校正的图像，选取地面控制点(GCP)，选取的原则包括：
(1)在图像放大2～3倍的条件下，选取覆盖全图像区域的控制点。
(2)对于平原地区图像，每景(约7580像元×7750像元)图像至少选取40个控制点；对于山区图像，每景图像至少选取50个控制点。目标图像地出南京附近，下点面有低山区、平原、水域、城镇等，选择45个控制点。
(3)在接近正交的道路等线状地物的交点、地物拐点上、图像清晰的明显地物点或固定的点状地物等位置选取控制点。不允许在高层建构筑物、河流等地物上选取控制点。
基于选取的控制点，利用三次多项式纠正模型分别对三个时间的多光谱OLI图像的几何变形进行纠正，满足后续数据处理实验要求。设置30米的采样间隔，采用最邻近法对几何校正后的图像进行重采样。
步骤三：对目标图像和参考图像进行云、阴影等检测；
基于图3所示流程，对于L1级别OLI数据进行厚云及其阴影检测，利用商业软件分别生成三个时间的OLI图像的云和云阴影掩膜C(0),C(-1),C(-2)，其中云和云阴影像元为0，晴空像元值为1。为提高线性回归建模(LRM)算法精度，降低漏判云或阴影像元的影响，云检测后生成掩膜文件，对云和阴影区域进行5个像元缓冲，生成最终的云和云阴影掩膜文件。
步骤四：目标图像中云和云阴影目标像元的估算；
将选定的五个1000*1000像元的研究区a-e的目标图像、参考图像和云掩膜分别都裁剪出。每个研究区都进行如下操作完成云、云阴影检测。将目标图像(p120r38_20160428)定义为T(0)，而参考图像(p120r38_20160329)设为T(-1)，参考图像(p120r38_20160412)设为T(-2)。为估算目标图像中所有可见光和近红外波段云像元的DN值，我们将选择OLI图像的2,3,4和5波段进行实验。T(0,i,k)，T(-1,i,k)，T(-2,i,k)分别表示目标图像T(0)和参考图像T(-1)，T(-2)的第i个像元在波段k的DN值，以此类推。估算云覆盖像元的DN值以去除云和云阴影像元方法的步骤如下(图1-2)：
(1)读取图像T(0),T(-1),T(-2)根据图3所示流程进行云、云阴影等信息提取，生成对应的云掩膜图像为C(0),C(-1),C(-2),所有像元都被划分为晴空像元(设为1)和厚云、云阴影像元(设为0)。
(2)确定目标图像中的云覆盖像元。流程如下：读取目标图像的云掩膜图像C(0,i)来确定目标图像T(0)的第i个像元是晴空还是云覆盖，若值为1，那么该像元是晴空像元。然后继续读取云掩膜图像下一个像元(i+1)的值C(0,i+1)，确定是晴空还是云覆盖。继续推进该过程直到找到一个目标图像中的云覆盖像元C(0,i)。
(3)确定参考图像相应的晴空像元C(0,i)。流程如下：对j＝-1，读C(-1,i)值来确定该像元是否为晴空。如果C(-1,i)＝0，则为云覆盖像元，就进入下一个参考图像，读取相应像元C(-2,i)云覆盖情况。
(4)在目标图像T(0)的云覆盖像元C(0,i)周围找相邻晴空像元，在参考图像T(j)中围绕C(j,i)找到与之相应的晴空像元；通过以下流程实现：分别围绕C(0,i)和C(j,i)设置2个足够大的窗口，窗口边长像元数为奇数，如5*5，7*7，...，301*301，最大窗口设置为301,并检查两个窗口中每个像元对的值,在目标图像T(0)的云覆盖像元C(0,i)周围尽可能最小的窗口下找光谱差值SDi,k小于阈值n的晴空点。
(5)可运用公式(2)计算波段k中的相邻像元与像元C(j,i)之间的光谱差，调整窗口大小，直到在两个尽可能小的窗口中得到10个相应的符合阈值条件的晴空像元；如果超出最大窗口设置范围301仍然不能找到相应的10个像元，将阈值n逐渐从5，10，15，20，30，40，50扩大，并继续逐窗口大小寻找,直到找到符合条件的点为止。
(6)建立波段k目标图像数据序列TIk和参考图像数据序列RIk。使用步骤(5)的10个像元，分别获取它们在目标图像T(0)和参考图像T(j)的像元值来建立数据序列TIk和RIk。
(7)建立数据序列TIk和RIk之间的光谱回归关系。可通过最小二乘回归分析方法来完成，参见方程(5)YTI,k＝a+bXRI,k，YTI,k为因变量，XRI,k则为独立自变量。
(8)估算目标图像中云覆盖像元的DN值。用参考图像T(j)中C(i,j)＝1的晴空像元对应的DN值，应用公式(6)Nq,k＝＝a+bNp,k，可实现目标图像中波段k云覆盖像元的DN值的估算。
(9)循环步骤(4)到(8)估算目标图像其他波段云覆盖像元的DN值。
(10)重复上述步骤(3)到(9)，分别对OLI图像的2,3,4,5波段的DN值估算，直到去除目标图像中的所有的云像元。
(11)重复上述步骤(2)到(10)，直到研究区的云和云阴影全部去除，研究区晴空图像全部完成。
Claims (5)
1.一种可见光和近红外波段厚云及阴影去除方法，其特征在于：该方法的步骤是：
S1、采集同一研究区不同时相的卫星图像数据，定义一个有云和阴影覆盖的图像为目标图像T(0)，针对目标图像区域，选择若干不同时相的同种传感器影像作为参考图像，记T(-1)，T(-2)，…，T(j)；对目标图像和参考图像进行几何校正，确保它们相对应的像元在空间上处于相同的几何位置；
S2、读取图像T(0)，T(-1)，T(-2)，…，T(j)进行云检测，对应的云掩膜图像设为C(0)，C(-1)，C(-2)，…，C(j)，所有像元都被划分为晴空像元和云像元，晴空像元设为1，云像元设为0；
S3、确定目标图像T(0)中的云覆盖像元；
读取云掩膜图像C(0,i)来确定目标图像T(0)的第i个像元是晴空还是云覆盖，若值为1，那么该像元是晴空像元p；然后继续读取云掩膜图像下一个像元(i+1)的值C(0,i+1)，确定是晴空还是云覆盖；继续推进该过程直到找到一个云覆盖像元C(0,i)；
S4、确定参考图像相应的晴空像元；
对j＝-1，读C(-1,i)值来确定该像元是否为晴空；如果C(-1,i)＝0，则为云覆盖像元，就进入下一个参考图像j＝j-1，读取相应像元C(j,i)的云覆盖情况，以此类推，直到找到一个晴空像元C(j,i)＝1；
S5、在参考图像T(j)波段k围绕C(j,i)找光谱最相似的10个相邻像元；
定义SDi,k是参考图像波段k中周围晴空像元i与像元p的光谱差，公式如下：
SDi,k＝|Di,k-Dp,k| (1)
SDi,k和Dp,k是像元i和像元p在波段k的值；
在目标图像T(0)的云覆盖像元C(0,i)周围找相邻晴空像元，在参考图像T(j)中围绕C(j,i)找到与之相应的晴空像元；分别围绕C(0,i)和C(j,i)设置2个足够大的窗口，窗口边长像元数为奇数，如5*5，7*7，9*9，...，301*301，...，最大窗口设置为m,并检查两个窗口中每个像元对的值,在目标图像T(0)的云覆盖像元C(0,i)周围尽可能最小的窗口下找光谱差值SDi,k小于阈值n的晴空点，并将其光谱差值进行排序，取前10个光谱最相似的相邻晴空像元；调整窗口大小，直到在两个尽可能小的窗口中得到10个相应的符合阈值条件的晴空像元；如果超出最大窗口设置范围m仍然不能找到相应的10个像元，将阈值n逐渐从5，10，15，20，30，40，50…100…扩大，并继续逐窗口大小寻找,直到找到符合条件的点为止；
S6、建立波段k目标图像数据序列TIk和参考图像数据序列RIk；使用步骤S5的10个像元，分别获取它们在目标图像T(0)和参考图像T(j)的像元值来建立数据序列TIk和RIk；
TIk＝(TI1,k,TI2,k,TI3,k,TI4,k,TI5,k,TI6,k,TI7,k,TI8,k,TI9,k,TI10,k) (2)
TIk是目标图像种子像元在波段k的数据序列，TI1,k,TI2,k,TI3,k,TI4,k,TI5,k,TI6,k,TI7,k,TI8,k,TI9,k,TI10,k分别是种子像元1，2，3，4…9，10的像元值；
RIk＝(RI1,k,RI2,k,RI3,k,RI4,k,RI5,k,RI6,k,RI7,k,RI8,k,RI9,k,RI10,k) (3)
RIk为参考图像种子像元在波段k的数据序列，RI1,k,RI2,k,RI3,k,RI4,k,RI5,k,RI6,k,RI7,k,RI8,k,RI9,k,RI10,k分别为种子像元1，2，3，4…9，10的像元值；
S7、通过最小二乘回归分析方法，建立数据序列TIk和RIk之间的光谱回归关系；
YTI,k＝a+bXRI,k (4)
其中，YTI,k是因变量，对应于目标图像波段k的像元值，XRI,k是自变量，对应于参考图像波段k的DN值，a和b分别是方程的回归系数；
S8、估算目标图像中云覆盖像元的值；用参考图像T(j)中C(i,j)＝1的晴空像元对应的像元值，实现目标图像中波段k云覆盖像元值的估算；
Dq,k＝＝a+bDp,k (5)
Dq,k是在目标图像云覆盖像元q在波段k的估算像元值，Dp,k是参考图像中波段k相应的晴空像元的真实像元值；
S9、循环步骤S5到步骤S8，估算目标图像其他波段云覆盖像元值；
S10、重复上述步骤S3到步骤S9，对所有云覆盖像元值进行估算，直到去除目标图像中的所有的云像元。
2.如权利要求1所述的可见光和近红外波段厚云及阴影去除方法，其特征在于：在步骤S1中，参考图像的时间越接近目标图像越好，且保证所有参考图像的晴空像元的合集能够全部覆盖目标图像中的云、阴影和雪的分布范围。
3.如权利要求1所述的可见光和近红外波段厚云及阴影去除方法，其特征在于：在步骤S1中，对目标图像与参考图像的进行几何校正时误差小于1个像元。
4.如权利要求1所述的可见光和近红外波段厚云及阴影去除方法，其特征在于：所述目标图像与参考图像采集于Landsat 8 OLI的L1级卫星图像数据。
5.如权利要求1所述的可见光和近红外波段厚云及阴影去除方法，其特征在于：利用目标云覆盖像元周围附近晴空像元的建立用于云像元值估算的线性回归方程。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811119642.8A CN109308688B (zh) | 2018-09-25 | 2018-09-25 | 一种可见光和近红外波段厚云及阴影去除方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811119642.8A CN109308688B (zh) | 2018-09-25 | 2018-09-25 | 一种可见光和近红外波段厚云及阴影去除方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109308688A CN109308688A (zh) | 2019-02-05 |
| CN109308688B true CN109308688B (zh) | 2021-06-25 |
Family
ID=65224545
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811119642.8A Active CN109308688B (zh) | 2018-09-25 | 2018-09-25 | 一种可见光和近红外波段厚云及阴影去除方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109308688B (zh) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111583132B (zh) * | 2020-04-20 | 2023-05-02 | 国家卫星气象中心(国家空间天气监测预警中心) | 一种去除遥感影像异常条带式噪声的方法、装置、设备及介质 |
| CN111898503B (zh) * | 2020-07-20 | 2021-02-26 | 中国农业科学院农业资源与农业区划研究所 | 基于云覆盖遥感影像和深度学习的作物识别方法和系统 |
| CN112102180B (zh) * | 2020-08-21 | 2022-10-11 | 电子科技大学 | 一种基于Landsat光学遥感图像的云识别方法 |
| CN113160183B (zh) * | 2021-04-26 | 2022-06-17 | 山东深蓝智谱数字科技有限公司 | 一种高光谱数据处理方法、设备及介质 |
| CN115082452B (zh) * | 2022-07-26 | 2022-11-04 | 北京数慧时空信息技术有限公司 | 基于云及阴影的遥感影像质量定量评价方法 |
| CN117408949A (zh) * | 2023-09-20 | 2024-01-16 | 宁波大学 | 一种季节性动态阈值的云及云阴影检测方法及装置 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101894382A (zh) * | 2010-07-23 | 2010-11-24 | 同济大学 | 一种集成LiDAR点云的卫星立体影像阴影计算方法 |
| CN102385694A (zh) * | 2010-09-06 | 2012-03-21 | 邬明权 | 一种基于地块的作物品种高光谱识别方法 |
| CN103778618A (zh) * | 2013-11-04 | 2014-05-07 | 国家电网公司 | 一种可见光图像和红外图像的融合方法 |
| CN105046242A (zh) * | 2015-08-24 | 2015-11-11 | 山东省农业可持续发展研究所 | 一种基于Landsat 8影像二维特征空间的芦笋种植面积提取方法 |
| CN105678777A (zh) * | 2016-01-12 | 2016-06-15 | 武汉大学 | 一种多特征联合的光学卫星影像云与云阴影检测方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20140240477A1 (en) * | 2013-02-26 | 2014-08-28 | Qualcomm Incorporated | Multi-spectral imaging system for shadow detection and attenuation |
-
2018
- 2018-09-25 CN CN201811119642.8A patent/CN109308688B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101894382A (zh) * | 2010-07-23 | 2010-11-24 | 同济大学 | 一种集成LiDAR点云的卫星立体影像阴影计算方法 |
| CN102385694A (zh) * | 2010-09-06 | 2012-03-21 | 邬明权 | 一种基于地块的作物品种高光谱识别方法 |
| CN103778618A (zh) * | 2013-11-04 | 2014-05-07 | 国家电网公司 | 一种可见光图像和红外图像的融合方法 |
| CN105046242A (zh) * | 2015-08-24 | 2015-11-11 | 山东省农业可持续发展研究所 | 一种基于Landsat 8影像二维特征空间的芦笋种植面积提取方法 |
| CN105678777A (zh) * | 2016-01-12 | 2016-06-15 | 武汉大学 | 一种多特征联合的光学卫星影像云与云阴影检测方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109308688A (zh) | 2019-02-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109308688B (zh) | 一种可见光和近红外波段厚云及阴影去除方法 | |
| Chen et al. | A practical approach to reconstruct high-quality Landsat NDVI time-series data by gap filling and the Savitzky–Golay filter | |
| Kuc et al. | Sentinel-2 imagery for mapping and monitoring imperviousness in urban areas | |
| CN112836610B (zh) | 一种基于遥感数据的土地利用变化与碳储量定量估算方法 | |
| EP3002729B1 (en) | Automated geospatial image mosaic generation with radiometric normalization | |
| Makarau et al. | Haze detection and removal in remotely sensed multispectral imagery | |
| Ma et al. | Coupling urbanization analyses for studying urban thermal environment and its interplay with biophysical parameters based on TM/ETM+ imagery | |
| Ishola et al. | Dynamics of surface urban biophysical compositions and its impact on land surface thermal field | |
| CN111738165A (zh) | 一种从高分辨率无人机可见光遥感影像中提取单株植物冠层的方法 | |
| WO2015195484A1 (en) | Automated geospatial image mosaic generation with multiple zoom level support | |
| CN113033670A (zh) | 一种基于Sentinel-2A/B数据的水稻种植面积提取方法 | |
| CN110991430A (zh) | 基于遥感图像的地物识别及覆盖率计算方法及系统 | |
| Militino et al. | Filling missing data and smoothing altered data in satellite imagery with a spatial functional procedure | |
| Morsy et al. | Impact of land use/land cover on land surface temperature and its relationship with spectral indices in Dakahlia Governorate, Egypt | |
| Mohamed et al. | Change detection techniques using optical remote sensing: a survey | |
| Liu et al. | Thick cloud removal under land cover changes using multisource satellite imagery and a spatiotemporal attention network | |
| Latif et al. | Preprocessing of low-resolution time series contaminated by clouds and shadows | |
| Sui et al. | Processing of multitemporal data and change detection | |
| Falcone et al. | Mapping impervious surface type and sub‐pixel abundance using hyperion hyperspectral imagery | |
| Seo et al. | Local-based iterative histogram matching for relative radiometric normalization | |
| Colditz | Time series generation and classification of MODIS data for land cover mapping | |
| CN110427961B (zh) | 一种基于规则和样本融合的建筑信息提取方法及系统 | |
| Danoedoro et al. | Preliminary study on the use of digital surface models for estimating vegetation cover density in a mountainous area | |
| Pérez-Goya et al. | RGISTools: Downloading, Customizing, and Processing Time Series of Remote Sensing Data in R | |
| Sohn et al. | Shadow-effect correction in aerial color imagery |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |