CN107437237B - 一种区域无云影像合成方法 - Google Patents

Abstract

本发明提供了一种直接合成目标区域无云影像的方法，包括步骤：根据待合成影像范围得到最终合成影像覆盖范围，并创建一个同样大小的模板影像和合成多光谱影像，并将像元值初始化；依次对合成多光谱影像每个像素依据无云影像影优先、相似性高影像优先的规则进行填充，并将影像编号记录在模板影像。根据模板影像记录的信息对合成后的无云影像以羽化宽度矩形范围内各影像编号出现的次数为权重对每一个像素做羽化，最终直接得到具有辐射一致的区域无云影像。本发明适用于有云覆盖的数字正射影像直接合成目标区域的无云影像，解决厚云遮挡不能直接获取目标区域无云影像的问题。

Description

一种区域无云影像合成方法

技术领域

本发明属于摄影测量与遥感影像处理领域，特别涉及一种区域无云影像合成方法。

背景技术

遥感影像被广泛的应用于陆地观测领域，然而成像时厚云的存在导致覆盖区域信息缺失，影像质量退化，影响影像的解译精度，因此厚云去除是遥感影像处理及遥感制图的重要步骤。

用于视频图像生成全景图的流形方法可以尽可能的保证合成影像内部的辐射一致性，但是由于卫星影像不具有视频图像在时空上的连续性，该方法难以适应。因此本方法针对卫星数据改进了流形方法用于直接合成区域无云影像。同时，对于异构景观该方法在合成后影像内部存在明显的拼接缝，并且由于云层细碎且形状不规则导致图斑与图斑之间的拼接缝异常复杂，传统的拼接缝消除方法难以适应。因此本发明方法在改进流形方法合成卫星数据无云影像的基础上，提出一种基于面积的拼接缝消除方法，解决合成无云影像后留下的复杂拼接缝问题，最终得到具有良好视觉一致性的区域无云影像。

发明内容

为了获得良好视觉一致性的区域无云影像，解决合成无云影像后留下的复杂拼接缝问题，本发明所采用的技术方案是：一种区域无云影像合成方法，包括以下步骤：

步骤1，确定每幅待合成影像I_i及其对应的云检测结果影像I'_i，其中i＝1,2,...,n，所述云检测结果影像I'_i为二值影像；

步骤2，根据每幅待合成影像I_i的左上角和右下角的地理坐标，确定最终合成结果的覆盖范围，并创建一个同样大小的模板影像以及合成影像，同时将像元值初始化为0；

步骤3，根据合成影像左上角第一个像素的地理位置，在待合成影像中选出在该地理位置具有唯一灰度值覆盖的影像，以影像编号和灰度值分别填充到模板影像、合成影像对应位置；

步骤4，对于下一个待填充的位置点，根据地理位置在待合成影像中寻找对该位置有灰度值覆盖的影像序列，并根据对应的云检测结果影像将影像序列分为有云覆盖和无云覆盖两类；从无云覆盖的影像中挑选在待填充点处周边像素最相似的影像，将影像编号填充到模板影像，将灰度值填充到合成影像；

步骤5，重复步骤是4，直至合成影像上每一个像元都被填充，得到合成的目标区域的无云影像和记录了合成影像每个像元来源的模板影像；

步骤6，设定羽化宽度，根据模板影像记录的信息对合成后的无云影像，以羽化宽度矩形范围内各影像编号出现的次数为权重对每一个像素做羽化；

所述步骤6的实现方式如下，

设拼接缝羽化的宽度为w，从行列号为(w/2，w/2)开始到(Rows-w/2,Cols-w/2)结束，对每一个像素做羽化处理，其中Rows、Cols为模板影像或者合成影像的行数和列数；

步骤6.1，以待羽化点(i，j)为中心，以羽化宽度为边长确定矩形区域，统计出现的像元来源的影像编号i₁、i₂、…、i_q及其对应的像素个数n_i1、n_i2、….、n_iq；

步骤6.2，以步骤6.1中的影像编号确定影像序列，根据在待羽化点是否有数据覆盖以及是否为云层剔除部分影像，保留在该点无云覆盖的有效影像Ij₁、Ij₂、...、Ij_r，其中j₁、j₂、...、j_r为影像编号i₁、i₂、…、i_q的子集；

步骤6.3，根据有效影像各自在羽化矩形区域内出现的像素个数n_j1、n_j2、….、n_jr，设置权重加权得到羽化后的像素值，加权羽化计算方式为：

其中，H为创建的合成影像在(i，j)处羽化后的灰度值，I_j1为影像Ij₁在对应地理位置的灰度值；

步骤6.4，对合成影像每一个波段都重复步骤6.3的方式处理；

步骤7，重复步骤6，对每一个像素做拼接缝羽化，得到经过拼接缝消除的目标区域的无云影像。

进一步的，所述步骤2的实现方式如下，

步骤2.1，对影像I_i，i＝1,2,...,n的x坐标排序，得到x方向的最大值和最小值X_max、X_min，对I_i，i＝1,2,...,n的y坐标排序，得到y方向的最大值和最小值Y_max、Y_min，以此值组合为模板影像和合成影像矩形角点坐标(X_min，Y_min)，(X_max，Y_min)，(X_max，Y_max)，(X_min，X_max)；

步骤2.2，根据模板影像和合成影像矩形角点坐标值，所述模板影像和合成影像的地理范围以及分辨率完全一致，得到模板影像和合成影像的行数和列数，其计算公式为：

其中，Cols表示模板影像列数，X_max、X_min、Res_x分别表示模板影像中X方向坐标最大值、X方向坐标最小值、X方向的分辨率，其中Res_x的值与待合成影像I₁在X方向分辨率值相等；Rows表示模板影像行数，Y_max、Y_min、Res_y分别表示模板影像中Y方向坐标最大值、Y方向坐标最小值、Y方向的分辨率，其中Res_y的值与待合成影像I₁在Y方向分辨率值相等；

然后创建单波段模板影像，每个像素1个字节存储，以0为初始值，写入数据到文件；合成影像的波段数以及每个像素需要的字节数与输入的待合成影像I₁相同。

进一步的，所述步骤3中行列号转为地理位置的计算方式如下，

x＝col×Res_x+X_min (3)

y＝row×Res_y+Y_min (4)

其中，x、y分别对应转换后的X方向地理坐标值和Y方向的地理坐标值，col、row为模板影像中待转换像素点的列号和行号。

进一步的，所述步骤4的实现方式如下，

步骤4.1，假设现在模板影像和合成影像已经填充到行列号为(i，j)位置，对于待填充位置(i，j+1)，根据行列号以及模板影像左上角地理位置坐标和空间分辨率，根据公式(3)、(4)得到(i，j+1)处的地理坐标(x_j+1，y_i)；根据该地理坐标比对I₁、I₂、...、I_n中每一幅影像的矩形地理范围，如果该点在Ii确定的矩形地理范围内，将该影像编号i放入到数组C中，得到覆盖(x_j+1，y_i)点的影像序列Ii₁、Ii₂、...、Ii_m；

步骤4.2，对于Ii₁、Ii₂、...、Ii_m对应的云检测结果影像I'i₁、I'i₂、...、I'i_m，将地理坐标(x_j+1，y_i)转换到每一幅云检测结果影像对应的行列号，转换计算的公式如下：

其中，表示影像Ii_m左上角x轴方向地理坐标、y轴方向地理坐标；根据云检测结果影像I'i_m在的灰度值，将影像Ii₁、Ii₂、...、Ii_m分为云覆盖和无云覆盖两类；

步骤4.3，以模板影像(i，j+1)为中心，取(i，j)、(i-1，j)、(i-1，j+1)这三个像素的行列号坐标，根据公式(3)、(4)、(5)、(6)可以得到Ii₁、Ii₂、...、Ii_m中每幅影像对应像素的行列号坐标，取出坐标位置处的灰度值并组成向量，即每一幅影像取出三个像素位置的多个灰度值组成一个向量；依次对Ii₁、Ii₂、...、Ii_m中每幅影像灰度值组成的向量与合成影像在(i，j)、(i-1，j)、(i-1，j+1)位置处灰度值组成的向量计算相关系数，其中相关系数的计算公式为：

其中，表示创建的合成影像在(i，j)、(i-1，j)、(i-1，j+1)处对应的多个波段灰度值组成的向量；表示Ii₁、Ii₂、...、Ii_m中某一幅无云影像与对应地理位置的多个波段灰度值组成的向量；表示与之间的协方差，表示的方差；取相关系数最高的影像Ii_p在(x_j+1，y_i)处灰度值填充到合成影像(i，j+1)处，并将模板影像(i，j+1)的值填充为i_p，其中Ii_p∈{Ii₁、Ii₂、...、Ii_m}。

进一步的，所述步骤1所确定的云检测结果是通过ENVI软件以监督分类方式获得的。

与现有技术相比，本发明的有益效果：本发明改进了针对视频图像合成全景图的流形方法，将其用于卫星影像合成目标区域影像。并针对卫星影像之间存在辐射差异，导致合成后影像存在极其复杂拼接缝现象的问题，进一步通过统计羽化宽度矩形范围内像素各自所占有比例设置权重对合成影像羽化，即一种基于面积的拼接缝消除方法来解决传统拼接缝羽化方法不能适应拼接缝极其复杂不规则的问题，得到最终具有良好视觉一致性的无云影像。

附图说明

图1是本发明实施实例基于面积的拼接缝羽化消除过程示意图，I₁、I₂、I₃、I₄、I₅表示编号为1、2、3、4、5的影像，p点为待羽化的像素点，小矩形待统计出现编号的范围区域，大矩形为创建的模板影像。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本实例针对有云层覆盖的数字正射影像，提供了直接合成目标区域的无云影像，包括以下步骤：

步骤1：确定每幅待合成影像及其对应的云检测结果。为叙述方便，待合成的n幅影像分别用I₁、I₂、...、I_n表示，对应的云检测结果影像为I'₁、I'₂、...、I'_n。其中I'₁、I'₂、...、I'_n为二值影像，以“0”表示背景和云层，以“1”表示非背景、非云层的有效值。云检测可以通过商业软件ENVI通过选取样本监督分类得到。

步骤2：根据I_i，i＝1,2,...,n的左上角和右下角的地理坐标，得到最终合成结果的覆盖范围，并创建一个同样大小的模板影像以及合成影像，同时将像元值初始化为0。具体包括以下子步骤：

步骤2.1：对影像I_i，i＝1,2,...,n的x坐标排序，得到x方向的最大值和最小值X_max、X_min，对I_i，i＝1,2,...,n的y坐标排序，得到y方向的最大值和最小值Y_max、Y_min，以此值组合为模板影像和合成影像矩形角点坐标(X_min，Y_min)，(X_max，Y_min)，(X_max，Y_max)，(X_min，X_max)。

步骤2.2：根据模板影像和合成影像(两者地理范围以及分辨率完全一致)矩形角点坐标值，得到模板影像和合成影像的行数和列数，其计算公式为：

其中Cols表示模板影像列数，X_max、X_min、Res_x分别表示模板影像中X方向坐标最大值、X方向坐标最小值、X方向的分辨率(合成影像中相同)，其中Res_x的值与待合成影像I₁在X方向分辨率值相等。Rows表示模板影像行数，Y_max、Y_min、Res_y分别表示模板影像中Y方向坐标最大值、Y方向坐标最小值、Y方向的分辨率(合成影像中相同)，其中Res_y的值与待合成影像I₁在Y方向分辨率值相等。然后创建模板影像，单波段即可，每个像素1个字节存储，以0为初始值，写入数据到文件。合成影像的波段数以及每个像素需要的字节数与输入的待合成影像I₁相同。例如，若影像I₁有5个波段，每个波段以2个字节存储，则创建的合成影像有5个波段，每个像素2个字节存储。

步骤3：根据合成影像行列号为(0,0)(左上角第一个像素)的地理位置，在待合成影像I₁、I₂、...、I_n中，选出在该地理位置有灰度值覆盖的影像I_i，以数字编号i填充到模板影像(0,0)位置上，以I_i在该位置覆盖的多个波段灰度值填充到合成影像(0,0)位置。在该地理位置一般只有一幅影像的灰度值有覆盖。其中行列号转为地理位置的计算公式为：

x＝col×Res_x+X_min (3)

y＝row×Res_y+Y_min (4)

其中，x、y分别对应转换后的X方向地理坐标值和Y方向的地理坐标值，col、row为模板影像中待转换像素点的列号和行号。例如计算模板影像中行列号坐标为(0,0)的像素点地理坐标时，分别带入col＝0、row＝0得到x＝X_min、y＝Y_min。即行列号坐标为(0,0)的像素点对应的地理坐标为(X_min，Y_min)。X_min、Y_min为步骤2.1中得到的值。

步骤4：假设现在模板影像和合成影像已经填充到行列号为(i，j)位置。对于下一个待填充的位置(i，j+1)，由公式(3)、(4)计算，得到对应地理位置(x_j+1，y_i)。计算过程如下：代入col＝j+1、row＝i得到x＝(j+1)×Res_x+X_min、y＝i×Res_y+Y_min。那么像素点(i，j+1)的对应的地理坐标为((j+1)×Res_x+X_min，i×Res_y+Y_min)。为简化书写，记x_j+1＝x＝(j+1)×Res_x+X_min，y_i＝y＝i×Res_y+Y_min，则行列号坐标为(i，j+1)的像素点对应的地理坐标为(x_j+1，y_i)。对待合成影像I₁、I₂、...、I_n中每一幅影像，判断点(x_j+1，y_i)是否在影像的矩形区域范围内，得到覆盖该点的影像序列Ii₁、Ii₂、...、Ii_m，其中对应的云检测结果影像为I'i₁、I'i₂、...、I'i_m。根据对应的云检测结果影像可以将影像Ii₁、Ii₂、...、Ii_m分为两类，一类影像在(x_j+1，y_i)处有云覆盖，一类影像在该位置无云覆盖。从无云覆盖的影像中挑选在(x_j+1，y_i)处周边像素最相似的影像Ii_p，将i_p填充到模板影像，将影像Ii_p的多个波段灰度值填充到合成影像。如果不存在无云影像覆盖该位置，就以同样的方式在云层覆盖的影像中寻找。具体包括以下子步骤：

步骤4.1：对于待填充位置(i，j+1)，根据行列号以及模板影像左上角地理位置坐标和空间分辨率，根据公式(3)、(4)得到(i，j+1)处的地理坐标(x_j+1，y_i)。根据该地理坐标比对I₁、I₂、...、I_n中每一幅影像的矩形地理范围，如果该点在Ii确定的矩形地理范围内，将该影像编号i放入到数组C中，得到覆盖(x_j+1，y_i)点的影像序列Ii₁、Ii₂、...、Ii_m。

步骤4.2：对于Ii₁、Ii₂、...、Ii_m对应的云检测结果影像I'i₁、I'i₂、...、I'i_m，将地理坐标(x_j+1，y_i)转换到每一幅云检测结果影像对应的行列号，转换计算的公式如下：

其中，表示影像Ii_m左上角x轴方向地理坐标、y轴方向地理坐标。取云检测结果影像I'i_m在的灰度值，如果灰度值为“0”，则将影像Ii_m划分到有云覆盖一类的影像中去；如果灰度值为“1”，则将影像Ii_m划分到无云覆盖一类的影像中去。将Ii₁、Ii₂、...、Ii_m中的每一幅影像都做计算判断，把影像Ii₁、Ii₂、...、Ii_m分为云覆盖和无云覆盖两类。

步骤4.3：以模板影像(i，j+1)为中心，取(i，j)、(i-1，j)、(i-1，j+1)这三个像素的行列号坐标，根据公式(3)、(4)、(5)、(6)可以得到Ii₁、Ii₂、...、Ii_m中每幅影像对应像素的行列号坐标，取出坐标位置处的灰度值并组成向量，即每一幅影像取出三个像素位置的多个灰度值组成一个向量。依次对Ii₁、Ii₂、...、Ii_m中每幅影像灰度值组成的向量与合成影像在(i，j)、(i-1，j)、(i-1，j+1)位置处灰度值组成的向量计算相关系数，其中相关系数的计算公式为：

其中表示创建的合成影像在(i，j)、(i-1，j)、(i-1，j+1)处对应的多个波段灰度值组成的向量。组成向量的方法如下：若合成影像含有三个波段，其三个像素位置处多个波段灰度值组成的向量可以是其中表示合成影像在(i，j)位置处第一波段的灰度值；表示Ii₁、Ii₂、...、Ii_m中某一幅无云影像与对应地理位置的多个波段灰度值组成的向量，该向量的构成方式与相同。表示与之间的协方差，表示的方差。取相关系数最高的影像Ii_p(其中Ii_p∈{Ii₁、Ii₂、...、Ii_m})在(x_j+1，y_i)处灰度值填充到合成影像(i，j+1)处，并将模板影像(i，j+1)的值填充为i_p。如果不存在无云影像覆盖该位置就在云层覆盖的影像中以同样的方式寻找相关系数最大的云覆盖影像，将灰度值和影像编号分别填充到合成影像和模板影像。

步骤5：重复步骤是4，直至合成影像上每一个像元都被填充，得到合成的目标区域无云影像和记录了合成影像每个像元来源的模板影像。下面对合成的无云影像做拼接缝消除。

步骤6:设拼接缝羽化的宽度为w，则从行列号为(w/2，w/2)开始到(Rows-w/2,Cols-w/2)结束，对每一个像素做羽化处理。其中Rows、Cols为模板影像或者合成影像的行数和列数。具体包括以下子步骤：

步骤6.1：对待羽化的中心点(i，j)，其中w/2<＝i<＝Rows-w/2，w/2<＝j<＝Cols-w/2。根据模板影像统计从左上角(i-w/2,j-w/2)到(i+w/2,j+w/2)的矩形范围内，出现的像元来源影像的编号i₁、i₂、…、i_q及其对应的像素个数n_i1、n_i2、….、n_iq。其中n_i1+n_i2+….+n_iq＝w*w。如图1所示，在待羽化的像素点P有出现的来源影像有I₁、I₃、I₄、I₅，对应的编号为1、3、4、5。其中I₁在以p点为中心矩形框内像素个数为n₁，I₃、I₄、I₅各自的像素个数为n₃、n₄、n₅。

步骤6.2：对影像Ii₁、Ii₂、...、Ii_q序列。如果存在某幅影像没有覆盖(i，j)，即(i，j)的地理位置不在该影像的范围内，则将该影像从影像序列Ii₁、Ii₂、...、Ii_q中移除。同时对每一幅影像对应的云检测结果影像I'i₁、I'i₂、...、I'i_q判断在模板(i，j)处对应的地理位置是否为云层。如果是云层，则将将该影像从影像序列Ii₁、Ii₂、...、Ii_q中移除。移除后序列影像还剩下Ij₁、Ij₂、...、Ij_r，其中在图1示例中，影像I₃在p点处有云层覆盖，将I₃从影像序列I₁、I₃、I₄、I₅中剔除，还剩下影像I₁、I₄、I₅。

步骤6.3：根据剩下序列影像Ij₁、Ij₂、…、Ij_r及其在羽化宽度矩形范围内统计的对应像素个数n_j1、n_j2、….、n_jr，设置权重将剩下的序列影像在(i，j)处加权得到羽化后的像素值。加权羽化计算方式为：

其中H为创建的合成影像在(i，j)处羽化后的灰度值，I_j1为影像Ij₁在对应地理位置的灰度值。对于图一中p点单个波段羽化后的灰度值计算为：

步骤6.4:对合成影像在(i，j)处每一个波段都采用步骤6.3的计算方式得到合成影像在(i，j)处羽化后的多个波段灰度值。

步骤7：重复步骤6，对合成影像从行列号(w/2，w/2)到(Rows-w/2,Cols-w/2)中的每一个像素做拼接缝羽化，得到经过拼接缝消除的目标区域的无云影像。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (5)

1.一种区域无云影像合成方法，其特征在于，包括以下步骤：

步骤1，确定每幅待合成影像I_i及其对应的云检测结果影像I'_i，其中i＝1,2,...,n，所述云检测结果影像I'_i为二值影像；

步骤2，根据每幅待合成影像I_i的左上角和右下角的地理坐标，确定最终合成结果的覆盖范围，并创建一个同样大小的模板影像以及合成影像，同时将像元值初始化为0；

步骤3，根据合成影像左上角第一个像素的地理位置，在待合成影像中选出在该地理位置具有唯一灰度值覆盖的影像，以影像编号和灰度值分别填充到模板影像、合成影像对应位置；

步骤4，对于下一个待填充的位置点，根据地理位置在待合成影像中寻找对该位置有灰度值覆盖的影像序列，并根据对应的云检测结果影像将影像序列分为有云覆盖和无云覆盖两类；从无云覆盖的影像中挑选在待填充点处周边像素最相似的影像，将影像编号填充到模板影像，将灰度值填充到合成影像；

步骤5，重复步骤是4，直至合成影像上每一个像元都被填充，得到合成的目标区域的无云影像和记录了合成影像每个像元来源的模板影像；

步骤6，设定羽化宽度，根据模板影像记录的信息对合成后的无云影像，以羽化宽度矩形范围内各影像编号出现的次数为权重对每一个像素做羽化；

所述步骤6的实现方式如下，

设拼接缝羽化的宽度为w，从行列号为(w/2，w/2)开始到(Rows-w/2,Cols-w/2)结束，对每一个像素做羽化处理，其中Rows、Cols为模板影像或者合成影像的行数和列数；

步骤6.1，以待羽化点(i，j)为中心，以羽化宽度为边长确定矩形区域，统计出现的像元来源的影像编号i₁、i₂、…、i_q及其对应的像素个数n_i1、n_i2、….、n_iq；

步骤6.2，以步骤6.1中的影像编号确定影像序列，根据在待羽化点是否有数据覆盖以及是否为云层剔除部分影像，保留在该点无云覆盖的有效影像Ij₁、Ij₂、...、Ij_r，其中j₁、j₂、...、j_r为影像编号i₁、i₂、…、i_q的子集；

步骤6.3，根据有效影像各自在羽化矩形区域内出现的像素个数n_j1、n_j2、….、n_jr，设置权重加权得到羽化后的像素值，加权羽化计算方式为：

其中，H为创建的合成影像在(i，j)处羽化后的灰度值，I_j1为影像Ij₁在对应地理位置的灰度值；

步骤6.4，对合成影像每一个波段都重复步骤6.3的方式处理；

步骤7，重复步骤6，对每一个像素做拼接缝羽化，得到经过拼接缝消除的目标区域的无云影像。

2.根据权利要求1所述的无云影像合成方法，其特征在于：所述步骤2的实现方式如下，

步骤2.1，对影像I_i，i＝1,2,...,n的x坐标排序，得到x方向的最大值和最小值X_max、X_min，对I_i，i＝1,2,...,n的y坐标排序，得到y方向的最大值和最小值Y_max、Y_min，以此值组合为模板影像和合成影像矩形角点坐标(X_min，Y_min)，(X_max，Y_min)，(X_max，Y_max)，(X_min，X_max)；

步骤2.2，根据模板影像和合成影像矩形角点坐标值，所述模板影像和合成影像的地理范围以及分辨率完全一致，得到模板影像和合成影像的行数和列数，其计算公式为：

其中，Cols表示模板影像列数，X_max、X_min、Res_x分别表示模板影像中X方向坐标最大值、X方向坐标最小值、X方向的分辨率，其中Res_x的值与待合成影像I₁在X方向分辨率值相等；Rows表示模板影像行数，Y_max、Y_min、Res_y分别表示模板影像中Y方向坐标最大值、Y方向坐标最小值、Y方向的分辨率，其中Res_y的值与待合成影像I₁在Y方向分辨率值相等；

然后创建单波段模板影像，每个像素1个字节存储，以0为初始值，写入数据到文件；合成影像的波段数以及每个像素需要的字节数与输入的待合成影像I₁相同。

3.根据权利要求2所述的无云影像合成方法，其特征在于，所述步骤3中行列号转为地理位置的计算方式如下，

x＝col×Res_x+X_min (3)

y＝row×Res_y+Y_min (4)

其中，x、y分别对应转换后的X方向地理坐标值和Y方向的地理坐标值，col、row为模板影像中待转换像素点的列号和行号。

4.根据权利要求3所述的无云影像合成方法，其特征在于：所述步骤4的实现方式如下，

步骤4.1，假设现在模板影像和合成影像已经填充到行列号为(i，j)位置，对于待填充位置(i，j+1)，根据行列号以及模板影像左上角地理位置坐标和空间分辨率，根据公式(3)、(4)得到(i，j+1)处的地理坐标(x_j+1，y_i)；根据该地理坐标比对I₁、I₂、...、I_n中每一幅影像的矩形地理范围，如果该点在Ii确定的矩形地理范围内，将该影像编号i放入到数组C中，得到覆盖(x_j+1，y_i)点的影像序列Ii₁、Ii₂、...、Ii_m；

步骤4.2，对于Ii₁、Ii₂、...、Ii_m对应的云检测结果影像I'i₁、I'i₂、...、I'i_m，将地理坐标(x_j+1，y_i)转换到每一幅云检测结果影像对应的行列号，转换计算的公式如下：

其中，表示影像Ii_m左上角x轴方向地理坐标、y轴方向地理坐标；根据云检测结果影像I'i_m在的灰度值，将影像Ii₁、Ii₂、...、Ii_m分为云覆盖和无云覆盖两类；

步骤4.3，以模板影像(i，j+1)为中心，取(i，j)、(i-1，j)、(i-1，j+1)这三个像素的行列号坐标，根据公式(3)、(4)、(5)、(6)可以得到Ii₁、Ii₂、...、Ii_m中每幅影像对应像素的行列号坐标，取出坐标位置处的灰度值并组成向量，即每一幅影像取出三个像素位置的多个灰度值组成一个向量；依次对Ii₁、Ii₂、...、Ii_m中每幅影像灰度值组成的向量与合成影像在(i，j)、(i-1，j)、(i-1，j+1)位置处灰度值组成的向量计算相关系数，其中相关系数的计算公式为：

其中，表示创建的合成影像在(i，j)、(i-1，j)、(i-1，j+1)处对应的多个波段灰度值组成的向量；表示Ii₁、Ii₂、...、Ii_m中某一幅无云影像与对应地理位置的多个波段灰度值组成的向量；表示与之间的协方差，表示的方差；取相关系数最高的影像Ii_p在(x_j+1，y_i)处灰度值填充到合成影像(i，j+1)处，并将模板影像(i，j+1)的值填充为i_p，其中Ii_p∈{Ii₁、Ii₂、...、Ii_m}。

5.根据权利要求1所述的无云影像合成方法，其特征在于：所述步骤1所确定的云检测结果是通过ENVI软件以监督分类方式获得的。