A natural resource right registration map spot linkage updating method and system

Abstract

The application provides a natural resource right registration plot linkage updating method and system, which comprises the following steps: acquiring reference period and monitoring period remote sensing images and right registration plots in the same geographical range, calculating the difference between the stable ground feature spectrum mean values in the two images, obtaining the system offset vector and correcting the each pixel spectrum change vector, constructing the band contribution matrix based on the ground feature spectrum knowledge base, weighting the corrected change vector, calculating the weighted change intensity and change direction angle, establishing the type transfer rule library, defining the reference direction angle for the preset plot type conversion path, constructing the tolerance function, determining the effective change pixel when the pixel weighted change intensity exceeds the initial type intensity threshold and the difference between the change direction angle and the reference direction angle of the conversion path is less than the direction angle tolerance determined by the tolerance function, aggregating all the effective change pixels, identifying the original right registration plot area as the to-be-updated plot if the adjacent cluster area exceeds the preset threshold, and generating the update list.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of linkage update, and particularly relates to a natural resource right registration plot linkage update method and system. BACKGROUND

[0002] In order to ensure the accuracy of the right registration data, it is necessary to establish a reliable plot data update mechanism to reflect the actual changes in time. At present, the change detection using multi-temporal remote sensing images is the main technical means to realize plot update. The traditional method mostly adopts the "classification comparison" strategy, that is, the images of the two periods before and after are respectively classified by supervised or unsupervised classification, and then the change information is extracted by superimposing the classification result maps. However, any error generated in the classification process will be accumulated and transmitted to the change detection result, which is easy to produce a large number of "pseudo change" plots, resulting in low reliability of the update. Therefore, the change detection method of directly comparing multi-temporal images is proposed, such as the change vector analysis (CVA) method. By directly calculating the change vector of the image element in the multi-dimensional spectral space, and using the vector length and direction angle to represent the change intensity and change direction respectively, the classification process is avoided. The traditional CVA method usually treats all spectral bands equally, ignoring the fact that the sensitivity and contribution of different bands to the conversion of specific ground object types are different, which may lead to the submersion and misjudgment of key change information. Lack of self-adjusting ability to the spectral variability of complex ground objects leads to low recognition accuracy of the change type, and it is difficult to realize the intelligent linkage from remote sensing monitoring change to right data update. SUMMARY

[0003] The application provides a natural resource right registration plot linkage update method, which comprises the following steps:

[0004] obtaining the reference period and monitoring period remote sensing images of the same geographical range and the corresponding right registration plots;

[0005] calculating the difference between the stable ground object spectral mean values in the reference period and monitoring period images to obtain a system offset vector, and subtracting the system offset vector from the spectral change vector of each image element from the reference period to the monitoring period to obtain a background corrected change vector; constructing a band contribution matrix based on a ground object spectral knowledge base, weighting the background corrected change vector, and calculating the weighted change intensity and change direction angle;

[0006] A type transition rule library is established, a reference direction angle is defined for each preset type transition path of the image patch, and a tolerance function with a weighted change intensity as input and a direction angle tolerance as output is established; for any pixel, when the weighted change intensity of the pixel exceeds the intensity threshold corresponding to the initial type of the image patch in which the pixel is located, and the absolute value of the difference between the change direction angle and the reference direction angle of at least one transition path starting from the initial type is less than the direction angle tolerance determined by the tolerance function of the path and according to the weighted change intensity of the pixel, the pixel is determined as a valid change pixel;

[0007] All valid change pixels are aggregated, if the area of a cluster formed by geographically adjacent valid change pixels exceeds a preset area threshold, the original legal right registered image patch area corresponding to the cluster is identified as an image patch to be updated, and an update list is generated.

[0008] Optionally, the difference between the mean values of the spectra of stable ground objects in the reference period and the monitoring period is calculated to obtain a system offset vector, including:

[0009] The normalized difference vegetation index method is used to screen, in the reference period and the monitoring period, pixels with stable normalized difference vegetation index values and a change amount less than 0.05 as stable ground objects;

[0010] The arithmetic mean values of the spectral values of all stable ground objects in each band in the reference period and the monitoring period are calculated respectively to obtain a reference period spectral mean vector and a monitoring period spectral mean vector;

[0011] The monitoring period spectral mean vector is subtracted from the reference period spectral mean vector to obtain a system offset vector.

[0012] Optionally, the band contribution degree matrix is constructed based on the ground object spectral knowledge base, including:

[0013] For a specific transition path from an initial type A to a target type B, the standard spectral reflectance of the two types of ground objects A and B is extracted from the ground object spectral knowledge base;

[0014] The absolute values of the reflectance differences of the two types of ground objects in each spectral band i are calculated ;

[0015] The difference values of all bands are normalized by the maximum value to obtain the contribution degree of each band under the transition path, thereby constructing the band contribution degree matrix.

[0016] Optionally, the background corrected change vector is weighted to calculate the weighted change intensity and the change direction angle, including:

[0017] For each preset type conversion path, multiply the change component of each band in the background corrected change vector by the band contribution degree corresponding to the conversion path to obtain a weighted change vector under the path;

[0018] Calculate the Euclidean norm of the weighted change vector as the weighted change intensity under the path;

[0019] Select the weighted change components of the near-infrared band and the red band as two-dimensional coordinates, and use the two-parameter inverse tangent function to calculate the change direction angle under the path.

[0020] Optionally, the establishment of the tolerance function with the weighted change intensity as input and the direction angle tolerance as output comprises:

[0021] A negative exponential function is used as the tolerance function, and the specific form is where M is the weighted change intensity, T(M) is the direction angle tolerance, A, k, and C are normal numbers fitted by sample data, so that the greater the weighted change intensity, the smaller the direction angle tolerance.

[0022] Optionally, for any pixel, when the pixel weighted change intensity exceeds the intensity threshold corresponding to the initial type of the pixel in the plot, and the absolute value of the difference between the change direction angle and the reference direction angle of at least one conversion path starting from the initial type is less than the direction angle tolerance determined by the tolerance function of the path and according to the weighted change intensity of the pixel, the pixel is determined as an effective change pixel, comprising:

[0023] For a conversion path in which an initial type pixel of forest land is converted into construction land, calculate the weighted change intensity and the change direction angle of the pixel under the conversion path;

[0024] If the weighted change intensity under the path exceeds the preset weighted change intensity threshold, and the absolute value of the difference between the change direction angle under the path and the reference direction angle of the forest land to construction land conversion path is less than the direction angle tolerance calculated by substituting the weighted change intensity under the path into the tolerance function, the pixel is determined as an effective change pixel.

[0025] Optionally, the aggregation of all effective change pixels comprises:

[0026] A binary matrix of the same size as the remote sensing image is created, and the positions of the effective change pixels are assigned a value of 1, and the remaining positions are assigned a value of 0;

[0027] An eight-neighbor connected component labeling algorithm is used to scan the binary matrix, and the pixels assigned a value of 1 that are spatially connected to each other are divided into the same cluster.

[0028] Optionally, the preset area threshold comprises:

[0029] When the initial type of the patch to be updated is cultivated land, the preset area threshold is 400㎡;

[0030] When the initial type is woodland, the preset area threshold is 600㎡.

[0031] Furthermore, this invention also relates to a system for linking and updating map features for natural resource ownership registration, comprising the following modules:

[0032] The acquisition module is used to acquire remote sensing images of the baseline and monitoring periods and the corresponding land rights registration patches within the same geographical area;

[0033] The calculation module is used to calculate the difference between the mean spectral values ​​of stable ground features in the images of the reference period and the monitoring period to obtain the system offset vector, and to subtract the system offset vector from the spectral change vector of each pixel from the reference period to the monitoring period to obtain the background correction change vector; based on the ground feature spectral knowledge base, a band contribution matrix is ​​constructed, and the background correction change vector is weighted to calculate the weighted change intensity and change direction angle.

[0034] The determination module is used to establish a type transfer rule base, define a reference direction angle for each preset patch type transformation path, and establish a tolerance function with weighted change intensity as input and direction angle tolerance as output. For any pixel, when the pixel's weighted change intensity exceeds the intensity threshold corresponding to the initial type of the patch, and the absolute value of the difference between the change direction angle and the reference direction angle of at least one transformation path starting from the initial type is less than the direction angle tolerance determined by the tolerance function of the path and based on the pixel's weighted change intensity, then the pixel is determined as a valid changed pixel.

[0035] The generation module is used to aggregate all valid changed pixels. If the cluster area formed by geographically adjacent valid changed pixels exceeds a preset area threshold, the original land registration area corresponding to the cluster is identified as a land patch to be updated, and an update list is generated.

[0036] Preferably, the step of calculating the difference between the mean spectral values ​​of stable ground features in the baseline and monitoring period images to obtain the system offset vector includes:

[0037] Using the normalized vegetation index (NVC) method, pixels with stable NVC values ​​and variations of less than 0.05 were selected as stable ground features in the baseline and monitoring period images.

[0038] The arithmetic mean of the spectral values ​​of all stable ground features in each band during the baseline period and the monitoring period were calculated respectively to obtain the spectral mean vector of the baseline period and the spectral mean vector of the monitoring period.

[0039] Subtracting the reference period spectrum mean vector from the monitoring period spectrum mean vector, a system offset vector is obtained.

[0040] Preferably, the band contribution degree matrix is constructed based on the ground object spectrum knowledge base, comprising:

[0041] For a specific conversion path from an initial type A to a target type B, the standard spectral reflectance of the A and B types of ground objects is extracted from the ground object spectrum knowledge base;

[0042] The absolute value of the reflectance difference of the two types of ground objects at each spectral band i is calculated ;

[0043] The difference values of all bands are normalized by maximum value to obtain the contribution degree of each band under the conversion path, thereby constructing the band contribution degree matrix.

[0044] Preferably, the background corrected change vector is weighted to calculate the weighted change intensity and change direction angle, comprising:

[0045] For each preset polygon type conversion path, the change component of each band in the background corrected change vector is multiplied by the band contribution degree corresponding to the conversion path to obtain the weighted change vector under the path;

[0046] The Euclidean norm of the weighted change vector is calculated as the weighted change intensity under the path;

[0047] The weighted change components of the near-infrared band and the red band are selected as two-dimensional coordinates, and the change direction angle under the path is calculated using the two-parameter arctangent function.

[0048] Preferably, the tolerance function is established with the weighted change intensity as input and the direction angle tolerance as output, comprising:

[0049] A negative exponential function is used as the tolerance function, and the specific form is , wherein M is the weighted change intensity, T(M) is the direction angle tolerance, A, k, and C are normal numbers fitted by sample data, so that the greater the weighted change intensity, the smaller the direction angle tolerance.

[0050] Preferably, for any pixel, when the pixel weighted change intensity exceeds the intensity threshold value corresponding to the initial type of the polygon, and the absolute value of the difference between the change direction angle and the reference direction angle of at least one conversion path starting from the initial type is less than the direction angle tolerance determined by the tolerance function of the path and according to the weighted change intensity of the pixel, the pixel is determined as an effective change pixel, comprising:

[0051] a transformation path of a pixel element of an initial type of forest land to construction land, calculating a weighted change intensity and a change direction angle of the pixel element under the transformation path;

[0052] If the weighted change intensity under the path exceeds a preset weighted change intensity threshold, and an absolute value of a difference between the change direction angle under the path and a reference direction angle of the transformation path from forest land to construction land is less than a direction angle tolerance calculated by substituting the weighted change intensity under the path into a tolerance function, the pixel element is determined as an effective change pixel element.

[0053] Preferably, the aggregation of all effective change pixel elements comprises:

[0054] A binary matrix with the same size as the remote sensing image is created, and the positions of the effective change pixel elements are assigned a value of 1, and the remaining positions are assigned a value of 0;

[0055] An eight-neighborhood connected component labeling algorithm is used to scan the binary matrix, and the pixel elements assigned a value of 1 that are spatially connected to each other are divided into the same cluster.

[0056] Preferably, the preset area threshold comprises:

[0057] When the initial type of the to-be-updated graph patch is cultivated land, the preset area threshold is 400 square meters;

[0058] When the initial type is forest land, the preset area threshold is 600 square meters.

[0059] The present application performs background correction on the image by calculating the system offset vector, eliminates the pseudo change information caused by factors such as temporal difference, sensor difference and non-ground object change, and improves the authenticity of change detection; uses the band contribution matrix based on the ground object spectrum knowledge base to weight the change vector, highlights the spectral information sensitive to the transformation of specific ground object types, and enhances the recognition ability of key change types; establishes a graph patch type transfer rule library and a reference direction angle, integrates prior knowledge of ground object transformation into the change discrimination process, so that the determination of the change type has a clear physical basis; the matching tolerance of the change direction is associated with the change intensity, and a more reasonable discrimination standard is set for changes of different intensity, reducing the confusion and misjudgment of the change type. By aggregating the effective change pixel elements and performing area constraint, it is ensured that the to-be-updated graph patch recognized has practical significance, and the reliability of the natural resource right registration graph patch linkage update is improved as a whole. BRIEF DESCRIPTION OF DRAWINGS

[0060] Figure 1 a flowchart for the first embodiment;

[0061] Figure 2 a tolerance function diagram;

[0062] Figure 3 To effectively change the pixel aggregation schematic diagram. DETAILED DESCRIPTION

[0063] In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present description. However, the present description can be practiced without the specific details, other than in the examples described herein, and it can be apparent to those skilled in the art that the present description can be practiced without such details. In other instances, well-known methods, procedures, components, and networks have not been described in detail so as not to unnecessarily obscure aspects of the present description.

[0064] The terminology used in one or more embodiments of the present description is for the purpose of describing particular embodiments only and is not intended to be limiting of one or more embodiments of the present description. As used in one or more embodiments of the present description and the accompanying claims, the singular forms "a," "an," and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms "comprises" and / or "comprising," when used in one or more embodiments of the present description, specify the presence of stated features, integers, steps, operations, elements, and / or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and / or groups thereof.

[0065] It will be understood that, although the terms first, second, etc. can be used herein to describe various information, these terms are not intended to denote a temporal or chronological order. Rather, these terms are used only to distinguish one from another. For example, without departing from the scope of one or more embodiments of the present description, first can be termed second, and similarly, second can be termed first. The word "if' as used herein, meaning "when" or "upon the occurrence of," can be understood differently, depending on the context.

[0066] In the first embodiment, the present application proposes a natural resource right registration plot linkage updating method, such as Figure 1 , comprising the following steps:

[0067] S1, acquiring reference period and monitoring period remote sensing images of the same geographical range and corresponding right registration plot;

[0068] The reference period image is, for example, a domestic high-resolution second satellite multispectral image acquired in summer 2020. The monitoring period image is a domestic high-resolution second satellite multispectral image of the same area acquired in summer 2022, ensuring consistency in imaging phase and sensor type to reduce differences. The two images are radiometrically calibrated, atmospherically corrected, and strictly geometrically corrected to ensure accurate pixel-to-pixel registration under the same coordinate system. The corresponding right registration plot is a vector format national space planning or natural resource registration data matched with the reference period image phase, including the boundaries and initial land cover type attributes of each plot, such as forest land, grassland, water area, and construction land.

[0069] S2, calculating the difference between the mean values of the stable feature spectrum in the reference period and the monitoring period image to obtain a system offset vector, and subtracting the system offset vector from the spectral change vector of each pixel from the reference period to the monitoring period to obtain a background corrected change vector; based on the feature spectrum knowledge base, a band contribution matrix is constructed, and the background corrected change vector is weighted to calculate the weighted change intensity and change direction angle;

[0070] Select the features in the image that are stable in nature and do not change substantially as pseudo-invariant features, such as airport runways, large building roofs, and clear deep water bodies. Extract the spectral values of all pixels in the pseudo-invariant feature region in the reference period and the monitoring period for each band. For each band, calculate the average of the spectral values of all pseudo-invariant feature pixels in the reference period and the average of the spectral values in the monitoring period. Subtract the average of the corresponding band in the reference period from the average of each band in the monitoring period to obtain an N-dimensional vector, where N is the number of bands. This vector is the system offset vector. For each pixel in the image, subtract the spectral values of the corresponding bands in the reference period image from the spectral values of the monitoring period image to obtain an original spectral change vector. Subtract the system offset vector calculated above from the original spectral change vector of this pixel to obtain a background corrected change vector, which more purely reflects the spectral response difference caused by the change of the feature itself. The feature spectrum knowledge base is a pre-collected and sorted standard spectrum curve data of various typical features, such as vegetation, water, bare soil, and buildings. According to the data, analyze the spectral response characteristics in the transformation process of specific feature types, for example, when forest land is converted into construction land, the reflectivity in the near-infrared band will decrease sharply, while the reflectivity in the visible red band will increase. Based on this, for each type of feature transformation, a weight vector is constructed, in which the most sensitive band to the transformation is given a high weight, and the insensitive or easily disturbed band is given a low weight. All weight vectors together constitute a band contribution matrix. For any pixel, multiply the background corrected change vector by the transformation path weight vector corresponding to the initial feature type in the band contribution matrix to obtain a weighted change vector. The weighted change intensity is the Euclidean norm of the weighted change vector, i.e. the square root of the sum of the squares of each component. The weighted change direction angle is the direction of the weighted change vector in the multi-dimensional spectral space, which can be calculated by the components, for example, the inverse tangent of the quotient of two components in a two-dimensional space.

[0071] S3, establishing a type transition rule library, defining a reference direction angle for each preset type transition path of the image patch, and establishing a tolerance function with the weighted change intensity as the input and the direction angle tolerance as the output; for any pixel, when the weighted change intensity of the pixel exceeds the intensity threshold corresponding to the initial type of the image patch, and the absolute value of the difference between the change direction angle and the reference direction angle of at least one transition path starting from the initial type is less than the direction angle tolerance determined by the tolerance function of the path and the weighted change intensity of the pixel, the pixel is determined as an effective change pixel;

[0072] Specifically, the type transition rule library is a set of logical rules established according to natural laws and urban and rural development plans, which is used to define which transitions between ground object types are reasonable, for example, forest land can be converted to grassland or construction land, but construction land usually will not be converted to forest land in a short period of time. For each reasonable transition path, such as from forest land to construction land, a standard weighted change vector in an ideal state is calculated using a spectral knowledge base or a large number of real change samples, and the direction of the vector is the reference direction angle of the transition path. A tolerance function is established, which takes the weighted change intensity as the independent variable and the direction angle tolerance as the dependent variable. Preferably, a reverse proportional function is used, that is, the greater the change intensity, the more intense and clear the change, and the closer the direction angle should be to the reference direction angle, so the tolerance is smaller; on the contrary, the change intensity is smaller, and there may be more noise, so the deviation range of the direction angle from the reference direction angle is allowed to be larger, that is, the tolerance is larger.

[0073] The weighted change intensity is calculated and compared with the preset forest land type change intensity threshold. If it is less than the threshold, it is considered that the change is not prominent, and it is determined as unchanged. If it is greater than the threshold, the direction angle judgment is entered. According to the type transition rule library, there are two reasonable transition paths starting from forest land, to grassland and to construction land. The absolute value of the difference between the weighted change direction angle of the pixel and the reference direction angle of the forest land to grassland path and the reference direction angle of the forest land to construction land path is calculated. At the same time, the weighted change intensity of the pixel is substituted into the tolerance function of the two paths of forest land to grassland and forest land to construction land respectively, to obtain two corresponding direction angle tolerance values. If the difference between the direction angle of the pixel and the reference direction angle of the forest land to construction land is less than the corresponding tolerance value, even if the difference with the reference direction angle of the forest land to grassland is larger, the pixel is still determined as an effective change pixel, and is preliminarily marked as converted from forest land to construction land.

[0074] S4, aggregating all effective change pixels, if the area of a cluster formed by geographically adjacent effective change pixels exceeds a preset area threshold, the original registered image patch area corresponding to the cluster is identified as a to-be-updated image patch, and an update list is generated.

[0075] All the valid change pixels determined above are generated into a binary image, where a valid change pixel is 1 and an unchanged pixel is 0. Using a connected domain analysis algorithm in image processing, the valid change pixels that are spatially connected or adjacent to each other are merged into individual clusters or patches. The total number of pixels contained in each cluster is calculated and multiplied by the actual geographic area represented by a single pixel to obtain the area of each change cluster. A minimum update area threshold is set, for example, 500 m2. All change clusters are traversed, and the clusters with an area smaller than the threshold are removed to eliminate isolated noise points or insignificant small changes. For the change clusters with an area exceeding the threshold, spatial position is superimposed and analyzed with the original right registration patch vector data, and any original patch that has a spatial overlap or inclusion relationship with the large-area change cluster is identified as a to-be-updated patch. The unique identification code, location information, etc. of all to-be-updated patches are summarized to obtain an update list, which is submitted to subsequent manual verification or right information change processes.

[0076] In an optional embodiment, the difference between the mean values of the stable ground object spectra in the reference period and the monitoring period images is calculated to obtain a system offset vector, including:

[0077] Using the normalized vegetation index method, the pixels with stable normalized vegetation index values and a change of less than 0.05 in the reference period and monitoring period images are selected as stable ground objects;

[0078] The arithmetic mean values of the spectral values of all stable ground objects in each band in the reference period and the monitoring period are calculated to obtain a reference period spectral mean vector and a monitoring period spectral mean vector;

[0079] The monitoring period spectral mean vector is subtracted from the reference period spectral mean vector to obtain a system offset vector.

[0080] Specifically, the reference period image of 2020 and the monitoring period image of 2021 are selected, and the normalized vegetation index of each pixel in the two images is calculated using the near-infrared band and the red light band. The normalized vegetation index values in the two images are compared pixel by pixel, for example, the index of a certain pixel in 2020 is 0.76, and in 2021 it is 0.78, the change is 0.02, which is less than the preset threshold of 0.05, so the pixel is a stable ground object. Through this method, all 50,000 stable ground object pixels in the image are selected. For the 50,000 stable ground objects, their values in each spectral band in the reference period and the monitoring period are counted respectively. For example, in the blue band of the reference period image, the values of the 50,000 pixels are added and divided by 50,000 to obtain the mean value of the blue band, which is 1050; the same operation is performed on all other bands such as green, red, etc. to form a multi-dimensional vector, that is, the reference period spectral mean vector, which is in the form of [1050, 1100, 900, 4500, 2500, 1500]. Similarly, the monitoring period spectral mean vector is calculated, which is in the form of [1070, 1125, 930, 4550, 2540, 1530]. Each component of the monitoring period spectral mean vector is subtracted from the corresponding component of the reference period spectral mean vector, and the result is the system offset vector [20, 25, 30, 50, 40, 30], which represents the radiation difference caused by factors such as sensor aging or atmospheric condition difference.

[0081] In an optional embodiment, the band contribution degree matrix is constructed based on the ground object spectral knowledge base, comprising:

[0082] For a specific conversion path from initial type A to target type B, the standard spectral reflectance of A and B ground objects is extracted from the ground object spectral knowledge base;

[0083] The absolute value of the reflectance difference of the two types of ground objects in each spectral band i is calculated ;

[0084] The difference values of all bands are normalized by the maximum value to obtain the contribution degree of each band under the conversion path, thereby constructing the band contribution degree matrix.

[0085] Specifically, taking the conversion of forest land to construction land as an example, the standard spectral reflectance curves of forest land and construction land are queried from the knowledge base such as the spectral library of the Geological Survey Bureau. Assuming that the remote sensing image contains four bands of blue, green, red and near-infrared, the standard reflectance of the forest land obtained by the query is 0.05, 0.08, 0.06 and 0.50 respectively in each band, and the standard reflectance of the construction land is 0.20, 0.22, 0.25 and 0.30 respectively. The absolute value of the reflectance difference of the two types of ground objects in each band is calculated, the difference value of the blue band is |0.05-0.20|, the result is 0.15; the difference value of the green band is |0.08-0.22|, the result is 0.14; the difference value of the red band is |0.06-0.25|, the result is 0.19; the difference value of the near-infrared band is |0.50-0.30|, the result is 0.20. A difference value vector [0.15, 0.14, 0.19, 0.20] is obtained. The maximum value in the vector is 0.20, and each element in the vector is divided by the maximum value. The normalized band contribution degree vector is calculated, the blue band is 0.75, the green band is 0.70, the red band is 0.95, and the near-infrared band is 1.00. The vector [0.75, 0.70, 0.95, 1.00] is the band contribution degree of the conversion path of forest land to construction land. Repeat the process for all preset conversion paths, and combine the obtained multiple contribution degree vectors to obtain a complete band contribution degree matrix.

[0086] In an optional embodiment, the background corrected change vector is weighted, and a weighted change intensity and a change direction angle are calculated, comprising:

[0087] For each preset conversion path of the plot type, the change component of each band in the background corrected change vector is multiplied by the band contribution degree corresponding to the conversion path to obtain a weighted change vector under the path;

[0088] The Euclidean norm of the weighted change vector is calculated as the weighted change intensity under the path;

[0089] The weighted change components of the near-infrared band and the red band are selected as two-dimensional coordinates, and the change direction angle under the path is calculated by using the two-parameter arctangent function.

[0090] Suppose a pixel has a change vector [30, 35, 80, -200] after background correction, corresponding to the changes in blue, green, red and near-infrared bands respectively. Now we need to determine whether the pixel has changed from forest land to construction land, and we know the band contribution degree vector of the transformation path is [0.75, 0.70, 0.95, 1.00]. Multiply each component of the change vector by the corresponding component of the contribution degree vector to obtain the weighted change vector, that is, [22.5, 24.5, 76, -200]. Calculate the Euclidean norm of the weighted change vector to obtain the weighted change intensity, the calculation result is approximately 216.5. The value 216.5 is the weighted change intensity of the pixel under the transformation path from forest land to construction land. Select the component 76 of the red light band in the weighted change vector as the X coordinate and the component -200 of the near-infrared band as the Y coordinate, and use the two-parameter inverse tangent function to calculate the change direction angle under the path, for example, -69.2°, indicating the direction of the change in the spectral space.

[0091] In an optional embodiment, the establishment of the tolerance function with the weighted change intensity as input and the direction angle tolerance as output comprises:

[0092] A negative exponential function is used as the tolerance function, and the specific form is where M is the weighted change intensity, T(M) is the direction angle tolerance, A, k, C are normal numbers fitted by sample data, so that the greater the weighted change intensity, the smaller the direction angle tolerance.

[0093] Specifically, from historical images and field survey data, a large number of sample pixels known to have changed from forest land to construction land are selected, as well as some pixels affected by noise but not having real changes. Calculate the weighted change intensity M and the difference between the change direction angle and the reference direction angle of each sample pixel. For example, the intensity M of a real change pixel is 200, and the angle difference is 3°; the intensity M of another noise pixel is 50, and the angle difference is 25°. After collecting hundreds of pairs of intensity M and angle difference data, plot the data points in the coordinate system with M as the horizontal axis and the angle difference as the vertical axis. Use mathematical methods such as nonlinear least squares to fit the data points to a negative exponential function. Through fitting, the specific parameters of the function can be determined, for example, A = 30, k = 0.02, C = 5, as shown in Figure 2 When the change intensity M of a pixel is very small, for example, 20, the direction angle tolerance T(20) is approximately 25.1°, allowing a larger deviation in the change direction; when the change intensity M is very large, for example, 200, the tolerance T(200) is sharply reduced to about 5.5°, requiring the change direction to be very close to the reference direction.

[0094] In an alternative embodiment, for any pixel, when the pixel weight change intensity exceeds the intensity threshold value corresponding to the initial type of the pixel's region, and the absolute value of the difference between the change direction angle and the reference direction angle of at least one conversion path from the initial type is less than the direction angle tolerance determined by the tolerance function and the weight change intensity of the pixel, the pixel is determined to be an effective change pixel, including:

[0095] For a conversion path in which a pixel of an initial type of forest land is converted to construction land, the weight change intensity and the change direction angle of the pixel under the conversion path are calculated.

[0096] If the weight change intensity under the path exceeds the preset weight change intensity threshold value, and the absolute value of the difference between the change direction angle under the path and the reference direction angle of the conversion path from forest land to construction land is less than the direction angle tolerance calculated by substituting the weight change intensity under the path into the tolerance function, the pixel is determined to be an effective change pixel.

[0097] Specifically, for the conversion path from forest land to construction land, two judgment criteria are preset, one is the minimum threshold of change intensity, which is 60, and the other is the function T(M) for calculating the angle tolerance, and the reference direction angle of the conversion path is -70°. Assuming that the weight change intensity M of the pixel is 216.5 and the change direction angle is -73° after calculation.

[0098] The weight change intensity 216.5 of the pixel is compared with the preset threshold value 60. Because 216.5>60, the intensity condition is met. First, the direction angle tolerance corresponding to the pixel is calculated, and the intensity value 216.5 is substituted into the tolerance function to obtain T(216.5)≈5.39°. The absolute value of the difference between the change direction angle of the pixel and the reference direction angle is calculated, and the result is 3°. Comparing the angle difference 3° with the calculated tolerance 5.39°, the pixel also meets the direction angle condition. Since both the intensity and the direction angle conditions are met, the pixel is an effective change pixel from forest land to construction land.

[0099] In an alternative embodiment, the aggregation of all effective change pixels includes:

[0100] A binary matrix of the same size as the remote sensing image is created, and the positions of the effective change pixels are assigned a value of 1, and the remaining positions are assigned a value of 0.

[0101] An eight-neighborhood connected component labeling algorithm is used to scan the binary matrix, and the pixels assigned a value of 1 that are spatially connected to each other are divided into the same cluster.

[0102] Suppose the processed remote sensing image region is 1000 pixels x 1000 pixels in size. After pixel-by-pixel change determination, a matrix of the same size, i.e. 1000 x 1000, is created and all elements are initialized to 0. For each pixel determined to be a valid change, the value in the corresponding spatial position of the binary matrix is changed from 0 to 1. For example, if the pixel at coordinates (15, 28) is a valid change pixel, the element value in the 15th row and 28th column of the matrix is set to 1. After this step, a graph containing only 0 and 1 is obtained, where 1 represents the position of the valid change pixel.

[0103] An eight-neighbor connected component labeling algorithm is applied to the above binary graph, starting from the top left corner of the image and scanning row by row. When the first pixel with a value of 1 is encountered, it is assigned a new label, for example, label 1, and placed in a pending queue. The eight neighbor pixels of the pixel, i.e. the pixels above, below, left, right and the four diagonal pixels, are checked. If a neighbor pixel also has a value of 1 and has not been labeled, it is assigned the same label 1 and also placed in the pending queue. The process is repeated until all pixels with a value of 1 that are connected to the original pixel are labeled with label 1. The image is scanned continuously to find the next unlabeled pixel with a value of 1, and a new label 2 is used to start a new round of labeling process. All pixels with a value of 1 are assigned to a label, and pixels with the same label form a spatially connected cluster, i.e. a change patch.

[0104] In an optional embodiment, the preset area threshold includes:

[0105] When the initial type of the to-be-updated patch is cultivated land, the preset area threshold is 400 m2;

[0106] When the initial type is forest land, the preset area threshold is 600 m2.

[0107] After aggregation to form multiple change patch clusters, the patches need to be screened according to a preset area threshold to eliminate excessive fine changes. The threshold refers to the minimum mapping area requirement of different land types in technical regulations or technical standards. For example, for an image with a spatial resolution of 10 m, each pixel represents an actual area of 100 m2. Suppose three change patches are obtained through the aggregation algorithm, patch one is composed of 5 pixels and the initial land type is cultivated land; patch two is composed of 7 pixels and the initial land type is forest land; and patch three is composed of 3 pixels and the initial land type is also cultivated land.

[0108] The actual area of each plot is calculated, the area of plot one is 500 m2, the area of plot two is 700 m2, and the area of plot three is 300 m2. The area is compared with the area threshold value of the corresponding class. When the initial type is cultivated land, the preset area threshold value is 400 m2. The area of plot one is 500 m2>400 m2, so it is retained as an effective change plot. The area of plot three is 300 m2<400 m2, so it is invalid change. When the initial type is forest land, the preset area threshold value is 600 m2. The area of plot two is 700 m2>600 m2, so it is also retained. Only the change plot that meets the minimum area standard will be confirmed as the land use update plot, as shown in Figure 3 .

[0109] In a second embodiment, a natural resource right registration plot linkage update system is provided, comprising the following modules:

[0110] The acquisition module is configured to acquire reference period and monitoring period remote sensing images and corresponding right registration plots in the same geographical range.

[0111] The calculation module is configured to calculate the difference between the stable ground object spectrum mean values in the reference period and the monitoring period images to obtain a system offset vector, and subtract the system offset vector from the spectrum change vector of each pixel from the reference period to the monitoring period to obtain a background corrected change vector. A band contribution matrix is constructed based on the ground object spectrum knowledge base, and the background corrected change vector is weighted to calculate the weighted change intensity and change direction angle.

[0112] The determination module is configured to establish a type transfer rule library, define a reference direction angle for each preset plot type conversion path, and establish a tolerance function with the weighted change intensity as the input and the direction angle tolerance as the output. For any pixel, when the pixel weighted change intensity exceeds the intensity threshold value corresponding to the initial type of the plot, and the absolute value of the difference between the change direction angle and the reference direction angle of at least one conversion path starting from the initial type is less than the direction angle tolerance determined by the tolerance function through the path and according to the weighted change intensity of the pixel, the pixel is determined as an effective change pixel.

[0113] The generation module is configured to aggregate all effective change pixels. If the area of a cluster formed by geographically adjacent effective change pixels exceeds a preset area threshold value, the original right registration plot area corresponding to the cluster is identified as a to-be-updated plot, and an update list is generated.

[0114] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0115] The above description is merely an embodiment of this application and is not intended to limit this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principle of this application should be included within the scope of the claims of this application.

Claims

1. A natural resource right registration plot linkage updating method, characterized in that, The method comprises the following steps: obtaining reference period and monitoring period remote sensing images of the same geographical range and corresponding right registration plot; calculating the difference between the mean spectrum of stable ground objects in the reference period and the monitoring period images to obtain a system offset vector, and subtracting the system offset vector from the spectral change vector of each pixel from the reference period to the monitoring period to obtain a background corrected change vector; constructing a band contribution matrix based on a ground object spectrum knowledge base, weighting the background corrected change vector, and calculating the weighted change intensity and change direction angle; establishing a type transfer rule base, defining a reference direction angle for each preset plot type conversion path, and establishing a tolerance function with the weighted change intensity as the input and the direction angle tolerance as the output; for any pixel, when the pixel weighted change intensity exceeds the intensity threshold corresponding to the initial type of the plot, and the absolute value of the difference between the change direction angle and the reference direction angle of at least one conversion path starting from the initial type is less than the direction angle tolerance determined by the tolerance function through the path and according to the weighted change intensity of the pixel, the pixel is determined as an effective change pixel; aggregating all effective change pixels, and if the area of a cluster formed by geographically adjacent effective change pixels exceeds a preset area threshold, the original right registration plot area corresponding to the cluster is identified as a plot to be updated, and an update list is generated; the band contribution matrix is constructed based on the ground object spectrum knowledge base, comprising: for a specific conversion path from initial type A to target type B, the standard spectral reflectance of A and B types of ground objects is extracted from the ground object spectrum knowledge base; The absolute value of the difference in reflectivity of the two types of ground objects at each spectral band i is calculated ; the difference values of all bands are normalized by maximum value to obtain the contribution degree of each band under the conversion path, thereby constructing the band contribution matrix; the weighting of the background corrected change vector, the calculation of the weighted change intensity and the change direction angle, comprising: for each preset plot type conversion path, the change components of each band in the background corrected change vector are multiplied by the band contribution degree corresponding to the conversion path to obtain the weighted change vector under the path; the Euclidean norm of the weighted change vector is calculated as the weighted change intensity under the path; the weighted change components of the near-infrared band and the red band are selected as two-dimensional coordinates, and the change direction angle under the path is calculated by using the two-parameter inverse tangent function.

2. The method of claim 1, wherein, the calculation of the difference between the mean spectrum of stable ground objects in the reference period and the monitoring period images to obtain a system offset vector, comprising: using the normalized difference vegetation index method, selecting the pixels with stable normalized difference vegetation index values and a change of less than 0.05 as stable ground objects in the reference period and the monitoring period images; respectively calculating the arithmetic mean of the spectral values of all stable ground objects in each band in the reference period and the monitoring period to obtain the reference period spectral mean vector and the monitoring period spectral mean vector; subtracting the reference period spectral mean vector from the monitoring period spectral mean vector to obtain the system offset vector.

3. The method of claim 1, wherein, the establishment of a tolerance function with the weighted change intensity as the input and the direction angle tolerance as the output, comprising: The tolerance function is a negative exponential function, which is where M is the weighted variation intensity, T(M) is the directional angle tolerance, A, k, C are normal numbers fitted by sample data, and the greater the weighted variation intensity, the smaller the directional angle tolerance.

4. The method of claim 1, wherein, The pixel is determined as an effective change pixel when the weighted change intensity of the pixel exceeds the intensity threshold corresponding to the initial type of the region in which the pixel is located, and the absolute value of the difference between the change direction angle and the reference direction angle of at least one conversion path starting from the initial type is less than the direction angle tolerance determined by the tolerance function of the path and according to the weighted change intensity of the pixel, and the pixel is determined as an effective change pixel, including: For a conversion path in which an initial type of a pixel is woodland and the pixel is converted into construction land, the weighted change intensity and the change direction angle of the pixel under the conversion path are calculated; If the weighted change intensity under the path exceeds the preset weighted change intensity threshold, and the absolute value of the difference between the change direction angle under the path and the reference direction angle of the conversion path from woodland to construction land is less than the direction angle tolerance calculated by substituting the weighted change intensity under the path into the tolerance function, the pixel is determined as an effective change pixel.

5. The method of claim 1, wherein, The effective change pixels are aggregated, including: A binary matrix with the same size as the remote sensing image is created, the positions of the effective change pixels are assigned as 1, and the remaining positions are assigned as 0; An eight-neighbor connected component labeling algorithm is used to scan the binary matrix, and the pixels with the value of 1 that are spatially connected to each other are divided into the same cluster.

6. The method of claim 1, wherein, The preset area threshold includes: When the initial type of the region to be updated is farmland, the preset area threshold is 400 square meters; When the initial type is woodland, the preset area threshold is 600 square meters.

7. A natural resource right registration plot linkage updating system, characterized in that, The method includes the following modules: An acquisition module is configured to acquire reference period and monitoring period remote sensing images and corresponding right-registered region polygons in the same geographical range; A calculation module is configured to calculate the difference between the mean values of stable ground object spectrums in the reference period and the monitoring period images to obtain a system offset vector, subtract the system offset vector from the spectral change vector of each pixel from the reference period to the monitoring period to obtain a background-corrected change vector, construct a band contribution matrix based on a ground object spectrum knowledge base, weight the background-corrected change vector, and calculate the weighted change intensity and the change direction angle; A determination module is configured to establish a type transition rule base, define a reference direction angle for each preset region type conversion path, and establish a tolerance function with the weighted change intensity as the input and the direction angle tolerance as the output; for any pixel, when the weighted change intensity of the pixel exceeds the intensity threshold corresponding to the initial type of the region in which the pixel is located, and the absolute value of the difference between the change direction angle and the reference direction angle of at least one conversion path starting from the initial type is less than the direction angle tolerance determined by the tolerance function of the path and according to the weighted change intensity of the pixel, the pixel is determined as an effective change pixel; A generation module is configured to aggregate all the effective change pixels, and if the area of a cluster formed by geographically adjacent effective change pixels exceeds a preset area threshold, the original right-registered region polygon corresponding to the cluster is identified as a region to be updated, and an update list is generated; The band contribution matrix is constructed based on the ground object spectrum knowledge base, including: For a specific conversion path from an initial type A to a target type B, the standard spectral reflectance of the A and B types of ground objects is extracted from the ground object spectrum knowledge base; The absolute value of the difference in reflectivity of the two types of ground objects at each spectral band i is calculated ; The difference values of all bands are subjected to maximum normalization to obtain the contribution degrees of each band under the conversion path, thereby constructing a band contribution degree matrix; The background corrected change vector is weighted to obtain a weighted change intensity and a change direction angle, including: For each preset patch type conversion path, the change component of each band in the background corrected change vector is multiplied by the band contribution degree corresponding to the conversion path to obtain a weighted change vector under the path; The Euclidean norm of the weighted change vector is calculated as the weighted change intensity under the path; The weighted change components of the near-infrared band and the red band are selected as two-dimensional coordinates, and the change direction angle under the path is calculated by using a two-parameter inverse tangent function.

8. The system of claim 7, wherein, The difference between the mean values of the stable ground object spectra in the reference period and the monitoring period images is calculated to obtain a system offset vector, including: The normalized difference vegetation index method is used to screen, in the reference period and the monitoring period images, the pixels with stable normalized difference vegetation index values and a change amount less than 0.05 as stable ground objects; The arithmetic mean values of the spectral values of all stable ground objects in each band in the reference period and the monitoring period are calculated to obtain a reference period spectral mean vector and a monitoring period spectral mean vector; The monitoring period spectral mean vector is subtracted from the reference period spectral mean vector to obtain a system offset vector.

Images