Method for identifying and extracting urban built-up area based on VIIRS-DNB night light data

Abstract

The application discloses a kind of urban built-up area identification extraction method based on VIIRS-DNB night light data, its method includes: obtaining the multi-source remote sensing time series data of study area;Night light data is carried out logarithmic transformation processing and is combined with NDVI vegetation index to construct night light vegetation correction index;The grid in geographic information system is divided into N1 tiles, the effective pixel of effective tile is extracted and is sorted in descending order according to night light vegetation correction index to obtain brightness distribution curve, the inflection point of brightness distribution curve is identified using Kneedle algorithm and the optimal threshold of grid is obtained;In each grid, its corresponding optimal threshold is used to identify, segment and extract city, finally the urban identification extraction result of study area is obtained.The application realizes the adaptive optimal threshold selection of multi-scale processing strategy, realizes the accurate, continuous extraction of long time series urban built-up area of large-scale study area, provides high-quality basic data support for global urbanization monitoring, urban space analysis and ecological environment evaluation.

Description

Technical Field

[0001] This invention relates to the field of ecological remote sensing monitoring, and in particular to a method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data. Background Technology

[0002] Accurate mapping of urban built-up areas is a crucial foundation for analyzing global urbanization, assessing the intensity of human activity, and understanding the evolution of socio-economic spatial patterns. With the continued advancement of global urbanization, urban spatial expansion not only profoundly alters regional land use structures but also exerts a long-term and far-reaching impact on the ecological environment, resource allocation, and population migration. Therefore, constructing temporally consistent and spatially comparable urban regional data products at a global scale has become a core research issue in urban geography and remote sensing. Currently, there is no unified definition of urban regions in academia; they are often equated with city categories in impermeable surface or land cover data, leading to significant differences in urban data products under different definitions. Nighttime light (NTL) data directly reflects the intensity of human nighttime lighting and socio-economic activities, serving as a representative remote sensing indicator for characterizing urban spatial distribution and human activity features. Compared to traditional urban identification methods that rely on daytime optical imagery, nighttime light data can reveal urban spatial structure from a functional activity perspective, offering unique advantages in identifying actual human activity areas. Early studies primarily utilized data from the DMSP-OLS (Defense Meteorological Satellite Program-OLS) for urban boundary extraction, urban expansion monitoring, and light pollution research. However, this data suffers from issues such as missing cross-sensor radiometric calibration, severe pixel saturation, and low spatial resolution. Furthermore, the data time series only extends to 2013, making it insufficient for continuous monitoring of the rapid urbanization process in recent years. The release of the Visible Infrared Imaging Radiometer Suite (VIIRS) Day / Night Band (DNB) in 2012 provided a new data foundation for nighttime light research. VIIRS-DNB data significantly outperforms DMSP-OLS products in radiometric calibration accuracy, dynamic range, and spatial resolution. It effectively alleviates saturation issues in high-brightness areas, more accurately characterizes brightness gradients and low-to-medium brightness areas within cities, and greatly improves the identification capabilities of small and medium-sized cities and urban fringe areas. With the accumulation of long-term continuous data, it has become a core data source for global urban expansion research.

[0003] At the technical level, threshold-based city extraction methods are widely used for extracting urban built-up areas from nighttime light data due to their simplicity. Threshold selection is crucial to the accuracy of the extraction results. However, due to significant differences in socio-economic development levels, energy usage patterns, and background environments across different regions, a single global threshold is difficult to apply to urban mapping at large regional or even global scales. Recent research has gradually shifted towards regional adaptive thresholding methods, which address spatial heterogeneity by dividing regions and calculating locally optimal thresholds. However, how to reasonably divide regions and stably determine the optimal threshold remains a pressing technical challenge.

[0004] Meanwhile, machine learning and deep learning methods have made some progress in the application of city identification. These methods can improve classification accuracy by fusing multi-source data, but they usually rely on a large number of training samples and high computational costs, facing the dual challenges of transferability and stability in long-term global applications. Other studies have attempted to improve city identification results by fusing auxiliary information such as points of interest (POIs), statistical data, or kernel density estimation. However, such data suffers from uneven spatial coverage and insufficient temporal consistency globally, limiting their application in long-term global urbanization monitoring. It is worth noting that relying solely on nighttime light intensity is easily affected by natural surface background interference. Low-brightness built-up areas, urban-rural transition zones, and agricultural activity areas often exhibit radiation characteristics similar to cities, easily leading to misjudgments in the extraction results. Therefore, recent research has increasingly emphasized constraining the city identification process from the perspective of multi-dimensional surface features. The National Distance Vegetation Index (NDVI) can stably characterize surface vegetation cover, and urban areas generally exhibit low vegetation cover. Combining the intensity of human activity reflected by nighttime light with the ecological background information characterized by NDVI can identify urban space from two complementary dimensions: functional activity and surface cover, effectively reducing the extraction uncertainty caused by single nighttime light data.

[0005] Despite significant progress in extracting urban built-up areas based on nighttime light data, several key challenges remain: First, existing methods struggle to simultaneously address regional differences and methodological consistency on a global scale; second, most studies focus on identifying static urban boundaries in small regions, neglecting the stability and continuity of long-term urban expansion processes; and third, further exploration is needed to develop an automated and generalizable method for extracting global urban built-up areas without relying on a large number of training samples. Summary of the Invention

[0006] The purpose of this invention is to provide a method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data. Using VIIRS-DNB nighttime light data as the core, it introduces the NDVI vegetation index to constrain the surface background, combines water body masking to remove non-urban interference, constructs a nighttime vegetation correction index to enhance the difference between urban and non-urban pixels, employs a multi-scale processing strategy of block gridding and a grid identification and automatic data switching method for arid regions, and automatically determines the optimal threshold for regional adaptation based on an inflection point detection algorithm. The extraction results are then optimized through temporal and spatial post-processing to ensure the temporal continuity and spatial accuracy of urban built-up area extraction, ultimately forming a globally standardized urban built-up area data product.

[0007] The objective of this invention is achieved through the following technical solution:

[0008] A method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data, the method comprising:

[0009] S1. Obtain nighttime light data, NDVI vegetation index, and water mask data of the study area in time sequence, and perform spatiotemporal registration to obtain multi-source remote sensing time series data.

[0010] S2. Perform logarithmic transformation on the nighttime light data and combine it with the NDVI vegetation index to construct a nighttime light vegetation correction index.

[0011] S3. Construct a gridded geographic map system with a spatial reference frame. Use multi-source remote sensing time-series data to obtain the nighttime vegetation correction index at the study time and express it in the geographic map system. Divide the grid in the geographic map system into N1 tiles, extract the effective pixels of the effective tiles, and sort them in descending order according to the nighttime vegetation correction index to obtain the brightness distribution curve. Use the Kneedle algorithm to identify the inflection point of the brightness distribution curve and use it as the candidate threshold of the effective tile. Take P1 or more candidate thresholds in the grid and calculate the median as the optimal threshold of the grid. In each grid, use its corresponding optimal threshold to perform city identification, segmentation and extraction, and finally obtain the city identification and extraction results of the study area.

[0012] To better implement this invention, in method S3, invalid pixels with dark backgrounds, including oceans, deserts, or large-scale water bodies, are removed from the tiles of the grid in the geographic map system. The pixels retained after removing invalid pixels are considered valid pixels. If the number of valid pixels in a tile is below a valid threshold, the tile is determined to be an invalid tile without urban features; otherwise, it is a valid tile. Furthermore, before entering the local threshold calculation, this invention introduces an adaptive identification mechanism for arid grids. The 90th percentile of the corresponding NDVI pixels within the grid is used as a macroscopic representation of the vegetation lushness of the region. If this value is below a preset aridity threshold (0.15), the grid is determined to be an extremely arid area, indicating that the vegetation constraint of the night-light vegetation correction index is locally ineffective. For such arid grids, the extraction data source is automatically switched from the noctilucent vegetation correction index to a pre-processed VIIRS-DNB image to avoid excessive amplification of the scattered light from the sandy background by the logarithmic transformation; for regular grids that are determined to be non-arid, the noctilucent vegetation correction index is continued to be used to maintain high accuracy in distinguishing between urban and rural boundaries.

[0013] Preferably, in method S3, if the number of candidate thresholds within a grid is less than P1, an extended search strategy for spatial neighborhood grids is introduced, as follows: the optimal threshold of the grid is marked as missing, and surrounding grids are searched in order of spatial distance from the grid to the farthest until the number of candidate thresholds collected from surrounding grids reaches more than P1. The median of the more than P1 candidate thresholds is taken as the optimal threshold of the grid. If sufficient candidate thresholds from neighborhood grids cannot be collected within the spatial distance of the preset search radius, the P2% quantile of non-zero pixels in the study area or globally is used as the extreme backoff threshold, and the extreme backoff threshold is used as the optimal threshold of the grid.

[0014] Preferably, in method S3, the night-light vegetation correction index in the brightness distribution curve is sorted and ranked by pixel, and the brightness value of pixel ranking k is... Calculate the curvature of pixel k in the brightness distribution curve. Brightness value The nighttime vegetation correction index for pixel k; the pixel ranking k and brightness value. Normalization transformations are performed on each curve, and a parametric curve is constructed. The curvature expression of the parametric curve is as follows:

[0015] ,in Let k be the normalized value of the pixel ranking k, i.e., the normalized ranking. For pixels The parameter, i.e., the brightness value The normalized value, for The first derivative, for The first derivative of , For pixels Brightness change rate The first derivative, For pixels Brightness change rate The second derivative;

[0016] curvature The maximum value is taken as the inflection point, and the night-light vegetation correction index of the pixel corresponding to the inflection point is used as the candidate threshold for the tile.

[0017] Preferably, in method S1, the NDVI vegetation index source includes MODIS NDVI data, the vegetation index is calculated from MODIS NDVI data, and the vegetation index is normalized to obtain the normalized vegetation index; the nighttime light data source includes VIIRS-DNB nighttime light data, an extreme value threshold is set, and the upper limit of the VIIRS-DNB nighttime light data is truncated using the extreme value threshold; only the pixel peak value is retained for the water body mask data, and pixel-level masking is performed in the multi-source remote sensing time series data using the water body mask data with retained pixel peak values; the data in the multi-source remote sensing time series data are spatiotemporally registered and aligned according to spatial distribution and time order.

[0018] Preferably, the expression for the bioluminescent vegetation correction index is as follows:

[0019] ,in Normalized Difference Vegetation Index (NDVI) This refers to the radiance data from nighttime light data.

[0020] Preferably, in method S3, the grid scale in the geographic map system is as follows: the grid side length is a rectangular or square grid with latitude and longitude of 3 to 6°; an overlap zone of 0.1 to 0.3° is set between adjacent grids in the geographic map system; the number of pixels in the grid is at least 256×256 pixels, and an overlap zone of no less than 64 pixels is set between adjacent tiles in the grid; a multi-scale processing strategy of three levels—grid, tile, and pixel—is formed in the geographic map system, and each grid adaptively obtains the corresponding optimal threshold according to the multi-scale processing strategy.

[0021] Preferably, method S3 further includes correction processing for the time and spatial dimensions, as follows:

[0022] S1. Obtain all pixel data of the study area within year N1, perform time-series majority filtering using a sliding time window, and retain only pixels that have been identified as city pixels at least three times.

[0023] S2. Use a time overlay strategy to prioritize the retention of city pixels from the previous year and pass them on to subsequent years, and use candidate pixels to supplement data for missing pixels.

[0024] S3. Use windows to detect and merge city cells that are connected across blocks and adjacent isolated cells, while removing isolated and invalid city patches.

[0025] Preferably, in the geographic map system, the grid, tile, and pixel levels are subjected to smooth fusion processing using overlapping feathering weights during the image stitching stage.

[0026] Preferably, land cover products and population distribution data of the study area are used to form a validation dataset. The validation dataset is then used to select specific areas for accuracy assessment based on spatial distribution characteristics, temporal expansion patterns, and correlation with socioeconomic data. The accuracy assessment includes the Kappa coefficient and overall accuracy. The Kappa coefficient is expressed as follows: ,in To correctly classify the proportion of pixels, The expected consistency rate; the overall precision 0A expression is as follows: , The number of city pixels that are correctly classified. The number of non-urban pixels misclassified as cities. The number of non-urban pixels misclassified as cities. The number of non-city pixels that are correctly classified.

[0027] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0028] (1) This invention utilizes multi-source remote sensing time-series data to construct a nighttime vegetation correction index and perform candidate and optimal threshold processing strategies at three levels (grid, tile, and pixel) in a geographic map system. Within each grid, the corresponding optimal threshold is used for city identification, segmentation, and extraction, ultimately obtaining the city identification and extraction results for the study area. This enables accurate and continuous extraction of long-term urban built-up areas in large-scale study areas, providing high-quality basic data support for global urbanization process monitoring, urban spatial pattern analysis, and ecological environment impact assessment.

[0029] (2) This invention takes VIIRS-DNB nighttime light data as the core, introduces NDVI vegetation index to constrain the surface background, combines water body mask to remove non-urban interference, constructs nighttime vegetation correction index to enhance the difference between urban and non-urban pixels, and through block gridding multi-scale processing strategy and arid area grid identification and automatic data switching method, realizes automatic determination of regional adaptive optimal threshold based on inflection point detection algorithm, and then optimizes the extraction results through temporal and spatial post-processing to ensure the temporal continuity and spatial accuracy of urban built-up area extraction, and finally forms a globally standardized urban built-up area data product.

[0030] (3) Based on the hybrid strategy of inflection point detection and block gridding, this invention first estimates and filters several candidate thresholds at the tile scale, then obtains stable local thresholds at a suitable grid scale, and finally generates continuous global urban built-up area classification results through feathering and voting fusion, which reduces the spillover effect of night light signals, enhances the difference in pixel data between urban and non-urban pixels for accurate extraction, and improves the accuracy of global urban built-up area extraction. Attached Figure Description

[0031] Figure 1 This is a flowchart of the method of the present invention;

[0032] Figure 2 This is a schematic diagram illustrating the principle in the embodiment;

[0033] Figure 3 This is a schematic diagram illustrating the principles of grid division and threshold calculation in the embodiment;

[0034] Figure 4 This is an example of the extraction results for some cities as the study area in the embodiment;

[0035] Figure 5 This is a graph showing the linear fitting results of the total urban area in the example;

[0036] Figure 6 This is a comparison image showing the city extraction results from satellite imagery for some cities in the example study area. Detailed Implementation

[0037] The present invention will be further described in detail below with reference to embodiments:

[0038] Example

[0039] like Figure 1 As shown, a method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data is described, the method comprising:

[0040] S1. Obtain nighttime light data, NDVI vegetation index, and water mask data of the study area in chronological order, and perform spatiotemporal registration to obtain multi-source remote sensing time-series data. Due to the significant spatial heterogeneity of urban spectra and background conditions globally, the simple single global threshold of traditional methods is usually difficult to simultaneously and correctly identify large cities and small towns, and vegetation-covered areas and urban built-up areas. Therefore, the block processing strategy based on local optimal thresholds (i.e., the optimal threshold of the grid and the candidate threshold of the tile) is more suitable for large-scale urban extraction problems. However, how to stably and accurately determine the "optimal threshold" of each local unit in an automated process remains a key challenge for achieving global urban monitoring. This invention is based on VIIRS-DNB nighttime light composite products, combined with MODIS NDVI and global water masking. Based on a hybrid strategy of knee detection and block gridding, several candidate thresholds are first estimated and screened at the tile scale. Then, stable local thresholds are obtained at an appropriate grid scale. Finally, continuous urban built-up area classification results for the study area (e.g., a specific geographic region or globally) are generated through feathering and voting fusion. To reduce the spillover effect of nighttime light signals, the Nighttime Light Vegetation Correction Index (NVRI) is proposed and used as an input index to enhance the difference between urban and non-urban pixels, thereby improving the accuracy of global urban built-up area extraction. In this embodiment, the nighttime light data source includes VIIRS-DNB nighttime light data. An extreme value threshold is set, and the VIIRS-DNB nighttime light data is truncated to an upper limit using this threshold. Considering that there may still be extreme high values ​​in the VIIRS-DNB composite products interfering with threshold estimation, this embodiment selects values ​​greater than 250 nW·cm from the annual "average_masked" imagery. -2 ·sr -1 The upper bound of the pixel is truncated to 250 nW·cm -2 ·sr -1 To limit the observation dynamic range to 0–250 nW·cm -2 ·sr -1 The range is adjusted to reduce the impact of rare high values ​​and improve the robustness of threshold estimation; a global water mask is used to perform pixel-level masking on VIIRS imagery, and the nightlight value of water pixels is set to 0 to avoid misjudgment of city extraction caused by reflected light from coastal or surface water.

[0041] like Figure 1 , Figure 2As shown, the NDVI vegetation index in this embodiment comes from MODIS NDVI data (preferably the MOD13A1 dataset in this embodiment). The vegetation index is calculated from the MODIS NDVI data, and then normalized to obtain the normalized vegetation index. This embodiment uses the MOD13A1 dataset as an example. The vegetation index can reliably characterize the distribution of surface vegetation, making it easier to distinguish between vegetation cover and real artificial lighting sources in addition to brightness information, thereby suppressing halo spillover and false detection of non-urban land. To construct a vegetation index baseline that is comparable to the study area (e.g., a certain geographical region or the world), several years of MOD13A1 datasets were obtained on the Google Earth Engine platform. Quality band masking was applied to the original images, pixel peak values ​​were calculated according to the season, missing data was filled, and the images were cropped and stitched together in blocks to finally generate a high-resolution (500-meter resolution) vegetation index data product that is consistent year by year and spatially complete. Combining this vegetation index data product with VIIRS nighttime light emission provides strong constraints on threshold selection and post-processing, resulting in significant improvements in spatial boundary, noise suppression, and cross-regional consistency in the final urban built-up area extraction. The vegetation index calculation formula is as follows: NIR and Red represent the actual reflectance of the Earth's surface in the near-infrared and red bands, respectively, derived from MODIS NDVI data. The vegetation index was normalized to obtain the Normalized Difference Vegetation Index (NDVI). Only pixel peak values ​​were retained for water body masking data. Pixel-level masking was performed on the multi-source remote sensing time-series data using this water body masking data (using water body masking for pixel-level masking of nighttime light data effectively eliminates non-urban interference). Data in the multi-source remote sensing time-series data underwent spatiotemporal registration and alignment according to spatial distribution and temporal order.

[0042] S2. Perform logarithmic transformation on the nighttime light data and combine it with the NDVI vegetation index to construct a nighttime light vegetation correction index; the expression for the nighttime light vegetation correction index is as follows:

[0043] ,in Normalized Difference Vegetation Index (NDVI) This refers to radiance data from nighttime light data. This is the nighttime light data after logarithmic transformation. This embodiment introduces the Normalized Difference Vegetation Index (NDVI) to amplify the differences between urban-rural boundaries and between high-vegetation patches and built-up area pixels within cities, thus facilitating more accurate extraction of urban boundaries. Logarithmic transformation (i.e., logarithmic transformation...) , To effectively enhance the distinction between built-up and unbuilt-up areas using a natural logarithm function with base e, this embodiment employs the Night Visible Vegetation Correction Index (NVRI). By combining the NDVI with the logarithmically transformed night vision signal, the accuracy of extracting urban built-up areas from night vision data is improved. When pixel vegetation cover is high... The lower NVRI weakens the normalized vegetation index (NVRI) of that pixel; although NVRI shows excellent urban-rural differentiation in most vegetated areas, in extremely arid or large-scale desert regions (such as the Middle East and North Africa), the surface NDVI is generally close to or even below 0. In this case, the formula... When the multiplier approaches 1, the physical constraint of vegetation becomes essentially ineffective, and NVRI degenerates into a pure logarithmic image of nighttime light. Combined with the high reflectivity of bare sandy soil in extremely arid areas and the interference from dense oil and gas field fire points, a severe "desert light spillover" effect can easily occur, leading to a significant overestimation of urban areas. To address this, this invention introduces an adaptive identification and dynamic data source switching mechanism for arid grids in the subsequent threshold calculation module. 1-NDVI approaching 1 indicates that the surface is almost devoid of vegetation (bare land, built-up land, paved roads); ln(1+VIIRS) approaching 1 indicates that the nighttime light intensity is at a moderate to high level, sufficient to be identified as an artificial light source; ln(1+VIIRS) approaching 0 indicates areas without artificial lighting, such as bare land and desert.

[0044] S3. Construct a gridded geographic map system that includes a spatial reference frame. Utilize multi-source remote sensing time-series data to obtain the nighttime bioluminescent vegetation correction index for the study time and express it in the geographic map system, such as... Figure 1 , Figure 2As shown, the grid in the geographic map system is divided into N1 tiles. In some embodiments, the grid scale in the geographic map system is as follows: the grid side length is a rectangular or square grid with a latitude and longitude of 3 to 6° (in this embodiment, a regular latitude and longitude grid of 4°×4° is used to characterize the regional differences in the relationship between nighttime light and vegetation in different geographic units. This grid scale can take into account global spatial heterogeneity and ensure sufficient sample size of urban nighttime light pixels in the region so as to accurately calculate the optimal grid threshold); an overlap zone of 0.1 to 0.3° latitude and longitude is set between adjacent grids in the geographic map system. In this embodiment, the overlap zone is 0.1°, which can effectively suppress the seam effect caused by the threshold abrupt change at the grid boundary. Setting a 0.1° overlap zone between adjacent grids ensures the spatial continuity and smooth transition of the threshold; before entering the local threshold calculation, this invention introduces an adaptive recognition mechanism for drought grids. The 90th percentile of the corresponding NDVI pixels within the grid is used as a macroscopic representation of the vegetation density in the region. If this value is lower than the preset drought threshold (0.15), the grid is determined to be an extremely arid area, indicating that the vegetation constraint of the night-photometric vegetation correction index is locally ineffective. For such arid grids, the data source is automatically switched from the night-photometric vegetation correction index back to a pre-processed VIIRS-DNB image to avoid excessive amplification of background scattered light from sandy areas by logarithmic transformation. For regular grids determined to be non-arid, the night-photometric vegetation correction index is continued to maintain high-precision distinction between urban and rural boundaries. After determining the optimal data source (NVRI or VIIRS-DNB) for each 4°×4° grid, the grid is further divided into 256×256 pixel local tiles (in this embodiment, an example is an overlap zone of 64 adjacent tiles). The multi-scale processing strategy preserves spatial details while improving the robustness of threshold estimation to noise and outliers. A multi-scale processing strategy is formed at three levels: grid, tile, and pixel, within the geographic map system. Each grid adaptively obtains its corresponding optimal threshold according to the multi-scale processing strategy. In some embodiments, invalid pixels with dark backgrounds, including oceans, deserts, or large-scale water bodies, are removed from the tiles of the grid in the geographic map system. Invalid pixels can be determined by the following methods: NDVI vegetation index equal to or close to 0, or pixels corresponding to the pixel-level mask of water body mask data, or a combination of other data sources. The pixels retained after removing invalid pixels are considered valid pixels. If the number of valid pixels in a tile is less than the valid threshold (e.g., less than 100), the tile is considered an invalid tile without urban features; otherwise, it is a valid tile. This invention extracts the effective pixels of effective tiles and sorts them in descending order according to the night-light vegetation correction index to obtain brightness distribution curves, such as... Figures 2 to 3As shown, the Kneedle algorithm is used to identify the inflection point of the brightness distribution curve and use it as a candidate threshold for effective tiles. The median of P1 or more candidate thresholds (P1 is 3 in this embodiment) within the grid is taken as the optimal threshold of the grid. Within each grid, the corresponding optimal threshold is used for city identification, segmentation and extraction, and finally the city identification and extraction results of the study area are obtained.

[0045] In some embodiments, if the number of candidate thresholds within a grid is less than P1 (P1 is 3 in this embodiment), an extended search strategy using spatial neighborhood grids is introduced. The method is as follows: the optimal threshold of the grid is marked as missing. The surrounding grids are searched in ascending order of spatial distance, centered on the grid, until the accumulated candidate thresholds of the surrounding grids reach P1 or more (P1 is 3 in this embodiment). The median of the searched P1 or more candidate thresholds is taken as the optimal threshold of the grid. If sufficient candidate thresholds from neighboring grids cannot be collected within a preset search radius, the P2% quantile of the non-zero pixels in the study area or globally (P2 is 90, i.e., the 90th quantile in this embodiment) is used as the extreme backoff threshold, which is then used as the optimal threshold for the grid. In the tile-level processing stage, only valid pixels are retained, and dark background-dominated areas, such as oceans and large-scale deserts, are removed. Areas judged as invalid tiles are skipped to avoid interfering with the threshold estimation results. In the threshold selection stage, the NVRI values ​​within each tile are first statistically analyzed to determine candidate thresholds: effective pixels within the tile are sorted in descending order of NVRI values, and the Kneedle algorithm is used to automatically identify the inflection points of the brightness distribution curve, which are then used as candidate thresholds for that tile. Simultaneously, the 90th percentile of global non-zero pixels is used as the lower limit of the threshold to filter out invalid candidate thresholds generated by tiles dominated by non-urban nighttime light pixels. Subsequently, the candidate thresholds for each tile are summarized and statistically analyzed at the grid scale to determine the final threshold for each grid. During the image stitching stage, overlap feathering weights are used to achieve smooth fusion. Existing technologies typically rely on high-resolution land cover products or annual statistical data to estimate this threshold; however, such data preparation is costly and updates are slow, making it difficult to meet the needs of large-scale, long-term automated monitoring. Other methods use contour lines or empirical rules for threshold statistics, but these often require manual setting of empirical parameters and lack generalization ability and repeatability. To overcome these limitations, the method of this invention starts from statistical features such as pixel brightness and vegetation distribution, automatically identifies the distribution inflection point within each spatial unit, and generates a region-adaptive threshold accordingly. This strategy retains the advantages of thresholding methods, such as simplicity and interpretability, and improves cross-regional consistency and responsiveness to urban dynamics through local adaptation and spatial smoothing. Therefore, it is more suitable for extracting urban built-up areas on a global scale and over long time series.

[0046] In some embodiments, the intensity of luminous pixels exhibits a long-tailed spatial distribution, with a small number of strong pixels coexisting with a large number of weak pixels. For example... Figure 2 , Figure 3 As shown, the Kneedle algorithm uses a geometric method to find the inflection point of the sorting curve, thereby automatically providing an optimal threshold that is highly consistent with the semantics of the urban background. Specifically, it treats the sorting curve as a parameterized function, calculates its curvature, and takes the pixel value corresponding to the point of maximum curvature as the threshold (this method can adaptively identify the "bright-dark" boundary without manually setting a local threshold). The specific method is as follows: The night-light vegetation correction index in the brightness distribution curve is sorted by pixel, and the brightness value of pixel rank k is... Calculate the curvature of pixel k in the brightness distribution curve. Brightness value The nighttime vegetation correction index for pixel k; the pixel ranking k and brightness value. Normalization transformations are performed on each curve, and a parametric curve is constructed. The curvature expression of the parametric curve is as follows:

[0047] ,in Let k be the normalized value of the pixel ranking k, i.e., the normalized ranking. N is the total number of valid pixels; For pixels The parameter, i.e., the brightness value The normalized value, for The first derivative, for The first derivative of , For pixels Brightness change rate The first derivative, For pixels Brightness change rate The second derivative of .

[0048] curvature The maximum value is taken as the inflection point, i.e., take... The night-luminescent vegetation correction index of the pixel corresponding to the inflection point is used as the candidate threshold for the tile.

[0049] This embodiment uses a square grid with side lengths of 4°×4° latitude and longitude as an example. The median threshold of all valid tiles within each 4°×4° grid is taken as the initial optimal threshold for that grid. If the number of valid tiles with thresholds in any grid is insufficient, the threshold for that grid is marked as missing. To alleviate the threshold instability caused by local data sparsity, a spatial neighborhood-based extended search strategy is introduced. Specifically, for each target grid, surrounding grids are searched in order of spatial distance from nearest to farthest, accumulating valid thresholds until at least three candidate thresholds are collected. The median of these candidate thresholds is taken as the final threshold for the corresponding grid. If sufficient neighborhood thresholds cannot be collected within a preset search radius, the 95th percentile of global non-zero pixels is used as the extreme backoff threshold to avoid missing potential city pixels. This multiple backoff mechanism preserves the local adaptability of the thresholds while improving the spatial consistency and stability of the thresholds.

[0050] In some embodiments, such as Figure 2 As shown, this invention also includes correction processing for the time and spatial dimensions (to obtain reasonable and stable urban built-up area extraction results, correction processing is performed on the initial binarized extraction results; utilizing the temporal information of the nighttime light time series is key to improving the accuracy and consistency of the urban built-up area extraction results), the method is as follows:

[0051] S1. Obtain all pixel data for the study area within year N1 (5 years in the example), and perform temporal majority filtering using a sliding time window, retaining only pixels that have been identified as city pixels at least three times. Within a continuous 5-year sliding window, only pixels identified as city pixels at least 3 times are retained, effectively suppressing single-year noise and transient anomalous observations in nighttime light data.

[0052] S2. Using a time overlay strategy, city pixels from the previous year are prioritized for retention and passed on to subsequent years. Missing pixels are supplemented with candidate pixels to ensure the temporal continuity and spatial integrity of the extraction results.

[0053] S3. Using window detection, connectable city pixels across blocks and adjacent isolated pixels are merged to eliminate seam artifacts. Simultaneously, isolated, ineffective city patches are removed (e.g., isolated connected components with an area less than 8 pixels are removed; these patches are small and considered noise patches rather than effective urban built-up area patches). This embodiment's extraction results maintain temporal consistency and optimize spatial boundary rationality through local fusion and small patch filtering, improving the overall robustness of the results. The final result is the annual extraction of urban built-up areas based on VIIRS nighttime light imagery. See [link to relevant documentation]. Figure 4 The horizontal bar for the Nighttime Light Visible Vegetation Correction Index (NDVI) represents the extraction results of urban built-up area identification and segmentation based on the NDVI. Figure 4This embodiment also provides a comparison chart of the Normalized Difference Vegetation Index (NDVI), the original nighttime light data (VIIRS), the population distribution data (LandScan), and the results extracted by this invention. Compared to the original VIIRS, the Nighttime Light Correction Index (NDVI) filters out weak light noise in non-urban areas (such as remote rural areas and natural areas), making urban boundaries clearer and more continuous; the urban outlines of Beijing and London are more regular, and the coastal / lakeside urban forms of Chicago and Sydney are more prominent. See also Figure 4 This embodiment verifies the reliability of the Nighttime Light Visible Vegetation Correction Index (NDVI) in extracting urban built-up areas: population agglomeration areas highly match the urban areas after light-visible vegetation correction; Beijing exhibits a clear single-center high-density core area; London and Sydney show multi-center population distribution; and Chicago shows population agglomeration along the lake. Preferably, in the geographic map system, the grid, tile, and pixel levels are smoothed and fused using overlap feathering weights during the image stitching stage. This invention utilizes the coefficient of determination R... 2 The root mean square error (RMSE) was used to evaluate the consistency of urban area, and the results are as follows: Figure 5 As shown. This embodiment can use urban built-up area datasets generated from high-resolution remote sensing data (such as MODIS MCD12Q1, MGUP) and LandScan population distribution data to conduct verification analysis on the extraction results, calculate accuracy evaluation indicators such as overall accuracy (OA) and Kappa coefficient, and comprehensively evaluate the practical applicability of this method.

[0054] This embodiment uses cities on major continents around the world as the study area. See [link / reference] Figure 6 This embodiment uses the Night Light Vegetation Correction Index (NDVI) as an example to illustrate the city extraction results (specifically the NTL Urban column) for urban built-up area identification and segmentation in satellite imagery of some cities as the study area. Based on nighttime light brightness, this invention can effectively distinguish between cities with human activity and non-cities without human activity. Therefore, in the four cities, the red areas are strictly locked within the core built-up area, filtering out a large amount of peripheral noise identified by MODIS (based on the NDVI principle, it is easy to identify bare land, grassland, and even some farmland around the city as cities, resulting in unclear boundaries of the built-up area).

[0055] In some embodiments, to verify the reliability of the urban product generated by this invention, land cover product and population distribution data are used to conduct multi-faceted verification of the extraction results. The performance of the urban product in terms of spatial structure representation and socio-economic correlation is evaluated from dimensions such as spatial distribution characteristics, temporal expansion patterns, and correlation with socio-economic data. Land cover product and population distribution data of the study area are obtained to form a verification dataset. Using this dataset, specific areas are selected for accuracy evaluation based on spatial distribution characteristics, temporal expansion patterns, and correlation with socio-economic data. Accuracy evaluation includes the Kappa coefficient and overall accuracy. The Kappa coefficient expression is as follows: ,in To correctly classify the proportion of pixels, The expected consistency rate, which is the sum of the products of the actual number of pixels in each category and the predicted number of pixels in the corresponding category, divided by the square of the total number of pixels, is expressed as: The overall precision 0A is expressed as follows: , The number of city pixels that are correctly classified (true positives). The number of non-urban pixels that were misclassified as cities (false positives). The number of non-urban pixels that were misclassified as cities (false positives). The number of non-urban pixels that are correctly classified (true negatives).

[0056] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data, characterized in that: The methods include: S1. Obtain nighttime light data, NDVI vegetation index, and water mask data of the study area in time sequence, and perform spatiotemporal registration to obtain multi-source remote sensing time series data. S2. Perform logarithmic transformation on the nighttime light data and combine it with the NDVI vegetation index to construct a nighttime light vegetation correction index. S3. Construct a gridded geographic map system with a spatial reference frame. Use multi-source remote sensing time-series data to obtain the nighttime vegetation correction index at the study time and express it in the geographic map system. Divide the grid in the geographic map system into N1 tiles, extract the effective pixels of the effective tiles, and sort them in descending order according to the nighttime vegetation correction index to obtain the brightness distribution curve. Use the Kneedle algorithm to identify the inflection point of the brightness distribution curve and use it as the candidate threshold of the effective tile. Take P1 or more candidate thresholds in the grid and calculate the median as the optimal threshold of the grid. In each grid, use its corresponding optimal threshold to perform city identification, segmentation and extraction, and finally obtain the city identification and extraction results of the study area.

2. The method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data according to claim 1, characterized in that: In method S3, invalid pixels with dark backgrounds, including oceans, deserts, or large-scale water bodies, are removed from the tiles of the grid in the geographic map system. The pixels retained after removing invalid pixels are valid pixels. If the number of valid pixels in a tile is lower than the valid threshold, the tile is determined to be an invalid tile without urban features; otherwise, it is a valid tile.

3. The method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data according to claim 1, characterized in that: In method S3, if the number of candidate thresholds within a grid is less than P1, an extended search strategy using spatial neighborhood grids is introduced, as follows: the optimal threshold of the grid is marked as missing, and surrounding grids are searched in order of spatial distance from the grid to the farthest until the number of candidate thresholds collected from surrounding grids reaches more than P1. The median of the more than P1 candidate thresholds is taken as the optimal threshold of the grid. If sufficient candidate thresholds from neighborhood grids cannot be collected within the spatial distance of the preset search radius, the P2% quantile of non-zero pixels in the study area or globally is used as the extreme backoff threshold, and the extreme backoff threshold is taken as the optimal threshold of the grid.

4. The method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data according to claim 1, characterized in that: In method S3, the night-light vegetation correction index in the brightness distribution curve is ranked by pixel, and the brightness value of pixel rank k is... Calculate the curvature of pixel k in the brightness distribution curve. Brightness value The nighttime vegetation correction index for pixel k; the pixel ranking k and brightness value. Normalization transformations are performed on each curve, and a parametric curve is constructed. The curvature expression of the parametric curve is as follows: ,in Let k be the normalized value of the pixel ranking k, i.e., the normalized ranking. For pixels The parameter, i.e., the brightness value The normalized value, for The first derivative, for The first derivative of , For pixels Brightness change rate The first derivative, For pixels Brightness change rate The second derivative; Curvature The maximum value is taken as the inflection point, and the night-light vegetation correction index of the pixel corresponding to the inflection point is used as the candidate threshold for the tile.

5. The method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data according to claim 1, characterized in that: In method S1, the NDVI vegetation index source includes MODIS NDVI data. The vegetation index is calculated from the MODIS NDVI data, and the vegetation index is normalized to obtain the normalized vegetation index. The nighttime light data source includes VIIRS-DNB nighttime light data. An extreme value threshold is set, and the upper limit of the VIIRS-DNB nighttime light data is truncated using the extreme value threshold. Only the pixel peak value is retained in the water body mask data. Pixel-level masking is performed on the water body mask data with retained pixel peak values ​​in multi-source remote sensing time series data. Data in multi-source remote sensing time series data are processed for spatiotemporal registration and alignment according to spatial distribution and temporal order.

6. The method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data according to claim 1, characterized in that: The expression for the bioluminescent vegetation correction index is as follows: ,in Normalized Difference Vegetation Index (NDVI) This refers to the radiance data from nighttime light data.

7. The method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data according to claim 1, characterized in that: In method S3, the grid scale in the geographic map system is as follows: the grid side length is a rectangular or square grid with latitude and longitude of 3 to 6°; an overlap zone of 0.1 to 0.3° is set between adjacent grids in the geographic map system; the number of pixels in the grid is at least 256×256 pixels, and an overlap zone of no less than 64 pixels is set between adjacent tiles in the grid; a multi-scale processing strategy of three levels, grid, tile, and pixel, is formed in the geographic map system, and each grid adaptively obtains the corresponding optimal threshold according to the multi-scale processing strategy.

8. The method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data according to claim 1 or 7, characterized in that: Method S3 also includes corrections for the time and space dimensions, as follows: S1. Obtain all pixel data of the study area within year N1, perform time-series majority filtering using a sliding time window, and retain only pixels that have been identified as city pixels at least three times. S2. Use a time overlay strategy to prioritize the retention of city pixels from the previous year and pass them on to subsequent years, and use candidate pixels to supplement data for missing pixels. S3. Use windows to detect and merge city cells that are connected across blocks and adjacent isolated cells, while removing isolated and invalid city patches.

9. The method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data according to claim 7, characterized in that: In the geographic mapping system, the grid, tile, and pixel levels are smoothly blended during the image stitching stage using overlapping feathering weights.

10. The method for identifying and extracting urban built-up areas based on VIIRS-DNB nighttime light data according to claim 1, characterized in that: Land cover products and population distribution data for the study area were collected to form a validation dataset. Using this dataset, specific regions were selected for accuracy assessment based on spatial distribution characteristics, temporal expansion patterns, and correlation with socioeconomic data. Accuracy assessment included the Kappa coefficient and overall accuracy. The Kappa coefficient expression is as follows: ,in To correctly classify the proportion of pixels, The expected consistency rate; the overall precision 0A expression is as follows: , The number of city pixels that are correctly classified. The number of non-urban pixels misclassified as cities. The number of non-urban pixels misclassified as cities. The number of non-city pixels that are correctly classified.

Images