Industrial activity intensity grid quantification method based on machine learning

Abstract

The invention relates to the technical field of city analysis, in particular to an industrial activity intensity grid quantification method based on machine learning, which comprises the following steps of: fusing night light grid data, population grid data, building volume grid data and normalized building indexes on the same grid scale and the mutual combination relationship among the data to obtain a regional feature vector; and based on the regional feature vectors and the statistical data, performing quantitative inversion on regional industrial activity intensity grids by using a machine learning method, and realizing conversion from a regional industrial total amount to a grid scale through an industrial potential weight to obtain an industrial intensity grid map. According to the method, high-cost data such as enterprise catalogues, fine land use and high-resolution noctilucence are not needed, the adaptability to data sources is high, and the method is convenient to popularize and apply in a large-scale distribution range.

Description

Technical Field

[0001] This invention belongs to the field of urban analysis technology, and in particular relates to a raster quantification method for industrial activity intensity based on machine learning. Background Technology

[0002] Traditional monitoring of industrial activity intensity relies heavily on statistical yearbook data and data reported by enterprises. Statistical yearbook data assesses regional industrial development levels using indicators such as GDP, industrial added value, and output value of large-scale industrial enterprises published by statistical departments at various levels. However, this type of data has a significant time lag (typically 1-2 years) and lacks spatial granularity, failing to reflect spatial heterogeneity within cities. Enterprise-reported data depends on enterprises periodically reporting information such as output value, profits, and electricity consumption, which suffers from inconsistent statistical standards, untimely or incomplete reporting, and difficulty in ensuring data continuity and reliability at regional scales and over long time series.

[0003] Building upon this foundation, existing research combines remote sensing and multi-source spatial data to indirectly characterize the intensity of industrial activities, forming an emerging technological path centered on remote sensing monitoring. However, most remote sensing monitoring methods rely on optical imagery such as Landsat and MODIS for land use classification or simple linear regression modeling of nighttime light and GDP. On the one hand, relying on land use classification limits the accuracy of identifying industrial land and has a long update cycle; on the other hand, a single linear regression of nighttime light and GDP cannot accurately depict the complex nonlinear relationships between multiple sources of information such as population, buildings, and nighttime light, resulting in poor generalization across different regions and time periods. In methods that heavily rely on auxiliary data, introducing high-precision auxiliary information such as enterprise statistics, Points of Interest (POIs), refined land use classification, tax data, and electricity data can significantly improve accuracy, but these methods are difficult to promote and apply in areas with scarce data, data confidentiality, or untimely updates.

[0004] In summary, the above methods mainly have the following prominent problems: (1) Most methods can only provide an assessment of the intensity of industrial activities at the administrative unit level, and cannot be refined to the city or grid level; (2) High-resolution nighttime light, POI and enterprise statistical data are difficult to obtain, prone to missing data, and have high statistical costs; (3) It relies heavily on a single data source or linear assumptions, and is insufficient in mining the nonlinear relationships of multi-source features; (4) Most of the work focuses on estimating the intensity of activities at a certain point in time, and lacks dynamic monitoring of the intensity of industrial activities. Summary of the Invention

[0005] In view of this, the present invention aims to provide a machine learning-based raster quantification method for industrial activity intensity. By using machine learning, a nonlinear mapping relationship is constructed between multi-source features of "night light - population - building volume - building area" and regional industrial output value, thereby quantifying the total regional industrial output to the raster level and realizing spatialized industrial activity intensity classification.

[0006] To achieve the above objectives, the technical solution created by this invention is implemented as follows: A machine learning-based grid quantization method for industrial activity intensity includes: S1: Obtain nighttime light raster data, population raster data, building volume raster data, and normalized building index, as well as the GDP data of secondary industries in each city; obtain the normalized raster data of each of the four types of data; S2: Based on nighttime light grid data and population grid data, construct per capita nighttime light intensity grid data; based on nighttime light grid data and building volume grid data, construct per unit building volume nighttime light intensity grid data; based on building volume grid data and population grid data, construct per capita building volume grid data. S3: Construct a regional feature set using the four types of normalized raster data obtained in step S1 and the three types of raster data obtained in step S2; extract the corresponding type of features from the elements in the regional feature set to obtain the regional feature vector; S4: Using the regional feature vector obtained in step S3 as input and the secondary industry GDP data obtained in step S1 as output, machine learning is used to train the regression model to obtain the GDP prediction model. S5: Following steps S1 to S3, process the raster data of the area to be analyzed to obtain the feature vector to be analyzed; input the feature vector to be analyzed into the GDP prediction model obtained in step S4 to obtain the predicted total industrial output. S6: Using the feature vector to be analyzed obtained in step S5, determine the industrial potential weight of each grid, and allocate the predicted total industrial output obtained in step S5 according to the industrial potential weight to obtain the industrial intensity grid.

[0007] Furthermore, step S1 also includes: sampling the nighttime light raster data, population raster data, building volume raster data, and normalized building index to the same resolution, and processing them using a uniform mask.

[0008] Furthermore, in step S2: The per capita nighttime light intensity grid data is as follows: ; in, N represents the average nighttime light intensity per person for grid cell i at time t. i,tP represents the nighttime radiance value of raster cell i in the nighttime light raster data at time t. i,t This represents the population count of raster i at time t in the population raster data; The grid data for nighttime light intensity per unit building volume is as follows: ; in, BV represents the nighttime luminescence intensity per unit building volume of grid cell i at time t. i,t This represents the building volume index of raster i at time t in the building volume raster data; The per capita building volume raster data is as follows: ; in, This represents the per capita building volume of grid cell i at time t; When per capita night light intensity grid data Luminous intensity per unit building volume and per capita building volume When the denominator is 0, it directly affects the per capita nighttime light intensity raster data. Luminous intensity per unit building volume and per capita building volume It is 0.

[0009] Furthermore, in step S3: the four types of data obtained in step S1 are summed, and the three types of data obtained in step S2 are averaged, thus completing the feature extraction of elements in the regional feature set and obtaining a regional feature vector composed of seven types of data.

[0010] Furthermore, in step S4, the random forest regression algorithm is used to train the regression model to obtain the GDP prediction model.

[0011] Furthermore, step S6 includes: performing a weighted summation of each feature in the feature vector to be analyzed corresponding to region j to obtain the industrial activity potential index; and calculating the sum of the potential indices of all grids within region j. The industrial potential weight W corresponding to grid i in the industrial intensity grid is obtained by the following formula. i,t : ; Among them, S i,t This represents the industrial activity potential index corresponding to grid i; the predicted total industrial output is allocated according to the industrial potential weights to obtain the industrial intensity grid.

[0012] Furthermore, the methods also include: S7: Based on the industrial intensity grid obtained in step S6, construct the industrial growth rate grid, and classify the grids of the industrial growth rate grid to form an industrial activity intensity distribution map.

[0013] Furthermore, the industrial growth rate of each grid cell in the industrial growth rate grid is obtained using the following formula: ; in, G represents the industrial growth rate at time t of grid i. i,t This represents the industrial intensity of grid i at time t in the industrial intensity grid. δ represents the industrial intensity of grid i in the industrial intensity grid at reference time t0, and δ represents the smoothing term.

[0014] Furthermore, the industrial growth rate grid is classified using the natural discontinuity method.

[0015] Compared with the prior art, the present invention can achieve the following beneficial effects: (1) In the machine learning-based raster quantification method for industrial activity intensity described in this invention, the machine learning method can automatically capture the nonlinear relationship and interaction effect between multi-source data features such as nighttime light raster data, population raster data, building volume raster data and normalized building index. Compared with traditional linear regression or empirical threshold methods, it has higher model accuracy and generalization ability. Through regional-raster two-level modeling and inversion, this invention can obtain the raster-level industrial intensity distribution while maintaining macro-statistical consistency, and characterize the spatial heterogeneity within the region and industrial cluster. This invention has strong compatibility with data sources and can select nighttime light, population or building data and statistical data indicators with different time and spatial resolutions according to the data conditions of different regions. It can also expand more related data (such as road density data), which is convenient for promotion and application in different regions and at different time scales. (2) In the machine learning-based grid quantification method for industrial activity intensity described in this invention, industrial growth and growth rate indicators are constructed on the basis of grid-level industrial intensity, and the natural discontinuity method is combined to classify the levels, thereby realizing dynamic monitoring of industrial activity intensity and providing decision-making basis for large-scale fine industrial distribution planning and industrial intensity change detection. Attached Figure Description

[0016] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments and descriptions of the invention are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings: Figure 1 This is a schematic flowchart of the machine learning-based raster quantization method for industrial activity intensity described in an embodiment of the present invention. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and do not constitute a limitation thereof. It should be noted that, unless otherwise specified, the embodiments and features described herein can be combined with each other.

[0018] The invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0019] like Figure 1 As shown in the embodiments of the present invention, the machine learning-based raster quantization method for industrial activity intensity includes: S1: Obtain nighttime light raster data, population raster data, building volume raster data, and normalized building index, as well as the GDP data of the secondary industry of each city; obtain the normalized raster data of each of the four types of data.

[0020] In some embodiments, step S1 further includes: sampling the nighttime light raster data, population raster data, building volume raster data, and normalized building index to the same resolution, and processing them using a uniform mask.

[0021] In this embodiment of the invention, the nighttime light grid data is derived from the monthly nighttime light composite product of NPP-VIIRS DNB. Extreme value and noise removal processing is performed on values ​​less than 0 and values ​​greater than those exceeding the maximum economic activity range in the original nighttime light grid data. An average composite is then performed based on the time series to obtain the nighttime radiance value N of grid cell i at time t. i,t ; The processed luminous radiance value N i,t Normalization using the following formula yields the first feature in the region feature set: ; in, N represents the normalized nighttime radiance value of grid cell i at time t. min N represents the minimum luminous intensity value. max This indicates the maximum luminous radiance value.

[0022] The population raster data comes from the LandScan population raster product, and the population raster has been resampled to the same spatial resolution as the nighttime light raster data; the processed population raster data P i,t Normalize the following formula: ; in, P represents the normalized population size of raster cell i at time t. i,t P represents the population count of raster i at time t in population raster data. min P represents the minimum population size. max Indicates the maximum population size; The normalized population raster data is spatially superimposed with the normalized nighttime radiance values, and the pixels that have valid values ​​in both are retained. The resulting feature is used as the second feature in the regional feature set.

[0023] Building volume raster data is derived from the Global Human Settlement Layer (GHSL) database. The building volume raster is resampled to a resolution consistent with the nighttime light raster and population raster data; the processed building volume raster (BV) is then analyzed. i,t Normalize the following formula: ; in, BV represents the normalized building volume index of grid cell i at time t. i,t BV represents the building volume index of raster i at time t in building volume raster data. min BV represents the minimum building volume index. max This indicates the maximum building volume index.

[0024] The normalized building volume raster data is subjected to the same masking process as the normalized nighttime light raster data and population raster data. The resulting feature is used as the third feature in the regional feature set.

[0025] The Normalized Difference Building Index (NDBI) is calculated based on MODIS data using the MOD09A1 V6.1 product, as shown in the following formula: ; Among them, NDBI i,t This represents the normalized building index of grid cell i at time t. This represents the surface reflectance of grid i at time t in the shortwave infrared band. This represents the surface reflectance of raster i at time t in the near-infrared band. Surface reflectance is obtained by cloud masking and temporal composite of images with low cloud cover throughout the year. The rasterized Normalized Building Index (NDBI) is resampled to a resolution consistent with the nighttime light raster data and population raster data; the processed NDBI... i,t Normalize the following formula: ; in, NDBI represents the normalized building index of raster cell i at time t. min The minimum normalized building index (NDBI) represents the minimum normalized building index. max This represents the maximum normalized building index. The normalized building index is then subjected to a masking process consistent with the normalized nighttime light raster data and population raster data. The resulting feature is used as the fourth feature in the regional feature set.

[0026] In this embodiment of the invention, the nighttime light raster data, population raster data, building volume raster data, and normalized building index required can all be obtained and calculated through publicly available data platforms, eliminating the need for extensive on-site enterprise surveys and cumbersome data statistics, thus significantly reducing data acquisition and update costs. To further enhance the ability to express the structural relationship between nighttime light, population, and buildings, this invention combines the normalized raster data of the nighttime light raster data, population raster data, and building volume raster data obtained in step S1 to construct three combined proportional indices as the latter three features in the regional feature set, used to characterize per capita nighttime light, nighttime light per building unit, and the combined intensity of nighttime light and building intensity.

[0027] S2: Based on nighttime light raster data and population raster data, construct per capita nighttime light intensity raster data; based on nighttime light raster data and building volume raster data, construct per unit building volume nighttime light intensity raster data; based on nighttime light raster data and normalized building index, construct per capita building volume raster data.

[0028] In some embodiments, the per capita nighttime light intensity raster data is as follows: ; in, This represents the average nighttime light intensity per person in grid cell i at time t. This feature is used to identify industrial, logistics, and other land uses with high nighttime light intensity and low population density. This feature is the fifth feature in the regional feature set. The grid data for nighttime light intensity per unit building volume is as follows: ; in, This represents the nighttime light intensity per unit building volume of grid cell i at time t. This feature helps to distinguish between manufacturing parks and low-intensity warehouses or idle land. This feature is the sixth feature in the set of regional features. The per capita building volume raster data is as follows: ; in, This represents the per capita building volume of grid cell i at time t. This feature is used to identify residential areas with high population density and large building volume. This feature is the seventh feature in the regional feature set.

[0029] In this embodiment of the invention, three combined proportional indicators are spatially superimposed, retaining elements where each indicator has a valid value; furthermore, when the three combined proportional indicators (i.e., per capita nighttime light intensity raster data) are... Luminous intensity per unit building volume and per capita building volume When the denominator of the fraction is 0, the per capita nighttime light intensity raster data is directly obtained. Luminous intensity per unit building volume and per capita building volume It is 0.

[0030] After steps S1 and S2, a set of regional features with uniform spatial resolution, uniform mask, and containing 7 types of features is obtained, which is used for subsequent region aggregation and modeling.

[0031] S3: Construct a regional feature set using the four types of normalized raster data obtained in step S1 and the three types of raster data obtained in step S2; extract the corresponding type of features from the elements in the regional feature set to obtain the regional feature vector.

[0032] In some embodiments, the feature extraction process in step S3 includes: The four types of data obtained in step S1 are summed, and the three types of data obtained in step S2 are averaged. Feature extraction is performed on the elements in the regional feature set, resulting in a regional feature vector composed of seven types of data.

[0033] In this embodiment of the invention, the region feature vector is specifically represented as follows: ; X j,t This represents the region feature vector of the j-th region at time t. This represents the total value of the k-th type of data at time t in the j-th region.

[0034] Region feature vector X j,t The first four elements are obtained from the following formula: ; Among them, Ω j This represents the total grid of the j-th region.

[0035] Region feature vector X j,t The last three elements are obtained from the following formula: .

[0036] S4: Using the regional feature vector obtained in step S3 as input and the secondary industry GDP data obtained in step S1 as output, machine learning is used to train a regression model to obtain a GDP prediction model. In some embodiments, a random forest regression algorithm is used to train the regression model to obtain a GDP prediction model. Random forest is an ensemble learning method that achieves nonlinear regression or classification by constructing multiple decision trees and integrating their results. This invention uses a random forest regression model to establish a nonlinear mapping relationship between multi-source remote sensing features and regional industrial activity indicators, which has advantages such as robustness to outliers, insensitivity to variable scale, and the ability to characterize high-dimensional nonlinear relationships.

[0037] In this embodiment of the invention, the training process in step S4 is as follows: Regional statistical industrial indicators Y for j regions across multiple years or time periods j,t For each combination (j,t), calculate the corresponding regional feature vector X. j,t ; Using the region feature vector X j,t As input for training, the corresponding regional statistical industrial index Y j,t As the output of training, construct the training set; The random forest regression algorithm is called using Python, and parameters such as the number of trees, maximum depth, and minimum number of samples are set. The hyperparameters are automatically tuned through cross-validation or grid search to obtain a trained GDP prediction model. Using the coefficient of determination R 2 Indicators for evaluating model performance: ; in, This represents the projected total industrial output obtained from the GDP forecasting model. Y represents the statistical industrial indicators for all regions. j,t The mean of R; 2 A higher value indicates that the model can better characterize the relationship between the features of "nighttime light – population grid – building volume grid – normalized difference building index" and the intensity of regional industrial activity.

[0038] S5: Following steps S1 to S3, process the raster data of the area to be analyzed to obtain the feature vector to be analyzed; input the feature vector to be analyzed into the GDP prediction model obtained in step S4 to obtain the predicted total industrial output.

[0039] After obtaining the predicted total industrial output, this invention quantifies the regional gross industrial output by industrial potential weights using a grid, thereby reversing the grid distribution of industrial activity intensity.

[0040] S6: Using the feature vector to be analyzed obtained in step S5, determine the industrial potential weight of each grid, and allocate the predicted total industrial output obtained in step S5 according to the industrial potential weight to obtain the industrial intensity grid.

[0041] In some embodiments, step S6 includes: The industrial activity potential index is obtained by weighted summing of the features in the feature vector to be analyzed for region j, as shown in the following formula: ; Among them, S i,t ω represents the industrial activity potential index of grid cell i at time t. k This represents the weight of the k-th data type. This indicates that the raster data of the region to be analyzed is processed according to steps S1 to S3 to obtain the k-th type of data, with a weight ω. k >0 can be set based on the importance of machine learning features or experience, such as having 7 equal weights; Calculate the sum of the potential indices of all grid cells within region j. As shown in the following formula: ; The industrial potential weight W corresponding to grid i in the industrial intensity grid is obtained by the following formula. i,t : ; Among them, S i,t This represents the industrial activity potential index corresponding to grid i; The predicted total industrial output is allocated according to industrial potential weights to obtain an industrial intensity grid, as shown in the following formula: ; Among them, G i,t This represents the industrial intensity of grid cell i at time t in the industrial intensity grid. This yields an industrial intensity grid map, enabling spatial inversion of industrial intensity from administrative units to the grid scale.

[0042] In order to dynamically monitor changes in the intensity of industrial activities, this invention further constructs an industrial growth rate grid and uses the natural discontinuity method to achieve automatic classification.

[0043] In some embodiments, the method further includes: S7: Based on the industrial intensity grid obtained in step S6, construct the industrial growth rate grid, and classify the grids of the industrial growth rate grid to form an industrial activity intensity distribution map. This distribution map can intuitively display the state and spatial distribution of industrial activity intensity at the grid scale.

[0044] Specifically, the industrial growth rate of each cell in the industrial growth rate grid is obtained using the following formula: ; in, This represents the industrial growth rate at time t of grid i. δ represents the industrial intensity of grid i in the industrial intensity grid at reference time t0, and δ represents the smoothing term; The natural discontinuity method is used to classify the industrial growth rate grid. The natural discontinuity method is a statistical classification method that divides continuous variables into several levels by minimizing within-class variance and maximizing between-class variance. This invention uses the natural discontinuity method to classify industrial activity intensity or industrial growth rate, automatically identifying levels such as "extremely high," "high," "medium," "low," and "extremely low," thereby achieving spatial classification of industrial prosperity.

[0045] It should be understood that the various forms of processes shown above can be used to reorder, add, or delete steps. For example, the steps described in this invention disclosure can be executed in parallel, sequentially, or in different orders, as long as the desired result of the technical solution disclosed in this invention can be achieved, and this is not limited herein.

[0046] The specific embodiments described above do not constitute a limitation on the scope of protection of this invention. Those skilled in the art should understand that various modifications, combinations, sub-combinations, and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this invention should be included within the scope of protection of this invention.

Claims

1. A machine learning-based grid quantization method for industrial activity intensity, characterized in that, include: S1: Acquire nighttime light raster data, population raster data, building volume raster data and normalized building index, as well as the GDP data of secondary industry in each city; Obtain the normalized raster data for each of the four data types; S2: Based on nighttime light grid data and population grid data, construct per capita nighttime light intensity grid data; based on nighttime light grid data and building volume grid data, construct per unit building volume nighttime light intensity grid data; based on building volume grid data and population grid data, construct per capita building volume grid data. S3: Construct a regional feature set using the four types of normalized raster data obtained in step S1 and the three types of raster data obtained in step S2; extract the corresponding type of features from the elements in the regional feature set to obtain the regional feature vector; S4: Using the regional feature vector obtained in step S3 as input and the secondary industry GDP data obtained in step S1 as output, machine learning is used to train the regression model to obtain the GDP prediction model. S5: Following steps S1 to S3, process the raster data of the area to be analyzed to obtain the feature vector to be analyzed; input the feature vector to be analyzed into the GDP prediction model obtained in step S4 to obtain the predicted total industrial output. S6: Using the feature vector to be analyzed obtained in step S5, determine the industrial potential weight of each grid, and allocate the predicted total industrial output obtained in step S5 according to the industrial potential weight to obtain the industrial intensity grid.

2. The machine learning-based grid quantization method for industrial activity intensity according to claim 1, characterized in that, Step S1 also includes: sampling the nighttime light raster data, population raster data, building volume raster data, and normalized building index to the same resolution, and processing them using a uniform mask.

3. The machine learning-based raster quantization method for industrial activity intensity according to claim 1, characterized in that, In step S2: The per capita nighttime light intensity grid data is as follows: ； in, N represents the average nighttime light intensity per person for grid cell i at time t. i,t P represents the nighttime radiance value of raster cell i in the nighttime light raster data at time t. i,t This represents the population count of raster i at time t in the population raster data; The grid data for nighttime light intensity per unit building volume is as follows: ； in, BV represents the nighttime luminescence intensity per unit building volume of grid cell i at time t. i,t This represents the building volume index of raster i at time t in the building volume raster data; The per capita building volume raster data is as follows: ； in, This represents the per capita building volume of grid cell i at time t; When per capita night light intensity grid data Luminous intensity per unit building volume and per capita building volume When the denominator is 0, it directly affects the per capita nighttime light intensity raster data. Luminous intensity per unit building volume and per capita building volume It is 0.

4. The machine learning-based grid quantization method for industrial activity intensity according to claim 1, characterized in that, In step S3: The four types of data obtained in step S1 are summed, and the three types of data obtained in step S2 are averaged. Feature extraction is performed on the elements in the regional feature set, resulting in a regional feature vector composed of seven types of data.

5. The machine learning-based grid quantization method for industrial activity intensity according to claim 1, characterized in that, In step S4, the random forest regression algorithm is used to train the regression model to obtain the GDP prediction model.

6. The machine learning-based grid quantization method for industrial activity intensity according to claim 1, characterized in that, Step S6 includes: The industrial activity potential index is obtained by weighted summing of each feature in the feature vector to be analyzed corresponding to region j. Calculate the sum of the potential indices of all grid cells within region j. ; The industrial potential weight W corresponding to grid i in the industrial intensity grid is obtained by the following formula. i,t : ； Among them, S i,t This represents the industrial activity potential index corresponding to grid i; The predicted total industrial output is allocated according to the industrial potential weights to obtain an industrial intensity grid.

7. The machine learning-based raster quantization method for industrial activity intensity according to claim 1, characterized in that, The method also includes: S7: Based on the industrial intensity grid obtained in step S6, construct the industrial growth rate grid, and classify the grids of the industrial growth rate grid to form an industrial activity intensity distribution map.

8. The machine learning-based raster quantization method for industrial activity intensity according to claim 7, characterized in that, The industrial growth rate of each cell in the industrial growth rate grid is obtained by the following formula: ； in, G represents the industrial growth rate at time t of grid i. i,t This represents the industrial intensity of grid i at time t in the industrial intensity grid. δ represents the industrial intensity of grid i in the industrial intensity grid at reference time t0, and δ represents the smoothing term.

9. The machine learning-based grid quantization method for industrial activity intensity according to claim 7, characterized in that, The industrial growth rate grid is classified using the natural discontinuity method.

Images