Introduction

Remote sensing imagery is a crucial data source for studying surface land use and environmental evolution. Currently, the longest-running global land use datasets and their derived data (such as wetland datasets1, water body datasets2, etc.) are primarily based on imagery from the Landsat series of satellites. However, constrained by satellite launch dates, global studies on large-scale surface processes often commence from the 1980s3. Additionally, the spatial resolution of early datasets is relatively low (60 m)4, making them less suitable for the demands of long-term geospatial analysis. Therefore, extracting feature information from graphical data preceding current remote sensing imagery, such as archived maps5 or early satellite images, is crucial for extending the timeline of surface change research.

Keyhole imagery refers to images obtained by a series of spy satellites launched by the United States starting in the 1960s6. Between 1995 and 2020, the U.S. government progressively declassified the images from the KH-1 to KH-9 satellites, covering most global surface areas from 1960 to 19807. Declassified Keyhole imagery could be divided into free and non-free data, the latter costing only about $30 per scene8. Due to several advantages of Keyhole imagery, such as early imaging dates, high ground resolution, extensive coverage, and low cost, this imagery holds significant research value globally in fields such as resource monitoring, disaster response, and climate change. Compared to the earliest Landsat data from the 1980s, Keyhole imagery (1960s) can extend remote sensing-based research by nearly two decades9,10.

The declassification of Keyhole satellite imagery has facilitated early surface changes research that were previously difficult to undertake11,12. Since 2020, the number of published papers on Keyhole imagery has nearly equaled the total from 2000 to 202013. Given the effectiveness of glacier identification with single-band panchromatic imagery, glacier area changes14, volume melting15, and responses to climate change16 have been investigated using 5 ~ 10 m resolution grayscale images from Keyhole satellites (KH-4, KH-9) or stereo image pairs. Nister et al.17applied 1.8 m resolution KH-4B imagery to analyze the conversion from farmland to forest in Budapest since the 1960s. Their research revealed that land abandonment in the study area emerged in the 1980s, earlier than the traditionally perceived period of the Soviet Union’s collapse in the early 1990s. This finding enriched the understanding of the drivers behind land use changes in Eastern Europe. Keyhole imagery has documented the land cover characteristics before urban expansion in developing countries18. Therefore, the high-resolution Keyhole imagery could detect ancient hidden cultural relics or tectonic geomorphology that may be completely lost in the modern landscape obscured by current land uses and land covers. Keyhole imagery could also be used to investigate long-term evolution of various degraded erosion areas, such as the gully erosion in the black soil region of Northeast China (KH-4A, 2.7 m resolution)19,20and the loess region of Northwest China (KH-4B, 1.8 m resolution)21, as well as the karst topography in Guizhou (KH-4A, KH-4B)22. Although the previously mentioned studies have highlighted the utility of Keyhole satellite imagery, they primarily focused on analyzing specific events or small-scale areas involving single or multiple scenes; however, research is scarce on the national or regional scales’ application of Keyhole imagery.

Applications of remote sensing imagery at national or regional scales require complete area coverage by satellite images. The frequently utilized ETM, SPOT, and other imagery feature regular revisit periods23, ensuring a relatively uniform spatial distribution of image coverage. However, Keyhole satellites were film-return satellites, and their lifetimes were short, ranging from several to dozens of days24. Additionally, the ground resolution of different types of Keyhole imagery varies from sub-meter to tens of meters. Hence, the spatial coverage frequency and resolution of Keyhole satellites might exhibit notable spatial variability, potentially affecting the selection and application of Keyhole imagery data at national or regional scales25. Due to the irregularity of above factors across different Keyhole images, detailed description of these factors’ distribution and heterogeneity become crucial to provide references for subsequent processes, such as feature recognition and interpretation. Those preprocessing steps are often neglected in published research, which typically focus on extraction methods26. Research on the archived maps demonstrated that preliminary analysis using metadata could provide valuable insights for image processing27. Similar to archived maps that include metadata, the Keyhole imagery provided by the United States Geological Survey (USGS) comes with the vector boundary files containing metadata, including bounding coordinates, imagery date, and resolution information. Hence, the metadata of Keyhole imagery offer great potential for automated analysis of distribution characteristics of Keyhole imagery.

Since early Keyhole satellites were primarily applied to monitor military activities in Asia and Eastern Europe, the imagery coverage in China is expected to be high. Therefore, this paper selected China as the case study area and analyzed the distribution characteristics of Keyhole imagery using quadrangle boundaries, resolution and imagery date. All images were categorized by resolution range into meter-level, five-meter-level, and ten-meter-level. Geostatistical and statistical methods investigated the spatial coverage distribution characteristics and spatial heterogeneity of meter-level, five-meter-level, and ten-meter-level Keyhole satellite imagery in China at national and provincial scales. Based on the assumption that image applications require triple coverage to ensure the selection of qualifying images (such as excluding cloud-covered images) and considering areas covered less than three times, the supplementary imagery cost for triple coverage imagery of the three different resolution categories in China was calculated. Given the large data volumes of Keyhole imagery and the regional extensibility of the spatial distribution analysis, this research automates the data analysis process and provides the necessary codes to promote the extensive application of Keyhole imagery beyond this case area.

Materials and methods

Data sources

The bounding vector files for all freely available Keyhole satellite images of the China region were downloaded from the United States Geological Survey (USGS) website28. The downloaded vector files were then transformed into the Krasovsky_1940_Albers projection. The attribute table of the vector files delineated the metadata of the image. Combined with the overall descriptive files of the Keyhole imagery29, the resolution of each image series was determined (Table 1). Based on the characteristics of the satellite system, the image series were named KH-1, KH-2, KH-3, KH-4, KH-4A, KH-4B, KH-5, KH-6, KH-7, KH-9L, and KH-9H, covering a range of ground resolutions from sub-meter to ten meters. Details on the resolution, acquisition time, single scene area, and number of images of different satellite series were provided in Table 1. China’s national and provincial boundaries data were obtained from the National Geomatics Center of China30.

Table 1 Metadata and quantities of different series of Keyhole images in China.
Full size table

Due to significant differences in resolution among the different series, this paper categorized Keyhole imagery into three classes: meter-level (C1), five-meter-level (C2), and ten-meter-level (C3), to summarize the distribution characteristics of Keyhole imagery. Since the Keyhole imagery provided by the United States Geological Survey were replicas of the original films, the actual resolution might be lower than the described resolution12. Therefore, KH-1 and KH-9H imagery with a resolution of 1 m or less were categorized as meter-level images (C1); KH-4A and KH-4B imagery, with resolutions higher than 1 m and less than 5 m, were classified as five-meter-level images (C2); KH-2, KH-3, and KH-4 imagery, with resolutions higher than 5 m and less than 10 m, were categorized as ten-meter-level images (C3). The total numbers of images for C1, C2, and C3 were 1,344, 11,237, and 5,040, respectively. To investigate the overall spatial distribution characteristics of Keyhole imagery, all images were compiled and classified as C0. The KH-5, with the lowest resolution of 140 m, was not included in the classification.

Data preparation

Imagery classification using attribute table

USGS restored the declassified Keyhole imageries in three imagery collections, namely declass1, declass2 and declass3. Each collection contained several types of satellite images. Since the attribute table of each collection described satellite information using mission number instead of mission designator (e.g. KH-4 or others), mapping between imagery resolution and mission number needs to be established to classify the images based on their resolution. Utilizing the correspondence between mission number and designator as well as imagery resolution described in the Declassified Intelligence Satellite Photographs report31, the imagery resolution and designator information were added to the attribute tables. Then all imageries were classified by resolution as C1, C2 or C3 (Table 1).

Generation of position-time pointset using fishnet tool

Here we took C1 imagery as an example to describe the data analysis process, and the approach applied for other imageries were similar. Appling the Fishnet function in ArcMap 10.2, a grid of points spaced at 10 km intervals within China was generated, with the Krasovsky_1940_Albers projection as the coordinate system. Each grid point could be covered by multiple C1 images, meaning that each grid point was associated with multiple time records. Therefore, the Spatial-Join function was applied to connect the grid points with boundaries of C1 images. The Join Operation parameter was set as JOIN_ONE_TO_MANY, and only the imagery time field from the attribute table of the image boundary was retained. The resulting shapefile was a Position-Time (PT for short) point set, where multiple PTs could exist at each grid ___location, with each PT recording the time information of a specific image at that ___location.

Data analysis

General distribution of keyhole imagery

We firstly visualized the spatial and temporal distribution characteristics of Keyhole imagery in China, allowing users to immediately gain a general understanding of the imageries’ distribution characteristics. For spatial distribution, all C1 images were overlaid with the research area to visually represent the coverage extent of C1 imagery. The temporal distribution referred to the variation of average imagery coverage area and frequency over the years. The total number of non-zero grid points and the total number of non-zero PT points were summarized annually from 1960 to 1984. The ratio of the total number of non-zero grid points in a given year to the total number of grid points represent the C1 imagery coverage proportion for that year. The ratio of the total number of non-zero PT points to the total number of non-zero grid points for that year represent the average coverage frequency of C1 imagery.

Spatial distribution of imagery coverage frequency at national scale

The coverage percentage of each imagery was calculated based on the number of points with non-zero coverage counts to total points in China. Classical statistical analysis was employed to analyze the mean, standard deviation, coefficient of variation, skewness, and kurtosis for C0, C1, C2, and C3 imagery with non-zero coverage. Linear regression analysis was conducted based on coverage frequency of C0, C1, C2, and C3 imagery at the same point, to assess the correlation among their coverage frequencies.

Since earth observations based on remote sensing require at least one image coverage, and more coverages increase the selection flexibility of images, the number of coverages at each point for each image resolution category was classified into six levels: 0, 1, 2, 3–5, 5–20, and 20 to the maximum. A coverage level of 0 indicated that no Keyhole imagery was available for the point. The coverage level of 1 indicated only one available Keyhole image for the point, meaning there was no alternative options. The coverage level of 2 indicated two available Keyhole images, meaning a choice was available. The coverage level of 3 to 5 indicated that multiple options were available, and coverage levels from 5 to 20 and 20 to the maximum indicated a substantial selection of Keyhole images were available. Frequency histograms for each coverage level of each resolution imagery were also generated.

Utilizing the number of coverages at each point, inverse distance weighting interpolation was applied to generate spatial distribution maps of the coverage frequency for each resolution category across China.

Spatial autocorrelation methods were applied to analyze image clustering to verify the spatial heterogeneity of the imagery. The Global Moran’s I index32 was used to assess the overall spatial characteristics of Keyhole imagery coverage frequency across the region. It determined whether the coverage frequency in a specific area was related to the coverage frequencies in adjacent areas. The Global Moran’s I value typically ranges between −1 and 1. A value greater than 0 indicates a significant positive spatial correlation, a value of 0 indicates no spatial autocorrelation and a value less than 0 indicates a significant negative spatial correlation. The Z-score measures the degree of spatial autocorrelation. A positive and significant Z-score indicates positive spatial autocorrelation, meaning that similar values (either high or low) are spatially clustered. Conversely, a negative and significant Z-score indicates negative spatial autocorrelation, meaning that similar values tend to be dispersed. A Z-score of 0 suggests that the values are randomly distributed with no spatial dependence. The P-value indicates the statistical significance of these findings. The Local Moran’s I index was applied to analyze the spatial distribution patterns of coverage frequencies at each point, measuring the degree of local spatial association between each area and its adjacent regions.

Provincial spatial distribution of keyhole imagery

Combining provincial boundaries, the spatial join method was employed to determine each province’s coverage frequency of each image category. Based on the number of fishnet points with non-zero coverage counts for C0, C1, C2, and C3 images, the coverage rate of each image type, the average coverage frequency, and the coefficient of variation (CV) were calculated for each province.

Calculation of triple coverage costs for images of different resolutions

Multiple image coverages were necessary for remote sensing applications at the regional scale to fulfill the requirements for image selection. Each resolution category of Keyhole imagery included images from different satellite series with varying parameters. Therefore, selecting a specific type of image as a priority supplementary imagery for uncovered areas could avoid the impact of differing image parameters and facilitate uniform data processing. Since each type of non-free Keyhole satellite image costs $30 per scene, a larger single-scene coverage area resulted in a lower cost per unit area. Therefore, this paper selected the single scene coverage area and the total number of images as indicators to determine the priority supplementary imagery among C1, C2, and C3 imageries.

Since this paper assumed that triple coverage was needed to construct the Keyhole imagery dataset for all areas, it was necessary to calculate the areas with 0, 1, and 2 times coverage for C1, C2, and C3 imageries to estimate the costs of acquiring supplementary images. If the image coverage within a specific area was once, then two additional imageries needed to be purchased. Due to a specified degree of overlap between images needed to facilitate adjacent image georeferencing, this study assumed an overlap rate of one-half for the supplementary images. This meant that three images needed to be acquired within the coverage area of two images, resulting in an overlap factor of 1.5. Therefore, the formula to calculate the acquisition cost of supplementary imagery for a certain coverage frequency of a specific type of imagery was as follows:

$${S}_{[i,j]}= {U}_{[i]} \times {A}_{\left[i,j\right]}\times {N}_{[i,j]}\times \text{ f}$$
(1)

Here, i represents a specific resolution of imagery, with values ranging from 1,2 and 3. The value '1' represents meter-level imagery (C1); the value '2' represents five-meter-level imagery (C2); and the value '3' represents ten-meter-level imagery (C3). j represents the coverage frequency of the imagery, with values of 0, 1, and 2. S[i, j] represents the cost of supplementary images for the area with j times coverage of Ci class imagery, expressed in ten thousand dollars. U[i] represents the unit area cost of Ci imagery, expressed in dollars per square kilometer. A[i, j] represents the area covered j times by Ci imagery, expressed in ten thousand square kilometers. N[i, j] represents the number of images that need to be acquired for areas with j times of coverage of Ci class imagery, with a value of 3—j. f represents the overlap factor between adjacent images, with a value of 1.5.

The supplementary imagery acquisition cost, Si​ for Ci imagery, is the sum of the acquisition costs for areas with 0, 1, and 2 times of coverage:

$${S}_{i}=\sum_{j}{S}_{[i,j]}$$
(2)

Results

Distribution of keyhole images

Spatial distribution

Overlaying C1, C2, and C3 class imagery with the regional boundaries of China revealed that each of these imagery classes could only cover part of the research area (Fig. 1 b, c, and d). Even C0, which combined the three types of imagery, could not cover the entire region (Fig. 1 a). Coverage gaps were primarily in southwestern Xinjiang, central Tibet, and parts of Sichuan. Additionally, the primary coverage areas of the different types of Keyhole imagery were predominantly land, with very few images covering oceanic regions. Therefore, the analysis in the following part focused on land areas’ imagery coverage.

Fig. 1
figure 1

Distribution of Keyhole imagery in China based on resolution classification (a, C0, all images; b, C1, meter-level images; c, C2, five-meter-level images; d, C3, ten-meter-level images).

Full size image

Yearly distribution

The sub-meter C1 imagery was available for the years 1964 to 1967 and 1971 to 1984, with broader coverage in the latter period (Fig. 2). The coverage was greatest in 1972 and 1973, at 27% and 22%, respectively. The C2 imagery covered the period from 1963 to 1972. Annual coverage consistently exceeded 25% between 1964 and 1970, with peak value of 45% at 1965. The ten-meter C3 imagery spanned from 1960 to 1964 and from 1973 to 1981. High coverage occurred in 1962 and from 1973 to 1975, ranging from 25 and 30%. The temporal gaps of in C1 and C3 imagery distribution were due to each resolution imagery being deprived from multiple satellites. The average annual coverage frequency for the three types of imagery mostly ranged between 1.0 and 1.5, indicating that the time interval in a given area was at least one year. Considering three resolution imageries, imagery was available from 1960 to 1984 with an average rate of 25%. The imagery coverage exceeded 30% almost each year between 1962 and 1974, with the maximum value being 64% in 1964. Although the overall annual coverage frequency considering three types of imagery together also mostly ranged between 1.0 and 1.5, it still exceeded the average coverage frequency of each resolution imagery. Research about muti-satellites effect on revisit intervals showed that Landsat-8, Sentinl-2A and Sentinl-2B together would provide a global median average revisit interval of 2.9 days, comparing with 16, 10 and 10 days for these three satellites, respectively33. Additionally, the Keyhole imagery discussed in this research referred exclusively to the free imagery, not the entire dataset that included non-free imagery. Across China, there were 334,779 images in total across China, while only 19,901 of these were freely available, constituting just 6% of the entire dataset.

Fig. 2
figure 2

Yearly distribution of Keyhole imagery average coverage area and frequency in China from 1960 to 1984.

Full size image

Statistical characteristics of keyhole imagery spatial coverage

All keyhole imagery (C0) covered 99% of China’s area, with an average coverage of 9.9 times (see Table 2). The average coverage frequencies of the three resolution categories (C1, C2, and C3) were relatively close, at 4.9, 4.5, and 3.6 times respectively. However, there were significant differences in the coverage areas of the three categories of imagery. Low-resolution imagery covered a larger area, with C2 and C3 covering 95% and 76%, respectively; meanwhile, high-resolution imagery had a lower coverage area, with C1 covering only 58%. Although the coverage area of C1 was lower, the average coverage frequency in the areas where images were available was the highest, at 4.9 times.

Table 2 Classical statistical analysis of Keyhole imagery coverage.
Full size table

This research applied standard deviation and coefficient of variation, which could reflect coverage frequencies’ absolute and relative dispersion, to analyze Keyhole imagery’s spatial coverage variability characteristics. Overall, the coefficient of variation for the coverage of the four categories of imagery ranged from 0.8 to 1.4, indicating substantial variability. The standard deviation of C0 was the highest at 9.1, while its coefficient of variation was relatively low at 0.9. This showed that although the absolute fluctuations of C0 image coverage frequency were significant, the dispersion relative to the average was comparatively small. For C1 imagery with the lowest coverage area, the coefficient of variation was the highest at 1.4, and the standard deviation was 6.8 times. This indicated that meter-level imagery with the highest resolution category exhibited the highest level of spatial variability. The variability levels for C2 and C3 were similar, around 0.8.

Due to the high coefficient of variation in image coverage frequency, skewness was further applied to characterize the asymmetry of coverage frequency distribution, and kurtosis was used to determine the heaviness of the distribution tails. All categories of imagery had positive skewness values, with the skewness for meter-level imagery being the highest at 3.8. Positive skewness values indicated a distribution of coverage frequencies skewed towards lower values, meaning that more points exhibited coverage frequencies below the average. In contrast, fewer points displayed coverage frequencies above the average. Consequently, most areas in the study region had a low imagery coverage frequency, with only a few areas having a higher coverage frequency. Additionally, the kurtosis values for all categories of imagery were positive and relatively high, indicating that the data distributions had long tails and sharp peaks. C2 and C1 had the highest kurtosis values, at 19.3 and 18.1 respectively. A positive skewness coefficient and a high positive kurtosis coefficient indicated that although most areas had a low coverage frequency, there were a few areas with coverage frequencies significantly higher than the average, resulting in a heavy-tailed distribution. In summary, a few areas were monitored intensely in all four imagery classes, while most areas exhibited coverage frequencies below the average.

Pearson correlation analysis based on the on the coverage frequency of four types of imagery at the same locations indicated that the p-values for all categories reached 0.0, demonstrating statistically significant correlations. The correlation coefficients between C0 and C1, C2 were 0.80 and 0.81 (Table 3), respectively, indicating strong positive correlations. This suggests that the increased coverage frequencies of meter-level and five-meter-level imagery significantly increased overall image coverage frequencies. Relative to C1 and C2, the correlation between C0 and C3 remained positive, while the strength of the correlation was reduced (0.68). The correlation coefficients between C1 and C2, C3 were 0.40 and 0.24, respectively, indicating a lower positive correlation. This showed that although there was a positive relationship between the coverage frequencies of meter-level imagery and those of five-meter and ten-meter imagery, the correlation was low, especially with ten-meter imagery. The correlation coefficient between C2 and C3 was 0.52, indicating a moderate positive correlation. This showed that although ten-meter and five-meter imagery were correlated, the relationship was weaker than the overall imagery correlations with meter-level or five-meter imagery.

Table 3 Correlation analysis of Keyhole imagery across different resolution categories.
Full size table

Spatial distribution characteristics of imagery coverages

Spatial interpolation analysis

Frequency histograms based on image coverage classification are shown in Fig. 3. For C1, C2, and C3, the proportions of areas with only one coverage were 15%, 13%, and 17%, respectively, and areas with two coverages displayed similar proportions as 13%, 18%, and 16%. Areas with 3 to 5 coverages exhibited the highest proportions, accounting for 16%, 37%, and 29% for C1, C2, and C3, respectively. In the coverage frequency range of 6 to 20 times, which could offer a wide selection of images, the coverage areas accounted for 11%, 27%, and 14% in the C1, C2, and C3, respectively. The areas with more than 20 coverages were relatively limited, with the proportions being 2.1%, 0.5%, and 0.0%, respectively. Among C1, C2, and C3, C1 had the highest resolution and the highest proportion of areas covered more than 20 times, indicating that certain regions were frequently covered by high-resolution imagery. Overall, the highest proportion of areas covered by Keyhole imagery (C0) in China fell within the 6 to 20 times coverage range, accounting for 57%. Areas with 3 to 5 times coverage were 23%, with only 12% having less than 3 times coverage.

Fig. 3
figure 3

Spatial and frequency distribution of Keyhole imagery coverage based on resolution classification (a, C0, all Keyhole images; b, C1, meter-level resolution images; c, C2, five-meter-level resolution images; d, C3, ten-meter-level resolution images).

Full size image

The spatial distribution of image coverage frequencies (Fig. 3) indicated that all four types of imagery were concentrated in the Jiangsu-Zhejiang area, southwestern Heilongjiang, the Beijing-Tianjin area, southern Henan, central Xinjiang, and western Xizang. This concentration might be related to the primary military application of Keyhole imagery. Additionally, Keyhole Imagery from neighboring regions in western China might contribute to the increased coverage rates along the border of the west. Image coverage rates were relatively low in southwestern China, northern Xinjiang, Inner Mongolia, Qinghai, and northern Xizang. Uncovered areas for all categories of imagery were concentrated in central and southern Xinjiang and northern Xizang.

The highest-resolution C1 imagery (Fig. 3b) within specific imagery categories was primarily concentrated in eastern and central China. The highest coverage was observed in the east of coastal areas, southwestern Heilongjiang, and parts of central China. Conversely, northwestern regions of China generally lacked imagery with at least two times coverage. The spatial distribution of C2 imagery (Fig. 3c) was relatively uniform across the study area. High-coverage areas (with more than 20 coverages) and low-coverage areas (with fewer than three coverages) were distributed in a scattered pattern, and there was a lack of imagery for northern Xizang and southern Xinjiang. C3 imagery was primarily distributed in the northeastern, eastern, southern coastal, and western border regions of China. However, there was a lack of imagery for the contiguous region extending from southeastern Xinjiang to Fujian and northern Inner Mongolia.

Spatial autocorrelation analysis

Spatial autocorrelation analysis was performed on Keyhole imagery coverage frequency at different resolutions. The study utilized Moran’s I index and was applied to the network of points spaced 10 km apart throughout China. Results indicated that the Moran’s I indices for all categories of imagery were significantly above zero, and all passed the significance test at the 1% level (Table 4). This demonstrated that the coverage frequencies were highly autocorrelated spatially, meaning similar coverage frequencies tended to cluster together in space.

Table 4 Global Moran’s I index of Keyhole imagery coverage frequency.
Full size table

The global Moran’s I index could not reveal the spatial clustering characteristics of specific areas, so further investigation into local clustering characteristics was needed. Figure 4 illustrates the coverage frequencies of the four types of Keyhole imagery in China, revealing the following characteristics of the regional spatial associations:

  1. 1.

    All four categories of imagery exhibited only high-high and low-low clusters, where the coverage frequency was significantly similar to surrounding areas. No high-low or low–high clusters were observed, where values significantly differed from surrounding areas. The occurrence of only high-high and low-low clusters with no high-low or low–high types suggested that the spatial distribution of imagery data across different resolutions be non-random and exhibit distinct clustering patterns. This could be attributed to the satellite monitoring areas being strongly influenced by human control and the lack of spatial continuity between different satellite monitoring areas. Although each satellite captured images continuously, the area covered by a single Keyhole satellite was limited. However, the mission area of Keyhole satellites was pre-set, resulting in the frequent monitoring of hotpot regions, while adjacent areas received significantly fewer observations. This pre-determined targeting likely contributed to the pronounced discontinuity in regional image coverage frequencies. The absence of high-low and low–high clustering types in the coverage frequencies of different kinds of Keyhole imagery indicated a lack of significant gaps in coverage frequency. Excluding the high-high and low-low clusters, satellite imagery coverage was relatively uniform across regions with low spatial variability.

  2. 2.

    In different resolution imagery, the areas classified as high-high clusters were larger than those in low-low clusters. In C0, C1, and C2, the proportions of high-high clusters were 12%, 11%, and 9%, while those of low-low clusters were 0%, 16%, and 9%, respectively. For C3, the proportion of high-high clusters was 15%, lower than that of low-low clusters at 22%.

  3. 3.

    Compared to low-low clusters, high-high clusters were more spatially concentrated. The high-high clustering areas were mainly distributed in eastern and northeastern China and along the western border. High-high areas of C2 in the southwestern part of Heilongjiang, which bordered Inner Mongolia, and in the center and northern part of Xinjiang were displayed as linear distribution. This linear distribution reflected the monitoring area of a single satellite and repeated monitoring of that specific region. This further explained the presence of only high-high and low-low clustering types as discussed earlier in (1). Low-low clusters primarily appeared around regions without imagery, such as most low-low regions in C3 and northern Xizang in C2. This indicated that the average value of imagery coverage frequency in these low-low areas was relatively low, close to the zero value of regions without any imagery.

Fig. 4
figure 4

LISA cluster map of Keyhole imagery coverage based on resolution classification (a, C0, all Keyhole images; b, C1, meter-level resolution images; c, C2, five-meter-level resolution images; d, C3, ten-meter-level resolution images).

Full size image

Characteristics of provinicial spatial distribution of imagery coverage frequency

Using provincial boundaries, the average coverage rate and spatial variability of Keyhole imagery for each province in China at different resolutions were statistically analyzed. C1 imagery exhibited the lowest coverage area proportion, yet the coverage frequency varied significantly between provinces (Fig. 5). Twenty-eight provinces had a C1 coverage area proportion greater than 50%. The average coverage frequency in areas having imagery between provinces ranged from a minimum of 1.1 times (Taiwan) to a maximum of 36.2 times (Shanghai). Fifteen provinces had more than 50% of their areas with an average coverage frequency exceeding three times. Twenty-six provinces had a coverage frequency coefficient of variation greater than 0.5. Additionally, C1 imagery also showed characteristics of concentrated coverage in hotpot areas within provinces. Although the coverage rates of C1 in Inner Mongolia and Xinjiang were only 50% and 68%, the average coverage frequencies in areas having imagery were 3.0 and 2.7 times, respectively, with high coefficients of variation of 1.7 and 1.3. This indicated that certain areas in Inner Mongolia and Xinjiang were monitored more frequently, while other regions within the province were relatively neglected.

Fig. 5
figure 5

Mean and variability of provincial coverage frequency of Keyhole imagery based on resolution classification.

Full size image

C2 imagery demonstrated high coverage frequencies in most provinces with generally low coefficients of variation. The coverage proportion in all provinces exceeded 75%. Twenty-five provinces had over 50% of their areas with coverage frequencies exceeding three times, and 11 provinces demonstrated a coefficient of variation in coverage frequency greater than 0.5. The average coverage frequency in areas having imagery between provinces ranged from a minimum of 2.0 times (Guangxi) to a maximum of 14.8 times (Shanghai). In Zhejiang, Hebei, and Jiangsu, over 30% of the area had a coverage frequency exceeding ten times.

Although the coverage rate of C3 imagery was higher than that of C1, the average coverage frequency for C3 was lower than that of C1. For C3, 22 provinces exhibited an area coverage proportion greater than 75%, while only two had coverage areas less than 50% (Chongqing and Hunan). For average coverage frequency across provinces for areas having imagery, the minimum value was 1.1 times (Chongqing), and the maximum was 9.2 times (Shanghai). Some provinces had limited coverage areas, yet the regions having imagery exhibited multiple coverages. In Sichuan and Gansu, the C3 coverage areas were 68% and 64%, respectively, while the average coverage frequency in these regions was 2.9 and 2.5 times.

For C0 imagery including C1, C2, and C3, Keyhole imagery showed a generally high area coverage rate across all provinces in China. In Xinjiang and Xizang, the keyhole imagery coverage rates were 98% and 95%, respectively, while all other provinces had a coverage rate of 100%. In all provinces, over 75% of the area had an average coverage frequency exceeding three times. The provinces where over 20% of the area had an average coverage frequency of less than two times were Xizang, Qinghai, and Xinjiang.

Purchase cost of three coverage of keyhole imagery in China

Since each resolution category included multiple types of imagery, and each type had varying orbital altitudes for different launch missions, the coverage area within each category also varied. Combining the average single-scene coverage area and the number of images of different satellite series within each category imagery as detailed in Table 1, KH-9H and KH-9L were selected as the priority supplementary imagery for C1 and C3, respectively (Table 5). The single-scene coverage area for KH-9H and KH-9L imagery was 8000 and 29,150 km2, respectively, with the cost per unit area being 6.0 and 1.0 per 1 k square kilometer (30$ per scene). Considering that the single-scene coverage areas and the number of images for KH-4A and KH-4B were closely comparable, both were chosen as the priority supplementary imagery of C2. The single-scene coverage area for KH-4A and KH-4B imagery was around 5000 km2, respectively, with the cost per unit area being 6.0 $ per square kilometer.

Table 5 Acquisition costs for triple coverage of Keyhole imagery across different resolution categories.
Full size table

For C1 meter-level imagery, the area with zero coverage across China was 3.99 million km2. The cost of acquiring three supplementary imageries with an overlap factor of 1.5 was 62 thousand dollars. The area with one coverage across China was 1.39 million km2, and the cost of acquiring two supplementary imageries was 22 thousand dollars. The area with two coverages across China was 1.26 million km2, and the cost of acquiring one supplementary imagery was 19 thousand dollars. Therefore, the total cost for triple coverages of C1 meter-level imagery across China was 130 thousand dollars. Similarly, the costs of acquiring supplementary imageries for C2 with zero, one, and two coverage areas were 14, 38, and 51 thousand dollars, respectively, with a total cost of 103 thousand dollars s. The costs of acquiring supplementary imageries for C3 with zero, one, and two coverage areas were 10, 7, and 7 thousand dollars, respectively, with a total cost of 24 thousand dollars.

Discussion

Rich data selection might influence imagery coverage characteristics

This research focused on freely downloadable Keyhole imagery, while the USGS also provided Keyhole imagery that required payment to access. According to the metadata of archived Keyhole imagery from USGS, the number of paid images within China was approximately 16 times that of the free images. The average coverage frequency of free data in the research area had reached 9.9 times. In contrast, the paid data far exceeded this free data, further demonstrating that Keyhole imagery extensively documented the landscapes during the 1960s and 1980s. The price for Keyhole imagery paid data was $30 per scene, and the cost per 1 k square kilometer for meter-level and five-meter-level imageries were 3.8 and 6.0 $, respectively. These costs were substantially lower than those for historical aerial imagery34. Therefore, Keyhole imagery could provide significant data support for extending the duration of existing studies on terrestrial landscape change and related research. On the other hand, the spatial coverage frequency and characteristics of paid imagery, which were not included in this study, might differ from those of the free imagery. Therefore, investigating the spatial distribution of paid imagery could provide valuable insight to supplement regions with limited coverage of free imagery, thereby further enhancing the application potential of Keyhole imagery.

The analysis of the spatial distribution of coverage frequencies in Keyhole imagery presented in this research did not account for the variability in imagery acquisition times, which might impact the potential application of the data. The earliest imagery dates from December 1960, while the latest is October 1984. Additionally, the period for the same series of Keyhole images could extend over several years. Thus, the extensive period could result in variations in acquisition times within the same region, potentially causing inaccuracies in earth observation, such as environmental change detection35. On the other hand, this research determined priority supplementary imagery based on the coverage area of individual scenes and associated costs. The advantage of these two metrics was that they provided many imagery options, facilitating the selection of supplementary imagery. However, these two metrics also did not account for imagery acquisition times, which might affect the temporal consistency of triple coverage across different regions. As the supplementary imagery for C1 and C3, KH-9H and Kh-9L could only be applied for research on terrestrial changes from the 1970s to the present, which was a decade later than the supplementary imagery from KH-4A or KH-4B for C2. With an average coverage frequency of 9.9 times in China, Keyhole imagery provided multiple images for the exact ___location at different time points, indicating a specific temporal span for the imagery of the exact ___location. These repeated observations could offer an opportunity to investigate the terrestrial landscape changes and related research during the 1960s and 1970s. Frank Paul36 utilized the abundant declassified satellite images (KH-4 and KH-9) covering the period 1961 to 1980 to decipher the surge history of 27 glaciers since the 1950s. Therefore, investigating the temporal spans and corresponding spatial distribution characteristics of Keyhole imagery at different resolutions could enhance the description of Keyhole imagery data in the study area and provide data support for deepening understanding of landscape pattern evolution processes.

Limitations of Keyhole imagery and future work

The potential of extensive Keyhole images has not been fully utilized due to the complexities of handling panoramic imaging geometry, film distortions, and the limited availability of metadata required for georeferencing. Since terrestrial surfaces from the 1960s and 1970s have generally undergone significant changes to today, finding sufficient ground control points for image correction represents the primary challenge in utilizing Keyhole imagery. Although most research utilizing Keyhole imagery relied on manually intensive identification of Ground Control Points (GCPs) to estimate unknown camera positions and orientation, researchers are also attempting to develop methods for automatically and rapidly identifying ground control points in Keyhole imagery. Sajid Ghuffar et al.37 designed an automated processing pipeline to generate GCPs between KH-4 scenes and modern satellite imagery of different spatial resolutions, such as the images from Planet Labs nano-satellites and Landsat-7. Amaury Dehecq et al.12 developed an automated workflow to generate Digital Elevation Models and orthophotos from scanned KH-9 mapping camera stereo images. Additionally, since Keyhole imagery is single-band images, classification methods based on particular information were difficult to apply, reducing the efficiency of Keyhole imagery utilization. Daniel et al.38 tested an integrated deep learning and object-based image analysis approach on panchromatic KH-4B to map debris-covered ice in the 1970s. Further research is needed for other land feature automated identification to facilitate the widespread application of Keyhole imagery.

Two-thirds of the earth’s surface is covered by clouds at all times39, and cloud coverage in mages from Landsat 5 and 7 could reach up to 40%40. Therefore, cloud detection is a fundamental requirement for the application of high-resolution remote sensing imagery. Since the metadata provided by USGS for Keyhole imagery does not include indicators describing cloud cover, this study did not account for cloud coverage in the spatial distribution analysis of the imagery. The presence of cloud cover might influence the spatial distribution characteristics of Keyhole imagery and further affect the acquiring cost of the triple-coverage dataset. Most current Keyhole imagery applications are based on manually selecting images without cloud coverage, while studies on cloud cover detection in Keyhole imagery have seldom been reported. Due to the lack of sufficient spectral information from Keyhole’s single-band imagery, it is difficult to obtain accurate cloud detection results using only the spectral features or separate clouds from particular bright ground objects (such as snow, coastlines, and buildings)41. The CDnet, a deep learning-based cloud detection neural network, incorporated feature pyramids and edge refinement modules to capture clear and detailed cloud boundaries42,43. This method could be applied for cloud detection in high-resolution images lacking spectral information. Since the original Keyhole imagery was divided into four parts and none had been georeferenced, merging images and georeferencing prior to cloud detection could consume significant processing resources. Utilizing Keyhole thumbnails could potentially avoid the abovementioned processes, thereby enhancing the efficiency of obtaining cloud coverage information.

Conclusions

This research calculated the coverage frequency for Keyhole imagery across different resolutions utilizing a 10 km spaced network of points across China to analyze the spatial distribution and heterogeneity of the imagery. The main conclusions of the research were as follows:

  1. 1.

    keyhole imagery of all resolution categories exhibited high area coverage rates and frequencies in China, suggesting substantial potential for widespread application. Keyhole imagery covered 99% of China’s area, with an average coverage frequency of 9.9 times. Additionally, regions with multiple image options (coverage frequency exceeding three times) comprised 88% of the research area. The coverage rates for five-meter and ten-meter resolution imagery were 95% and 76%, respectively, with coverage frequencies of 4.5 times and 3.6 times. Since the declassified Keyhole imageries were captured during the 1960s and 1970s, research into China’s territorial surface and environmental changes since the 1960s could be enabled by combining this data with other remote sensing imagery like Landsat, extending existing national or regional research by two decades.

  2. 2.

    Although Keyhole imagery has high coverage rates and frequencies, there is significant spatial heterogeneity in the distribution of coverage frequencies. The coverage areas for meter-level, five-meter-level, and ten-meter-level resolution imagery were 58%, 95%, and 76, respectively, with coefficients of variation of 1.4, 0.9, and 0.9. All these three imageries exhibited significant spatial clustering characteristics. The areas of high-high clustering were primarily Beijing-Tianjin, Heilongjiang, Xinjiang, and southeastern China, and low-low clustering areas were mainly found in the Northwest region. The diversity in resolutions intensified the variability in the spatial distribution of coverage frequencies of Keyhole imagery, and the correlations between coverage frequencies of different images were weak. Therefore, selecting suitable imagery for other research areas requires examining detailed data characteristics, such as resolution and coverage frequency.

  3. 3.

    Considering resolution, coverage area, and coverage frequency, the five-meter resolution imagery could be most suitable for national or regional scales application of Keyhole imagery. The five-meter resolution imageries include the KH-4A and KH-4B. The acquisition periods for these two imageries were from September 1967 to May 1972 and August 1963 to September 1969, respectively. These images covered 95% of China, with an average coverage frequency of 4.5 times. As the Keyhole images provided by USGS were scanned copies from the original films, the actual ground resolution might be lower than the stated resolution from USGS31. Therefore, KH-4A and KH-4B imageries with stated resolutions of 1.8 m and 2.7 m were classified as five-meter resolution images. However, the actual ground resolution was expected to be better than the five-meter level.

  4. 4.

    Combining the actual coverage area per image and the number of images, KH-9H, KH-4A or KH-4B, and KH-9L could be priority supplementary images for C1, C2, and C3 imageries. Assuming that selecting suitable images for remote sensing applications requires three times coverages, the image acquisition costs for areas not achieving three times coverage in China for C1, C2, and C3 are 103, 103, and 23 thousand dollars, respectively.

In conclusion, this research automatically reorganized the large data amount Keyhole imagery that could not analyzed manually using the metadata, to facilitate the spatial distribution and cost estimation analysis. Such automation can be further applied to investigate the temporal distribution of Keyhole imagery (e.g. revisit span), and to extend imagery characteristics beyond this case area. The high coverage rates and frequencies of Keyhole imagery demonstrates its significant potential for extending the research period on land surface using satellite images. The five-meter resolution imagery offers the most suitable coverage for China, covering almost the entire area with an average frequency of 4.5 times. However, spatial variation in coverage frequencies among different resolutions requires careful selection of imagery to meet certain research needs. The KH-9H, KH-4A or KH-4B, and KH-9L are recommended as supplementary data when the coverage frequency of meter-level, five-meter-level, and ten-meter-level imagery is insufficient.