Large-Area Full-Coverage Remote Sensing Image Collection Filtering Algorithm for Individual Demands

: Remote sensing is the main technical means for urban researchers and planners to effectively observe targeted urban areas. Generally, it is difﬁcult for only one image to cover a whole urban area and one image cannot support the demands of urban planning tasks for spatial statistical analysis of a whole city. Therefore, people often artiﬁcially ﬁnd multiple images with complementary regions in an urban area on the premise of meeting the basic requirements for resolution, cloudiness, and timeliness. However, with the rapid increase of remote sensing satellites and data in recent years, time-consuming and low performance manual ﬁlter results have become more and more unac-ceptable. Therefore, the issue of efﬁciently and automatically selecting an optimal image collection from massive image data to meet individual demands of whole urban observation has become an urgent problem. To solve this problem, this paper proposes a large-area full-coverage remote sensing image collection ﬁltering algorithm for individual demands (LFCF-ID). This algorithm achieves a new image ﬁltering mode and solves the difﬁcult problem of selecting a full-coverage remote sensing image collection from a vast amount of data. Additionally, this is the ﬁrst study to achieve full-coverage image ﬁltering that considers user preferences concerning spatial resolution, timeliness, and cloud percentage. The algorithm ﬁrst quantitatively models demand indicators, such as cloudiness, timeliness, resolution, and coverage, and then coarsely ﬁlters the image collection according to the ranking of model scores to meet the different needs of different users for images. Then, relying on map gridding, the image collection is genetically optimized for individuals using a genetic algorithm (GA), which can quickly remove redundant images from the image collection to produce the ﬁnal ﬁltering result according to the ﬁtness score. The proposed method is compared with manual ﬁltering and greedy retrieval to verify its computing speed and ﬁltering effect. The experiments show that the proposed method has great speed advantages over traditional methods and exceeds the results of manual ﬁltering in terms of ﬁltering effect.


Introduction
To meet the needs of urban planning, land surveying, and other applications for large-scale regional observation, remote sensing with high spatial, temporal, and spectral resolutions is rapidly developing [1] and several countries have established relatively adequate satellites for ground observation; these satellites include the Landsat series of satellites and the moderate-resolution imaging spectroradiometer (MODIS) sensors on

Related Work
Current research on image filtering focuses mainly on retrieving image slices from massive images that are similar in land cover to input images; this research includes studies by Akshara [10], Liu [11], and Li [12]. The fundamental purpose of the abovementioned remote sensing image retrieval work was to analyze and retrieve images containing similar scenes. However, the above methods cannot solve the problem presented by users that filter images according to coverage, resolution, timeliness, etc.
Other scholars have conducted some research on remote sensing data filtering based on image attributes. Li proposed an image pyramid-oriented spatial indexing algorithm based on a linear quadtree; this algorithm improved the efficiency of coding and retrieving single-view images [13]. Xie designed a remote sensing image catalog data description framework and, based on this framework, designed an efficient catalog data storage organization scheme and positioning retrieval, qualitative retrieval, and combined retrieval algorithms to solve the problem of achieving fast and accurate image localization under distributed storage [14]. However, the studies by these scholars focused on parallelism and indexing mainly at the data level. Area coverage and filtering of image collections to meet individual demands still require manual intervention. There is still no solution for meeting these key issues of full coverage and the satisfaction of demands. Egenhofer studied the sketch spatial data retrieval method [15]; Lee proposed a visual query based on topological relationships in GIS [16]; and Shekhar, Volker, Gaede, and others studied information retrieval related to spatial information [17,18].
To date, few studies have focused on satisfying both full coverage and individual demands in remote sensing data filtering. He proposed a single-phase full-coverage filtering algorithm for remote sensing images [19]. Although this algorithm could complete full coverage of a given area, the algorithm did not consider factors such as cloudiness and resolution. Thus, it was difficult to obtain satisfactory results using this algorithm. Zuo started from remote sensing tiled data and conducted data filtering in units of tiles. A full-coverage retrieval model was designed, but this model still required manual review and interactive filtering of the results, which could not avoid the problem of low efficiency when the search area or the amount of data was large [5]. Li proposed a data set filtering model for optimal remote sensing image area coverage [20]. This model normalized the user-defined cloud range, time range, sensor type, and resolution. The score of each image on a regular grid was calculated, and the image with the highest score on each grid was selected as the filtering result. The filtering results of this model did not consider the image overlap ratio, so the filtered image collections were very repetitive.
Currently, the public can obtain remote sensing images for free from some websites (e.g., USGS Earth Explorer, the European Space Agency (ESA)'s sentinel mission, and the National Aeronautics and Space Administration (NASA)'s Reverb); these websites can also simply filter images according to cloud cover, latitude and longitude, and other conditions. However, unlike the proposed method, inputting cloud cover and latitude and longitude in the website for retrieval will enable hundreds, or even more, qualified images to be obtained, and images will greatly overlap. Users who want to achieve full coverage of very large areas by selecting multiple remote sensing images still cannot avoid the problem of selecting spatially complementary images from many images that the website filters. Even if several complementary images are selected to fully cover the area, there is no guarantee that, compared with other highly overlapping images, the image at this time is of relatively good quality.
In summary, there is currently no mature image collection filtering method that can both achieve full-coverage and satisfy the different demands of users.

Proposed Method
The workflow of the LFCF-ID proposed in this study is shown in Figure 1. This framework includes four main steps: obtaining the image collection to be filtered; fast, coarse filtering of the image collection; executing a genetic filtering algorithm; and result optimization. Generally, attributes such as cloudiness, location, and time vary greatly among remote sensing images. Most of the images that are available cannot satisfy users' specific demands. Thus, in the first step, it is necessary to obtain an image collection that is worth filtering by removing many redundant images according to the user's specific preset area, time interval, and minimum thresholds for cloudiness and resolution. The output is used as the image collection to be filtered to ensure that any image in the collection can meet the basic needs of users. Then, fast, coarse filtering of the image collection is performed. In this step, a score is designed, mainly in grid units, that can represent how well the image collections satisfy the user's demands. Images with the top k scores in each grid are selected as the result of coarse filtering. Next, the GA is used to further filter the coarse filtering results according to the fitness score, which can provide a comprehensive expression of the user's demand satisfaction and coverage, and the optimal full-coverage image collection result filtered to meet the individual demands is obtained. Finally, the final image collection result is obtained by result optimization, which is conducted to slightly change the optimal full-coverage image collection result by adding or subtracting a few images.
The blue area in Figure 1 is the area preset by the user, and the yellow rectangles are the boundaries of each remote sensing image. The figure shows that with the implementation of each step, the filtered image collection becomes increasingly optimal. Considering that the image collection filtering algorithm needs to meet the individual demands of different users, mathematical modeling of the requirements is essential. The user's demands are reflected mainly in terms of cloudiness, resolution, and timeliness, and this paper uses a weighted triplet to model the user demands.
W gsd , W cloud , and W time are the input weights preset by the user. W gsd indicates the user's demand weight for the spatial resolution of the image. A larger value of W gsd indicates that the user prefers images with a higher resolution in the filtering result. W cloud indicates the user's demand for image cloudiness. A larger value of W cloud indicates that the user prefers images with fewer clouds in the filtering result. The weight of W time indicates the user's demand for the timeliness of the image. The larger the value, the more the user hopes to obtain images in the filtering results that are close to the deadline selected by the user themselves.

Obtaining the Image Collection to Be Filtered
In this method, we first need to predefine several required parameters in the filtering workflow, as shown in Table 1: Table 1. Definition of the preset parameters.

Parameter Definition
Target_Area Specific area preset by users. The goal of the LFCF-ID is to filter and obtain the full-coverage image collection of this area T start ,T end The time interval is preset by users to obtain the images to be filtered GSD min The lowest image resolution that the user can accept in the filter result Cloud max The highest image cloudiness that the user can accept in the filtered result Coverage min The lowest coverage that the user can accept W gsd The user's demand weight for the spatial resolution of the image W cloud The user's demand weight for the image cloudiness W time The user's demand weight for the timeliness of the image In this step, all the images available are used as the input data, which can be thousands of scenes. We obtain the images that can satisfy the preset parameters Target_Area, T start , T end , GSD min , and Cloud max as the initial image collection, which will be filtered in the next steps. The initial image collection to be filtered is denoted by I_ori = {I 1 , I 2 , . . .}. I_ori is used as the input of the fast, coarse filtering step, which performs a further filtering operation according to W gsd , W cloud , and W time .

Fast, Coarse Filtering of the Image Collection
In this step, we conduct a coarse filtering of I_ori by designing a series of standard scores consisting of W gsd , W cloud , and W time . A subset of I_ori = {I 1 , I 2 , . . .} can be obtained as the output, which will be used as the input of the GA in the next step.

Creation of Global Grids
In this paper, coarse filtering and genetic filtering algorithms need to calculate several scores to rank images by using the grid as the computing unit. Large grids may cause fineness errors due to excessively large benchmarks, and small grids will be too fragmented and thus can result in many calculations, so the grid size should be as reasonable as possible. In the experiment, we divide the grid according to different standards such as 0.01 • × 0.01 • , 0.05 • × 0.05 • , 0.1 • × 0.1 • , 0.3 • × 0.3 • , and 0.5 • × 0.5 • . We find that among these standards, 0.1 • × 0.1 • works best. Therefore, we use 0.1 • × 0.1 • in all experiments.

Coarse Filtering Based on Grids
First, we calculate the grids denoted by g = {g 1 , g 2 , . . . , g i , . . . , g n }(i = 1, 2, . . . , n), which can cover the Target_Area preset by the user. Then, we search for the images denoted by I_ori i = {I_ori ∩ g i }, which cover g i in I_ori. Because I_ori i contains many images, only a few images belong to the optimal image collection that meets the user's demands. It is inappropriate to use all images in I_ori as the input of the GA, as doing so can cause data redundancy. Data redundancy may make it difficult for the GA to converge in the next step and also reduces the calculation efficiency of the GA. Thus, to reduce the number of images to be filtered by the GA, an image scoring formula S I is designed to calculate the satisfaction of the user's demands by each image in this paper. S I is calculated as where S gsd is the resolution score of the remote sensing image, S cloud is the cloudiness score of the remote sensing image, and S time is the timeliness score of the remote sensing image. S gsd is calculated as follows: where gsd max is the highest resolution of all images in I_ori, gsd min is the lowest resolution of all images in I_ori, and gsd present is the resolution of the image to be calculated. S cloud is calculated as follows: where C max represents the theoretical maximum value of the image cloudiness (this value is actually 100%), C min represents the theoretical minimum value of the image cloudiness (this value is actually 0%), and C present represents the cloudiness of the image to be calculated. S time is calculated as follows: where T start and T end are defined in Table 1 and T present represents the image shooting, which is calculated in days.
The formula for S I can be used to calculate the score of each image under the user's demands, which are modeled by Model user = W gsd , W cloud , W time . The best k images are selected as the coarse filtering result of each grid g i by ranking the S I values of the images in I_ori i ; the best k images of all grids are then merged together as the final coarse filtering result of the Target_Area; this result is denoted by where n is the number of grids that can fully cover the Target_Area. The coarse filtering method above can quickly remove many images, thereby producing a result that can satisfy the preset parameters Target_Area, T start , T end , GSD min , and Cloud max , but is not suitable enough for the final result of the LFCF-ID. After coarse filtering, the number of images in the image collection is reduced to the same magnitude as the number of grids n, thus possibly greatly improving the efficiency of the GA algorithm.

Further Filtering by the Genetic Algorithm
The ultimate purpose of the filtering task in this paper is to obtain an image collection containing the fewest images that can exactly fully cover the Target_Area specified by the user while best satisfying the user's demands. Since each image usually covers more than one grid, we select the k best images in each grid as the output of the coarse filtering; this selection will result in a high level of overlap between images. Many images are unnecessary and the image collection is still redundant. Therefore, it is time to find an optimization method to select the best combination of images as the optimal image collection that can meet the full-coverage requirement and the lowest repeated coverage requirement.
As a well-known search and optimization method, the GA has been successfully applied to intelligent optimization problems in various remote sensing applications, such as image classification [21][22][23], image segmentation [24], feature extraction [25][26][27], and quantitative inversion [28,29]. The quasi code for a GA is Algorithm 1.

Algorithm 1 Genetic Algorithm.
Input: Pc: possibility of cross Pm: possibility of mutation m: the number of genes in a population Output: optimal population 1: initialize population including m genes 2: calculate the fitness score for each gene 3: repeat 4: select m genes from population by using the roulette method 5: if(random(0,1) < Pc) { select two genes randomly cross between the two selected genes } 6: if(random(0,1) < Pm) { select one gene randomly mutation for this selected gene } 7: calculate fitness score for each gene 8: until(reaches stop condition) Inspired by the GA, this study uses a grid as the calculation unit to calculate the fitness score, which can model the image collection's state of coverage and satisfaction of personalized demands. The GA is used to continuously optimize the fitness score to obtain the optimal filtering results. In this algorithm, the input form and fitness functions should be defined according to the problem. The unit of the input is called the chromosome in the GA; this input can be a solution to the problem. In this paper, the chromosome refers to an image collection. The chromosome is composed of many genes that can be represented in binary, and refers to the images in the image collection. At the beginning of the GA, a population is formed by randomly generating multiple possible chromosomes. Biological selection, crossover, and mutation operations of the population are used to simulate the biological genetic evolution and achieve the next-generation population according to a preset fitness function. The higher the chromosome's fitness score, which is calculated by the fitness function, the higher the probability of saving the chromosome for the next generation. Therefore, the fitness score of the chromosome can improve generation-bygeneration. Eventually, the chromosome with the best fitness score is obtained as the final optimization result.

Population Initialization
First, we need to define the meaning of a chromosome in the GA according to the problem to be solved in this paper and generate an initial population. Since the final result of the LFCF-ID is an image collection, we define the chromosome as an image collection, which is denoted by G j .
The variable m, which denotes the number of images in the coarse filtering result, is set to the length of G j . The h-th gene A h = 0 indicates that image collection G j does not include the h-th image in the coarse filtering result, and A h = 1 indicates that image collection G j includes the h-th image in the coarse filtering result. The coding method for an image collection is shown in Figure 2. After randomly generating multiple binary codes of length m, we can form the initialization population. In this paper, we set P, which denotes the number of image collections in the population, to an empirical value of 10.

Fitness Score Function
The second important process of the GA is to design a reasonable fitness function to describe the satisfaction of each image collection G j in terms of coverage and user demands, which are denoted by S G j . In this paper, we use a grid as the computing unit to calculate the fitness score of the image collection.
To obtain S G j , we first need to define S g i , which denotes the fitness function of the image collection on each grid g i . Since usually more than one image can cover grid g i in G j in most cases, we calculate the S I (Equation (2)) of all images that can cover grid g i in G j and take the highest score as S g i . The calculation formula for S g i is as follows: where I j ∩ g i = ∅ indicates that the intersection of the j-th image and the i-th grid is not empty. Then, we need to define the fitness function S G j of the image collection. Considering that the fitness score is related to more than the satisfaction of the user's demands, the coverage also has a certain impact. We define the parameter OL to represent the coverage of image collection G j over the Target_Area. OL is calculated as follows: where cr denotes the coverage rate, n GridsO f Covered denotes the number of grids covered by the image, ct i denotes the number of times the i-th grid is covered by the image, and n Grids denotes the number of grids in the administrative area. A value of less than 1 for OL indicates that the Target_Area cannot be completely covered by G j ; a value of 1 for OL indicates that the Target_Area is completely covered by G j ; and a value larger than 1 for OL indicates the number of layers that repeatedly cover the Target_Area. Commonsensically, when OL is equal to 1, G j can cover only the Target_Area without redundancy, and S G j should be the largest at this time. When OL is greater than 1, S G j will decrease by a multiple. When OL is less than 1, G j cannot completely cover the Target_Area; users cannot tolerate the lack of coverage. Considering the above characteristics, S G j is calculated as follows: The above formula shows that when OL = 1, S G j is equal to adding 1 to the scores of all grids S g i . Bias 1 is added to ensure that the score for full coverage is greater than the score for incomplete coverage.
After the population is initialized by the GA, reproduction, crossover, and mutation are carried out to simulate the inheritance process to obtain the best gene.
We iterate the selection, crossover, and mutation steps until convergence is reached. The condition for convergence is set as the highest fitness score of the image collection not changing for fifty iterations. In Figure 3, the red line means that although the preset threshold is not reached, as the generation increases, the fitness score changes by less than 0.1 after 50 iterations. This means that the algorithm has converged and the current optimal image collection can be the final result.

Filtering Result Optimization
Through the above steps, the optimal image collection in the last population can be obtained; this collection is denoted by G best . Although the GA can guarantee evolution to a better population, there are some small flaws, which are generally caused by one or two redundant images. Therefore, we need to optimize the structure of G best . The optimization involves traversing each gene A j in G best . When A j = 1, A j will be set to 0 to generate a new chromosome G best new . If S G best is then increased, the new chromosome G best new will be replaced with G best until there is no further increase in S G best . G best is then the final image collection filtering result.

Experimental Region
Different levels of government urban planning departments need to observe regions of different scales. Provincial governments need to observe and conduct spatial statistical analysis on the whole province, and urban planning departments at the city level are concerned with whole cities. In order to verify the accuracy and efficiency of the LFCF-ID for different levels of administrative regions, the experiment used multiple administrative regions in China as experimental regions, including Shijiazhuang (medium city), Beijing (large city), and Hebei (province), which are shown in Figure 4.

Data Source
To effectively verify the ability of the LFCF-ID to filter image collections, the data used in the experiment should be as diverse as possible. Therefore, data from multiple satellites was selected; the main satellites selected include GaoFen, ZiYuan, and Jilin, which use multispectral, synthetic-aperture radar (SAR), hyperspectral, and other sensor types with spatial resolutions from sub-meter to hundreds of meters and widths from 5 km to 720 km. Detailed information is shown in Table 2.

Data Source
To effectively verify the ability of the LFCF-ID to filter image collections, the data used in the experiment should be as diverse as possible. Therefore, data from multiple satellites was selected; the main satellites selected include GaoFen, ZiYuan, and Jilin, which use multispectral, synthetic-aperture radar (SAR), hyperspectral, and other sensor types with spatial resolutions from sub-meter to hundreds of meters and widths from 5 km to 720 km. Detailed information is shown in Table 2. In this paper, we selected 111,414 images from the above satellites within the time range of 2008 to 2018. We calculated the statistics of the data from the satellites and sensors; these statistics are shown in Figure 5.

Image Data Distribution of Different Types and Regions
To more intuitively describe the distribution of data, we statistically distributed the data in three aspects: time, resolution, and cloudiness. Figure 6a shows the statistics of the images in terms of four levels: cloudless, low cloudiness, medium cloudiness, and high cloudiness. The number of cloudless images is the largest, and the number of images decreases with increasing cloudiness. Figure 6b shows the statistics of the number of images at different spatial resolutions. Although there is a certain difference in the amount of data in each section, there are enough images in each section. Figure 6c shows the statistics of the amount of image data for each year from 2008 to 2018 in units of years. This image collection covers all years and is concentrated from 2015-2017. To verify the applicability of the proposed method in different regions, this paper selected three areas of different sizes as experimental areas: Shijiazhuang, Beijing, and Hebei; additionally, we statistically analyzed the data distributions of different areas. Figure 6d shows that the number of images is proportional to the area of the region, and the image data for each area are sufficient. In summary, the distribution of the data used in the experiment in different types and regions is sufficient and is thus able to support the verification of this study.

Experiment and Analysis
To fully verify the accuracy, robustness, and efficiency of the proposed method, multiple experiments were designed from different perspectives. First, the top k images were selected for each grid during coarse filtering; this selection can directly affect the efficiency of the algorithm by affecting the length of the chromosome. This paper compared experiments with different k values. Then, under the optimal k value, experiments were performed in different regions, with different numbers of images to be filtered, and different weights of users' demands. Designing experiments using different regions can verify the robustness of the algorithm. Experiments using different numbers of images to be filtered can directly verify the speed of the algorithm. Designing experiments with different demands can verify whether the proposed method can meet the users' individual demands. In addition, this article also designed experiments to compare this method with manual methods and other algorithms to verify the superiority of the proposed method.

a. Experimental results for different k values
We take the data from Shijiazhuang in 2018 as the input and set W gsd , W cloud , W time to {0.333, 0.333, 0.333}. We set the k value from 1 to 6. The experimental results are shown in Figure 7 and Table 3. The black box in Figure 7 indicates the boundary of each remote sensing image.   Figure 7 shows that regardless of the k value, the full-coverage requirements of the filtering results are basically met, but the composition of the image collection is different. Table 3 shows that when k = 1, S G best can reach only 1.661959, which is far lower than the results obtained under other k values. This result occurs because the coarse filtering in the algorithm when k = 1 can be understood as equivalent to the greedy solution, which roughly finds a single image that best meets the user's demands for each grid, and the filtering results can be obtained by combining the results of all grids. Since the size of each remote sensing image is larger than the grid, greedy conditions will make the OL of the image collection very large; this can result in a low fitness score.
When k gradually increases from 1 to 3, the total score continues to increase, and the image number of G best gradually decreases. This result indicates that better image collection can be achieved when the number of images participating in the GA gradually increases.
When k continues to increase from 3 to 6, the score gradually decreases. This result indicates that when the value of k increases to a certain level, the quality of the newly added images will decrease compared to when k is small, and the improvement is not obvious. When k further increases, the length of G i will also increase exponentially, thus affecting the results of the GA by leading to a local optimum and making it difficult to find the optimal solution. For example, when k = 6 in Table 3, the length of G i reaches 163 and S G best is reduced to 1.691373.

b. Experimental results for different regions
To verify that the proposed method has robustness for different areas with irregular shapes, this study performed experiments in Shijiazhuang, Beijing, and Hebei within a fixed time interval under the premise of setting W gsd , W cloud , W time to {0.333, 0.333, 0.333}. The filtering results are shown in Figure 8.  Figure 8 shows that when the proposed method is used for filtering in different regions, although the area and shape may differ greatly, good coverage can be obtained for all regions. Table 4 shows that the OL values of the image collection after filtering in different regions are 1.320225, 1.243119, and 1.432177, all of which are below 2, and there are no repeated overlays caused by image redundancy. To verify the robustness of our algorithm, we performed experiments with 50 iterations in Beijing, Hebei, and Shijiazhuang under the same input (k set to 3, the cloud weight set to 0.3, the timeliness weight set to 0.35, and the resolution weight set to 0.35). We produced two boxplots of fitness scores and OL. The boxplot is a statistical graph used to display the dispersion of a dataset and can therefore be used to verify the robustness of the algorithm. The boxplots are shown in Figure 9.
In the boxplot, the middle line of the box is the median of the data and thus represents the average level of the sample data. The upper and lower limits of the box are the upper and lower quartiles of the data, thus indicating that the box contains 50% of the data, so the broadband of the box reflects the degree of data fluctuation to a certain extent. There are additional lines above and below the box; these lines represent the maximum and minimum values. Sometimes, some points are outside of the upper and lower limits, and these points are outliers. As shown by the boxplots of the fitness scores, the distance between the upper quartile and the lower quartile of Hebei is within 0.07, while the distances of Beijing and Shijiazhuang are within 0.05, thus indicating that the fluctuation of the fitness score in the 50 experiments was very small. As the OL boxplot shows, the distance between the upper quartile and the lower quartile of Hebei is within 0.18, while the distances of Beijing and Shijiazhuang are within 0.1, thus indicating that the fluctuation of OL in the 50 experiments was very small. The distance between the maximum and minimum values in the two boxplots is also within the acceptable range, and there are no outliers. In summary, this algorithm is very robust.

d. Experimental results for different image numbers of I_ori
To verify the efficiency of this method when processing images of different orders of magnitude, we performed experiments under the premise of fixing other conditions and controlling the number of images processed by setting time intervals of different lengths. The results are shown in Figure 10 and Table 5.  Figure 10 shows that the proposed method exhibits good filtering results for input images of all orders of magnitude. In detail, Table 5 shows that with the gradual increase in the amount of I_ori, the OL of I_ori increases from 14.93578 to 500.2661, and the OL of G best decreases to approximately 2. The fitness score of G best increases steadily from 1.5579 to 1.750479. This result shows that the proposed method can obtain better results regardless of whether the amount of data is in the hundreds, thousands, or tens of thousands. An increasing amount of input data provides more choices for the filtering process and produces better results. Table 5 shows that the time consumed by the proposed method remains stable when the order of magnitude of the input images increases; the time consumed by traditional methods should increase exponentially. Different types of urban study and planning tasks have different demands for remote sensing images. For example, in the task of urban road planning, planners need to master the distribution of all roads in the city, which requires high resolution of remote sensing images. In the urban real-time monitoring task, it is necessary to carry out high-frequency real-time monitoring of key areas such as illegal buildings in the urban area, and high timeliness of remote sensing images is required. In order to effectively verify whether the proposed method can meet the demands of urban study and planning tasks, multiple experiments were designed from high to low in terms of timeliness, resolution, and cloudiness.
To quantitatively evaluate the cloudiness, resolution, and timeliness of the image collection, three measures are defined: the average grid cloudiness is denoted by Aver_Cloud, the average grid resolution is denoted by Aver_GSD, and the average grid timeliness is denoted by Aver_Time. We obtain the g i cloudiness score, which is denoted by S cloud g i ; the g i resolution score, which is denoted by S gsd g i ; and the g i timeliness score, which is denoted by calculating S gsd (Equation (4)), S cloud (Equation (5)), and S time (Equation (6)) of the image with the highest score in grid g i . We obtain Aver_Cloud, Aver_GSD, and Aver_Time by calculating the average S cloud g i , S gsd g i , and S time g i for all the grids. The formulas are as follows: where n is the number of grids covered by the Target_Area. In Figure 11, W gsd in the second line of images gradually decreases from left to right from 1 to 0.1. The color of the images in the filtered image collection gradually changes from all green to partially red and partially pink until blue appears, thus indicating that the image resolution in the filtered results is gradually decreasing. In the task of urban road planning, planners can set W gsd to a larger value such as the first figure in the second row. Planners can obtain an image collection with a resolution better than 3 m to achieve full coverage of the Shijiazhuang area in which the road and other small-scale landcover can be seen. W time in the third line of the image gradually decreases from left to right from 1 to 0.1. This figure shows that the color of the images in the filtered image collection gradually changes from all green to partially red, and from having little blue to partial blue appearing, thus indicating that the time of the images in the filtered result is far from the last date. In the task of urban real-time monitoring, people can set W time to a larger value such as the first figure in the third row. Then, an image collection within 10 days can be obtained to observe buildings and other illegal landcover.
W cloud in the first row of images gradually decreases from left to right from 1 to 0.1. The color of the images in the filtered image collection gradually changes from all green to partially red and partially pink until blue appears, thus indicating that the cloudiness of the images in the filtered result is gradually increasing. For the task of a large-scale landcover survey, people can set W cloud to a larger value such as the first figure in the first row. Then, an image collection without any clouds can be obtained to avoid situations in which land is blocked. Table 6 provides detailed information for each image in Figure 12. Figure 12 is a visualization of the changes in Aver_Cloud, Aver_GSD, and Aver_Time according to W gsd , W cloud , W time ; these changes are also shown in Table 6. As shown in Figure 12a, as the cloudiness increases, the Aver_Cloud score gradually increases from 0.965056 to 1, and the scores of Aver_GSD and Aver_Time gradually decrease from 0.961545 and 0.946545 to 0.708643 and 0.225515, respectively. In Figure 12 b, as the grid resolution increases, the Aver_GSD score gradually increases from 0.737714 to 0.986794, and the Aver_Cloud and Aver_GSD scores gradually decrease from 0.999555 and 0.986116 to 0.958371 and 0.675474, respectively. In Figure 12 c, as time increases, the Aver_Time score gradually increases from 0.869792 to 0.985105, and the Aver_Cloud and Aver_GSD scores gradually decrease from 0.994994 and 0.977812 to 0.99764 and 0.629191, respectively. In summary, the LFCF-ID can filter image collections according to preset weights to meet individual demands.

Comparison with Other Methods
There is currently no other study similar to this paper, so the algorithm in this paper is only compared with manual filtering and greedy methods. The manual filtering method refers to the process by which a person selects from several images to obtain the results that he or she feels meet the demands. The greedy method refers to traversing all the images of each grid to obtain the best image and combining the best images of all the grids as the filtering result. The experimental results are shown in Table 7.    Table 7 shows that in terms of time consumption, the greedy LFCF-ID (k = 1) and the LFCF-ID (k = 3) outperform the manual method. The greedy algorithm lacks optimization of the filtering results and directly uses the results of the grid retrieval as the output, so this algorithm takes the least time. However, the fitness score of the greedy method is poor compared to the fitness scores of the manual, LFCF-ID (k = 1), and LFCF-ID (k = 3) methods, thus indicating that the filtering effect of the greedy method is very poor. The LFCF-ID (k = 1) can be considered genetic optimization based on greedy retrieval. The LFCF-ID (k = 1) is slower than the greedy algorithm, but the fitness score is greatly improved. Compared with the LFCF-ID (k = 1), the LFCF-ID (k = 3) allows more images to be subsequently added to genetic optimization. Due to the increase in gene length, the time is slightly increased, and the fitness score is higher than that of other methods and is even better than that of the manual method. In summary, the comparison indicates that the LFCF-ID is 10 times faster than the manual method and exhibits a filtering effect that is better than the filtering effects of all existing methods.

Conclusions
This paper proposed a new remote sensing image filtering method that can support people to maximize their preferences in urban observation for study and planning tasks, on the premise of ensuring full coverage of a region. We first designed a coarse filtering strategy to reduce the dimensionality of the data; this reduction can save the most useful images and save time for subsequent optimization calculations. Then, a grid was used as the basic unit to calculate the fitness score of the image collections in terms of resolution, cloudiness, and coverage; the score was used as the index to evaluate the performance of the image filtering algorithm. Finally, the evaluation index was optimized by genetic iteration to optimize the image collections. We designed different sets of experiments to evaluate the performance of the LFCF-ID; the experimental results showed that the LFCF-ID could quickly achieve full-coverage filtering of images and simultaneously minimize repeated coverage as much as possible with different regions, different data volumes, and different demand preferences. The LFCF-ID showed great potential in solving how to automatically and quickly obtain full-coverage image collections of areas that meet the demands of preferences, in a context where people are facing an explosive increase in remote sensing image data. It is foreseeable that, in the future, the LFCF-ID will have wide engineering applications in urban study and planning.