Land change monitoring method and system based on unmanned measurement

Abstract

The invention relates to the technical field of land change monitoring, and discloses a land change monitoring method and system based on unmanned measurement, and the method comprises the steps: dividing a target land region into a plurality of grid units which do not coincide with each other, and generating land change monitoring indexes of the grid units based on multispectral remote sensing data; a multi-parameter trend identification algorithm is adopted to generate land change intensity of the grid units; and in combination with the space structure of the grid unit, smoothing the land change intensity of the grid unit by adopting a spatial filtering mode to form a visual change intensity monitoring graph. According to the method, a target land area is divided into grid units, land change monitoring indexes are calculated, a sequential sequence is constructed, change intensity is evaluated by using a multi-parameter trend identification algorithm, a visual change intensity monitoring graph is realized in combination with spatial filtering, dominant change parameters are identified, and land change types are classified and identified. And the spatial precision and the time sequence stability of land change identification are improved.

Description

Technical Field

[0001] The present invention relates to the field of image generation technology, in particular to the field of land change monitoring, and specifically to a land change monitoring method and system based on unmanned measurement. Background Art

[0002] With the acceleration of urbanization, the continuous adjustment of agricultural structure and the increasing demand for refined management of land and resources, changes in land use types are becoming increasingly frequent. Land change monitoring has become an important issue that needs to be urgently addressed in the fields of agriculture, ecology, planning, and environment. In particular, in typical application scenarios such as farmland protection, illegal land occupation supervision, identification of urban-rural boundary expansion, and farmland degradation monitoring, monitoring methods are required to have not only high timeliness and high resolution, but also key technical advantages such as high spatial positioning accuracy, strong change recognition capabilities, and accurate type discrimination. Traditional manual field survey methods are time-consuming and labor-intensive, with long cycles and delayed data updates, and can no longer meet the field needs of rapid response and dynamic updates.

[0003] In recent years, remote sensing technology has become increasingly widely used in land change monitoring. Change detection techniques, particularly those based on multispectral imagery, have played a significant role in large-scale regional land use renewal. However, most existing technologies still rely on satellite remote sensing data. These technologies are limited by factors such as long imaging cycles, limited resolution, and a fixed viewing angle. Consequently, they suffer from inaccuracy and delayed response in areas with complex terrain and drastic or frequent local changes, such as mountainous farmland and urban-rural fringe areas.

[0004] Unmanned aerial vehicle (UAV) remote sensing platforms are ideal for ground-based remote sensing due to their flexibility, high resolution, and low-altitude flight without cloud interference. Equipped with multispectral cameras, inertial measurement units (IMUs), and GPS positioning systems, UAVs can collect high-precision, multi-temporal remote sensing data of ground features, enabling refined measurements in complex areas. UAV-based data offers high spatial accuracy, strong timeliness, and low reliance on ground control points, making it particularly suitable for rapid identification of sensitive local changes and boundary change monitoring, which requires high structural accuracy.

[0005] In existing research, methods for identifying farmland changes based on satellite remote sensing have been widely used. For example, patent publication CN118608939B proposes an agricultural drought monitoring method based on multi-source remote sensing data. By analyzing the land assessment coefficients of target farmland subregions, this method accurately identifies cracked areas of farmland and assists in guiding crop planting choices, improving the spatiotemporal efficiency and data interpretation capabilities of drought monitoring. This technology shows promising application prospects in agricultural drought scenarios. However, this method inherently relies on fixed-phase satellite imagery and specific index differences, requiring low geometric accuracy for images. It also fails to consider the interference of remote sensing imagery with attitude disturbances and terrain undulations on change identification. Furthermore, this method constructs land assessment indicators based solely on the spectral characteristics of land objects and empirical factors, lacking systematic temporal modeling capabilities and spatial change trend analysis mechanisms. This makes it difficult to apply to more complex and fine-grained land change scenarios (such as slow degradation, illegal land occupation, and boundary erosion).

[0006] To address this issue, this application proposes a land change monitoring method and system based on unmanned measurement. By collecting multispectral remote sensing image data, the target land area is divided into multiple grid units, and the land change monitoring index is periodically calculated and processed, and the land change intensity is quantified. Summary of the Invention

[0007] The present invention provides a land change monitoring method and system based on unmanned measurement, which introduces UAV IMU positioning data and GPS positioning data, combines with digital elevation model, and uses ray projection and intersection solution method to perform precise geometric correction on the original pixel coordinates. Based on vegetation index and water body index, a variety of fusion monitoring indices are constructed to monitor the intensity of land change from multiple dimensions. Multi-index collaborative scoring and a multi-parameter trend recognition algorithm combined with index change trends are introduced to generate the land change intensity of grid units, thereby quantifying the change amplitude of surface cover in the target land area and characterizing its long-term evolution behavior. It has stronger physical continuity and temporal interpretability, and realizes land change monitoring based on unmanned measurement.

[0008] To achieve the above objectives, the present invention provides a land change monitoring method based on unmanned measurement, comprising the following steps:

[0009] S1: Use a drone equipped with a multispectral camera to collect multispectral remote sensing image data of the target land area, combine the drone's attitude information to perform geometric correction on the pixel coordinates in the multispectral remote sensing image data, and correct the pixel coordinates to the ground coordinates of the target land area to obtain multispectral remote sensing data of the ground coordinates;

[0010] S2: Divide the target land area into multiple non-overlapping grid cells and generate a land change monitoring index for each grid cell based on multispectral remote sensing data;

[0011] S3: Construct a time series of land change monitoring indices for grid cells and use a multi-parameter trend identification algorithm to generate the land change intensity for grid cells;

[0012] S4: Based on the spatial structure of the grid cells, spatial filtering is used to smooth the land change intensity of the grid cells to form a visual change intensity monitoring map of the target land area. The changes in the dominant change parameters in the land change monitoring index are extracted to classify and identify land change behaviors.

[0013] As a further improvement method of the present invention:

[0014] Optionally, a drone equipped with a multispectral camera is used to collect multispectral remote sensing image data of the target land area, and pixel coordinates in the multispectral remote sensing image data are geometrically corrected in combination with the drone attitude information, including:

[0015] The drone is equipped with a high-precision GPS and IMU positioning system, as well as a multispectral camera;

[0016] Based on the scope 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 drone 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 drone, and the IMU positioning data is the attitude information of the drone.

[0017] The process of geometric correction of pixel coordinates is as follows:

[0018] Get the pixel coordinates (u, v) and the intrinsic parameter matrix K of the multispectral camera, and convert the pixel coordinates (u, v) into the normalized coordinates (u in the camera coordinate system) 1 ,v 1 ,1):

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

[0020] Where T represents transpose;

[0021] Combining the position coordinates and attitude information of the drone, we can get the ray equation of the pixel coordinates (u, v) in the world coordinate system:

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

[0023] Where P(u,v;λ) represents the ray equation of the pixel coordinate (u,v) in the world coordinate system, λ represents the ray iteration parameter, λ>0, and R(u,v) represents the posture information of the UAV when shooting the multispectral remote sensing image associated with the pixel coordinate (u,v).

[0024] ((u * ,v * ,h * ) represents the position coordinates of the UAV when capturing the multispectral remote sensing image associated with the pixel coordinates (u, v);

[0025] Obtaining a digital elevation model of the target land area, the digital elevation model consisting of ground coordinates of the target land area 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 coordinates ((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 a starting point of the ground coordinates ((x, y)) and a length of z(x, y);

[0026] Gradually increase the ray iteration parameter λ until there is an intersection between the ray equation P(u, v; λ) and an arbitrary elevation function, and use the ground coordinates associated with the elevation function at the intersection as the geometric correction result of the pixel coordinates ((u, v)), and obtain the multispectral remote sensing data of the ground coordinates;

[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 for a grid cell is generated based on 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 unit 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 unit is calculated based on the basic monitoring index, including:

[0031] The land disturbance response index of the grid cell is calculated as follows:

[0032]

[0033] Among them, LDRI t (n,m) represents the land disturbance response index of grid cell B(n,m) in the tth period, NDVI t (n,m) represents the vegetation index of grid cell B(n,m) in the tth period, NDVI t-1 (n,m) represents the vegetation index of grid cell B(n,m) in the t-1th period, Red t (n,m),NIR t (n,m),SWIR t (n,m) represent the reflectance averages of the red band, near infrared band, and short infrared band of the grid cell B(n,m) in the tth cycle, respectively. The tth cycle is the tth period of the UAV’s periodic collection of multispectral remote sensing image data.

[0034] σ(NDVI1(n,m):NDVI t (n,m)) represents the sequence (NDVI1(n,m),NDVI2((n,m),...,NDVI t The standard deviation of ((n,m)), NDVI1(n,m), and NDVI2(n,m) represent the vegetation index of grid cell B((n,m)) in the first and second cycles respectively;

[0035] Grid cell B(n,m) represents the grid cell at the nth row and mth 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;

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

[0037]

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

[0039] The calculation method of the multi-band difference ratio vector index of the grid cell is:

[0040]

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

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

[0043] Optionally, constructing a time series of land change monitoring indices for grid cells includes:

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

[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] Among them, Q 1→t (n,m) represents the time series of land change monitoring index of grid unit B(n,m), Q t ((n,m) represents the land change monitoring index of grid unit B(n,m) in the tth cycle, Q1((n,m), Q2((n,m)) represent the land change monitoring index of grid unit B(n,m) in the first and second cycles respectively;

[0048] NDSI t (n,m) represents the soil index of grid cell B(n,m) in the tth period, grid cell B(n,m) represents the grid cell in the nth row and mth 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, vegetation-water interaction index and multi-band difference ratio vector index of grid cell B((n,m)) in the tth 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 unit, including:

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

[0052] Calculate the normalized land change monitoring index time series The mean and standard deviation of the land change monitoring index in the series are calculated, and the tail value of the land change monitoring index in the series is extracted. The tail 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] Calculate the degree to which the end value of the land change monitoring index deviates from the mean, and use it as the land change monitoring index in the normalized land change monitoring index time series. Changes in

[0054] Based on the dynamic weight of the land change monitoring index, the change items are weighted and fused to obtain the land change intensity of the grid unit, where the land change intensity of the grid unit B (n, m) in the tth period is H t (n,m).

[0055] Optionally, the land change intensity of the grid cells is smoothed using a spatial filtering method in combination with the spatial structure of the grid cells to form a visual change intensity monitoring map of the target land area, including:

[0056] Obtaining the neighboring grid cells of the grid cell to be smoothed and the land change intensity of the neighboring grid cells;

[0057] Combined with the land change intensity of the neighboring grid cells, the land change intensity of the smoothed grid cells is smoothed using spatial filtering, where the land change intensity H t The smoothing result of (n,m) is

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

[0059]

[0060] Among them, the visual change intensity monitoring map t It is a matrix with N rows and M columns. Represents a visual change intensity monitoring mapt The matrix element at row n and column m in ;

[0061] If the smoothed result of the land change intensity of a grid unit is higher than the preset change intensity threshold, the changes in the dominant change parameters in the land change monitoring index of the grid unit are extracted, and the land change behavior of the grid unit is classified and identified. The dominant change parameters in the land change monitoring index include the vegetation index, water body index, soil index, and land disturbance response index.

[0062] In order to solve the above problems, the present invention further provides a land change monitoring system based on unmanned measurement, which includes a data acquisition device and a change intensity monitoring module:

[0063] The data acquisition device is used to collect multispectral remote sensing image data of a target land area using an unmanned aerial vehicle equipped with a multispectral camera, and geometrically correct pixel coordinates in the multispectral remote sensing image data in combination with the attitude information of the unmanned aerial vehicle, and correct the pixel coordinates to the ground coordinates of the target land area to obtain multispectral remote sensing data of the ground coordinates;

[0064] The index monitoring module is used to divide the target land area into multiple non-overlapping grid cells and generate a land change monitoring index for the grid cells based on multispectral remote sensing data;

[0065] The change intensity monitoring module is used to construct a time series of land change monitoring indices for grid cells, use a multi-parameter trend recognition algorithm to generate the land change intensity of grid cells, and use a spatial filtering method to smooth the land change intensity of grid cells in combination with the spatial structure of the grid cells to form a visual change intensity monitoring map of the target land area, and extract the changes in the leading change parameters in the land change monitoring index to classify and identify land change behaviors;

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

[0067] In order to solve the above problem, the present invention further provides an electronic device, comprising:

[0068] a memory storing at least one instruction;

[0069] Communication interfaces to enable electronic equipment to communicate; and

[0070] The processor executes the instructions stored in the memory to implement the above-mentioned land change monitoring method based on unmanned measurement.

[0071] In order to solve the above problems, the present invention also provides a computer-readable storage medium, which stores at least one instruction. 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 existing technology, the present invention proposes a land change monitoring method and system based on unmanned measurement, which has the following beneficial effects:

[0073] First, using high-precision GPS / IMU positioning data recorded by drones, the sensor's external parameters (i.e., position and attitude) are acquired. Combined with a digital elevation model, ray projection of the remote sensing imagery is used to solve the ground intersection problem. This accurately maps pixels from the sensor's perspective to geospatial coordinates, enabling precise correction of image geometric distortion caused by terrain undulations. This process converts tilted and multi-angle remote sensing images into orthophotos with a uniform scale and accurate spatial reference. Compared with traditional approximate geometric correction or plane-hypothetical projection methods, this method significantly improves image positioning accuracy in complex terrain areas such as mountains, hills, and wetlands, avoiding mismatches and false change detections caused by terrain occlusion or perspective differences. Using attitude information to correct ray direction accurately determines pixel line-of-sight direction, eliminating tilt errors. Combined with a digital elevation model, the true intersection of pixels and the ground surface is accurately obtained, correcting for vertical terrain disturbances.

[0074] At the same time, this 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 integrating the temporal change difference of the vegetation index, the historical stability standard deviation and the multispectral reflectance ratio. The first half of the land disturbance response index calculation formula uses the normalized difference to measure the short-term vegetation index change amplitude, and uses the historical vegetation index standard deviation for normalization to suppress false responses under the background of natural fluctuations; the second half combines the near-infrared, red light and short-wave infrared reflectance ratio to enhance the disturbance amplification response to non-vegetation targets such as bare land, excavation areas, and construction areas, thereby enhancing the sensitivity to non-natural disturbances. The land disturbance response index can not only distinguish between drastic changes caused by human activities and natural seasonal changes, but also has good cross-temporal adaptability and regional migration capabilities. Compared with the traditional differential vegetation index or change vector method, the land disturbance response index has stronger stability and abnormal response significance in complex terrain, changeable climate or mixed surface type scenarios. It is particularly suitable for application scenarios such as farmland degradation monitoring, grassland reclamation judgment, and urban-rural expansion identification. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 A schematic flow chart of a land change monitoring method based on unmanned measurement provided by one embodiment of the present invention;

[0076] The purpose, features and advantages of the present invention will be further described with reference to the accompanying drawings and in conjunction with the embodiments. DETAILED DESCRIPTION

[0077] It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0078] The present embodiment provides a method for monitoring land changes based on unmanned measurement. The method can be executed by at least one of a server, a terminal, or other electronic device capable of executing the method provided by the present embodiment. In other words, the method can be executed by software or hardware installed on a terminal or server device, where the software can be a blockchain platform. The server can include, but is not limited to, a single server, a server cluster, a cloud server, or a cloud server cluster.

[0079] Reference Figure 1 , embodiment 1 of the present invention is:

[0080] S1: Use a drone equipped with a multispectral camera to collect multispectral remote sensing image data of the target land area. Combined with the drone's attitude information, the pixel coordinates in the multispectral remote sensing image data are geometrically corrected 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.

[0081] Multispectral remote sensing image data of the target land area is collected using a drone equipped with a multispectral camera. The pixel coordinates in the multispectral remote sensing image data are geometrically corrected based on the drone's attitude information, including:

[0082] The drone is equipped with a high-precision GPS and IMU positioning system, as well as a multispectral camera;

[0083] Based on the scope 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 drone 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 drone, and the IMU positioning data is the attitude information of the drone.

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

[0085] Get the pixel coordinates (u, v) and the intrinsic parameter matrix K of the multispectral camera, and convert the pixel coordinates (u, v) into the normalized coordinates (u in the camera coordinate system) 1 ,v 1 ,1):

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

[0087] Where T represents transpose;

[0088] Combining the position coordinates and attitude information of the drone, we can get the ray equation of the pixel coordinates (u, v) in the world coordinate system:

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

[0090] Where P(u,v;λ) represents the ray equation of the pixel coordinate (u,v) in the world coordinate system, λ represents the ray iteration parameter, λ>0, and R(u,v) represents the posture information of the UAV when shooting the multispectral remote sensing image associated with the pixel coordinate (u,v).

[0091] (u * ,v * ,h * ) represents the position coordinates of the UAV when capturing the multispectral remote sensing image associated with the pixel coordinates (u, v);

[0092] Obtain a digital elevation model of the target land area, the digital elevation model consisting of ground coordinates of the target land area and elevations associated with the ground coordinates, convert the elevations associated with the ground coordinates into an elevation function, wherein the elevation associated with the ground coordinates (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 a starting point at the ground coordinates (x, y) and a length of z(x, y);

[0093] Gradually increase the ray iteration parameter λ until there is an intersection between the ray equation P(u, v; λ) and an arbitrary elevation function, and use the ground coordinates associated with the elevation function at the intersection as the geometric correction result of the pixel coordinates ((u, v)), and obtain the multispectral remote sensing data of the ground coordinates;

[0094] 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.

[0095] S2: Divide the target land area into multiple non-overlapping grid cells and generate a land change monitoring index for each grid cell based on multispectral remote sensing data.

[0096] Generate land change monitoring index for grid cells based on multispectral remote sensing data, including:

[0097] 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 unit 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.

[0098] Specifically, the basic monitoring index of the grid unit is calculated as follows:

[0099]

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

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

[0102] The tth cycle is the tth period during which the UAV periodically collects multispectral remote sensing image data, and t represents the total number of cycles during which the UAV currently collects multispectral remote sensing image data.

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

[0104] The land disturbance response index of the grid cell is calculated as follows:

[0105]

[0106] Among them, LDRI t (n,m) represents the land disturbance response index of grid cell B(n,m) in the tth period, NDVI t (n,m) represents the vegetation index of grid cell B(n,m) in the tth period, NDVI t-1 (n,m) represents the vegetation index of grid cell B(n,m) in the t-1th period, Red t (n,m),NIR t (n,m),SWIR t (n,m) represent the reflectance averages of the red band, near infrared band, and short infrared band of the grid cell B(n,m) in the tth cycle, respectively. The tth cycle is the tth period of the UAV’s periodic collection of multispectral remote sensing image data.

[0107] σ(NDVI1(n,m):NDVI t (n,m)) represents the sequence (NDVI1(n,m),NDVI2((n,m),...,NDVI t The standard deviation of ((n,m)), NDVI1(n,m), and NDVI2(n,m) represent the vegetation index of grid cell B((n,m)) in the first and second cycles respectively;

[0108] Grid cell B(n,m) represents the grid cell at the nth row and mth 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;

[0109] Specifically, although natural variations (such as seasonal greening and post-rainfall greening) can also cause changes in the vegetation index, they are generally cyclical and continuous. The land disturbance response index, by dividing by the historical standard deviation, suppresses dramatic but predictable natural fluctuations. However, sudden human disturbances, such as excavation and reclamation, are often accompanied by a sharp drop in the vegetation index but have no historical basis for fluctuations. Therefore, the land disturbance response index is significantly amplified.

[0110] The enhanced term in the land disturbance response index is For bare land or construction areas, the reflectivity in the near-infrared band is low, while the reflectivity in the red and short-infrared bands is high, which significantly reduces the enhancement factor. This is especially evident when bare soil is exposed after vegetation removal, making the land disturbance response index more sensitive to non-green disturbances.

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

[0112]

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

[0114] As a preferred embodiment of the present invention, the vegetation-water interaction index is obtained by coupling the vegetation index and the water index, which describes the coupling relationship and dynamic changes between vegetation coverage and water distribution in space. It is suitable for remote sensing monitoring of intersecting areas such as wetlands, rice fields, and irrigated farmland. The vegetation-water interaction index is designed to use a product form to enhance the response intensity of areas where vegetation and water coexist, and introduces NDVI in the denominator. t (n,m)+NDWI t ((n,m)+∈ achieves a balance in normalization and enhances numerical stability in low-coverage areas, thereby effectively avoiding false identification caused by a single high index. Compared with traditional independent vegetation indices or water body indices, the vegetation-water interaction index is more sensitive to subtle changes in the coexistence of vegetation and water bodies. It can promptly reflect phenomena such as changes in the ratio of water surface to seedlings during the rice planting period, changes in the boundaries of returning farmland to wetlands projects, and seasonal water level fluctuations.

[0115] The calculation method of the multi-band difference ratio vector index of the grid cell is:

[0116]

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

[0118] They represent the multi-band difference ratio vector index F respectively. t 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);

[0119] In this embodiment, compared to single-band or traditional indices, the multi-band difference ratio vector index emphasizes relative trends between bands, effectively eliminating the effects of light intensity and sensor gain, and improving data consistency across time and devices. For example, the difference ratios between red and green light, near-infrared and red light, and short-wave infrared and near-infrared can reflect subtle changes in vegetation red edges, moisture content, and soil structure, helping to identify fine-grained changes in surface conditions. The dimensions of this multi-band difference ratio vector index jointly represent the gradient trend of the spectral curve, providing a more sensitive response to vegetation growth stages, water drying processes, and land disturbance levels. Compared to traditional vegetation, soil, and water indices, it offers greater dimensional representation and information gain.

[0120] S3: Construct a time series of land change monitoring indices for grid cells and use a multi-parameter trend identification algorithm to generate the land change intensity of grid cells.

[0121] Construct a time series of land change monitoring indices for grid cells, including:

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

[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] Among them, Q 1→t (n,m) represents the time series of land change monitoring index of grid unit B(n,m), Q t ((n,m) represents the land change monitoring index of grid unit B(n,m) in the tth cycle, Q1((n,m), Q2((n,m)) represent the land change monitoring index of grid unit B(n,m) in the first and second cycles respectively;

[0126] NDSI t(n,m) represents the soil index of grid cell B(n,m) in the tth period, grid cell B(n,m) represents the grid cell in the nth row and mth 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;

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

[0128] 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 unit, including:

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

[0130] Specifically,

[0131] in:

[0132] They are Q1(n,m),Q2(n,m),...,Q t Normalization result of (n,m);

[0133] NDSI t (n,m),LDRI t ((n,m),VWI t (n,m),F t Normalization result of (n,m);

[0134] NDSI t (n,m),LDRI t (n,m),VWI t (n,m),F t The normalization method of (n,m) is the minimum-maximum normalization method;

[0135] Calculate the normalized land change monitoring index time series The mean and standard deviation of the land change monitoring index in the series are calculated, and the tail value of the land change monitoring index in the series is extracted. The tail 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.

[0136] In this embodiment, the dynamic weight calculation method of the land disturbance response index in the land change monitoring index is:

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

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

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

[0140] Calculate the degree to which the end value of the land change monitoring index deviates from the mean, and use it as the land change monitoring index in the normalized land change monitoring index time series. Changes in

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

[0142]

[0143] in, Represents the normalized land change monitoring index time series Changes to the Earth Disturbance Response Index, Represents the normalized land change monitoring index time series The mean value of the land disturbance response index, Represents the normalized land change monitoring index time series the standard deviation of the soil disturbance response index;

[0144] Based on the dynamic weight of the land change monitoring index, the change items are weighted and fused to obtain the land change intensity of the grid unit, where the land change intensity of the grid unit B (n, m) in the tth period is H t(n,m).

[0145] S4: Based on the spatial structure of the grid cells, spatial filtering is used to smooth the land change intensity of the grid cells to form a visual change intensity monitoring map of the target land area. The changes in the dominant change parameters in the land change monitoring index are extracted to classify and identify land change behaviors.

[0146] Combined with the spatial structure of the grid cells, spatial filtering is used to smooth the land change intensity of the grid cells to form a visual change intensity monitoring map of the target land area, including:

[0147] Obtaining neighboring grid cells of the grid cell to be smoothed and the land change intensity of the neighboring grid cells; as an embodiment of the present invention, the neighboring grid cells of the grid cell to be smoothed are grid cells within a 3×3 grid cell area centered on the grid cell to be smoothed;

[0148] Combined with the land change intensity of the neighboring grid cells, the land change intensity of the grid cells to be smoothed is smoothed, where the smoothing formula is:

[0149]

[0150] in, Represents the land change intensity H of grid cell B(n,m) in the tth period t ((n,m) smoothing result, β represents the filter coefficient, set β to 0.8, represents the land change intensity of the i-th neighboring grid cell of grid cell B(n,m);

[0151] Grid cell B(n,m) represents the grid cell at the nth row and mth 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;

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

[0153]

[0154] Among them, the visual change intensity monitoring map t It is a matrix with N rows and M columns. Represents a visual change intensity monitoring map t The matrix element at row n and column m in ;

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

[0156] Specifically, if the vegetation index decreases rapidly and the land disturbance response index increases, the land change behavior of the grid unit is construction occupation, that is, the land is converted into construction land. If the vegetation index decreases slowly and the soil index increases, the land change behavior of the grid unit is land degradation. If the water body index continues to decrease, the land change behavior of the grid unit is water body shrinkage or water source loss. If both the water body index and the vegetation index increase, the land change behavior of the grid unit is green space expansion or waterside greening.

[0157] Example 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 used to collect multispectral remote sensing image data of a target land area using an unmanned aerial vehicle equipped with a multispectral camera, and geometrically correct pixel coordinates in the multispectral remote sensing image data in combination with the attitude information of the unmanned aerial vehicle, and correct the pixel coordinates to the ground coordinates of the target land area to obtain multispectral remote sensing data of the ground coordinates;

[0160] The index monitoring module is used to divide the target land area into multiple non-overlapping grid cells and generate a land change monitoring index for the grid cells based on multispectral remote sensing data;

[0161] The change intensity monitoring module is used to construct a time series of land change monitoring indices for grid cells, use a multi-parameter trend recognition algorithm to generate the land change intensity of grid cells, and use spatial filtering to smooth the land change intensity of grid cells in combination with the spatial structure of the grid cells to form a visual change intensity monitoring map of the target land area, extract changes in the dominant change parameters in the land change monitoring index, and classify and identify land change behaviors.

[0162] This is to realize the land change monitoring method based on unmanned measurement as described in Example 1.

[0163] It should be understood that the embodiment is for illustration only and the scope of the patent application is not limited to this structure.

[0164] It should be noted that the serial numbers of the above-mentioned embodiments of the present invention are for descriptive purposes only and do not represent the advantages or disadvantages of the embodiments. In addition, the terms "including", "comprising" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, device, article or method comprising a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, device, article or method. In the absence of further restrictions, an element defined by the sentence "including a ..." does not exclude the presence of other identical elements in the process, device, article or method comprising the element.

[0165] Through the description of the above embodiments, those skilled in the art can clearly understand that the above-mentioned embodiment methods can be implemented by means of software plus the necessary general hardware platform, and of course can also be implemented by hardware, but in many cases the former is a better embodiment. Based on this understanding, the technical solution of the present invention is essentially or the part that contributes to the prior art 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 number of instructions for enabling a terminal device (which can be a mobile phone, computer, server, or network device, etc.) to execute the methods described in each embodiment of the present invention.

[0166] The above are only preferred embodiments of the present invention and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made using the contents of the present invention description and drawings, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.

Claims

1. A land change monitoring method based on unmanned measurement, characterized in that: The method comprises: S1: Use a drone equipped with a multispectral camera to collect multispectral remote sensing image data of the target land area, combine the drone's attitude information to perform geometric correction on the pixel coordinates in the multispectral remote sensing image data, and correct the pixel coordinates to the ground coordinates of the target land area to obtain multispectral remote sensing data of the ground coordinates; S2: Divide the target land area into multiple non-overlapping grid cells and generate a land change monitoring index for each grid cell based on multispectral remote sensing data; S3: Construct a time series of land change monitoring indices for grid cells and use a multi-parameter trend identification algorithm to generate the land change intensity for grid cells; S4: Based on the spatial structure of the grid cells, spatial filtering is used to smooth the land change intensity of the grid cells to form a visual change intensity monitoring map of the target land area. The changes in the dominant change parameters in the land change monitoring index are extracted to classify and identify land change behaviors.

2. The land change monitoring method based on unmanned measurement according to claim 1, characterized in that: Multispectral remote sensing image data of the target land area is collected using a drone equipped with a multispectral camera. The pixel coordinates in the multispectral remote sensing image data are geometrically corrected based on 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 scope 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 drone 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 drone, and the IMU positioning data is the attitude information of the drone. The process of geometric correction of pixel coordinates is as follows: Get the pixel coordinates (u, v) and the intrinsic parameter matrix K of the multispectral camera, and convert the pixel coordinates (u, v) into the normalized coordinates (u in the camera coordinate system) 1 ,v 1 ,1): [u 1 ,v 1 ,1] T =K -1 ·[u,v,1]; Where T represents transpose; Combining the position coordinates and attitude information of the drone, we can get the ray equation of the pixel coordinates (u, v) in the world coordinate system: P(u,v;λ)=[u * ,v * ,h * ] T +λ·R(u,v)·[u 1 ,v 1 ,1] T ; Where P(u,v;λ) represents the ray equation of the pixel coordinate (u,v) in the world coordinate system, λ represents the ray iteration parameter, λ>0, and R(u,v) represents the posture information of the UAV when shooting the multispectral remote sensing image associated with the pixel coordinate (u,v). ((u * ,v * ,h * ) represents the position coordinates of the UAV when capturing the multispectral remote sensing image associated with the pixel coordinates (u, v); Obtaining a digital elevation model of the target land area, the digital elevation model consisting of ground coordinates of the target land area 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 coordinates ((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 a starting point of the ground coordinates ((x, y)) and a length of z(x, y); Gradually increase the ray iteration parameter λ until there is an intersection between the ray equation P(u, v; λ) and an arbitrary elevation function, and use the ground coordinates associated with the elevation function at the intersection as the geometric correction result of the pixel coordinates ((u, v)), and obtain the multispectral remote sensing data of the ground coordinates; 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.

3. The land change monitoring method based on unmanned measurement according to claim 2, characterized in that: Generate land change monitoring index for grid cells based on multispectral remote sensing data, 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, and the fusion monitoring index of the grid unit 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.

4. The land change monitoring method based on unmanned measurement according to claim 3, characterized in that: The fusion monitoring index of the grid unit is calculated based on the basic monitoring index, including: The land disturbance response index of the grid cell is calculated as follows: Among them, LDRI t (n,m) represents the land disturbance response index of grid cell B(n,m) in the tth period, NDVI t (n,m) represents the vegetation index of grid cell B(n,m) in the tth period, NDVI t-1 (n,m) represents the vegetation index of grid cell B(n,m) in the t-1th period, Red t (n,m),NIR t (n,m),SWIR t (n,m) represent the reflectance averages of the red band, near infrared band, and short infrared band of the grid cell B(n,m) in the tth cycle, respectively. The tth cycle is the tth period of the UAV’s periodic collection of multispectral remote sensing image data. σ(NDVI1(n,m):NDVI t (n,m)) represents the sequence (NDVI1(n,m),NDVI2((n,m),...,NDVI t The standard deviation of ((n,m)), NDVI1(n,m), and NDVI2(n,m) represent the vegetation index of grid cell B((n,m)) in the first and second cycles respectively; Grid cell B(n,m) represents the grid cell at the nth row and mth 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; The vegetation-water interaction index of the grid cell is calculated as follows: Among them, VWI t (n,m) represents the vegetation-water interaction index of grid cell B(n,m) in the tth period, NDWI t (n,m) represents the water index of grid cell B(n,m) in the tth period, ∈ represents the control parameter, and ∈ is set to 0.1; The calculation method of the multi-band difference ratio vector index of the grid cell is: Among them, F t (n,m) represents the multi-band difference ratio vector index of the grid cell B(n,m) in the tth period, They represent the multi-band difference ratio vector index F respectively. t The red-green band difference ratio, near-infrared-red band difference ratio, and short-infrared-near-infrared band difference ratio in (n,m).

5. The land change monitoring method based on unmanned measurement according to claim 1, characterized in that: Construct a time series of land change monitoring indices for grid cells, including: The vegetation index and water body index in the basic monitoring index are removed to form the land change monitoring index time series of the grid unit: Q 1→t (n,m)=[Q1(n,m),Q2(n,m),...,Q t (n,m)],n∈[1,N],m∈[[1,M]; Q t (n,m)=[NDSI t (n,m),LDRI t (n,m),VWI t (n,m), F t (n,m)]? Among them, Q 1→t (n,m) represents the time series of land change monitoring index of grid unit B(n,m), Q t ((n,m) represents the land change monitoring index of grid unit B(n,m) in the tth cycle, Q1((n,m), Q2((n,m)) represent the land change monitoring index of grid unit B(n,m) in the first and second cycles respectively; NDSI t (n,m) represents the soil index of grid cell B(n,m) in the tth period, grid cell B(n,m) represents the grid cell in the nth row and mth 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; LDRI t (n,m),VWI t (n,m),F t (n,m) represent the land disturbance response index, vegetation-water interaction index and multi-band difference ratio vector index of grid cell B((n,m)) in the tth period respectively.

6. The land change monitoring method based on unmanned measurement according to claim 5, characterized in that: 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 unit, including: The land change monitoring index time series Q 1→t The land change monitoring index in (n,m) is normalized to obtain the normalized land change monitoring index time series Calculate the normalized land change monitoring index time series The mean and standard deviation of the land change monitoring index in the series are calculated, and the tail value of the land change monitoring index in the series is extracted. The tail 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. Calculate the degree to which the end value of the land change monitoring index deviates from the mean, and use it as the land change monitoring index in the normalized land change monitoring index time series. Changes in Based on the dynamic weight of the land change monitoring index, the change items are weighted and fused to obtain the land change intensity of the grid unit, where the land change intensity of the grid unit B (n, m) in the tth period is H t (n,m).

7. The land change monitoring method based on unmanned measurement according to claim 6, characterized in that: Combined with the spatial structure of the grid cells, spatial filtering is used to smooth the land change intensity of the grid cells to form a visual change intensity monitoring map of the target land area, including: Obtaining the neighboring grid cells of the grid cell to be smoothed and the land change intensity of the neighboring grid cells; Combined with the land change intensity of the neighboring grid cells, the land change intensity of the smoothed grid cells is smoothed using spatial filtering, where the land change intensity H t The smoothing result of (n,m) is The visual change intensity monitoring map of the target land area in the tth period is Map t : Among them, the visual change intensity monitoring map t It is a matrix with N rows and M columns. Represents a visual change intensity monitoring map t The matrix element at row n and column m in ; If the smoothed result of the land change intensity of a grid unit is higher than the preset change intensity threshold, the changes in the dominant change parameters in the land change monitoring index of the grid unit are extracted, and the land change behavior of the grid unit is classified and identified. The dominant change parameters in the land change monitoring index include the vegetation index, water body index, soil index, and land disturbance response index.

8. A land change monitoring system based on unmanned measurement, 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 a target land area using an unmanned aerial vehicle equipped with a multispectral camera, and geometrically correct pixel coordinates in the multispectral remote sensing image data in combination with the attitude information of the unmanned aerial vehicle, and correct the pixel coordinates 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 cells and generate a land change monitoring index for the grid cells based on multispectral remote sensing data; The change intensity monitoring module is used to construct a time series of land change monitoring indices for grid cells, use a multi-parameter trend recognition algorithm to generate the land change intensity of grid cells, and use a spatial filtering method to smooth the land change intensity of grid cells in combination with the spatial structure of the grid cells to form a visual change intensity monitoring map of the target land area, and extract the changes in the leading change parameters in the land change monitoring index to classify and identify land change behaviors; To realize the land change monitoring method based on unmanned measurement as described in any one of claims 1-7.

Images