近年来，遥感卫星技术不断发展，卫星数量不断增加，遥感数据的数量呈几何级数增长，已经从TB级数据规模向数十PB级规模增长。遥感数据资料的丰富推动了遥感应用的发展，目前针对大尺度区域的变化分析需求日益增加，如全国一张图建设^[1]、“一带一路”经济发展监测^[2]、全球森林土地覆盖变化^[3]等。针对时间序列的大尺度区域影像的快速生成需求，就是在时空条件约束下，从海量遥感影像资料中快速筛选出最佳影像组合，在完全覆盖目标区域的前提下，保证组合中影像的成像时间尽可能接近感兴趣时间，影像的云量尽可能低，为后续的遥感应用提供数据服务。遥感技术在我国环境监测、农业生产、矿产资源勘探、林业调查等许多行业都有广泛的应用，面对大量的遥感应用需求，从海量数据中快速、准确地筛选出所需数据至关重要^[4]。

国内外很多遥感平台提供了数据服务，例如国外的USGS(United States Geological Survey)^[5]和GLCF(The Global Land Cover Facility)^[6]等，国内的遥感数据共享平台有中国资源卫星应用中心数据服务平台^[7]、地理空间数据云^[8]、中国遥感卫星地面站提供的对地观测数据共享平台^[9]等。这些平台提供了数据检索功能，用户可以通过设置目标区域、成像时间范围和最高云量阈值等进行查询，但是若条件设置得过于宽泛，导致结果中的影像互相重叠，数据冗余度较高，需进一步人工挑选较优的数据；若条件设置得过于严格，则会导致结果中的影像无法完全覆盖目标区域。目前一般采用系统检索和人工挑选结合的方式筛选数据，一方面，人工挑选准确度难以保证，影像间可能会存在人眼识别不了的缝隙，另一方面，人工挑选耗时耗力，严重影响了用户使用遥感数据的效率。

国内外学者对遥感数据检索方法进行了大量研究。熊明豪^[10]、原发杰^[11]开展了遥感数据检索相关的研究，但是其中的遥感数据以瓦片的形式进行存储和检索。王金杰^[12]为帮助用户更好地获取遥感数据信息，进行基于语义的遥感影像数据检索关键技术研究，提高了检索的效率和质量，但没有涉及区域覆盖的遥感数据检索问题。苗立新等^[13]针对各卫星数据查询系统互相独立的问题，实现了多源遥感数据覆盖范围快速查询技术，但是仅做了数据查询条件的初筛，没有做覆盖整个区域的最少影像的筛选。Demir和Bruzzone^[14]以及Sukhia等^[15]开展了基于内容的遥感影像检索方法研究，Bretschneider等^[16]提出基于图像光谱信息的检索查询技术，Chaudhuri等^[17]匹配策略进行遥感图像的区域检索，但是这些研究都没有涉及区域覆盖的遥感影像筛选问题。

另一些学者对遥感数据的进一步筛选进行了探讨。李峰等^[18]提出一种针对遥感影像的区域覆盖最优数据集的筛选模型，对目标区域网格拆分和筛选参数归一化处理，保留覆盖每个网格权重最大的影像，并对结果进行去重。保留覆盖每个网格权重最大的影像，难以保证筛选结果影像数量最少，且网格划分的大小影响算法性能。网格划分的边界和目标区域存在一定的误差，如果网格划分过大，会产生在目标区域边界上数据多选或者漏选的问题，如果网格划分过小，计算量呈指数上升，导致算法效率大大降低。杨曦和颜静^[19]对资源3号卫星影像最优区域覆盖进行了算法研究，利用遥感影像的内接矩形对目标区域进行连续切割得到影像集合，最后进行去重。但是在影像筛选过程中，没有考虑属性条件的对比，且采用影像的内接矩形作为有效范围，对影像的真实空间覆盖范围损失较大。何方圆等^[20]采用基于空间二次过滤的遥感数据单时相全覆盖检索算法，对候选数据集进行权重组合排序，依次将权重最大的影像加入结果集，没有考虑每次新加入的影像与未被覆盖目标区域相交面积的大小，导致筛选结果可能会出现多张权重大的影像堆叠在一起的情况，筛选结果冗余度较高。

对于目前遥感数据筛选效率低、准确度低等问题，本文提出一种遥感影像区域覆盖数据集筛选方法，将遥感影像与目标区域相交的部分作为遥感影像的有效范围，利用遥感影像的有效范围进行相互切割，得到覆盖目标区域的互不重叠的碎块，对每景影像包含的碎块数量、成像时间距离感兴趣时间的时间间隔、云量等参数归一化得到一个综合代价，通过代价的大小进行选择，直至已选择影像对目标区域的覆盖率达到最大，得到一组最优影像组合。采用Landsat8数据，选择青藏高原地区作为目标区域进行实验，实验结果表明，该算法能自动地筛选出接近感兴趣时间、数量少且云量优的遥感影像组合，为进一步提高算法效率，对算法耗时的部分进行改进，采用多个进程并发执行。

1 原理与方法 1.1 方法原理

区域覆盖的数据筛选问题可以看作一个集合覆盖问题(set covering problem，SCP)，SCP是运筹学中经典的组合优化问题。集合覆盖问题可以描述为：给定2个集合E和S，分别为元素的集合E和E的子集的集合S，求出S的子集C，使得C中所有集合的并等于E，同时使得|C|最小。

对于集合覆盖问题，算法思想就是按照某种规则选择若干列覆盖所有的行，即每行至少被一列覆盖，并且使得花费的代价最小^[21]。例如图 1中A为一个6×8矩阵，向量C为矩阵A的列向量代价值。向量x为矩阵A的一个集合覆盖，而本例的最小覆盖为x=(0，0，0，1，0，1，0，1)，代价为15。

影像的成像时间越接近用户的感兴趣时间，越能够直接地反映问题。云一直是影响遥感成像的关键因素，它严重影响了信号的传输，被云层遮挡的区域难以识别地物信息，因此云量是衡量影像质量的一个关键要素。本文综合考虑了成像时间和云量2个属性要素。

遥感影像区域覆盖数据集筛选方法设计思路如下：首先，通过空间查询获取与目标区域相交的所有影像数据集，并进行简单的条件筛选，作为初始数据集。然后，利用目标区域对初始数据集进行裁剪，得到这些影像在目标区域的有效范围。利用各景影像的有效范围进行互相切割，得到互不重叠的碎块，相当于整个目标范围被分割为碎块。每景影像包含的碎块数量相当于图中每列数值为1的数量，数值为1代表影像包含对应行的碎块。通过对影像包含的碎块数量、成像时间与感兴趣时间的时间差以及云量归一化计算得到一个综合代价，通过代价的大小进行选择。以用户感兴趣时间为中心，影像成像时间越接近中心时间，相对代价越小；影像包含的碎块数量越多，相对代价越小；云量越低，相对代价越小。影像选择的规则是，每个碎块至少被一景影像包含，对应上面的每行至少被一列覆盖，采用贪婪算法的思想，每次优先选择代价最小的影像，为了避免结果有冗余，每次选择后需重新计算影像包含的碎块数量。

1.2 方法流程

遥感影像区域覆盖数据集筛选方法流程如图 2所示。

首先对遥感影像数据进行预筛选。通过设定目标区域范围、成像时间和云量阈值等条件，将满足上述条件的遥感影像作为候选数据集A。

然后利用目标区域对候选数据集A中的影像进行裁切，获得各景影像在目标区域的有效范围，其中有效范围指的是影像和目标区域相交的部分。将裁剪后的各遥感影像和其他裁剪后的遥感影像进行相互切割处理，以将各遥感影像与目标区域的重叠部分拆分成多个碎块，将相互切割处理结果中的碎块进行去重后生成碎片化几何要素集合，作为候选集合M。对裁切结果和碎片化几何要素集合进行空间连接，得到每景影像包含的碎块数量等信息。

为保证整个区域被完全覆盖，优先选择目标区域的一部分仅被该影像覆盖，而不被其他影像覆盖的影像。即当部分区域数据量较少时，仅有一景影像覆盖，那么该景影像被选择，具体步骤如下：

1) 首先对每景影像求包含的碎块的集合a，然后求除这景影像外，剩余影像包含的碎块的并集b。

2) 若集合a不被集合b包含，那么集合a对应的影像被优先选择，得到优先选择的影像集合B。

3) 选择结束后，从候选集合M中剔除已被选择影像所包含的碎块，以免对后续的选择产生干扰。

对剩余的影像，计算每景影像所包含的候选集合M中的碎块的数量，记为count；以及影像的成像时间距离需求时间的时间差，以天数为单位，记为costDate，其中需求时间指的是想要获取的遥感影像的最佳时间；影像的云量参数记为cloud。然后对这3个参数进行最值归一化处理，由于影像包含的碎块数量和代价呈反比，因此对碎块数量求倒数之后再归一化到(0，100)区间范围，记为count′；将costDate的值归一化到(0，100)区间范围，记为costDate′；云量本身就在(0，100)区间，记为cloud′。根据遥感应用场景，设置归一化后的影像包含碎块数量、成像时间差以及云量的权重系数分别为w₁、w₂、w₃(0≤w₁, w₂, w₃≤1 and w₁+w₂+w₃=1)，则遥感影像的综合代价为

$ \text { Cost }=w_1 \times \text { count }^{\prime}+w_2 \times \text { costDate }{ }^{\prime}+w_3 \times \text { cloud }^{\prime} . $

(1)

遍历剩余影像，计算每景影像的综合代价，采用贪婪的思想，每一步选择综合代价最小的影像，加入到影像集合B，并从候选集合M中剔除已被选择影像所包含的碎块，然后在剩下的影像中再次选取Cost最小值对应的影像，迭代循环。

当已选择的影像能够将所有碎块完全包含，即影像能够将目标区域完全覆盖，选择结束，所得到的影像集合B即为最终结果。当目标区域的全部数据无法对目标区域完全覆盖时，那么对应地筛选结果实现的是对目标区域的最大化覆盖。其中，每次选择影像后，将该影像所包含的碎块从候选集合中剔除，是为保证每一次选择时，都是选择包含剩下集合中碎块数量、成像时间和云量综合最优的影像，即对未被覆盖的目标区域贡献率最大的影像。

1.3 Dask并行计算

在上述方法流程中，各遥感影像的有效范围相互切割的步骤可以使用Dask进行并行化改进。Dask是一款用于分析计算的并行计算库，它能够在单机和集群中进行分布式计算，以一种更方便简洁的方式处理大数据量，与Spark这些大数据处理框架相比较，Dask更加轻量化。Dask更侧重与其他框架，如Numpy、Pandas、Scikit-learning相结合，从而使其能更加方便地进行分布式计算。

Dask中存在3种最基本的数据结构，分别是Arrays、Dataframes以及Bags，允许使用多核处理器并行运行它们的操作。它们可以存储大于RAM的数据，这些集合类型中的每一个都能够使用在RAM和硬盘之间分区的数据，以及分布在集群中多个节点上的数据。其中dask Dataframe是基于pandas Dataframe改进的一个可以并行处理大数据量的数据结构，即使对于大于内存的数据也能够处理。Dask DataFrame在单台机器上可进行大于内存的计算，也可运行在集群中，将空间限制由内存扩展到硬盘。

在Dask分布式中，共有3种角色：Client端，Scheduler端以及Worker端。Client负责提交task给Scheduler，Scheduler负责将提交的task按照一定的策略分发给Worker，Worker进行实际的计算、数据存储。在此期间，Scheduler时刻关注着Worker的状态。

当问题不适用于dask.array、dask.dataframe以及dask.bag集合时，可以使用更简单的dask.delayed和future并行化自定义算法。其中Delayed为延时计算，只是构建了一个Task Graph，并没有进行实际的计算，只有调用compute时，才开始真正计算；Future是立即执行的，可以通过submit、map方法将一个Function提交给Scheduler，Scheduler会在后台对提交的任务进行处理并分发给worker进行实际的计算，当任务提交后，会返回一个指向任务运行结果的Key值，即为Future对象。我们可以跟踪其状态，也可以通过result和gather方法等待任务完成后将结果收集到本地。Dask的工作原理如图 3所示。

2 实验及分析 2.1 实验数据

由于Landsat8卫星的成像质量稳定，且云量评估误差较小，选择Landsat8卫星的遥感影像数据作为实验数据，选择青藏高原地区作为目标区域，成像时间跨度为2018年1月1日到2020年12月31日。

Landsat8卫星上携带陆地卫星成像仪(operational land imager，OLI)和热红外传感器(thermal infrared sensor，TIRS)，陆地成像仪OLI包括9个波段，空间分辨率为3 m，其中包括一个15 m的全色波段，成像宽幅为185 km×185 km。热红外传感器TIRS包括2个单独的热红外波段，分辨率为100 m。卫星每16 d可以实现一次全球覆盖。

青藏高原是中国最大、世界海拔最高的高原，南起喜马拉雅山脉南缘，北至帕米尔高原北缘、西昆仑山和祁连山北缘，西自兴都库什山脉和帕米尔高原北缘，东及东北部与秦岭山脉西段和黄土高原相接^[22]，平均海拔约4 32 m，总面积为308.34万km²。

2.2 实验环境

实验使用处理器Intel(R) Core(TM) i7-7700 CPU@3.60 GHz，配备16 GB内存、8个逻辑处理器，实验环境为Windows 10专业版操作系统，编程语言为python3.7，主要函数库为GeoPandas0.9.0、Dask等。

2.3 数据处理

为便于进行数据处理，采用矢量文件表达遥感影像的边界范围，相关的元数据信息存储在属性表里面。对2018、2019、2020年这3年的数据进行筛选，以每年的7月15日为中心时间，验证本方法的有效性。首先根据成像时间、云量条件过滤数据，其中云严重影响信号的传输，导致成像质量变差，云量高的数据对于一般的遥感应用用处不大，因此过滤掉云量高于20%的数据。

在各影像的有效范围进行相互切割时，由于GeoPandas是一个独立的开源函数库，相对于ArcGIS商业软件，成本低，且容易和dask并行计算库结合，选择采用GeoPandas库调用overlay函数中的union模块。该函数的参数必须为2个GeoDataFrame类型的数据，为实现每景影像都被其他影像切割的目的，采用每2个要素先进行拆分，再将结果两两结合进行拆分的方式处理，如图 4(a)所示，图中2个要素作为输入，得到3个要素，分别为要素1和要素2单独的部分和2个要素的公共部分。将这3个要素作为其中一个输入和要素3进行union，得到图 4(b)的结果，每个要素都和要素3分别进行了拆分。

由于影像之间的边界距离过近，GeoPandas在进行空间计算时保留精度过高，出现了非节点交叉错误。因此，将每2个要素进行union结果中的要素几何属性坐标精度降低，然后作为下一次union操作的输入，以此方法避免因产生非节点交叉错误而执行失败。具体流程如图 5所示。

由于影像之间进行相互切割是算法最为耗时的部分，为提升方法效率，使用dask中的future计算模块对该部分进行并行化改进。每一层级中，每2个要素集合进行union操作是相互独立的，因此采用多个worker，即多个进程并发执行各个子任务。

2.4 实验结果分析

根据不同的应用场景，可以自定义影像包含的碎块数量、成像时间差和云量参数的权重占比，其中影像包含的碎块数量参数决定了筛选结果中影像的数量。本实验采用的综合代价计算方式为

$ \text { Cost }=\frac{1}{2} \text { count }^{\prime}+\frac{1}{4} \operatorname{costDate}^{\prime}+\frac{1}{4} \text { cloud }^{\prime} . $

(2)

对2018—2020年这3年覆盖青藏高原地区的数据筛选结果如图 6所示。

将本方法和何方圆等^[20]提出的方法(方法2)进行对比，为保证对比的一致性，方法2采用“时间优先”、“云量低于20%”作为主要权重因子，同样对2018至2020年这3年覆盖青藏高原地区的数据进行筛选，两种方法筛选结果中的影像数量、平均云量和平均成像时间差对比如表 1所示。

本方法会优先选择云量低、最接近用户感兴趣时间的影像，当还没有达到对目标区域完全覆盖时，会不断做妥协，选择质量较差的影像。在对2018—2020年连续3年的筛选结果中，平均云量在8.3%左右，平均成像时间差在50 d左右。从图 6中可以看出在青藏高原的西部边缘、中部和东部地区，相邻影像之间的地物光谱信息差异较小，其他地区地物光谱信息差异较大，主要原因在于其他地区时间在7月15日前后云量低的数据较少，只能退而求其次选择成像时间较远的数据。从表 1可以看出，2020年青藏高原地区云量优的数据较充足，而图 6中2020年的筛选结果相比其他2年光谱信息差异最小，说明目标区域的数据越充足，筛选结果中相邻影像之间的差异越小。

由于选择时没有考虑每景影像对未覆盖区域的覆盖面积，方法2对2018—2020年这3年的筛选结果中的影像数量均比本方法筛选出的数量多，且平均云量均高于本方法。方法2虽然优先选择成像时间接近的影像，但是结果中的成像时间差均高于本方法。且该方法仅滤除了云量20%以上的数据，在选择时没有做到尽可能选择云量低的数据，而本方法综合考虑了影像对未覆盖区域的空间贡献度、成像时间和云量，且根据用户对不同条件的关注程度，可以自由设置参数的权重占比。

采用dask并行计算库，对2018年数据影像相互切割的处理步骤使用多个worker并行计算，处理时间随worker数量增加的变化如图 7，其中原始指的是一个worker，即没有并行。从图中可以得出，本方法随着worker数量的增加，数据处理效率得到了很大提升。采用2个worker相对于原始处理时间缩短约38.6%，4个worker相对于原始处理时间缩短58.2%，速度提升效果最明显，worker数量为6或8时，速度提升幅度变小，数据处理时间逐渐趋于稳定，原因是worker数量的增加带来了多个worker之间数据通信成本增加的问题，因此采用4个worker进行并行处理是最佳的方案。

采用4个worker对数据处理部分进行并行改进，对3年的数据union并行处理前和union并行处理后的数据筛选耗时进行对比。如图 8所示，2018年数据筛选耗时缩短55.5%，2019年数据筛选耗时缩短56.1%，2020年数据筛选耗时缩短55.1%，效率提升明显。

采用单机对青藏高原地区的一年的Landsat8数据进行数据筛选，耗时大约在6 min，人工挑选这些数据要消耗1~2 d的时间，本方法效率相对于人工提升了几十倍甚至上百倍，若目标区域进一步扩大，人工挑选难度更大，且均为重复性的工作，而本方法可以在集群环境下运行，适用于处理大区域的数据筛选工作。

3 结语

本文提出一种遥感影像区域覆盖数据集筛选方法，该方法通过各影像在目标区域的部分进行相互切割，将目标区域拆分为互不重叠的碎块，通过影像包含的碎块数量、成像时间距离感兴趣时间的时间差以及云量参数选择数据，根据不同的遥感应用场景，用户可以自定义每个参数的权重系数，实现需求定制化。与人工选取的方法相比，该方法可以快速、自动地筛选出一组最优影像组合，尽可能保证组合中影像的数量少、成像时间接近感兴趣时间及云量优。采用尽可能少的质量优的影像可以减少后续的镶嵌工作，更好地为各类遥感应用提供数据服务。且人工选取容易出现数据多选或漏选，本方法采用迭代循环的方式选取影像，当达到对目标区域的完整覆盖或最大化覆盖时终止循环，不会出现多选或漏选的问题。相对于格网切分的方式，这种方法是更为精确的空间计算，避免了网格在目标区域边界上难以完全拟合的问题。

通过遥感影像区域覆盖数据集筛选方法进行并行化改进，大幅度提升了运行效率，证明该方法适用于解决大区域的数据筛选问题。当目标区域进一步扩大或候选遥感数据量变多时，计算量也随之增加，可以考虑在多节点的分布式集群环境下进行计算。