一种对高分辨率光学遥感图像的自适应云区检测方法与流程

本发明涉及遥感图像处理技术领域，具体而言，涉及一种对高分辨率光学遥感图像的自适应云区检测方法。

背景技术：

光学遥感卫星传感器经常采用的可见近红外波段(0.38-0.90μm)无法穿透云层，在实际应用时难以获取云层覆盖区域内的地物信息，国际云气候学计划(internationalsatellitecloudclimatologyproject，isccp)根据其建立的isccp-fd(isccp-fluxdata)辐射数据集对遥感数据的云覆盖量估测结果显示，全球每年获取的光学遥感数据云覆盖约为66％[1]。云层遮挡密集的情况会严重影响遥感数据的判读、解译精度以及应用。因此，云区识别技术作为遥感图像数据应用分析前一个基础的重要步骤，既能提高影像数据的利用率，同时也节省了数据收集成本。

目前，中国高分辨率对地观测系统已发射并获取光学遥感数据的卫星包括高分一号(gf-1)、高分二号(gf-2)以及高分四号(gf-4)。高分系列光学卫星实现了全天候全天时对地观测，极大丰富了高分辨率国产空间数据集内容。gf-1多光谱数据成像波段包括四个波段：蓝波段(0.45-0.52μm)、绿波段(0.52-0.59μm)、红波段(0.63-0.69μm)以及近红外波段(0.77-0.89μm)。国际上目前使用的一些中低分辨率星载影像，如avhrr、modis、landsat系列和sentinel等数据，均包含一个或多个有利于识别云雾的热红外和水汽吸收波段[2-4]。在有限波段的局限下，对高分辨率多光谱数据进行云区判定，易造成与影像中其他具有高反射率的地表物体(如道路、房屋、裸地等)发生混淆。因此，如何对日益增多的高分辨率卫星遥感数据进行有效且高效的云区自动识别，仍是一个十分具有挑战且亟待解决的问题。

在云区识别方面，常见的云区识别方法通常包括基于多时相影像的云区识别和基于单幅影像的云区识别。基于多时相影像的云区识别方法，需要采用不同时相的无云数据作为参考，识别精度较高；基于单幅影像的云区识别方法无需参考数据，直接对影像本身进行处理即可生成云区掩膜数据，因此具有更强的适用性和实效性。目前，同态滤波法仍是最常用的针对单幅遥感影像的云区识别方法[5-6]，同态滤波通过高通滤波器在频域内抑制云区分布的低频区域，使得云区辐射信息被削弱，非云区辐射信息得到增强。影像经过同态滤波后，亮度较高的云层反射率降低，而原来相对较暗的无云区域反射率变高，经过与原始输入进行比对，便可得出云区的掩膜结果。但同时，同态滤波的去云效果极大地依赖高通滤波器的截止频率，前人多数的研究都对截止频率采用经验值，显然经验截止频率无法适应大批量云区层次复杂的遥感数据处理，且对不同厚度的云区，如厚云和薄云等，采用经验值无法达到最佳识别效果。

[1]：zhangy,rossowwb,lacisaa,etal.calculationofradiativefluxesfromthesurfacetotopofatmospherebasedonisccpandotherglobaldatasets:refinementsoftheradiativetransfermodelandtheinputdata[j].journalofgeophysicalresearch:atmospheres,2004,109(d19).

[2]：ackermansa,strabalaki,menzelwp,etal.discriminatingclearskyfromcloudswithmodis[j].journalofgeophysicalresearch:atmospheres,1998,103(d24):32141-32157.

[3]：derrienm,farkib,harangl,etal.automaticclouddetectionappliedtonoaa-11/avhrrimagery[j].remotesensingofenvironment,1993,46(3):246-267.

[4]：zhuz,woodcockce.object-basedcloudandcloudshadowdetectioninlandsatimagery[j].remotesensingofenvironment,2012,118:83-94.

[5]：shenh,lih,qiany,etal.aneffectivethincloudremovalprocedureforvisibleremotesensingimages[j].isprsjournalofphotogrammetryandremotesensing,2014,96:224-235.

[6]：lvh,wangy,sheny.anempiricalandradiativetransfermodelbasedalgorithmtoremovethincloudsinvisiblebands[j].remotesensingofenvironment,2016,179:183-195.

技术实现要素：

本发明提供一种对高分辨率光学遥感图像的自适应云区检测方法，以有效解决现有技术中同态滤波算法无法根据输入影像自适应调整滤波器截止频率的问题。

为达到上述目的，本发明提供了一种对高分辨率光学遥感图像的自适应云区检测方法，其中，高分辨率光学遥感图像为高分一号遥感图像，其包括以下步骤：

s1：对高分辨率光学遥感图像f，计算其云雾厚度图f1；

s2：根据f1的径向能量谱与同态滤波过程中应用的高通滤波器的截止频率f之间的关系，计算f的值；

s3：以f作为同态滤波过程中高通滤波器的截止频率，对f1进行同态滤波，得到滤波后的图像f2；

s4：计算f2中每一像元的白度指数并从中滤除白度指数大于一预设阈值的像元，得到图像f3；

s5：选用圆形结构元素，对f3先进行闭运算再进行开运算，以对其进行形态学优化，得到云区识别结果f4。

在本发明的一实施例中，步骤s1中，先进行以下计算：

其中，i(x,y)＝2ib(x,y)-0.95ig(x,y)，ib(x,y)、ig(x,y)分别为高分辨率光学遥感图像f的蓝色波段和绿色波段，i(x,y)为输入波段，w(x′,y′)为中心位于点(x′,y^′)的w×w的局部窗口，

之后再利用中值滤波器对htm(x0,y0)进行平滑，并采样到原始图像大小，得到云雾厚度图f1。

在本发明的一实施例中，所述中值滤波器采用3×3滤波窗口。

在本发明的一实施例中，w＝5。

在本发明的一实施例中，f1的径向能量谱r(r)与同态滤波过程中应用的高通滤波器的截止频率f之间的关系如下：

其中，r(u,v)、i(u,v)分别为对云雾厚度图f1进行傅里叶变换后得到的复数频谱的实部与虚部。

在本发明的一实施例中，步骤s5中选用直径为7个像元的圆形结构元素。

在本发明的一实施例中，步骤s4包括如下步骤：

s41：对f2中的可见光波段进行光谱密度积分，计算出云亮度fbr，

其中，λi为波长，λmax、λmin分别为可见光区域内的最大波长和最小波长，单位均为nm，bv为可见光的波段范围，f(λi)表示波长为λi波段的像元值；

s42：计算白度指数fwh，

其中，e(λi)＝|f(λi)-fbr|；

s43：将fwh＞β的像元对应的物体判断为非白色物体并将其从f2滤除，得到过滤后的图像f3。

在本发明的一实施例中，β＝0.1。

在本发明的一实施例中，步骤s2中，同态滤波过程中应用的高通滤波器为高斯差分高通滤波器，高斯差分高通滤波器的高频增益为1，低频增益为0.05。

本发明提供的对高分辨率光学遥感图像的自适应云区检测方法通过对大量高分一号(gf-1)遥感图像进行试验分析，得出了图像的径向能量谱与同态滤波过程中应用的高通滤波器的截止频率之间的量化关系，能够根据不同的输入图像自适应调整截止频率，实现gf-1遥感图像的批量云区识别。本发明能够对高分辨率光学遥感图像进行批量化高效的云区检测，总体识别精度达93.81％。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1a-图1d为选用圆形结构元素对云区进行形态学优化的具体实施例，图1a为7个像元的圆形结构元素，图1b为云区示例f3，图1c为f3经过闭运算处理后的结果，图1d为图1c再经过开运算处理的结果，即图像f4；

图2a-图2c为利用白度指数以及形态学运算对同态滤波结果f2进行优化的实例，图2a为原始输入图像f，图2b为同态滤波结果f2，图2c为经过后处理优化的图像f4；

图3a为8幅不同云量覆盖的高分辨率光学遥感图像(1024×1024)，8幅图像云量覆盖从左上至右下依次减少；

图3b为采用本发明的方法对图3a的云区检测结果，图3a中为图像f，图3b中为图像f4，图3a、图3b中位于同一位置的图像一一对应；

图4中第一列的前两幅图为gf-1卫星中的8m分辨率多光谱高分相机(pms)采集到的遥感图像f，后两幅图为16m分辨率多光谱宽幅相机(wfv)采集到的遥感图像f，对第一列原始影像采用本发明的方法处理后得到第二列的云区掩膜结果f4，第三列为从第一列的遥感影像中除去云区后的结果影像。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在具体实施例之前，先介绍一下本发明中应用到的同态滤波的原理：

在云雾区域，遥感数据成像模型包括两个部分：太阳辐射经云层反射部分以及太阳辐射经地物反射后再穿透云层的部分。用f(x,y)表示传感器接收到的信号，fi(x,y)为云层反射率，能够反映云的分布，fr(x,y)表示地物反射率，则云的成像模型可简化为：

f(x,y)＝fi(x,y)·fr(x,y)

在频率域中，表示云层分布的fi(x,y)集中于低频，反映影像细节内容的fr(x,y)分布在高频区域。同态滤波通过低频成分与高频成分的分离，并抑制云分布的低频分量并增强高频信息，实现空域的去云处理。

通过对上式两边取对数得：

in(f(x,y))＝in(fi(x,y))+in(fr(x,y))

再进行傅里叶变换进行空-频转换：

f{in(f(x,y))}＝f{in(fi(x,y))}+f{in(fr(x,y))}，简记为f(u,v)＝l(u,v)+h(u,v)

对频谱采用高通滤波器，在频域上抑制云层反射率的低频成分:

s(u,v)＝l(u,v)filter(u,v)+h(u,v)filter(u,v)

其中filter(u,v)为高通滤波器的传递函数，滤波函数采用高斯差分高通滤波器，

式中：γh为高频增益，γl为低频增益，通常设置γh≥1且γl<<1，本实施例中采用γh＝1，γl＝0.05；d(u^,v)代表点(u，v)到傅里叶变换中心点的距离；d0为高通滤波器截止频率。

将滤波后的结果进行傅立叶反变换：

s(x,y)＝in′(fi(x,y))+in′(fr(x,y))

再进行指数变换：

f′(x,y)＝fi′(x,y)+fr′(x,y)。

将变换结果f′(x,y)拉伸到原来输入数据的幅度：

其中maxi和mini分别为输入数据f(x,y)的最大值与最小值，maxo和mino分别为f′(x,y)的最大值与最小值，并由此得到云掩膜结果：

对于同态滤波过程中应用的高通滤波器的截止频率，说明如下：

截止频率取值不同会带来不一样的滤波结果，同态滤波的关键是要根据不同的输入确定合适高通滤波器的截止频率。若截止频率取值较低，高通滤波器无法完全抑制低频的亮度信息；若截止频率取值较高，地表物体的细节信息会被高通滤波器过度削减。

在高通滤波器的截止频率点处，滤波器会极大衰减截止频率前的低频信号，使输出信号强度衰减为输入信号的一半(-3db)，同时允许截止频率后的信号通过。截止频率的取值极大影响了同态滤波的结果，而它的选择与照度场和反射场对应的频谱能量有关，需要大量实践来确定。本发明从图像的照度场和反射场的频谱能量出发，通过分析照度分量的能量占比，从而选择合适的截止频率。

图像经过二维快速傅里叶变换后，将原点平移至频谱中心，此时频率原点值表示影像的平均灰度级，频率原点值越大表示图像平均灰度越高，反之图像平均灰度越低。傅里叶频谱中心向四周过渡分别表示频率由低频逐渐向高频变化。此时输入影像f(x,y)的频谱f(u,v)表示有：

f(u,v)＝|f(u,v)|e^jθ

其中，幅值谱表示为r(u,v)和i(u,v)分别为傅里叶变换复数频谱的实部与虚部，相角假设输入影像尺寸为m×n，则有u＝m/2和v＝n/2。假设光照条件是均匀的，则光照场的傅里叶频谱只有直流分量，随着光照不均匀度增加，开始出现反射变量，在傅里叶频谱上表示为除了直流分量，同时增加了谐波分量。此时，傅里叶变换径向频谱r(r)表示影像频域谐波圆环内各个方向的总能量之和，反映频域中谐波频率圆环内的能量大小，随着谐波频率逐渐趋于奈奎斯特频率，r(r)幅值越来越小并衰减趋于0。r(r)和其占总能量的百分比r(r)％表示为：

式中：θ＝0o～360o。云量较多时图像径向频谱整体能量水平越高，平均灰度值r[0]也越高，此时，体现云的照度分量高度集中于低频部分，因此截止频率设置在接近圆心的低频处较为合理；无云覆盖时图像径向频谱能量水平较低，谐波能量收敛趋近于0的速度更慢，此时，照度分量较弱，体现反射分量的谐波能量较强，则截止频率应设置得更高来充分抑制地物信息，避免造成欠滤波。

本案发明人经实践发现，图像的平均灰度能量占比r(r)％与截止频率呈现正相关，即云量覆盖较多(r(r)％低)其截止频率较低，当云量较少时(r(r)％高)，其截止频率较高。本案发明人经过对79幅实例数据的截止频率进行分析，建立了r(r)％与截止频率处径向频谱和的定量模型，从而可得出每景影像对应的截止频率值。本发明的回归模型精度较优，决定系数r²达0.9275，模型可靠，具体见以下详述。

本发明中的高分辨率光学遥感图像为高分一号(gf-1)遥感图像，gf-1卫星搭载有两台2m分辨率全色和8m分辨率多光谱高分相机(pms)，4台16m分辨率多光谱宽幅相机(wfv)。

本发明提供的对高分辨率光学遥感图像的自适应云区检测方法包括以下步骤：

s1：对高分辨率光学遥感图像f，计算其云雾厚度图f1；

于本发明中，假设在某场景中总是存在某些反射值极低的像元，其像元值接近于0。由于受到云雾反射的影响，使上述暗像元的像素值并未完全的“纯黑”，则认为该类暗像元的像元值客观上反映了场景中云雾的厚度。

在本实施例中，步骤s1中，先进行以下计算：

上式是通过设置局部窗口由图像从左到右从上到下寻找暗目标，从而得到图像的htm结果。其中，i(x,y)＝2ib(x,y)-0.95ig(x,y)，ib(x,y)、ig(x,y)分别为高分辨率光学遥感图像f的蓝色波段和绿色波段，i(x,y)为输入波段，通常采用极易受云雾影响的蓝波段作为输入波段i(x,y)，但由于蓝波段计算得到的htm结果会造成云雾厚度过检测，此处采用线性差值合成波段i(x,y)＝2ib(x,y)-0.95ig(x,y)，以该合成波段计算的htm结果能更好的抑制地面高反射地物的亮度，缓和过度检测的效果。w(x′,y′)为中心位于点(x′,y′)的w×w的局部窗口，本实施例中w＝5，即采用5×5的窗口。

之后再利用中值滤波器对htm(x0,y0)进行平滑，此处的中值滤波器较佳为采用3×3滤波窗口，并采样到原始图像大小，得到云雾厚度图f1。

s2：根据f1的径向能量谱与同态滤波过程中应用的高通滤波器的截止频率f之间的关系，计算f的值；

本实施例中，f1的径向能量谱r(r)与同态滤波过程中应用的高通滤波器的截止频率f之间的关系如下：

其中，r(u,v)、i(u,v)分别为对云雾厚度图f1进行傅里叶变换后得到的复数频谱的实部与虚部。

s3：以f作为同态滤波过程中高通滤波器的截止频率，对f1进行同态滤波，得到滤波后的图像f2；

通过同态滤波自适应处理，能将高亮反射的云区准确识别，由于其他非云区的高反射率地物同样存在被过度识别的风险，因此，下一步采用白度指数和形态学算法对云掩膜初步结果进行优化后处理。

s4：计算f2中每一像元的白度指数并从中滤除白度指数大于一预设阈值的像元，得到图像f3；

白度指数基于云最主要的特性，即呈现高亮白色并且可见光光谱曲线平坦。高亮白色表示云区的光谱曲线在可见光波段的值相对较高，通过对可见光波段进行光谱密度的积分，可得出云亮度。

本实施例中，步骤s4包括如下步骤：

s41：对f2中的可见光波段进行光谱密度积分，计算出云亮度fbr，

其中，λi为波长，λmax、λmin分别为可见光区域内的最大波长和最小波长，单位均为nm，bv为可见光的波段范围，f(λi)表示波长为λi波段的像元值；

可见光光谱曲线呈现“白谱”表示光谱曲线平坦，则光谱曲线的一阶微分值较低，由于存在噪声和量化误差，因此微分值的精度会受到不同程度的影响。将每个像元与fbr的误差的微分作为白度指数。

s42：计算白度指数fwh，

其中，e(λi)＝|f(λi)-fbr|；

e(λi)表示波长为λi波段的像元值与fbr的误差。显然，云区越白，e(λi)越小，白度指数越低，反之白度指数越高。利用白度指数对同态滤波结果进行优化，能将像元值较高但非白色的物体剔除，降低错误识别率。

s43：将fwh＞β的像元对应的物体判断为非白色物体并将其从f2滤除，得到过滤后的图像f3。

本实施例中将β的值选定为0.1，在其他实施例中，可以对β的值进行调整以使其适应实际图像。

s5：选用圆形结构元素，对f3先进行闭运算再进行开运算，以对其进行形态学优化，得到云区识别结果f4。

本实施例中选用直径为7个像元的圆形结构元素。

形态学滤波基本思想可理解为利用结构元素对二值图像进行滤波。形态学滤波效果与结构元素选择有关，由于云区通常呈现不规则状，因此本发明选用圆形结构元素获得平滑的云区轮廓，并采用先进行闭运算再进行开运算的策略对云区掩膜结果进行后处理。先进行闭运算能填充云区内的空隙同时不显著改变其面积，再进行开运算能消除云区边界离散噪声使轮廓平滑。经形态学处理后，较破碎的小斑块云点被去除，一些错误识别的高亮离散点也被滤除掉，且云区周围的轮廓显示更加分明，得到较优的后处理效果。

图1a-图1d为选用圆形结构元素对云区进行形态学优化的具体实施例，图1a为7个像元的圆形结构元素，图1b为云区示例f3，其中白色为云像元，黑色为背景像元，图1c为f3经过闭运算处理后的结果，图1d为图1c再经过开运算处理的结果，即图像f4；

图2a-图2c为利用白度指数以及形态学运算对同态滤波结果f2进行优化的实例，图2a为原始输入图像f，图2b为同态滤波结果f2，图2c为经过后处理优化的图像f4；

图3a为8幅不同云量覆盖的高分辨率光学遥感图像(1024×1024)，8幅图像云量覆盖从左上至右下依次减少；

图3b为采用本发明的方法对图3a的云区检测结果，图3a中为图像f，图3b中为图像f4，图3a、图3b中位于同一位置的图像一一对应；

图4中第一列的前两幅图为gf-1卫星中的8m分辨率多光谱高分相机(pms)采集到的遥感图像f，后两幅图为16m分辨率多光谱宽幅相机(wfv)采集到的遥感图像f，对第一列原始影像采用本发明的方法处理后得到第二列的云区掩膜结果f4，第三列为从第一列的遥感影像中除去云区后的结果影像。

通过对79景影像的云区识别结果进行人工解译，精度验证结果显示，本发明精度较高，总体精度为93.81％。其中，云区掩膜内的精度最高，达97.69％，云区掩膜缓冲区精度较低，为90.15％。从精度误差分析来看，云区边界轮廓处产生较多误差，靠近云中心点以及缓冲区远离点精度均较高。从云区中心逐渐向云区边界过渡时，识别的正确率逐渐下降；从云区边界向缓冲区边缘过渡时，正确率又逐渐升高。由于云区边界多为不规则形状，当云区边界像素点与地物对比不明显时，基于同态滤波的算法不能有效的准确识别。其次，本发明中采用的先闭后开的形态学滤波处理也会使得云区边界产生变化。总体看来，本发明能够相对准确的识别gf-1遥感数据的云区位置，获得较为精准的云区掩膜结果。

本发明提供的对高分辨率光学遥感图像的自适应云区检测方法通过对大量高分一号(gf-1)遥感图像进行试验分析，得出了图像的径向能量谱与同态滤波过程中应用的高通滤波器的截止频率之间的量化关系，能够根据不同的输入图像自适应调整截止频率，实现gf-1遥感图像的批量云区识别。本发明能够对高分辨率光学遥感图像进行批量化高效的云区检测，总体识别精度达93.81％。

本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。

本领域普通技术人员可以理解：实施例中的装置中的模块可以按照实施例描述分布于实施例的装置中，也可以进行相应变化位于不同于本实施例的一个或多个装置中。上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。