Method, equipment and medium for monitoring cultivated land in hilly and mountainous areas by integrating sky, ground and air

Abstract

The present invention discloses a method, equipment and medium for monitoring cultivated land in hilly mountainous areas with sky-ground integration. The method comprises the following steps: S1: obtaining satellite orbit data and laser point cloud data of hilly mountainous areas; S2: registering laser point cloud data using a registration method with multi-feature constraints; S3: dividing a monitoring study area on a map, and using SBAS-InSAR technology to perform time-series InSAR analysis on satellite orbit data within the monitoring study area; S4: fusing the analysis results of SBAS-InSAR technology with the laser point cloud data of the nearest point cloud point to obtain the deformation area and settlement degree of cultivated land in hilly mountainous areas within the monitoring study area. The present invention overcomes the problem of insufficient recognition accuracy in the prior art when performing deformation monitoring on cultivated land in hilly mountainous areas. The method fuses laser point cloud data with the processing results of SBAS-InSAR technology, fully applies laser point cloud data to complete the terrain information of the study area, and realizes high-precision monitoring and effective prediction of settlement of cultivated land areas in hilly mountainous areas.

Description

Technical Field

[0001] The present invention relates to the field of terrain monitoring technology, and in particular to a sky-ground integrated farmland monitoring method, equipment and medium in hilly and mountainous areas. Background Art

[0002] Arable land is an important resource and the most basic condition for agricultural production. Over time, the quality of arable land and its suitable distribution area will change due to terrain subsidence, especially in areas with more hilly and mountainous areas. The terrain changes are more complex than in plain areas. Therefore, it is necessary to continuously monitor the arable land in hilly and mountainous areas. Due to the complex terrain of hilly and mountainous areas, it poses a considerable challenge to the exploration and detection technology.

[0003] At present, the monitoring of terrain subsidence mostly uses satellite orbit data images for data analysis and comparison to estimate the approximate deformation rate and deformation trend in the selected area. This technology is mostly used to monitor changes in mining terrain and key area terrain. Due to the difference in the environment, the effect of directly applying this technology to hilly and mountainous areas is often not ideal, and there are fewer monitoring methods for hilly and mountainous terrain. Chinese patent publication number CN117310703A, publication date December 29, 2023, named "InSAR mountainous and hilly area deformation monitoring suitability evaluation method and system" discloses a method that can be used to monitor hilly and mountainous area deformation, grid the target area, and use SAR data, DEM data, and MDDG mathematical function model to obtain the InSAR suitability grade distribution map of the target area, forming a set of evaluation methods for hilly and mountainous area deformation. However, this method is important to solve the problem of terrain deformation evaluation based on existing InSAR analysis results, and the data used is limited to SAR data, DEM data and optical image information. Since the accuracy of the available satellite orbit data is generally low, the data distribution density is low in the target area, and the actual satellite orbit data cannot be fully and accurately covered in a large range. As a result, many areas need to use data points farther away from the surrounding area for deformation analysis because there is no available accurate data, which leads to the fuzziness of the judgment area of ​​terrain change, and the obtained terrain subsidence degree and subsidence rate are not accurate enough. Summary of the invention

[0004] The present invention overcomes the problem of insufficient recognition accuracy in the prior art when monitoring deformation of cultivated land in hilly and mountainous areas, and provides a sky-ground integrated method, equipment and medium for monitoring cultivated land in hilly and mountainous areas. The method integrates laser point cloud data with the processing results of SBAS-InSAR technology, fully utilizes laser point cloud data to complete the terrain information of the study area, and realizes high-precision monitoring and effective prediction of regional subsidence of cultivated land in hilly and mountainous areas.

[0005] In order to achieve the above object, the present invention adopts the following scheme:

[0006] The method for monitoring cultivated land in hilly and mountainous areas by integrating air, ground and hills includes the following steps:

[0007] S1: Acquire satellite orbit data and laser point cloud data of hilly and mountainous areas, wherein the laser point cloud data includes station-borne point cloud data collected by ground base stations and airborne point cloud data collected by aircraft equipment;

[0008] S2: Using a multi-feature constraint registration method to register the laser point cloud data and construct a common error model, all the laser point cloud data are integrated into a common geographic reference system, wherein the multi-feature constraint includes surface constraint, line constraint and point constraint;

[0009] S3: Divide the monitoring study area on the map and use SBAS-InSAR technology to perform time series InSAR analysis on satellite orbit data within the monitoring study area;

[0010] S4: Distribute the point cloud points of the laser point cloud data to the geocentric coordinate system on the map, fuse the analysis results of the SBAS-InSAR technology with the laser point cloud data of the nearest point cloud points, and obtain the deformation area and subsidence degree of cultivated land in the hilly and mountainous areas within the monitoring study area.

[0011] Preferably, the face constraint in step 2 is constructed as follows:

[0012] Based on the normal parallelism and center point coplanarity criteria, the surface constraint expression is obtained: Among them, R is the rotation matrix, t is the translation parameter, and F 0 and F 1 is the normal vector, M 0 and M 1 is the center point coordinate; the surface constraint in the point cloud to be registered is represented as S 0 (M 0 ,F 0 ) and S 1 (M 1 ,F 1 ), and the error equation is obtained by separating the related variables:

[0013]

[0014] Among them, V 21 is the normal error, V 22 is the center error, A 21 is the coefficient matrix composed of normal vectors, A 22 is the coefficient matrix composed of the center coordinates, r is the rotation correction parameter, B 22 is the distance constraint coefficient matrix, L 21 is the direction residual, L 22is the distance residual; the error equation is combined into V 2 =A 2 r+B 2 tL 2 , where A 2 =[A 21 T A 22 T ] T , B 2 =[0B 22 T ] T , L 2 =[L 21 T L 22 T ] T .

[0015] Preferably, the line constraint in step 2 is constructed as follows:

[0016] Take a point P(p x ,p y ,p z ) and direction quantity F(f x ,f y ,f z ), with the spatial line feature L 0 (P 0 ,F 0 ) and L 1 (P 1 ,F 1 ) is taken as the criterion, and the line constraint expression is obtained:

[0017]

[0018] Where Δx p , Δy p , Δz p Point P 1 To point P 0 The translation amount on the x, y, and z axes, where point P 1 Located on line L 0 superior.

[0019] Preferably, the unified rotation and translation error equation expression of the error model in step 2 is:

[0020]

[0021] Among them, A is the coefficient matrix composed of observations, L is the corresponding residual matrix, R is the rotation matrix parameter, B is the coefficient matrix composed of multi-constraint distance constraints, t is the translation parameter, W is the corresponding residual matrix, m, n, l are the number of point constraints, line constraints and surface constraints involved in the registration respectively.

[0022] Preferably, the laser point cloud data in step 1 also includes image feature point cloud data, and the image feature point cloud data is obtained by analyzing and extracting the hilly and mountainous area images taken by the UAV using the SFM algorithm.

[0023] Preferably, in step 4, the point cloud points of the laser point cloud data are distributed into the geocentric coordinate system using the following method:

[0024] The laser point cloud data is projected as the LOS vector in the UTM coordinate system, and the position of the PS point on the laser point cloud data is located in combination with the PS point. The UTM coordinates of the laser point cloud data are converted into the longitude and latitude of the geocentric coordinate system WGS84 using the Earth-centered Earth-fixed coordinate system (ECEF) conversion formula.

[0025] Preferably, the step S3 of performing time series InSAR analysis on satellite orbit data using SBAS-InSAR technology comprises the following steps:

[0026] N+1 SAR images of different periods in the monitoring study area are extracted from the satellite orbit data, and k differential interferograms are obtained by processing them using the set time and space thresholds. The differential interferogram phase is calculated using the following formula:

[0027]

[0028] in, t b Time, t a The phase at the moment, d(t b ,x,r)、d(t a ,x,r) is t b Time, t a Keep your sight line pointed to the accumulated settlement at all times;

[0029] Calculate the average rate of ground subsidence v t :

[0030]

[0031] Calculate the phase:

[0032]

[0033] Among them, E J is the main image acquisition time, S J The acquisition time of the secondary image, vk is the pixel deformation rate at time k, and A(j,k)=t k -t k-1 , where A is an m×n rank-deficient matrix, we get the matrix equation A v =δφ, φ is the phase of the permanent scatterer point PS in the image in the interference pattern. The average settlement rate phase value and ground settlement amount are calculated using this matrix equation.

[0034] Preferably, step S1 further comprises: storing the acquired satellite orbit data and laser point cloud data in a spatial database, wherein the spatial database operates as follows: abstracting the acquired multiple data and storing them respectively in corresponding local ER models, integrating the local ER models into a global ER model by an iterative refinement method, converting the global ER model into a UML structure, and constructing the logical structure of the data model in the database accordingly.

[0035] The present application also provides an electronic device, comprising: at least one processor, and a memory communicatively connected to the at least one processor, wherein the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to perform a method as described in any one of the above schemes.

[0036] The present application also provides a storage medium on which a computer program is stored. When the program is executed by a processor, the method described in any one of the above solutions is implemented.

[0037] The present invention includes at least the following beneficial effects: (1) combining satellite orbit data and laser point cloud data, and using the three-dimensional characteristics of laser point cloud data to mark complex hilly and mountainous terrain, thereby increasing the effectiveness of the overall data. At the same time, the high density of laser point cloud data is used to increase the coverage of sampling data in the study area, thereby effectively improving the position accuracy of settlement identification; (2) fusing the analysis results of SBAS-InSAR technology with laser point cloud data, and analyzing the terrain characteristics of the study area from multiple data from three angles of sky and ground, so that the monitoring results are closer to the changes in terrain; (3) a dedicated database structure is designed for the data characteristics of cultivated land in hilly and mountainous areas, which can more quickly organize the acquired data and improve the efficiency of identifying terrain changes. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 A method principle flow chart of the present invention;

[0039] Figure 2 It is the structure diagram of the monitoring object ER model of the spatial database architecture of the present invention;

[0040] Figure 3It is a sensor ER model structure diagram of the spatial database architecture of the present invention;

[0041] Figure 4 This is a structure diagram of the ER model of the laser radar point cloud data of the spatial database architecture of the present invention;

[0042] Figure 5 It is a typical defect ER model structure diagram of the spatial database architecture of the present invention. DETAILED DESCRIPTION

[0043] The present invention is further described in detail below in conjunction with the accompanying drawings so that those skilled in the art can implement the invention with reference to the description.

[0044] like Figure 1 As shown, the sky-ground integrated hilly and mountainous area farmland monitoring method provided by the present invention comprises the following steps:

[0045] S1: Obtain satellite orbit data and laser point cloud data in hilly and mountainous areas. The laser point cloud data includes station-borne point cloud data collected by ground base stations and airborne point cloud data collected by aircraft equipment. Satellite orbit data is usually obtained free of charge through the corresponding official data release website. Generally, data from multiple satellite orbit systems are obtained at the same time for combined use. The advantages of data from each system can be combined and the data update cycle of each system can be adapted to ensure that the latest satellite orbit data is obtained. As an example, the satellite orbit data obtained in real time include: Sentinel-1 satellite orbit data provided by ASF (Alaska Satellite Facility), DEM data obtained from multiple channels around the world, and ORB precise orbit data. The acquired laser point cloud data includes station-mounted point cloud data and airborne point cloud data. Station-mounted point cloud data is ground monitoring point cloud data, and airborne point cloud data is aerial monitoring point cloud data. TLS and ALS, two Lidar technologies, are usually used to collect three-dimensional spatial coordinates in hilly and mountainous areas. TLS is scanned at a low angle, while ALS has a high-angle overlooking scene. The low-angle and high-angle measurement methods have good complementary characteristics and produce all-round point cloud data. The above measurement methods have good time consistency and can collect data separately and complete them all within a few hours. Therefore, it can be regarded as no data time difference. The measured point cloud data will be uniformly registered as an overall point cloud in the future, which is timely, spatially accurate, and complete. The overall point cloud can support a variety of analysis methods, including ellipticity, central axis inclination, convexity and concavity, and surface unfolding flatness, forming a comprehensive and accurate analysis technology for cultivated land deformation in hilly and mountainous areas.

[0046] S2: Use a multi-feature constraint registration method to register laser point cloud data and build a common error model, integrate all laser point cloud data into a common geographic reference system, and the multi-feature constraints include surface constraints, line constraints, and point constraints; register laser point cloud data from different sources, analyze the point cloud data errors to build an error model, and record the feature information of all laser point cloud data in the form of a multi-level fused point cloud data map in the geographic reference system that integrates all point cloud data. The accuracy and robustness of laser point cloud data registration can be improved by combining different feature constraints. Feature registration can be divided into data matching at different times, different perspectives, and different sensors, and high-dimensional structural features can be restored from low-dimensional data content to achieve three-dimensional reconstruction of the target scene and stitching of remote sensing panoramic images.

[0047] S3: Divide the monitoring study area on the map, and use SBAS-InSAR technology to perform time-series InSAR analysis on satellite orbit data within the monitoring study area; when monitoring the deformation of cultivated land in hilly and mountainous areas, select the area to be analyzed on the map as required, and perform data analysis on this area. SBAS-InSAR technology is a time-series InSAR method based on multiple main images. It measures the slow deformation of the surface by analyzing the time series of multi-view SAR images. It can monitor deformation over a large range with high accuracy. Using the short baseline principle, a large amount of SAR data is combined into an interference subset with multiple main images. The interference baseline length in each subset is lower than the critical baseline value, and the time baseline is very short, which solves the problem of de-correlation in time and space, so that reliable monitoring results can be obtained with less data. SBAS-InSAR technology overcomes the shortcomings of the original InSAR method that the single main image leads to poor interference results, improves the coherence and computational efficiency of the interference pattern, reduces the required SAR image data, and is better suitable for data processing in hilly and mountainous areas.

[0048] S4: Distribute the point cloud points of the laser point cloud data to the geocentric coordinate system on the map, fuse the analysis results of the SBAS-InSAR technology with the laser point cloud data of the nearest point cloud point, and obtain the deformation area and settlement degree of the cultivated land in the hilly and mountainous areas in the monitoring study area. The object directly analyzed by the SBAS-InSAR technology is the satellite orbit data of a specified point in the geocentric coordinate system. By distributing the point cloud points of the laser point cloud data to the geocentric coordinate system on the map, when using the SBAS-InSAR technology to perform time series analysis on the specified point, the spatial data structure KD tree is constructed to quickly find the point cloud point closest to the point in the geocentric coordinate system. The KD tree is widely used in high-dimensional space search, such as nearest neighbor search and range search. It recursively divides the space into multiple hyper-rectangular areas, each node represents a hyper-rectangular area, and has corresponding segmentation dimensions and segmentation points. First, select a dimension and median point of the point cloud point as the root node, and then divide the point set into two parts according to the value of the point in the dimension. Recursively construct subtrees for the two parts, and construct a spatial data structure KD tree with point cloud point data. When searching for the closest point cloud point, based on the spatial partitioning characteristics of the KD tree search algorithm, we start from the root node and search downwards according to the comparison results between the specified point analyzed by SBAS-InSAR technology and the current node until we reach the leaf node. Then, by backtracking, we determine whether there are closer points in the unvisited branches, and by calculating the distance between the query point and the subtree boundary, we decide whether further search is needed, and finally obtain the point cloud point closest to the specified point. The SBAS-InSAR analysis results are fused with the data of the point cloud points. The fused objects may include the basic attributes of the point cloud data of the point cloud points and the prediction results obtained by using the feature information of the point cloud points for deformation analysis preprocessing. The SBAS-InSAR analysis results preliminarily obtain the time series change amount and change rate of each point, the overall terrain change range within the segmented area, the transfer direction of the settlement area and other information. The basic attributes of the point cloud points include the current height of each point, the terrain continuity and the three-dimensional contour of the terrain. The preprocessing results of the point cloud points include the three-dimensional change information such as the regional change range and deformation volume within the monitoring time period. The basic attributes and preprocessing results of the point cloud points are combined with the SBAS-InSAR analysis results to form a comprehensive fusion settlement prediction model. The fusion results are converted and calculated. During the calculation process, the SBAS-InSAR analysis results and similar information in the point cloud data are merged and compared. According to the importance and confidence of each data, weighted averaging operations are performed, and the edge information is fuzzy transition processed. Finally, the calculation results are mapped on the map of the monitoring study area as the terrain settlement degree and settlement rate of the point in the monitoring study area.

[0049] This method combines satellite orbit data and laser point cloud data, and uses the three-dimensional characteristics of laser point cloud data to annotate and model complex hilly and mountainous terrain, increasing the authenticity and effectiveness of the overall data. At the same time, the high density of laser point cloud data is used to increase the coverage of sampling data in the study area, enrich the data set used for calculation and analysis, and effectively improve the location accuracy of subsidence identification. By fusing the analysis results of SBAS-InSAR technology with laser point cloud data, the terrain characteristics of the study area are analyzed from multiple data from three angles of sky and ground, making the monitoring results closer to the changes in the terrain.

[0050] In another technical solution, the face constraint in step 2 is constructed as follows:

[0051] Based on the normal parallelism and center point coplanarity criteria, the surface constraint expression is obtained: Among them, R is the rotation matrix, t is the translation parameter, and F 0 and F 1 is the normal vector, M 0 and M 1 is the center point coordinate; the surface constraint in the point cloud to be registered is represented as S 0 (M 0 ,F 0 ) and S 1 (M 1 ,F 1 ), and the error equation is obtained by separating the related variables:

[0052]

[0053] Among them, V 21 is the normal error, V 22 is the center error, A 21 is the coefficient matrix composed of normal vectors, A 22 is the coefficient matrix composed of the center coordinates, r is the rotation correction parameter, B 22 is the distance constraint coefficient matrix, L 21 is the direction residual, L 22 is the distance residual; the error equation is combined into V 2 =A 2 r+B 2 tL 2 , where A 2 =[A 21 T A 22 T ] T , B 2 =[0B 22 T ] T , L 2 =[L21 T L 22 T ] T .

[0054] In another technical solution, the line constraint in step 2 is constructed as follows:

[0055] Take a point P(p x ,p y ,p z ) and direction quantity F(f x ,f y ,f z ), with the spatial line feature L 0 (P 0 ,F 0 ) and L 1 (P 1 ,F 1 ) is taken as the criterion, and the line constraint expression is obtained:

[0056]

[0057] Where Δx p , Δy p , Δz p Point P 1 To point P 0 The translation amount on the x, y, and z axes, where point P 1 Located on line L 0 The normal parallel formula used in the line constraint is the same as that in the above surface constraint formula. In the calculation, let P = RP, and the translation Δx can be obtained. p , Δy p , Δz p .

[0058] In another technical solution, the unified rotation and translation error equation of the error model in step 2 is expressed as:

[0059]

[0060] Among them, A is the coefficient matrix composed of observations, L is the corresponding residual matrix, R is the rotation matrix parameter, B is the coefficient matrix composed of multi-constraint distance constraints, t is the translation parameter, W is the corresponding residual matrix, m, n, l are the number of point constraints, line constraints and surface constraints involved in the registration respectively.

[0061] In another technical solution, the laser point cloud data in step 1 also includes image feature point cloud data, which is obtained by analyzing and extracting the hilly and mountainous area images taken by the drone using the SFM algorithm. Corresponding to the way of collecting laser point cloud data by radar scanning, image acquisition uses two photogrammetry methods: close-range photography and oblique photography, which can collect aerial images of the monitoring study area and generate dense point clouds. Aerial images can provide color information, where close-range photography is taken at a low angle, while oblique photography has a high-angle bird's-eye view, forming a perspective coordination to obtain comprehensive terrain information. By processing the aerial images and integrating them into the laser point cloud data, the point cloud data not only has three-dimensional information such as shape contours, but also has color information that can assist in judgment.

[0062] Step S1 also includes: storing the acquired satellite orbit data and laser point cloud data in a spatial database, wherein the spatial database operates as follows: abstracting the acquired multiple data and storing them in corresponding local ER models respectively, integrating the local ER models into a global ER model by an iterative refinement method, converting the global ER model into a UML structure, and constructing the logical structure of the data model in the database accordingly.

[0063] The ER model (entity-relationship model) is a common method for data modeling. It represents the relationship between entities in a graphical way. Entities are used to represent real objects. These entities have their own attributes, and relationships are used to represent the connection and interaction between entities. Figure 2-5 As shown in Figure 1, in the spatial database of cultivated land monitoring in hilly and mountainous areas, entities include monitoring objects, sensors, Sentinel-1 data, lidar point cloud data, typical defects, etc. Figure 2 As shown in , the attributes of the monitoring object entity include monitoring object ID, name, location, and remarks; Figure 3 As shown in , the attributes of the sensor entity include sensor ID, monitoring object ID, installation location, installation time, warning threshold, danger threshold, sampling period, use environment, and remarks; Figure 4 As shown in , the attributes of the laser radar point cloud data entity include time, three-dimensional point coordinates, color information, and reflection intensity; Figure 5 As shown in the figure, the attributes of a typical defect entity include defect ID, defect type, location, discovery time, severity, and repair status; the local ER model composed of other entities and their attributes in the spatial database architecture is constructed in the above manner. A separate local ER model is created based on each entity and its attributes. The attributes of each local ER model are terminal data items, which cannot be further differentiated into sub-attributes and cannot be related to other entities. Relationships can only occur between entities.

[0064] When the iterative refinement method is used to integrate local ER models into a global ER model, the local ER models of cultivated land monitoring objects and sensors in hilly and mountainous areas are integrated first, and then entities such as lidar monitoring data and typical defects are gradually integrated. Every time an entity or relationship is added, the model is re-evaluated to identify and resolve possible structural conflicts, naming conflicts, and attribute conflicts to optimize the redundancy of data and connections and ensure the gradual improvement and refinement of the model.

[0065] According to the relationship between the entities of the global ER model, the database tables are designed, including the lidar acquisition data table, accelerometer data table, static data table, typical defect table, monitoring object data table, and sensor data table. Considering the one-to-many relationship between the monitoring object and the sensor, one monitoring object can be associated with multiple sensors to monitor it. When designing the database, a primary key is defined for each entity to set the primary attribute to a non-null value, and a certain attribute of "one" is set as a foreign key at the "many" end. Therefore, the ID of the monitoring object is used as a unique identifier and also as a foreign key in the sensor entity table, point cloud data entity table, and typical defect entity table. At the same time, considering that different types of sensors have not only different name attributes, but also their own working parameters and other attributes, an independent entity is created for the sensor type. There is also a one-to-many relationship between the sensor type and the sensor itself, which is associated through the name of the sensor type. Based on the spatial database management system, the data of the hilly and mountainous area spatial database is visualized and integrated with data operation functions.

[0066] In another technical solution, the following method is used to distribute the point cloud points of the laser point cloud data to the geocentric coordinate system in step 4:

[0067] The laser point cloud data is projected as the LOS vector in the UTM coordinate system, and the position of the PS point on the laser point cloud data is located in combination with the PS point. The UTM coordinates of the laser point cloud data are converted into the longitude and latitude of the geocentric coordinate system WGS84 using the Earth-centered Earth-fixed coordinate system (ECEF) conversion formula.

[0068] The original point cloud data is projected in the UTM Zone 48N coordinate system. After point cloud preprocessing, the denoised point cloud is obtained. The point cloud data is converted into the WGS84 coordinate system using the Earth-Centered Earth-Fixed Coordinate System (ECEF) conversion formula. For the design of this application, the azimuth of the satellite is from north to east, that is, the east direction in the geographic coordinate system is 0 degrees. The line of sight direction is calculated based on the incident angle and the azimuth. The component of the LOS vector in the Northeast Sky (ENU, East-North-Up) coordinate system is calculated, and the dot product is used to project each point in the point cloud to the line of sight (LOS) direction. The projection distance or projection point is calculated by the point cloud product. For each point P = (X, Y, Z) in the point cloud, its projection in the LOS direction is calculated. The spatial data structure KD tree is used to quickly find the point cloud point closest to each SBAS-InSAR point and fuse their attributes.

[0069] In another technical solution, performing time series InSAR analysis on satellite orbit data using SBAS-InSAR technology in step S3 includes the following steps:

[0070] N+1 SAR images of different periods in the monitoring study area are extracted from the satellite orbit data, and k differential interferograms are obtained by processing them using the set time and space thresholds. The differential interferogram phase is calculated using the following formula:

[0071]

[0072] in, t b Time, t a The phase at the moment, d(t b ,x,r)、d(t a ,x,r) is t b Time, t a Keep your sight line pointed to the accumulated settlement at all times;

[0073] Calculate the average rate of ground subsidence v t :

[0074]

[0075] Calculate the phase:

[0076]

[0077] Among them, E J is the main image acquisition time, S J The acquisition time of the secondary image, v k is the pixel deformation rate at time k, and A(j,k)=t k -t k-1, where A is an m×n rank-deficient matrix, we get the matrix equation A v =δφ, φ is the phase of the permanent scatterer point PS in the image in the interference pattern. The average settlement rate phase value and ground settlement amount are obtained by using this matrix equation through singular value decomposition and least squares method.

[0078] The present invention also provides an electronic device, comprising: at least one processor, and a memory connected to the at least one processor in communication, wherein the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor, so that the at least one processor executes the sky-ground integrated hilly and mountainous area farmland monitoring method provided by any of the above schemes. The electronic device can be any terminal device including a mobile phone, a laptop computer, a desktop computer, a tablet computer, a PDA (Personal Digital Assistant), a POS (Point of Sales), a vehicle-mounted computer, etc.

[0079] When this electronic device is used for deformation detection of cultivated land in hilly and mountainous areas, it has all the beneficial effects of the method provided in the present application, and can achieve high-precision monitoring and effective prediction of regional settlement of cultivated land in hilly and mountainous areas.

[0080] The present invention also provides a storage medium on which a computer program is stored, characterized in that when the program is executed by a processor, the sky-ground integrated hilly and mountainous area farmland monitoring method provided by any of the above schemes is implemented.

[0081] Through the description of the above implementation mode, the technicians in the relevant field can clearly understand that the present invention can be implemented by means of software plus necessary general hardware, and of course, it can also be implemented by special hardware including special integrated circuits, special CPUs, special memories, special components, etc. In general, all functions completed by computer programs can be easily implemented by corresponding hardware, and the specific hardware structure used to implement the same function can also be various, such as analog circuits, digital circuits or special circuits. However, for the present invention, software program implementation is a better implementation mode in more cases. Based on such an 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 readable storage medium, such as a computer floppy disk, a U disk, a mobile hard disk, a read-only memory (ROM, Read-Only Memory), a random access memory (RAM, Random Access Memory), a disk or an optical disk, etc., including a number of instructions to enable a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in each embodiment of the present invention.

[0082] It should be noted that, although the above describes the steps in a specific order, it does not mean that the steps must be performed in the above specific order. In fact, some of these steps can be performed concurrently or even in a different order, as long as the required functions can be achieved. The number of devices and processing scales described here are used to simplify the description of the present invention, and the application, modification and variation of the present invention are obvious to those skilled in the art.

[0083] Although the embodiments of the present invention have been disclosed as above, they are not limited to the applications listed in the specification and the implementation modes, and they can be fully applied to various fields suitable for the present invention. For those familiar with the art, additional modifications can be easily implemented. Therefore, without departing from the general concept defined by the claims and the scope of equivalents, the present invention is not limited to the specific details and the illustrations shown and described herein.

Claims

1. A method for monitoring cultivated land in hilly and mountainous areas by integrating air, ground and soil, characterized in that: The following steps are involved: S1: Acquire satellite orbit data and laser point cloud data of hilly and mountainous areas, wherein the laser point cloud data includes station-borne point cloud data collected by ground base stations and airborne point cloud data collected by aircraft equipment; S2: Using a multi-feature constraint registration method to register the laser point cloud data and construct a common error model, all the laser point cloud data are integrated into a common geographic reference system, wherein the multi-feature constraint includes surface constraint, line constraint and point constraint; S3: Divide the monitoring study area on the map and use SBAS-InSAR technology to perform time series InSAR analysis on satellite orbit data within the monitoring study area; S4: Distribute the point cloud points of the laser point cloud data to the geocentric coordinate system on the map, fuse the analysis results of the SBAS-InSAR technology with the laser point cloud data of the nearest point cloud points, and obtain the deformation area and subsidence degree of cultivated land in the hilly and mountainous areas within the monitoring study area.

2. The method for monitoring cultivated land in hilly and mountainous areas by integrating sky, ground and hills according to claim 1, characterized in that: The face constraint described in step 2 is constructed as follows: Based on the normal parallelism and center point coplanarity criteria, the surface constraint expression is obtained: Among them, R is the rotation matrix, t is the translation parameter, F0 and F1 are normal vectors, M0 and M1 are the center point coordinates; the surface constraints in the point cloud to be registered are expressed as S0(M0, F0) and S1(M1, F1), and the error equation is obtained by separating the related variables: Among them, V 21 is the normal error, V 22 is the center error, A 21 is the coefficient matrix composed of normal vectors, A 22 is the coefficient matrix composed of the center coordinates, r is the rotation correction parameter, B 22 is the distance constraint coefficient matrix, L 21 is the direction residual, L 22 is the distance residual; the error equation is combined into V2 = A2r + B2t-L2, where A2 = [A 21 T A 22 T ] T , B2=[0B 22 T ] T , L2=[L 21 T L 22 T ] T .

3. The sky-ground integrated hilly and mountainous area farmland monitoring method according to claim 1 is characterized in that: The line constraint described in step 2 is constructed as follows: Take a point P(p x ,p y ,p z ) and direction quantity F(f x ,f y ,f z ), taking the parallelism of the normal vectors of the space line features L0(P0,F0) and L1(P1,F1) as the criterion, the line constraint expression is obtained: Where Δx p , Δy p , Δz p are the translation amounts from point P1 to point P0 on the x, y, and z axes, respectively, where point P1 is on line L0.

4. The sky-ground integrated farmland monitoring method in hilly and mountainous areas according to claim 1 is characterized in that: The unified rotation and translation error equations of the error model described in step 2 are expressed as: Among them, A is the coefficient matrix composed of observations, L is the corresponding residual matrix, R is the rotation matrix parameter, B is the coefficient matrix composed of multi-constraint distance constraints, t is the translation parameter, W is the corresponding residual matrix, m, n, l are the number of point constraints, line constraints and surface constraints involved in the registration respectively.

5. The sky-ground integrated hilly and mountainous area farmland monitoring method according to claim 1 is characterized in that: The laser point cloud data in step 1 also includes image feature point cloud data, which is obtained by analyzing and extracting the hilly and mountainous area images taken by the UAV using the SFM algorithm.

6. The sky-ground integrated hilly and mountainous area farmland monitoring method according to claim 1 is characterized in that: In step 4, the point cloud points of the laser point cloud data are distributed into the geocentric coordinate system using the following method: The laser point cloud data is projected as the LOS vector in the UTM coordinate system, and the position of the PS point on the laser point cloud data is located in combination with the PS point. The UTM coordinates of the laser point cloud data are converted into the longitude and latitude of the geocentric coordinate system WGS84 using the Earth-centered Earth-fixed coordinate system (ECEF) conversion formula.

7. The method for monitoring cultivated land in hilly and mountainous areas by integrating sky, ground and hills according to claim 1, characterized in that: In step S3, using SBAS-InSAR technology to perform time series InSAR analysis on satellite orbit data includes the following steps: N+1 SAR images of different periods in the monitoring study area are extracted from the satellite orbit data, and k differential interferograms are obtained by processing them using the set time and space thresholds. The differential interferogram phase is calculated using the following formula: in, t b Time, t a The phase at the moment, d(t b ,x,r)、d(t a ,x,r) is t b Time, t a Keep your sight line pointed to the accumulated settlement at all times; Calculate the average rate of ground subsidence v t : Calculate the phase: Among them, E J is the main image acquisition time, S J The acquisition time of the secondary image, v k is the pixel deformation rate at time k, and A(j,k)=t k -t k-1 , where A is an m×n rank-deficient matrix, we get the matrix equation A v =δφ, φ is the phase of the permanent scatterer point PS in the image in the interference pattern. The average settlement rate phase value and ground settlement amount are calculated using this matrix equation.

8. The sky-ground integrated farmland monitoring method in hilly and mountainous areas according to claim 1 is characterized in that: Step S1 also includes: storing the acquired satellite orbit data and laser point cloud data in a spatial database, wherein the spatial database operates as follows: abstracting the acquired multiple data and storing them in corresponding local ER models respectively, integrating the local ER models into a global ER model by an iterative refinement method, converting the global ER model into a UML structure, and constructing the logical structure of the data model in the database accordingly.

9. An electronic device, characterized in that include: At least one processor, and a memory communicatively connected to the at least one processor, wherein the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor so that the at least one processor executes the method according to any one of claims 1 to 8.

10. A storage medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the method described in any one of claims 1 to 8 is implemented.

Images