A land change monitoring method and system based on unmanned measurement

Abstract

The application relates to the technical field of land change monitoring, and discloses a land change monitoring method and system based on unmanned measurement, which comprises the following steps: dividing a target land area into a plurality of mutually non-overlapping grid units, generating a land change monitoring index of the grid units based on multispectral remote sensing data; adopting a multi-parameter trend identification algorithm to generate land change intensity of the grid units; combining the spatial structure of the grid units, adopting a spatial filtering mode to perform smoothing processing on the land change intensity of the grid units to form a visual change intensity monitoring graph. The target land area is divided into grid units, a land change monitoring index is calculated and a time sequence is constructed, a multi-parameter trend identification algorithm is used to evaluate the change intensity, a spatial filtering mode is combined to realize the visual change intensity monitoring graph, dominant change parameters are identified, and classification identification of land change types is performed, so that the spatial precision and time sequence stability of land change identification are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of image generation, especially to the field of land change monitoring, in particular to a land change monitoring method and system based on unmanned measurement. BACKGROUND

[0002] With the acceleration of urbanization, the continuous adjustment of agricultural structure and the increasing demand for fine management of land and resources, the types of land use change frequently, and land change monitoring has become an important problem that needs to be solved in the fields of agriculture, ecology, planning, environment, etc. Especially in typical application scenarios such as farmland protection, illegal land occupation supervision, urban-rural boundary expansion identification, farmland degradation monitoring, etc., the monitoring means are required to have high timeliness, high resolution, high spatial positioning accuracy, strong change recognition ability, and accurate type discrimination. Traditional manual field investigation is time-consuming and labor-intensive, with long cycle and lagging data update, which has been difficult to meet the on-site needs of rapid response and dynamic update.

[0003] In recent years, remote sensing technology has been widely applied in land change monitoring, especially the change detection technology based on multi-spectral images, which has played an important role in large-scale regional land use update. However, most of the existing technologies still rely on satellite remote sensing data, which are limited by long imaging cycle, limited resolution, fixed viewing angle, etc. In complex terrain, local changes are intense or frequent (such as mountainous farmland, urban-rural junction), there are problems such as insufficient accuracy and lagging response.

[0004] Unmanned aerial vehicle (UAV) remote sensing platform has become an ideal ground remote sensing measurement method due to its flexibility, high resolution, and low-altitude flight not affected by clouds. With multi-spectral cameras, inertial measurement units (IMUs), and GPS positioning systems, UAVs can collect multi-temporal, multi-angle, and high-precision remote sensing data, and achieve fine measurement in complex areas. The data based on UAV measurement has high spatial accuracy, strong timeliness, and low dependence on ground control points, which is especially suitable for rapid identification of local sensitive changes and boundary change monitoring tasks with high structural accuracy requirements.

[0005] In existing research, satellite remote sensing-based farmland change identification methods have been widely applied. For example, the public patent CN118608939B proposes an agricultural drought monitoring method based on multi-source remote sensing data. By analyzing the land evaluation coefficient of the target farmland sub-region, the accurate identification of the farmland cracking area is realized, and the crop planting selection is assisted and guided, which improves the spatio-temporal efficiency and data interpretation ability of drought monitoring. This technology shows good application prospects in the agricultural drought scene. However, its method essentially depends on fixed-phase satellite images and specific index differences, and the geometric accuracy requirement of the image is not high, without considering the disturbance of attitude disturbance and terrain undulation of remote sensing images on change identification. At the same time, this method only constructs a land evaluation index based on the spectral characteristics of ground objects and empirical factors, lacks systematic time series modeling ability and spatial change trend analysis mechanism, and is difficult to apply to more complex and fine-grained change scenarios in land change (such as slow degradation, illegal land occupation, boundary erosion, etc.).

[0006] To solve this problem, the present application provides a land change monitoring method and system based on unmanned measurement, which divides the target soil region into multiple grid cells by collecting multispectral remote sensing image data, periodically calculates and processes the land change monitoring index, and quantifies the land change intensity. SUMMARY

[0007] The present application provides a land change monitoring method and system based on unmanned measurement, which introduces IMU positioning data and GPS positioning data of unmanned aerial vehicles, combines digital elevation models, uses ray projection and intersection solving methods to accurately geometrically correct the original pixel coordinates, constructs multiple fusion monitoring indexes based on vegetation indexes and water indexes, monitors the land change intensity from multiple dimensions, introduces multi-index collaborative scoring and multi-parameter trend recognition algorithm combined with index change trend to generate the land change intensity of the grid cell, and then quantifies the change amplitude of the ground cover in the target soil region, describes its long-term evolution behavior, has stronger physical continuity and time interpretation, and realizes land change monitoring based on unmanned measurement.

[0008] To achieve the above purpose, the present application provides a land change monitoring method based on unmanned measurement, which comprises the following steps:

[0009] S1: Collecting multispectral remote sensing image data of the target soil region by using unmanned aerial vehicles equipped with multispectral cameras, and geometrically correcting the pixel coordinates in the multispectral remote sensing image data combined with the attitude information of the unmanned aerial vehicles, correcting the pixel coordinates to the ground coordinates of the target soil region, and obtaining the multispectral remote sensing data of the ground coordinates;

[0010] S2: Dividing the target soil region into multiple non-overlapping grid cells, and generating the land change monitoring index of the grid cell based on the multispectral remote sensing data;

[0011] S3: constructing a land change monitoring index time sequence of the grid cell, and generating a land change intensity of the grid cell by using a multi-parameter trend identification algorithm;

[0012] S4: combining the spatial structure of the grid cell, smoothing the land change intensity of the grid cell by using a spatial filtering method, constructing a visual change intensity monitoring map of the target land area, and extracting the change of a dominant change parameter in the land change monitoring index to classify and identify the land change behavior.

[0013] As a further improved method of the present application:

[0014] Optionally, multi-spectral remote sensing image data of the target land area is collected by using a multi-spectral camera-equipped unmanned aerial vehicle, and pixel coordinates in the multi-spectral remote sensing image data are geometrically corrected in combination with unmanned aerial vehicle attitude information, including:

[0015] The unmanned aerial vehicle is equipped with a high-precision GPS and IMU positioning system, and a multi-spectral camera;

[0016] According to the range of the target land area, a flight route plan covering the entire target land area is designed, and a periodic flight time interval is set, and the unmanned aerial vehicle periodically collects multi-spectral remote sensing image data, corresponding GPS positioning data and IMU positioning data according to the preset flight route plan, the GPS positioning data being position coordinates of the unmanned aerial vehicle, and the IMU positioning data being attitude information of the unmanned aerial vehicle;

[0017] The process of geometrically correcting the pixel coordinates is as follows:

[0018] Obtain pixel coordinates (u, v) and an intrinsic matrix K of the multi-spectral camera, and convert the pixel coordinates (u, v) into normalized coordinates (u 1 ,v 1 ,1) in a camera coordinate system:

[0019] [u 1 ,v 1 ,1] T =K -1 ·[u,v,1];

[0020] Wherein, T represents transposition;

[0021] In combination with the position coordinates and the attitude information of the unmanned aerial vehicle, a ray equation of the pixel coordinates (u, v) in a world coordinate system is obtained:

[0022] P(u,v;λ)=[u * ,v * ,h * ] T +λ·R(u,v)·[u 1v 1 ,1] T ;

[0023] wherein, P(u, v; l) represents a ray equation of the pixel coordinate (u, v) in the world coordinate system, l represents a ray iteration parameter, l > 0, and R(u, v) represents attitude information of the unmanned aerial vehicle when shooting a multispectral remote sensing image associated with the pixel coordinate (u, v),

[0024] ((u * ,v * ,h * ) represents a position coordinate of the unmanned aerial vehicle when shooting a multispectral remote sensing image associated with the pixel coordinate (u, v);

[0025] obtaining a digital elevation model of a target soil region, the digital elevation model being composed of ground coordinates of the target soil region and elevations associated with the ground coordinates, converting the elevations associated with the ground coordinates into an elevation function, wherein the elevation associated with the ground coordinate ((x, y) is z(x, y), and the elevation function converted from the elevation z(x, y) is Z = z((x, y), representing a straight line segment perpendicular to the ground with the ground coordinate ((x, y) as a starting point and a length of z(x, y);

[0026] gradually increasing the ray iteration parameter l until there is an intersection point between the ray equation P(u, v; l) and any elevation function, and taking the ground coordinate associated with the elevation function having the intersection point as a geometric correction result of the pixel coordinate (u, v), to obtain multispectral remote sensing data of the ground coordinate;

[0027] The spectral bands in the multispectral remote sensing data include a green light band, a red light band, a near-infrared band, and a short-infrared band.

[0028] Optionally, a land change monitoring index of a grid cell is generated based on the multispectral remote sensing data, including:

[0029] The land change monitoring index includes a basic monitoring index and a fusion monitoring index, the basic monitoring index includes a vegetation index, a soil index, and a water body index, and the fusion monitoring index of the grid cell is calculated based on the basic monitoring index, wherein the fusion monitoring index includes a land disturbance response index, a vegetation-water body interaction index, and a multi-band difference ratio vector index.

[0030] Optionally, the fusion monitoring index of the grid cell is calculated based on the basic monitoring index, including:

[0031] The land disturbance response index of the grid cell is calculated in the following manner:

[0032]

[0033] wherein, LDRI t (n,m) represents the land disturbance response index of the grid cell B(n,m) in the tth cycle, NDVI t (n,m) represents the vegetation index of the grid cell B(n,m) in the tth cycle, NDVI t-1 (n,m) represents the vegetation index of the grid cell B(n,m) in the (t-1)th cycle, Red t (n,m), NIR t (n,m), SWIR t (n,m) represents the reflectivity average of the red light band, near-infrared band and short-infrared band of the grid cell B(n,m) in the tth cycle, respectively, and the tth cycle is the cycle of the tth periodic collection of multispectral remote sensing image data by the unmanned aerial vehicle;

[0034] σ(NDVI1(n,m):NDVI t (n,m)) represents the standard deviation of the sequence (NDVI1(n,m), NDVI2((n,m),..., NDVI t ((n,m)), and NDVI1(n,m), NDVI2(n,m) represent the vegetation index of the grid cell B((n,m) in the 1st cycle and the 2nd cycle, respectively;

[0035] The grid cell B(n,m) represents the grid cell in the nth row and the mth column in the target soil region, N represents the number of grid cells in the horizontal direction of the target region, and M represents the number of grid cells in the vertical direction of the target region;

[0036] The vegetation-water body interaction index of the grid cell is calculated in the following manner:

[0037]

[0038] wherein, VWI t (n,m) represents the vegetation-water body interaction index of the grid cell B(n,m) in the tth cycle, NDWI t (n,m) represents the water body index of the grid cell B(n,m) in the tth cycle, and ∈ represents a control parameter, which is set to 0.1;

[0039] The multi-band difference ratio vector index of the grid cell is calculated in the following manner:

[0040]

[0041] wherein, F t (n,m) represents the multi-band difference ratio vector index of the grid cell B(n,m) in the tth cycle,

[0042] F (n,m) represents the land change monitoring index time series of the grid cell B(n,m), F ((n,m) represents the land change monitoring index of the grid cell B(n,m) in the t th period, F ((n,m) and F ((n,m) represent the land change monitoring index of the grid cell B(n,m) in the 1 st and 2 nd periods, respectively. t (n,m) represents the red-green band difference ratio, the near-infrared-red band difference ratio, and the short-infrared-near-infrared band difference ratio in (n,m).

[0043] Optionally, the land change monitoring index time series of the grid cell is constructed, including:

[0044] The vegetation index and the water body index in the basic monitoring index are removed to form the land change monitoring index time series of the grid cell:

[0045] Q 1→t (n,m) = [Q1(n,m), Q2(n,m),..., Q t (n,m)], n∈[1,N], m∈[1,M];

[0046] Q t (n,m) = [NDSI t (n,m), LDRI t (n,m), VWI t (n,m), F t (n,m)];

[0047] Wherein, Q 1→t (n,m) represents the land change monitoring index time series of the grid cell B(n,m), Q t ((n,m) represents the land change monitoring index of the grid cell B(n,m) in the t th period, F ((n,m) and F ((n,m) represent the land change monitoring index of the grid cell B(n,m) in the 1 st and 2 nd periods, respectively.

[0048] NDSI t (n,m) represents the soil index of the grid cell B(n,m) in the t th period, the grid cell B(n,m) represents the grid cell in the n th row and the m th column in the target land area, N represents the number of grid cells in the horizontal direction of the target area, and M represents the number of grid cells in the vertical direction of the target area.

[0049] LDRI t (n,m), VWI t (n,m), F t (n,m) represent the land disturbance response index, the vegetation-water body interaction index, and the multi-band difference ratio vector index of the grid cell B((n,m) in the t th period, respectively.

[0050] Optionally, a multi-parameter trend identification algorithm is used to identify the land change monitoring index time series, to generate the land change intensity of the grid cell, including:

[0051] The land change monitoring index time series Q 1→t (n,m) is normalized to obtain the normalized land change monitoring index time series

[0052] The normalized land change monitoring index time series The mean and standard deviation of the land change monitoring index in the sequence are calculated, and the last value of the land change monitoring index in the sequence is extracted. The last value is used to calculate the trend item of the land change monitoring index in the normalized land change monitoring index time series, and the trend item is converted into the dynamic weight of the land change monitoring index.

[0053] The degree of deviation of the last value of the land change monitoring index from the mean is calculated as the change item of the land change monitoring index in the normalized land change monitoring index time series .

[0054] The change item is weighted and fused based on the dynamic weight of the land change monitoring index, to obtain the land change intensity of the grid cell B(n,m), wherein the land change intensity of the grid cell B(n,m) in the tth cycle is H t (n,m).

[0055] Optionally, in combination with the spatial structure of the grid cell, a spatial filtering method is used to smooth the land change intensity of the grid cell, to form a visual change intensity monitoring map of the target soil region, including:

[0056] The neighborhood grid cells of the grid cell to be smoothed are obtained, as well as the land change intensity of the neighborhood grid cells.

[0057] In combination with the land change intensity of the neighborhood grid cells, a spatial filtering method is used to smooth the land change intensity of the grid cell to be smoothed, wherein the smoothing result of the land change intensity H t ((n,m) is

[0058] The visual change intensity monitoring map of the target soil region in the tth cycle is Map t .

[0059]

[0060] The visual change intensity monitoring map Map t is in the form of an N-row and M-column matrix, , and the visual change intensity monitoring map Mapt a matrix element in the nth row and the mth column;

[0061] If the smoothing result of the land change intensity of the grid cell is higher than the preset change intensity threshold, the change of the dominant change parameter in the land change monitoring index of the grid cell is extracted, and the land change behavior of the grid cell is classified and identified, wherein the dominant change parameter in the land change monitoring index includes a vegetation index, a water body index, a soil index and a land disturbance response index.

[0062] To solve the above problems, the application further provides a land change monitoring system based on unmanned measurement, which comprises a data acquisition device and a change intensity monitoring module.

[0063] The data acquisition device is used for acquiring multispectral remote sensing image data of a target soil area by using a drone equipped with a multispectral camera, and geometrically correcting pixel coordinates in the multispectral remote sensing image data in combination with attitude information of the drone to correct the pixel coordinates to ground coordinates of the target soil area, so as to obtain multispectral remote sensing data of the ground coordinates.

[0064] The index monitoring module is used for dividing the target soil area into a plurality of non-overlapping grid cells, and generating a land change monitoring index of the grid cells based on the multispectral remote sensing data.

[0065] The change intensity monitoring module is used for constructing a time sequence of the land change monitoring index of the grid cells, generating a land change intensity of the grid cells by using a multi-parameter trend identification algorithm, combining a spatial structure of the grid cells, and smoothing the land change intensity of the grid cells by using a spatial filtering method to construct a visual change intensity monitoring map of the target soil area, and extracting the change of the dominant change parameter in the land change monitoring index to classify and identify the land change behavior.

[0066] To achieve the land change monitoring method based on unmanned measurement as described above.

[0067] To solve the above problems, the application further provides an electronic device, which comprises:

[0068] a memory for storing at least one instruction;

[0069] a communication interface for realizing communication of the electronic device; and

[0070] a processor for executing the instruction stored in the memory to realize the land change monitoring method based on unmanned measurement as described above.

[0071] To solve the above problems, the application further provides a computer readable storage medium, wherein at least one instruction is stored in the computer readable storage medium, and the at least one instruction is executed by a processor in an electronic device to implement the above-mentioned land change monitoring method based on unmanned measurement.

[0072] Compared with the prior art, the application provides a land change monitoring method and system based on unmanned measurement, which has the following beneficial effects:

[0073] Firstly, the external parameters (i.e. position and attitude) of the sensor are obtained by using high-precision GPS / IMU positioning data recorded by the unmanned aerial vehicle, the ray projection and ground intersection of the remote sensing image are solved in combination with the digital elevation model, the pixel points under the visual angle of the sensor can be accurately mapped to the geographical spatial coordinates, and the geometric distortion of the image caused by the terrain undulation can be accurately corrected. Through the process, the remote sensing image shot at an angle or from multiple angles can be converted into an orthographic image with unified scale and accurate spatial reference. Compared with the traditional approximate geometric correction or plane assumption projection method, the method can significantly improve the image positioning accuracy in complex terrain areas such as mountains, hills and wetlands, and can avoid false matching and false change detection caused by terrain obstruction or visual angle difference. The use of attitude information to correct the ray direction can accurately obtain the pixel visual line direction and eliminate the tilt error. In combination with the digital elevation model, the true intersection of the pixel and the ground surface can be accurately obtained, and the terrain vertical disturbance can be corrected.

[0074] Meanwhile, the application introduces a land disturbance response index in the process of land change monitoring, the land disturbance response index effectively identifies the disturbance behavior of the land surface by fusing the time variation difference of the vegetation index, the historical stability standard deviation and the multispectral reflectance ratio, wherein the first half of the land disturbance response index calculation formula measures the short-term vegetation index change amplitude by using the normalized difference value, and the historical vegetation index standard deviation is used for normalization to suppress the false response in the natural fluctuation background. The second half combines the near-infrared, red and short-wave infrared reflectance ratios to strengthen the disturbance amplification response of non-vegetation targets such as bare land, excavation area and construction area, thereby enhancing the sensitivity to non-natural disturbance. The land disturbance response index can not only distinguish between the dramatic changes caused by human activities and the seasonal changes, but also has good cross-time adaptability and regional migration ability. Compared with the traditional difference vegetation index or change vector method, the land disturbance response index has stronger stability and significant abnormal response in complex terrain, variable climate or mixed ground type scenes, and is particularly suitable for application scenes such as farmland degradation monitoring, grassland reclamation discrimination and urban expansion identification. BRIEF DESCRIPTION OF DRAWINGS

[0075] Figure 1 A flowchart of a land change monitoring method based on unmanned measurement is provided for an embodiment of the application.

[0076] The objectives, functional characteristics and advantages of the present application will be further illustrated in conjunction with the embodiments with reference to the accompanying drawings. DETAILED DESCRIPTION

[0077] It should be understood that the specific embodiments described herein are merely illustrative of the present application and do not limit the present application.

[0078] Embodiments of the present application provide a land change monitoring method based on unmanned measurement. The execution subject of the land change monitoring method based on unmanned measurement includes but is not limited to at least one of electronic devices such as a server, a terminal, etc. which can be configured to execute the method provided by the embodiments of the present application. In other words, the land change monitoring method based on unmanned measurement can be executed by software or hardware installed in a terminal device or a server device, and the software can be a blockchain platform. The server includes but is not limited to a single server, a server cluster, a cloud server, or a cloud server cluster, etc.

[0079] Reference Figure 1 , Embodiment 1 of the present application is:

[0080] S1: Collecting multi-spectral remote sensing image data of a target land area by using a multi-spectral camera-equipped unmanned aerial vehicle, and performing geometric correction on pixel coordinates in the multi-spectral remote sensing image data in combination with attitude information of the unmanned aerial vehicle to correct the pixel coordinates to ground coordinates of the target land area, thereby obtaining multi-spectral remote sensing data of the ground coordinates.

[0081] Collecting multi-spectral remote sensing image data of a target land area by using a multi-spectral camera-equipped unmanned aerial vehicle, and performing geometric correction on pixel coordinates in the multi-spectral remote sensing image data in combination with attitude information of the unmanned aerial vehicle, including:

[0082] The unmanned aerial vehicle is equipped with a high-precision GPS and an IMU positioning system, as well as a multi-spectral camera;

[0083] According to the range of the target land area, a flight route plan covering the entire target land area is designed, and a periodic flight time interval is set. The unmanned aerial vehicle periodically collects multi-spectral remote sensing image data, corresponding GPS positioning data, and IMU positioning data according to the preset flight route plan. The GPS positioning data is the position coordinates of the unmanned aerial vehicle, and the IMU positioning data is the attitude information of the unmanned aerial vehicle;

[0084] The process of geometric correction on the pixel coordinates is as follows:

[0085] Obtaining pixel coordinates (u, v) and an intrinsic matrix K of the multi-spectral camera, and converting the pixel coordinates (u, v) into normalized coordinates (u 1 ,v 1 ,1) in a camera coordinate system;

[0086] [u 1 ,v 1 ,1] T =K -1 ·[u,v,1];

[0087] wherein T represents transpose;

[0088] Combined with the position coordinates and attitude information of the unmanned aerial vehicle, a ray equation of the pixel coordinates (u, v) in the world coordinate system is obtained:

[0089] P(u, v; l) = [u * ,v * ,h * ] T + l·R(u, v)·[u 1 ,v 1 ,1] T ;

[0090] wherein P(u, v; l) represents the ray equation of the pixel coordinates (u, v) in the world coordinate system, l represents a ray iteration parameter, l > 0, R(u, v) represents the attitude information of the unmanned aerial vehicle when shooting the multi-spectral remote sensing image associated with the pixel coordinates (u, v),

[0091] (u * ,v * ,h * ) represents the position coordinates of the unmanned aerial vehicle when shooting the multi-spectral remote sensing image associated with the pixel coordinates (u, v);

[0092] A digital elevation model of the target soil area is obtained, the digital elevation model is composed of ground coordinates of the target soil area and elevations associated with the ground coordinates, and the elevations associated with the ground coordinates are converted into an elevation function, wherein the elevation associated with the ground coordinates (x, y) is z(x, y), the converted elevation function of the elevation z(x, y) is Z = z((x, y), and a straight line segment perpendicular to the ground with the ground coordinates (x, y) as the starting point and the length of z(x, y) as the length;

[0093] The ray iteration parameter l is gradually increased until the ray equation P(u, v; l) and any elevation function have an intersection point, and the ground coordinates associated with the elevation function having the intersection point are taken as the geometric correction result of the pixel coordinates (u, v), and multi-spectral remote sensing data of the ground coordinates are obtained;

[0094] The spectral bands in the multi-spectral remote sensing data include green light bands, red light bands, near-infrared bands, and short-infrared bands.

[0095] S2: dividing the target soil region into a plurality of non-overlapping grid cells, and generating a land change monitoring index of the grid cells based on the multispectral remote sensing data.

[0096] The generating of the land change monitoring index of the grid cells based on the multispectral remote sensing data comprises:

[0097] The land change monitoring index comprises a basic monitoring index and a fusion monitoring index, the basic monitoring index comprises a vegetation index, a soil index and a water index, and the fusion monitoring index of the grid cells is calculated based on the basic monitoring index, wherein the fusion monitoring index comprises a land disturbance response index, a vegetation-water interaction index and a multi-band difference ratio vector index.

[0098] Specifically, the basic monitoring index of the grid cells is calculated in the following manner:

[0099]

[0100] n∈[1,N],m∈[[1,M];

[0101] wherein NDVI t ((n,m),NDSI t (n,m),NDWI t (n,m) represent the vegetation index, the soil index and the water index of the grid cell B(n,m) in the tth period respectively, the grid cell B(n,m) represents the grid cell in the nth row and the mth column in the target soil region, N represents the number of grid cells in the horizontal direction of the target region, M represents the number of grid cells in the vertical direction of the target region, Green t (n,m),Red t (n,m),NIR t ((n,m),SWIR t ((n,m) represent the average reflectivity of the green band, the red band, the near-infrared band and the short-infrared band of the grid cell B(n,m) in the tth period respectively, and the multispectral remote sensing data is the reflectivity of the ground coordinates in different spectral bands.

[0102] The tth period is the period of the tth periodic collection of multispectral remote sensing image data by the unmanned aerial vehicle, and t represents the total number of periods of the current collection of multispectral remote sensing image data by the unmanned aerial vehicle.

[0103] The fusion monitoring index of the grid cells is calculated based on the basic monitoring index, comprising:

[0104] The land disturbance response index of the grid cells is calculated in the following manner:

[0105]

[0106] wherein, LDRI t (n,m) represents the land disturbance response index of the grid cell B(n,m) in the tth cycle, NDVI t (n,m) represents the vegetation index of the grid cell B(n,m) in the tth cycle, NDVI t-1 (n,m) represents the vegetation index of the grid cell B(n,m) in the t-1th cycle, Red t (n,m), NIR t (n,m), SWIR t (n,m) represents the reflectivity average of the red light band, near-infrared band and short-infrared band of the grid cell B(n,m) in the tth cycle, respectively, and the tth cycle is the period of the tth periodic collection of multispectral remote sensing image data by the UAV;

[0107] σ(NDVI1(n,m):NDVI t (n,m)) represents the standard deviation of the sequence (NDVI1(n,m), NDVI2((n,m),..., NDVI t ((n,m)), and NDVI1(n,m), NDVI2(n,m) represent the vegetation index of the grid cell B((n,m) in the 1st cycle and the 2nd cycle, respectively;

[0108] The grid cell B(n,m) represents the grid cell in the nth row and the mth column in the target soil region, N represents the number of grid cells in the horizontal direction of the target region, and M represents the number of grid cells in the vertical direction of the target region;

[0109] Specifically, although natural changes (such as seasonal greening and post-rain green recovery) can also cause changes in the vegetation index, they usually have periodicity and continuity. The land disturbance response index is obtained by dividing by the historical standard deviation, which suppresses the severe but predictable natural fluctuations. For sudden human disturbances such as excavation and reclamation, the vegetation index usually decreases sharply but has no historical fluctuation basis, so the land disturbance response index is significantly amplified;

[0110] And in the enhancement term of the land disturbance response index, for bare land or construction sites, the reflectivity of the near-infrared band is low, and the reflectivity of the red light band and the short-infrared band is high, which makes the enhancement term decrease significantly, especially when the bare soil is exposed after the vegetation is removed, which makes the land disturbance response index have enhanced perception ability for non-green disturbance;

[0111] The vegetation-water body interaction index of the grid cell is calculated as follows:

[0112]

[0113] wherein, VWI t (n,m) represents the vegetation-water interaction index of the grid cell B(n,m) in the tth period, NDWI t (n,m) represents the water index of the grid cell B(n,m) in the tth period, and ∈ represents a control parameter, which is set to 0.1.

[0114] As a preferred embodiment of the present application, the vegetation-water interaction index is obtained by coupling the vegetation index and the water index, which depicts the coupling relationship and dynamic change of the vegetation coverage and water distribution in space, and is suitable for remote sensing monitoring of the interlaced areas such as wetlands, rice fields, and irrigated farmlands. The vegetation-water interaction index is designed in a product form to strengthen the response intensity of the area where vegetation and water coexist, and NDVI t (n,m) + NDWI t ((n,m) + ∈ realizes the balance of normalization and enhances the numerical stability in the low coverage area, thereby effectively avoiding the false identification caused by the high single index. Compared with the traditional independent vegetation index or water index, the vegetation-water interaction index is more sensitive to the slight change of the vegetation-water coexistence state, and can timely reflect phenomena such as the change of the water surface and seedling ratio in the rice planting period, the boundary change of the returning farmland to wetland project, and the seasonal water level fluctuation.

[0115] The multi-band difference ratio vector index of the grid cell is calculated in the following manner:

[0116]

[0117] wherein, F t (n,m) represents the multi-band difference ratio vector index of the grid cell B(n,m) in the tth period,

[0118] respectively represent the red-green band difference ratio, the near-infrared-red band difference ratio, and the short-infrared-near-infrared band difference ratio in the multi-band difference ratio vector index F t (n,m);

[0119] In the embodiment, compared with single band or traditional index, the multi-band difference ratio vector index emphasizes the relative change trend between bands, can effectively eliminate the influence of light intensity, sensor gain, etc., and improve the consistency of data at different times or devices. Taking red and green, near-infrared and red, and short-wave infrared and near-infrared as examples, the difference ratio can reflect the change of vegetation red edge, water change and small change of soil structure, which is helpful to identify fine-grained ground state change. The dimensions of the multi-band difference ratio vector index jointly represent the gradient change trend of the spectral curve, and have more sensitive response ability to vegetation growth stage, water drying process, land disturbance degree, etc. Compared with traditional vegetation index, soil index and water index, the multi-band difference ratio vector index has stronger dimension expression ability and information gain.

[0120] S3: constructing a land change monitoring index time series of the grid cell, and generating a land change intensity of the grid cell by using a multi-parameter trend identification algorithm.

[0121] Constructing a land change monitoring index time series of the grid cell comprises:

[0122] Removing the vegetation index and the water index in the basic monitoring index to constitute the land change monitoring index time series of the grid cell.

[0123] Q 1→t (n,m)=[Q1(n,m),Q2(n,m),...,Q t (n,m)],n∈[1,N],m∈[[1,M];

[0124] Q t (n,m)=[NDSI t (n,m),LDRI t (n,m),VWI t (n,m),F t (n,m)];

[0125] Wherein, Q 1→t (n,m) represents a land change monitoring index time series of the grid cell B(n,m), Q t ((n,m) represents a land change monitoring index of the grid cell B(n,m) in the t th period, Q1((n,m) and Q2((n,m) represent land change monitoring indexes of the grid cell B(n,m) in the 1 st and 2 nd periods, respectively.

[0126] NDSI t(n,m) represents the soil index of the grid cell B(n,m) in the t-th period, the grid cell B(n,m) represents the grid cell in the n-th row and the m-th column in the target soil region, N represents the number of grid cells in the horizontal direction of the target region, and M represents the number of grid cells in the vertical direction of the target region;

[0127] LDRI t (n,m),VWI t (n,m),F t (n,m) represent the land disturbance response index, the vegetation-water body interaction index, and the multi-band difference ratio vector index of the grid cell B((n,m) in the t-th period, respectively.

[0128] The multi-parameter trend identification algorithm is used to identify the land change monitoring index time series to generate the land change intensity of the grid cell, including:

[0129] The land change monitoring index in the land change monitoring index time series Q 1→t (n,m) is normalized to obtain a normalized land change monitoring index time series

[0130] Specifically,

[0131] Wherein:

[0132] Q1(n,m), Q2(n,m),..., Q t (n,m) are the normalized processing results of NDSI

[0133] NDSI t (n,m), LDRI t ((n,m), VWI t (n,m), F t (n,m) are the normalized processing results of NDSI

[0134] NDSI t (n,m), LDRI t (n,m), VWI t (n,m), F t (n,m) are normalized in a minimum-maximum normalization manner;

[0135] The normalized land change monitoring index time series The mean and the standard deviation of the land change monitoring index are calculated, and the last value of the land change monitoring index in the sequence is extracted, and the last value is used to calculate the trend item of the land change monitoring index in the normalized land change monitoring index time sequence, and the trend item is converted into the dynamic weight of the land disturbance response index in the land change monitoring index;

[0136] In the embodiment, the dynamic weight of the land disturbance response index in the land change monitoring index is calculated in the following manner:

[0137] Ω={NDSI, LDRI, VWI, F};

[0138] wherein, represents the dynamic weight of the land disturbance response index in the normalized land change monitoring index time sequence, Ω represents an index set composed of a soil index, a land disturbance response index, a vegetation-water interaction index, and a multi-band difference ratio vector index, c e Ω, c represents any index in the index set Ω, represents the normalized processing result of the index c of the grid cell B(n, m) in the tth period, exp(·) represents an exponential function with a natural constant as a base, and ||·||2 represents an L2 norm;

[0139] The higher the dynamic weight is, the higher the change range of the land change monitoring index in multiple periods is;

[0140] The degree of deviation of the last value of the land change monitoring index from the mean is calculated as a change item of the land change monitoring index in the normalized land change monitoring index time sequence ;

[0141] Specifically, the change item of the land disturbance response index is:

[0142]

[0143] wherein, represents the change item of the land disturbance response index in the normalized land change monitoring index time sequence, represents the mean of the land disturbance response index in the normalized land change monitoring index time sequence, represents the standard deviation of the land disturbance response index in the normalized land change monitoring index time sequence;

[0144] The change item is weighted and fused based on the dynamic weight of the land change monitoring index to obtain the land change intensity of the grid cell, wherein the land change intensity of the grid cell B(n, m) in the tth period is H t ​​​(n, m).

[0145] S4: combined with the spatial structure of the grid unit, the land change intensity of the grid unit is smoothed by using spatial filtering method, the visual change intensity monitoring map of the target land area is constructed, and the change of the dominant change parameter in the land change monitoring index is extracted to classify and identify the land change behavior.

[0146] Combined with the spatial structure of the grid unit, the land change intensity of the grid unit is smoothed by using spatial filtering method, the visual change intensity monitoring map of the target land area is constructed, and the change of the dominant change parameter in the land change monitoring index is extracted to classify and identify the land change behavior.

[0147] The neighborhood grid unit of the grid unit to be smoothed is obtained, and the land change intensity of the neighborhood grid unit is obtained; as an embodiment of the present application, the neighborhood grid unit of the grid unit to be smoothed is the grid unit in the 3*3 grid unit region centered on the grid unit to be smoothed;

[0148] Combined with the land change intensity of the neighborhood grid unit, the land change intensity of the grid unit to be smoothed is smoothed, and the smoothing formula is:

[0149]

[0150] Wherein, H B (n, m) represents the land change intensity of the grid unit B(n, m) in the tth period t ((n, m) is the smoothing result, β represents the filtering coefficient, and β is set to 0.8, H B (n, m) represents the land change intensity of the grid unit B(n, m) in the tth period

[0151] The grid unit B(n, m) represents the grid unit in the n th row and the m th column in the target land area, N represents the number of grid units in the horizontal direction of the target area, and M represents the number of grid units in the vertical direction of the target area;

[0152] The visual change intensity monitoring map of the target land area in the tth period is Map t :

[0153]

[0154] Wherein, the visual change intensity monitoring map Map t is in the form of N rows and M columns, H B (n, m) represents the land change intensity of the grid unit B(n, m) in the tth period t is the matrix element in the n th row and the m th column in the visual change intensity monitoring map Map t .

[0155] If the smoothing result of the land change intensity of the grid cell is higher than a preset change intensity threshold (set to 0.6), the change of the dominant change parameter in the land change monitoring index of the grid cell is extracted, and the land change behavior of the grid cell is classified and identified, wherein the dominant change parameter in the land change monitoring index includes a vegetation index, a water body index, a soil index, and a land disturbance response index.

[0156] Specifically, if the vegetation index decreases rapidly and the land disturbance response index rises, the land change behavior of the grid cell is construction occupation, that is, the land becomes construction land; if the vegetation index decreases slowly and the soil index rises, the land change behavior of the grid cell is land degradation; if the water body index continuously decreases, the land change behavior of the grid cell is water body shrinkage or water source loss; and if the water body index and the vegetation index both rise, the land change behavior of the grid cell is green space expansion or water bank greening.

[0157] Embodiment 2:

[0158] A land change monitoring system based on unmanned measurement includes a data acquisition device and a change intensity monitoring module.

[0159] The data acquisition device is configured to acquire multi-spectral remote sensing image data of a target soil region by using a multi-spectral camera-equipped unmanned aerial vehicle, and to perform geometric correction on pixel coordinates in the multi-spectral remote sensing image data in combination with attitude information of the unmanned aerial vehicle, so as to correct the pixel coordinates to ground coordinates of the target soil region and obtain multi-spectral remote sensing data of the ground coordinates.

[0160] The index monitoring module is configured to divide the target soil region into a plurality of non-overlapping grid cells, and to generate a land change monitoring index of the grid cells based on the multi-spectral remote sensing data.

[0161] The change intensity monitoring module is configured to construct a time sequence of the land change monitoring index of the grid cells, to generate a land change intensity of the grid cells by using a multi-parameter trend identification algorithm, to perform smoothing processing on the land change intensity of the grid cells by using a spatial filtering method in combination with a spatial structure of the grid cells, to construct a visual change intensity monitoring map of the target soil region, and to extract the change of the dominant change parameter in the land change monitoring index and classify and identify the land change behavior.

[0162] The land change monitoring system is configured to implement the land change monitoring method based on unmanned measurement as described in Embodiment 1.

[0163] It should be understood that the embodiments are only for illustration and are not limited in the scope of the patent application by the structure.

[0164] It should be noted that the above-mentioned embodiment serial numbers of the present application are only for description, not representing the advantages and disadvantages of the embodiments. Also, the terms "comprising", "containing" or any other variants thereof in this document are intended to cover the non-exclusive inclusion, so that the process, device, article or method including a series of elements not only includes those elements, but also includes other elements not explicitly listed or inherent to such process, device, article or method. Without more limitations, the element defined by the statement "comprising a" does not exclude the presence of other identical elements in the process, device, article or method including the element.

[0165] Through the above description of the embodiments, those skilled in the art can clearly understand that the above-mentioned embodiment methods can be realized by means of software and necessary general hardware platforms, of course, they can also be realized by hardware, but in many cases the former is a better embodiment. Based on such understanding, the technical solutions of the present application can be embodied in the form of a software product, which is stored in a storage medium (such as ROM / RAM, magnetic disk, optical disk) as described above, and includes a plurality of instructions for making a terminal device (which can be a mobile phone, computer, server, or network device, etc.) execute the methods described in various embodiments of the present application.

[0166] The above is only the preferred embodiment of the present application, and does not limit the patent scope of the present application, and any equivalent structure or equivalent flow transformation made by using the content of the present application specification and drawings, or directly or indirectly applied to other related technical fields, are also included in the patent protection scope of the present application.

Claims

1. A land change monitoring method based on unmanned surveying, characterized in that, The method includes: S1: Use a drone equipped with a multispectral camera to collect multispectral remote sensing image data of the target land area, and combine the drone's attitude information to perform geometric correction on the pixel coordinates in the multispectral remote sensing image data, correcting the pixel coordinates to the ground coordinates of the target land area, and obtaining multispectral remote sensing data of the ground coordinates. S2: Divide the target land area into multiple non-overlapping grid units and generate land change monitoring indices for the grid units based on multispectral remote sensing data; The land change monitoring index includes a land disturbance response index. This index effectively identifies land surface disturbance behavior by integrating the temporal variation differences of vegetation indices, historical stability standard deviations, and multispectral reflectance ratios. The first part of the land disturbance response index calculation formula measures the short-term vegetation index variation using normalized differences and normalizes it using historical vegetation index standard deviations to suppress spurious responses in the context of natural fluctuations. The second part combines near-infrared, red, and short-wave infrared reflectance ratios to amplify the response to disturbances to non-vegetated targets and enhance sensitivity to unnatural disturbances. S3: Construct a time series of land change monitoring indices for grid units, and use a multi-parameter trend recognition algorithm to generate the land change intensity of grid units; The time series of the land change monitoring index The land change monitoring index in the data was normalized to obtain the normalized land change monitoring index time series. ; Calculation of the time series of the normalized land change monitoring index The mean and standard deviation of the land change monitoring index are calculated, and the last value of the land change monitoring index in the series is extracted. The trend term of the land change monitoring index in the normalized land change monitoring index time series is calculated using the last value, and the trend term is converted into the dynamic weight of the land change monitoring index. The degree to which the last value of the land change monitoring index deviates from the mean is calculated and used as the land change monitoring index in the normalized land change monitoring index time series. Changes in the document; The land change intensity of the grid unit is obtained by weighting and fusing the change items based on the dynamic weights of the land change monitoring index. The land change intensity in the t-th period is ; grid cell This represents the grid cell in the nth row and mth column of the target land area; S4: Combining the spatial structure of the grid cells, the land change intensity of the grid cells is smoothed by spatial filtering to form a visual change intensity monitoring map of the target land area, and the changes of the dominant change parameters in the land change monitoring index are extracted to classify and identify land change behavior. The dominant change parameters include vegetation index, water index, soil index, and land disturbance response index.

2. The land change monitoring method based on unmanned measurement as described in claim 1, characterized in that, Multispectral remote sensing image data of the target land area is acquired using a drone equipped with a multispectral camera. Geometric correction of the pixel coordinates in the multispectral remote sensing image data is then performed using the drone's attitude information, including: The drone is equipped with a high-precision GPS and IMU positioning system, as well as a multispectral camera; Based on the target land area, a flight route plan covering the entire target land area is designed, and a periodic flight time interval is set. The UAV periodically collects multispectral remote sensing image data, as well as corresponding GPS positioning data and IMU positioning data according to the preset flight route plan. The GPS positioning data is the position coordinates of the UAV, and the IMU positioning data is the attitude information of the UAV. The process for geometrically correcting pixel coordinates is as follows: Get pixel coordinates And the intrinsic parameter matrix K of the multispectral camera, which will store the pixel coordinates. Convert to normalized coordinates in camera coordinate system : ； Where T represents transpose; By combining the drone's position coordinates and attitude information, pixel coordinates are obtained. Ray equation in world coordinate system: ； in, Represents pixel coordinates Ray equations in the world coordinate system Represents the ray iteration parameters. , Indicates the pixel coordinates captured by the drone. The attitude information of the associated multispectral remote sensing image. Indicates the pixel coordinates captured by the drone. The location coordinates of the associated multispectral remote sensing image; A digital elevation model (DEM) of the target land area is obtained. This DEM consists of the ground coordinates of the target land area and the associated elevations. The elevations associated with the ground coordinates are then converted into an elevation function. The ground coordinates... The associated elevation is Elevation The converted elevation function is , indicating ground coordinates Starting point, length is A straight line segment perpendicular to the ground; Gradually increase the ray iteration parameters Until the ray equation There are intersection points with any elevation function, and the ground coordinates associated with the elevation function that has intersection points are used as pixel coordinates. The geometric correction results yielded multispectral remote sensing data of the ground coordinates; The spectral bands in the multispectral remote sensing data include the green band, red band, near-infrared band, and short-infrared band.

3. The land change monitoring method based on unmanned measurement as described in claim 2, characterized in that, Land change monitoring indices are generated based on multispectral remote sensing data to form grid cells, including: The land change monitoring index includes a basic monitoring index and a fusion monitoring index. The basic monitoring index includes a vegetation index, a soil index, and a water body index. The fusion monitoring index of the grid unit is calculated based on the basic monitoring index. The fusion monitoring index includes a land disturbance response index, a vegetation-water interaction index, and a multi-band difference ratio vector index.

4. The land change monitoring method based on unmanned measurement as described in claim 1, characterized in that, The fusion monitoring index for grid cells is calculated based on the basic monitoring index, including: The land disturbance response index of a grid cell is calculated as follows: ； in, Represents grid cells The land disturbance response index in the t-th period, Represents grid cells The vegetation index in the t-th period, Represents grid cells The vegetation index in the (t-1)th period, The grid cells are represented in sequence. The average reflectance of the red band, near-infrared band and short-infrared band in the t-th period, where the t-th period is the period of the t-th periodic acquisition of multispectral remote sensing image data by the UAV. Represents a sequence standard deviation The grid cells are represented in sequence. Vegetation indices in the first and second cycles; Grid cells This represents the grid cell in the nth row and mth column of the target land area, where N represents the number of grid cells in the horizontal direction of the target area, and M represents the number of grid cells in the vertical direction of the target area. The vegetation-water interaction index of the grid cell is calculated as follows: ； in, Represents grid cells The vegetation-water interaction index in the t-th period Represents grid cells Water index in the t-th period Indicates control parameters, settings It is 0.1; The calculation method for the multi-band inter-band difference ratio vector exponent of the grid cell is as follows: ； ； ； ； in, Represents grid cells The multi-band inter-difference ratio vector exponent in the t-th period These represent the vector exponents of the inter-band difference ratio, respectively. The ratio of red-green band difference, the ratio of near-infrared band difference, and the ratio of short-infrared band difference are measured.

5. A land change monitoring method based on unmanned surveying as described in claim 1, characterized in that, Constructing a time series of land change monitoring indices for grid units, including: After removing the vegetation index and water body index from the basic monitoring indices, the land change monitoring index time series of the grid unit is constructed as follows: ； ； in, Represents grid cells The time series of land change monitoring index Represents grid cells The land change monitoring index in the t-th period, The grid cells are represented in sequence. Land change monitoring index in the first and second cycles; Represents grid cells Soil index in the t-th period, grid cell This represents the grid cell in the nth row and mth column of the target land area, where N represents the number of grid cells in the horizontal direction of the target area, and M represents the number of grid cells in the vertical direction of the target area. The grid cells are represented in sequence. The land disturbance response index, vegetation-water interaction index, and inter-band difference ratio vector index in the t-th period.

6. A land change monitoring method based on unmanned surveying as described in claim 1, characterized in that, Based on the spatial structure of the grid cells, spatial filtering is used to smooth the land change intensity of the grid cells, thus constructing a visual land change intensity monitoring map of the target land area, including: Obtain the neighboring grid cells of the grid cell to be smoothed, and the land change intensity of the neighboring grid cells; Based on the land change intensity of neighboring grid cells, a spatial filtering method is used to smooth the land change intensity of the grid cells to be smoothed, where the land change intensity... The smoothing result is ; The visualization change intensity monitoring map of the target land area in the t-th period is as follows. : ； Among them, the visual change intensity monitoring map It is an N-row, M-column matrix. Visualized change intensity monitoring map The matrix element in the nth row and mth column; If the smoothing result of the land change intensity of the grid cell is higher than the preset change intensity threshold, the changes of the dominant change parameters in the land change monitoring index of the grid cell are extracted to classify and identify the land change behavior of the grid cell.

7. A land change monitoring system based on unmanned surveying, characterized in that, The land change monitoring system based on unmanned measurement includes a data acquisition device and a change intensity monitoring module: The data acquisition device is used to collect multispectral remote sensing image data of the target land area using a drone equipped with a multispectral camera. The pixel coordinates in the multispectral remote sensing image data are geometrically corrected by combining the attitude information of the drone, and the pixel coordinates are corrected to the ground coordinates of the target land area to obtain multispectral remote sensing data of the ground coordinates. The index monitoring module is used to divide the target land area into multiple non-overlapping grid units and generate a land change monitoring index for the grid units based on multispectral remote sensing data. The land change intensity monitoring module is used to construct a time series of land change monitoring index for grid units, generate land change intensity for grid units using a multi-parameter trend recognition algorithm, and smooth the land change intensity of grid units by combining the spatial structure of grid units and using spatial filtering to form a visual land change intensity monitoring map of the target land area. It also extracts the changes of the dominant change parameters in the land change monitoring index and classifies and identifies land change behavior. To achieve the land change monitoring method based on unmanned measurement as described in any one of claims 1-6.

Images