A method for generating a sequence of urban active areas taking into account the temporal and spatial interference of night light remote sensing

Abstract

The present invention discloses a method for generating a sequence of urban vitality zones taking into account the temporal and spatial interference of night light remote sensing, comprising the following steps: obtaining the annual time series of night light remote sensing data and impermeable layer data, calculating the proportion of impermeable layers in a single pixel of the night light remote sensing data; establishing a calculation formula for the overflow brightness of eight neighborhood pixels to the central pixel, and removing the overflow brightness of the neighborhood pixels to the central pixel; calculating the brightness difference histogram of two adjacent periods of night light remote sensing data, distinguishing between fluctuation points and normal change points; performing linear fitting on continuous fluctuation points, and updating the brightness values ​​of the remaining fluctuation points; setting a threshold T, and obtaining a time series of the range of urban vitality zones. The present invention overcomes the boundary deviation and continuity problems caused by the previous extraction with a single threshold, removes the errors caused by the overflow effect and the interannual fluctuation effect, and realizes the extraction of urban vitality zones with higher precision, which helps to accurately understand the urban spatial structure and development process, and provides a decision-making basis for urban sustainable development management.

Description

Technical Field

[0001] The present invention relates to the field of satellite remote sensing technology, and in particular to a method for generating a sequence of urban active areas taking into account the spatiotemporal interference of night light remote sensing. Background Art

[0002] Urban vitality zones refer to areas in cities with frequent economic activities, dense populations, developed transportation, and prosperous commerce. They are important basic information for studying urban spatial structure and are also the basis for urban planning and policy making. Clearly defining the spatial scope of urban vitality zones is an important task. Different from traditional optical remote sensing, night light remote sensing uses satellite sensors to capture light radiation information on the earth's surface at night. It has the unique advantage of directly reflecting human activities and can effectively reflect the activity intensity, spatial distribution, and temporal variation characteristics of urban vitality zones. Using night light remote sensing data to scientifically delineate urban vitality zones has broad application prospects.

[0003] In the previous work of extracting urban vitality zones through night light remote sensing data, urban vitality zones were mainly identified through threshold segmentation and texture feature analysis. Existing methods focus on obtaining areas with high brightness and consistent texture, but ignore the impact of the characteristics of night light data on the extraction results. Night light remote sensing data has light spillover effects in the spatial dimension and interannual brightness fluctuation effects in the temporal dimension. If the spatial and temporal interference is not removed, the scope of the urban vitality zone obtained is too large and its continuity is difficult to guarantee. Summary of the invention

[0004] The present invention provides a method for generating a sequence of urban active areas taking into account the spatiotemporal interference of night light remote sensing, thereby solving the above-mentioned problems of the prior art.

[0005] A method for generating a sequence of urban vitality zones taking into account the spatiotemporal interference of night light remote sensing comprises the following steps:

[0006] S1, obtain the annual time series of night light remote sensing data and impervious layer data, align the spatial range of single-scene night light remote sensing data and impervious layer data, and calculate the proportion of impervious layer in a single pixel of night light remote sensing data;

[0007] S2, establishing a calculation formula for the overflow brightness of the eight-neighborhood pixels to the central pixel, screening the central pixel that is easily affected by the eight-neighborhood pixels, and removing the overflow brightness of the eight-neighborhood pixels to the central pixel according to the overflow brightness calculation formula;

[0008] S3, calculate the brightness difference histogram of two adjacent periods of night light remote sensing data, fit the histogram with Laplace distribution, calculate the cumulative probability density to determine the brightness fluctuation range, and distinguish the fluctuation point and the normal change point;

[0009] S4, performing linear fitting on the continuous fluctuation points, removing the fluctuation points that deviate too much from the fitting straight line, performing linear fitting again, and updating the brightness values ​​of the remaining fluctuation points;

[0010] S5, set a threshold T, regard pixels greater than T as urban vitality areas, extract the night light time series data year by year, and obtain the time series of the urban vitality area range.

[0011] In some specific embodiments, the specific steps of step S1 include:

[0012] S11, obtain annual synthetic night light remote sensing data for more than 5 years, and obtain annual time series of impermeable layer data. The pixel size of night light remote sensing data needs to be more than 10 times the pixel size of impermeable layer data;

[0013] S12, cutting and splicing the night light remote sensing data and the impervious layer data to generate a single image with the same starting point and a spatial range of 5°×5°, and performing subsequent processing based on the single image;

[0014] S13, for single-scene night light remote sensing data and impermeable layer data of the same year, count the area of ​​the impermeable layer within a single night light pixel, and calculate the area ratio of the impermeable layer within a single night light pixel.

[0015] In some specific embodiments, the specific steps in step S2 include:

[0016] S21, based on the proportion of the impermeable layer of the luminous pixel, screen the pixels that do not emit light or emit a small amount of light, and whose brightness is mainly affected by the surrounding pixel images. The screening condition is that the proportion of the impermeable layer is less than 5%. The brightness L of the pixel C that meets the screening condition is also the brightness overflow value L of the surrounding pixels. N ;

[0017] S22, for a single pixel c that is easily affected by surrounding pixels, take its eight neighboring pixels a0, a1, a2, a3, b0, b1, b2, b3, and divide the brightness of the surrounding pixels into A and B according to the different distances from the surrounding pixels to the pixel c, where L i is the light brightness of pixel i;

[0018] S23, add a weight factor q to the surrounding pixel brightness B, q∈[0,1], set the sum of the brightness values ​​of the pixels around pixel c S to S=A+qB, increase q from 0 to 0.05, and calculate the surrounding pixel brightness overflow value L of pixel c N The Pearson correlation coefficient of the surrounding pixel brightness sum S is obtained, the weight factor q when the maximum correlation coefficient is obtained, and the surrounding pixel brightness sum S is calculated;

[0019] S24, using the least squares method to calculate the surrounding brightness sum S and the surrounding pixel brightness overflow value L N Perform linear fitting and calculate k and b, which are expressed as L N =kS+b;

[0020] S25, screening pixels that simultaneously meet the conditions that the pixel brightness value L is less than 10, the proportion of the pixel impermeable layer is less than 20%, and the pixel is not the local brightness maximum value, that is, the pixels that meet 8L-AB<4;

[0021] S26, for the pixels that meet the conditions in step S25, calculate the sum of the surrounding pixels and A, B, and substitute them into the formula in step S24 to calculate the overflow brightness L of the surrounding pixels N , the actual brightness value of the pixel itself L S LL N To obtain.

[0022] In some specific embodiments, the specific steps of step S3 include:

[0023] S31, selecting pixels whose impermeable layer difference ratio between two consecutive years is less than 10%, and calculating the light brightness difference between the two consecutive years;

[0024] S32, using Laplace distribution to define the probability density function of the brightness difference ΔL, then the probability density function of the brightness difference ΔL is: The maximum likelihood estimation method is used to calculate the parameters, and the calculation formula of the parameter μ is: The calculation formula of parameter λ is Where N is the number of pixels, ΔLi is the brightness difference of the i-th pixel;

[0025] S33, calculate t = ln20 × λ, so that the cumulative probability density in the interval [μ-t, μ+t] is equal to 0.95;

[0026] S34, judging whether the brightness change is a fluctuation according to the following rules: if the difference in light brightness of a pixel in two consecutive years is within the interval [μ-t,μ+t], the brightness values ​​of the pixel in two years are marked as a fluctuation change, otherwise they are marked as a normal brightness change.

[0027] In some specific embodiments, the specific steps in step S4 include:

[0028] S41, for the consecutive years marked as fluctuating changes, the least square method is used to perform linear fitting on the brightness values ​​of each year, and the difference between the fitted value and the actual value and the average value a of the difference are calculated. If the difference between the fitted value and the actual value is greater than the average value a, the brightness of the year is regarded as a normal brightness change point, otherwise it is confirmed as a fluctuation point;

[0029] S42, for the consecutive years confirmed to be fluctuating, the least square method is used to perform linear fitting again, and the annual brightness is set as the linear fitting value.

[0030] The beneficial effects of the present invention are:

[0031] The present invention discloses a method for generating a sequence of urban vitality zones taking into account the temporal and spatial interference of night light remote sensing, comprising the following steps: obtaining the annual time series of night light remote sensing data and impermeable layer data, aligning the spatial range of single-scene night light remote sensing data and impermeable layer data, and calculating the proportion of impermeable layers in a single pixel of night light remote sensing data; establishing a calculation formula for the overflow brightness of eight neighborhood pixels to the central pixel, screening the central pixel that is easily affected by the neighborhood pixels, and removing the overflow brightness of the neighborhood pixels to the central pixel according to the overflow brightness calculation formula. Calculate the brightness difference histogram of two adjacent periods of night light remote sensing data, fit the histogram with Laplace distribution, calculate the cumulative probability density to determine the brightness fluctuation range, and distinguish between fluctuation points and normal change points; perform linear fitting on continuous fluctuation points, remove the fluctuation points that deviate too much from the fitting straight line, and perform linear fitting again, and update the brightness values ​​of the remaining fluctuation points; set a threshold T, regard pixels greater than T as urban vitality zones, extract the night light time series data year by year, and obtain the range time series of urban vitality zones. The present invention utilizes the influence of impermeable layers on night light remote sensing data, takes into account the errors caused by spillover effects and interannual fluctuation effects, overcomes the boundary deviation and continuity problems caused by the previous use of only a single threshold for extraction, removes the errors caused by spillover effects and interannual fluctuation effects, and achieves higher-precision extraction of urban vitality areas, which helps to accurately understand the urban spatial structure and development process, and provide a decision-making basis for urban sustainable development, urban planning, and urban management. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0033] Figure 1 It is a flow chart of a method for generating a sequence of urban active areas taking into account the spatiotemporal interference of night light remote sensing according to the present invention;

[0034] Figure 2 is a schematic diagram of the spillover effect of the present invention;

[0035] Figure 3 is a comparison diagram of the urban vitality zone range extraction results before and after removing the spillover effect of the present invention;

[0036] Figure 4is a histogram of brightness difference distribution between two adjacent years of the present invention;

[0037] Figure 5 It is a result graph of time series data of the scope of the urban vitality zone of the present invention;

[0038] Figure 6 is a diagram showing changes in the range of urban vitality zones from 2013 to 2023 according to an embodiment of the present invention;

[0039] Figure 7 is a method flow chart of another embodiment of the present invention. DETAILED DESCRIPTION

[0040] In order to make the purpose, technical solution and advantages of the present invention clearer, the present invention is further described in detail below in conjunction with the accompanying drawings. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention. It should be noted that the terms "including" and "having" and any variations thereof are intended to cover non-exclusive inclusions, so that the processes, methods, systems, products or devices of a series of units are not necessarily limited to those units, but may include other units that are not clearly listed or inherent to these processes, methods, products or devices.

[0041] Reference Figure 1 , Figure 2 , Figure 3 , Figure 4 , Figure 5 , Figure 6 and Figure 7 The method for generating a sequence of urban vitality zones taking into account the temporal and spatial interference of night light remote sensing is shown, comprising the following steps:

[0042] S1, obtain the annual time series of night light remote sensing data and impervious layer data, align the spatial range of single-scene night light remote sensing data and impervious layer data, and calculate the proportion of impervious layer in a single pixel of night light remote sensing data.

[0043] S2, establishing a calculation formula for the overflow brightness of the eight-neighborhood pixels to the central pixel, screening the central pixel that is easily affected by the eight-neighborhood pixels, and removing the overflow brightness of the eight-neighborhood pixels to the central pixel according to the calculation formula for the overflow brightness.

[0044] S3, calculate the brightness difference histogram of two adjacent periods of night light remote sensing data, fit the histogram with Laplace distribution, calculate the cumulative probability density to determine the brightness fluctuation range, and distinguish the fluctuation point and the normal change point.

[0045] S4, performing linear fitting on the continuous fluctuation points, removing the fluctuation points that deviate too much from the fitting straight line, performing linear fitting again, and updating the brightness values ​​of the remaining fluctuation points.

[0046] S5, set a threshold T, regard pixels greater than T as urban vitality areas, extract the night light time series data year by year, and obtain the time series of the urban vitality area range.

[0047] In this embodiment, Figure 6 It shows the changes in the scope of urban vitality zones from 2013 to 2023. Different colors represent the scope of newly added urban vitality zones in the corresponding years, intuitively expressing the expansion trend of the city.

[0048] It should be noted that in previous work on extracting urban vitality zones through night light remote sensing data, urban vitality zones were mainly identified through threshold segmentation and texture feature analysis. Existing methods focus on obtaining areas with high brightness and consistent texture, but ignore the impact of the characteristics of night light data on the extraction results. Night light remote sensing data has light spillover effects in the spatial dimension and interannual brightness fluctuation effects in the temporal dimension. If the spatial and temporal interference is not removed, the scope of the urban vitality zone obtained is too large and its continuity is difficult to guarantee.

[0049] In some specific embodiments, the specific steps of step S1 include:

[0050] S11, obtain annual synthetic night light remote sensing data for more than 5 years, and obtain annual time series of impermeable layer data. The pixel size of night light remote sensing data needs to be more than 10 times the pixel size of impermeable layer data;

[0051] S12, cutting and splicing the night light remote sensing data and the impervious layer data to generate a single image with the same starting point and a spatial range of 5°×5°, and performing subsequent processing based on the single image;

[0052] S13, for single-scene night light remote sensing data and impermeable layer data of the same year, count the area of ​​the impermeable layer within a single night light pixel, and calculate the area ratio of the impermeable layer within a single night light pixel.

[0053] In some specific embodiments, the specific steps in step S2 include:

[0054] S21, based on the proportion of the impermeable layer of the luminous pixel, screen the pixels that do not emit light or emit a small amount of light, and whose brightness is mainly affected by the surrounding pixel images. The screening condition is that the proportion of the impermeable layer is less than 5%. The brightness L of the pixel C that meets the screening condition is also the brightness overflow value L of the surrounding pixels. N ;

[0055] S22, for a single pixel c that is easily affected by surrounding pixels, take its eight neighboring pixels a0, a1, a2, a3, b0, b1, b2, b3, and divide the brightness of the surrounding pixels into A and B according to the different distances from the surrounding pixels to the pixel c, where L iis the light brightness of pixel i;

[0056] S23, add a weight factor q to the surrounding pixel brightness B, q∈[0,1], set the sum of the brightness values ​​of the pixels around pixel c S to S=A+qB, increase q from 0 to 0.05, and calculate the surrounding pixel brightness overflow value L of pixel c N The Pearson correlation coefficient of the surrounding pixel brightness sum S is obtained, the weight factor q when the maximum correlation coefficient is obtained, and the surrounding pixel brightness sum S is calculated;

[0057] S24, using the least squares method to calculate the surrounding brightness sum S and the surrounding pixel brightness overflow value L N Perform linear fitting and calculate k and b, which are expressed as L N =kS+b;

[0058] S25, screening must simultaneously meet the requirements that the pixel brightness value L is less than 10, the proportion of the pixel impermeable layer is less than 20%, and the pixel is not the local brightness maximum value, that is, it meets 8L-AB<4;

[0059] S26, for the pixels that meet the conditions in step S25, calculate the sum of the surrounding pixels and A, B, and substitute them into the formula in step S24 to calculate the overflow brightness L of the surrounding pixels N , the actual brightness value of the pixel itself L S LL N In this embodiment, Figure 3 The results of urban vitality zone range extraction before and after removing the spillover effect in this embodiment are shown, wherein Figure (a) is the result map of urban vitality zone extraction using original night light remote sensing data, and Figure (b) is the result map of urban vitality zone extraction after removing the spillover effect, indicating that removing the spillover effect achieves higher-precision extraction of urban vitality zone boundaries.

[0060] In some specific embodiments, the specific steps of step S3 include:

[0061] S31, selecting pixels whose impermeable layer difference between two consecutive years is less than 10%, and calculating the light brightness difference between the two consecutive years.

[0062] S32, using Laplace distribution to define the probability density function of the brightness difference ΔL, then the probability density function of the brightness difference ΔL is: The maximum likelihood estimation method is used to calculate the parameters, and the calculation formula of the parameter μ is: The calculation formula of parameter λ is Where N is the number of pixels, ΔLi is the brightness difference of the i-th pixel;

[0063] S33, calculate t = ln20 × λ, so that the cumulative probability density in the interval [μ-t, μ+t] is equal to 0.95;

[0064] S34, judging whether the brightness change is a fluctuation according to the following rules: if the difference in light brightness of a pixel in two consecutive years is within the interval [μ-t,μ+t], the brightness values ​​of the pixel in two years are marked as a fluctuation change, otherwise they are marked as a normal brightness change.

[0065] In some specific embodiments, the specific steps in step S4 include:

[0066] S41, for the consecutive years marked as fluctuating changes, the least square method is used to perform linear fitting on the brightness values ​​of each year, and the difference between the fitted value and the actual value and the average value a of the difference are calculated. If the difference between the fitted value and the actual value is greater than the average value a, the brightness of the year is regarded as a normal brightness change point, otherwise it is confirmed as a fluctuation point;

[0067] S42, for the consecutive years confirmed to be fluctuating, the least square method is used to perform linear fitting again, and the annual brightness is set as the linear fitting value.

[0068] In this embodiment, Figure 5 This is the result graph of the time series data of the urban vitality zone range after removing the interannual fluctuation effect in this embodiment, where the blue points represent fluctuation points, the green points represent normal brightness change points, and a new fitting curve is formed after updating the fluctuation points.

[0069] It should be noted in this embodiment that the least square method (also known as the least square method) is a mathematical optimization technique that finds the best function match for data by minimizing the sum of squares of errors. This method is widely used in many disciplines such as error estimation, uncertainty, system identification, prediction, and forecasting data processing.

[0070] Basic principle: When we study the relationship between two variables (x, y), we usually get a series of paired data points (x1, y1, x2, y2, ..., x m ,y m ). Plot these data points in the xy rectangular coordinate system. If we find that these points are roughly distributed near a straight line, we can assume that the equation of this line is Y=a0+a1X (where a0 and a1 are the coefficients to be determined). In order to determine the best fit of this line, that is, to determine the values ​​of a0 and a1, we apply the principle of least squares to minimize the sum of the squares of the vertical distances (i.e., residuals) from all data points to this line.

[0071] The present invention utilizes the influence of impermeable layers on night light remote sensing data, takes into account the errors caused by spillover effects and interannual fluctuation effects, overcomes the boundary deviation and continuity problems caused by the previous use of only a single threshold for extraction, removes the errors caused by spillover effects and interannual fluctuation effects, and achieves higher-precision extraction of urban vitality areas, which helps to accurately understand the urban spatial structure and development process, and provide a decision-making basis for urban sustainable development, urban planning, and urban management.

Claims

1. A method for generating a sequence of urban vitality zones taking into account the spatiotemporal interference of night light remote sensing, characterized in that: The following steps are involved: S1, obtaining the annual time series of night light remote sensing data and impermeable layer data, aligning the spatial ranges of the night light remote sensing data and the impermeable layer data of a single scene, and calculating the proportion of the impermeable layer in a single pixel of the night light remote sensing data; S2, establishing a calculation formula for overflow brightness of eight neighboring pixels to a central pixel, screening the central pixel that is easily affected by the eight neighboring pixels, and removing the overflow brightness of the eight neighboring pixels to the central pixel according to the calculation formula for overflow brightness; S3, calculating the brightness difference histogram of the night light remote sensing data of two adjacent periods, fitting the histogram with Laplace distribution, calculating the cumulative probability density to determine the brightness fluctuation range, and distinguishing the fluctuation point and the normal change point; S4, performing linear fitting on the continuous fluctuation points, removing the fluctuation points that deviate too much from the fitting straight line, and then performing the linear fitting again to update the brightness values ​​of the remaining fluctuation points; S5, set a threshold T, regard pixels greater than T as urban vitality areas, extract the night light time series data year by year, and obtain the time series of the urban vitality area range.

2. The method for generating a sequence of urban active areas taking into account the spatiotemporal interference of night light remote sensing according to claim 1 is characterized in that: The specific steps of step S1 include: S11, obtaining the night light remote sensing data synthesized annually for more than 5 years, and obtaining the annual time series of the impermeable layer data, wherein the pixel size of the night light remote sensing data needs to be more than 10 times the pixel size of the impermeable layer data; S12, cutting and splicing the night light remote sensing data and the impermeable layer data to generate a single image with the same starting point and a spatial range of 5°×5°, and performing subsequent processing based on the single image; S13, for the night light remote sensing data and the impermeable layer data of a single scene in the same year, count the area of ​​the impermeable layer falling within a single night light pixel, and calculate the area ratio of the impermeable layer within a single night light pixel.

3. The method for generating a sequence of urban active areas taking into account the spatiotemporal interference of night light remote sensing according to claim 1, characterized in that: The specific steps in step S2 include: S21, based on the proportion of the impermeable layer of the luminous pixel, screen the pixels that do not emit light or emit a small amount of light, and whose brightness is mainly affected by the surrounding pixel images. The screening condition is that the proportion of the impermeable layer is less than 5%. The brightness L of the pixel C that meets the screening condition is also the brightness overflow value L of the surrounding pixels. N ; S22, for a single pixel c that is easily affected by surrounding pixels, take its eight neighboring pixels a0, a1, a2, a3, b0, b1, b2, b3, and divide the brightness of the surrounding pixels into A and B according to the different distances from the surrounding pixels to the pixel c, where A=L a0 +L a1 +L a2 +L a3 , B=L b0 +L b1 +L b2 +L b3 , L i is the light brightness of pixel i; S23, add a weight factor q to the surrounding pixel brightness sum B, q∈[0,1], set the sum of the brightness values ​​of the pixels surrounding the pixel c S to S=A+qB, increase q from 0 by 0.05, and calculate the surrounding pixel brightness overflow value L of the pixel c N The Pearson correlation coefficient of the surrounding pixel brightness sum S is obtained, the weight factor q when the maximum correlation coefficient is obtained, and the surrounding pixel brightness sum S is calculated; S24, using the least squares method to calculate the surrounding brightness sum S and the surrounding pixel brightness overflow value L N Perform linear fitting and calculate k and b, which are expressed as L N =kS+b; S25, screening pixels that simultaneously meet the conditions that the pixel brightness value L is less than 10, the proportion of the pixel impermeable layer is less than 20%, and the pixel is not the local brightness maximum value, that is, the pixels that meet 8L-AB<4; S26, for the pixels that meet the conditions in step S25, calculate the sum of the surrounding pixels and A, B, and substitute them into the formula in step S24 to calculate the overflow brightness L of the surrounding pixels N , the actual brightness value of the pixel itself L S LL N To obtain.

4. The method for generating a sequence of urban active areas taking into account the spatiotemporal interference of night light remote sensing according to claim 1, characterized in that: The specific steps of step S3 include: S31, selecting pixels whose impermeable layer difference ratio between two consecutive years is less than 10%, and calculating the light brightness difference between the two consecutive years; S32, using Laplace distribution to define the probability density function of the brightness difference ΔL, then the probability density function of the brightness difference ΔL is: The maximum likelihood estimation method is used to calculate the parameters, and the calculation formula of the parameter μ is: The calculation formula of parameter λ is Where N is the number of pixels, ΔL i is the brightness difference of the i-th pixel; S33, calculate t = ln20 × λ, so that the cumulative probability density in the interval [μ-t, μ+t] is equal to 0.95; S34, judging whether the brightness change is a fluctuation according to the following rules: if the difference in light brightness of a pixel in two consecutive years is within the interval [μ-t,μ+t], the brightness values ​​of the pixel in two years are marked as a fluctuation change, otherwise they are marked as a normal brightness change.

5. The method for generating a sequence of urban active areas taking into account the spatiotemporal interference of night light remote sensing according to claim 1, characterized in that: The specific steps in step S4 include: S41, for the consecutive years marked as fluctuating changes, the least square method is used to perform linear fitting on the brightness values ​​of each year, and the difference between the fitted value and the actual value and the average value a of the difference are calculated. If the difference between the fitted value and the actual value is greater than the average value a, the brightness of the year is regarded as a normal brightness change point, otherwise it is confirmed as a fluctuation point; S42, for the consecutive years confirmed to be fluctuating, the least square method is used to perform linear fitting again, and the annual brightness is set as the linear fitting value.

Images