[发明专利]一种可见光和近红外波段厚云及阴影去除方法

摘要

本发明公开了一种可见光和近红外波段厚云及阴影去除方法，目标晴空图像的构建主要是应用多时相图像的时空信息对云和云阴影覆盖像元值的进行估算，包括以下步骤：1)目标图像与参考图像的几何校正；2)对目标图像和参考图像进行云和阴影检测；3)云和阴影覆盖像元值(DN值、反射率等)估算。本发明根据云分布的动态变化性、遥感平台传感器的成像周期性以及下垫面在短时间内无明显变化等特点，基于多时相图像的时空信息，利用线性回归数学模型的逐波段逐像元进行定量估算云下像元值，建立了一套实现逐像元逐波段自动化厚云及其阴影去除的流程方法，提高了云和云阴影去除精度，实现全自动厚云及阴影去除。

主权项

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)的像元值来建立数据序列RIk和TIk；RIk＝(RI1,k,RI2,k,RI3,k,RI4,k,RI5,k,RI6,k,RI7,k,RI8,k,RI9,k,RI10,k)    (2)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)    (3)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的像元值；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，对所有云覆盖像元值进行估算，直到去除目标图像中的所有的云像元。