Landslide deformation continuous monitoring method integrating unmanned aerial vehicle photogrammetry and GNSS

Abstract

The invention belongs to the technical field of satellite navigation and geological disaster monitoring, and particularly relates to a landslide deformation planar continuous monitoring method integrating unmanned aerial vehicle photogrammetry and GNSS, which uses an unmanned aerial vehicle to carry a camera to carry out image acquisition on a landslide area, and uses a base station in a GNSS monitoring system to provide a reference for an unmanned aerial vehicle airborne receiver. Unifying unmanned aerial vehicle image coordinates to a GNSS monitoring system based on a double-difference model, a cubic spline interpolation function and a coordinate compensation value of an unmanned aerial vehicle photographing record file, and performing aerial triangulation; a joint monitoring model is constructed based on the point cloud derived from the multi-stage unmanned aerial vehicle image data and the GNSS monitoring sequence; the combined monitoring model is used for calculating the elevation type variable of the deformation area in the monitoring time period and is combined with the DOM for display, the elevation deformation condition of the landslide surface-shaped monitoring area is visually displayed, continuous and visual monitoring of landslide deformation is achieved, and the monitoring precision is high.

Description

Technical Field

[0001] The present invention belongs to the technical field of satellite navigation and signal recognition, and particularly relates to a method for continuous monitoring of landslide deformation by integrating UAV photogrammetry and GNSS. Background Art

[0002] As a serious geological disaster, landslides pose a huge threat to the natural environment, residents' safety and transportation facilities. Therefore, efficient and accurate monitoring technologies for landslides have attracted much attention. GNSS (Global Navigation Satellite System) has been widely used in landslide monitoring due to its high precision and all-weather characteristics. However, its monitoring range is limited to point targets and it is difficult to reflect the dynamic changes of the overall landslide instability surface. Although UAV photogrammetry technology can obtain the entire planar terrain data, limited by the flight frequency and data processing delay, the monitoring accuracy and timeliness cannot meet the requirements of continuous monitoring.

[0003] In existing research, methods of combining GNSS and UAV technologies have been proposed. For example, a slope intelligent monitoring and early warning cloud platform based on multi-source data fusion (CN119169770A) uses intelligent optimization and fusion of different types of data including GNSS displacement data and UAV aerial photography slope topography data from existing independent monitoring systems in open-pit mines, and comprehensively analyzes and judges to provide technical support for the safety of open-pit mine slopes. However, this early warning platform actually mainly relies on inclinometers to monitor the displacement of each monitoring point in real time. Its monitoring range is limited to point targets and does not involve the data fusion of GNSS and UAVs, making it difficult to conduct planar real-time monitoring of the dynamic changes of the landslide instability surface.

[0004] Currently, the methods involving the combination of GNSS and UAV technologies are mainly used for the mutual verification and analysis of GNSS monitoring data and UAV monitoring results, and a systematic method for integrated monitoring has not been formed. There are the following problems in constructing an "air-ground" integrated monitoring system by integrating the two technical means: First, the GNSS monitoring system usually has an independent coordinate framework, and there is a systematic deviation from the coordinate system commonly used by UAVs, which poses a challenge to the unified analysis of the monitoring processes of the two. Second, it is usually only possible to obtain the terrain change amount within the flight time interval using a limited number of UAV-derived data, and it is difficult to obtain planar monitoring data for any date, thus making it impossible to accurately grasp the evolution process of the entire landslide surface. Summary of the Invention

[0005] Aiming at the problem that the existing technology cannot conduct real-time planar monitoring on the dynamic changes of the landslide instability surface, this paper proposes a continuous planar monitoring method for landslide deformation that integrates UAV photogrammetry and GNSS. The present invention uses a UAV equipped with a camera to collect images of the landslide area, uses the reference station in the GNSS monitoring system to provide a reference benchmark for the UAV airborne receiver, and uses the double-difference model formula in GNSS-RTK technology to make the UAV photogrammetry coordinate system consistent with the GNSS monitoring framework. Based on multi-period UAV image data and GNSS elevation monitoring sequences, a joint monitoring model is constructed. The elevation deformation variables of the deformation area during the monitoring period are calculated using the joint monitoring model and displayed in combination with DOM, intuitively showing the elevation deformation of the planar monitoring area of the landslide, realizing continuous planar monitoring of the elevation changes in the landslide area, and thus accurately grasping the evolution process of the entire landslide surface.

[0006] Based on the above objectives, the technical solution adopted by the present invention is as follows:

[0007] A continuous monitoring method for landslide deformation that integrates UAV photogrammetry and GNSS, comprising the following steps:

[0008] (1) In the GNSS monitoring system, based on the position information of the reference station and GNSS monitoring stations, GNSS monitoring station elevation monitoring sequence data is obtained based on the double-difference model;

[0009] (2) UAV photogrammetry obtains UAV airborne receiver data, exposure record file data, and image data;

[0010] (3) Using the reference station in the GNSS monitoring system to provide a reference benchmark for the UAV airborne receiver, based on the double-difference model, the coordinates of the UAV airborne receiver are obtained, and combined with the exposure record file data, the image center coordinates are obtained to achieve the unification of the coordinate frameworks of GNSS monitoring stations and UAVs;

[0011] (4) Under the unified coordinate framework, the image data and image center coordinates are used to obtain UAV-derived data through aerial triangulation, and then the deformation area is determined by DEM differential of multi-period UAV-derived data, and discrete points of different flight dates in the deformation area are obtained;

[0012] (5) Using the GNSS elevation monitoring sequence data obtained in step (1) and the discrete points of different flight dates in the deformation area obtained in step (4) to construct a joint measurement model;

[0013] (6) Based on the joint measurement model, planar monitoring data for any date is obtained, a planar area deformation data set is constructed, and the landslide monitoring results are visualized to achieve continuous and intuitive monitoring of landslide deformation.

[0014] The present invention uses the reference station in the GNSS monitoring system to provide a reference benchmark for the airborne receiver of the unmanned aerial vehicle (UAV). The double-difference model is used to calculate the position of the airborne receiver of the UAV. Then, cubic spline interpolation and the coordinate compensation value in the UAV photo record file are used to calculate the image center coordinates, ensuring the consistency of the UAV data with the coordinate system of the GNSS monitoring system. This is the key to realizing the fusion of GNSS and UAV photogrammetry data. Based on the image data and the image center coordinates, aerial triangulation is carried out by combining the collinearity equations to restore the surface coordinates and obtain UAV-derived data. Then, the deformation area is determined by DEM differential of multi-period UAV-derived data, and discrete points on different flight dates in the deformation area are obtained. A joint measurement model is constructed with the GNSS elevation monitoring sequence data. The joint measurement model constructed by fusing GNSS and UAV in the present invention can not only obtain continuous and high-precision point displacement information, but also facilitate the acquisition of the continuous change trend of the planar shape, effectively compensating for the problem of insufficient spatial sampling caused by the "point-like" monitoring mode of GNSS.

[0015] Preferably, the double-difference model in step (1) is as follows:

[0016]

[0017] In the formula: represents double-difference calculation, r represents the GNSS monitoring station, and b represents the reference station; and are respectively the pseudorange and carrier phase observation values of the GNSS monitoring station r to the satellite s; and are respectively the pseudorange and carrier phase observation values of the GNSS monitoring station r to the satellite k; is the geometric distance between the satellite s and the GNSS monitoring station r, is the geometric distance between the satellite k and the GNSS monitoring station r; λ is the carrier wavelength; and are the carrier phase ambiguities; and are respectively the pseudorange multipath error and the carrier phase multipath error; ε P and ε Φ are respectively the observation noises of the pseudorange and the carrier phase;

[0018] From the position coordinates of the reference station b, and simultaneously based on and in formulas (1) and (2), calculate the distance from the GNSS monitoring station to the reference station, and obtain the GNSS monitoring station coordinates by inverse coordinate calculation, and then obtain the elevation monitoring sequence of the GNSS monitoring station.

[0019] Preferably, the calculation formula of the double-difference model in step (3) is the same as that in step (1), where b represents the reference station in the GNSS monitoring system, which is the same as the reference station in the double-difference model in step (1), except that r represents the airborne receiver of the UAV; based on the coordinates of the reference station b, the distance from the airborne receiver of the UAV to the reference station is calculated by and in the double-difference model, and the coordinates of the airborne receiver of the UAV are obtained by inverse coordinate calculation.

[0020] Preferably, since the frequency of the airborne receiver of the UAV is not synchronized with the camera exposure interval frequency, a cubic spline interpolation function is used to calculate the coordinates of the airborne receiver of the UAV at the exposure moment in step (3), and then the image center coordinates are obtained based on the coordinates of the airborne receiver of the UAV at the exposure moment combined with the data of the exposure record file; the cubic spline interpolation function is:

[0021]

[0022] This function divides an interval [a, b] into n sub-intervals through n + 1 data points, and constructs a cubic function S0(x)~S n-1 (x) on each sub-interval, and all functions pass through the known nodes within the interval; in the formula, the known nodes are the coordinates of the airborne receiver after solution, h i is the distance between two adjacent points at the sampling interval of the airborne receiver, M i is the second derivative at the coordinate node of the airborne receiver, and S i (x) is the coordinate position of the receiver at the required exposure moment.

[0023] Preferably, the method for obtaining the image center coordinates from the coordinates of the airborne receiver of the UAV at the exposure moment combined with the data of the exposure record file is as follows:

[0024] Based on the exposure moment coordinate compensation value in the exposure record file data, the coordinates of the airborne receiver at the exposure moment in the rectangular coordinate system are reduced and corrected to obtain the image center coordinates at the exposure moment.

[0025] Preferably, based on the image center coordinates at the exposure moment, the coordinates corresponding to in the double-difference model in step (3) are interpolated and transformed into the initial position of the UAV aerial photograph; the initial position of the UAV aerial photograph is used as the weighted approximate value of (X S , Y S , Z S ) in Equation 3, and is triangulated aerially together with the image data, and the ground coordinates of the object photographed by the image are obtained using the collinearity equation of Equation 3, that is, the UAV-derived data is obtained; the collinearity equation is as follows: The collinearity equation is as follows:

[0026]

[0027] Where: x and y are the coordinates of the image point; f is the camera focal length; x0 and y0 are the offsets of the principal point; a i , b i , c i (i = 1, 2, 3) are the three angular elements, the pitch angle the roll angle ω, and the direction cosines formed by the yaw angle κ; (X A , Y A , Z A ), (X S , Y S , Z S ) are the coordinates of the ground feature point A and the photography center S in the ground photogrammetry coordinate system, respectively.

[0028] Preferably, in step (5), a combined measurement model is constructed with the obtained GNSS elevation monitoring sequence data and the discrete points in the deformation area, including the following steps:

[0029] The first step: Average the second-level elevation monitoring values of GNSS by date to obtain the elevation measurement for the day.

[0030] The second step: Perform equally spaced discretization on the monitoring area in each period of DEM.

[0031] The third step: Calculate the elevation ratio for each discretized point and the GNSS monitoring point on the flight date, and establish a ratio change function through interpolation.

[0032] The fourth step: Calculate the elevation ratio of each discretized point to GNSS on the target monitoring date, and calculate the elevation value of each point in the area based on the GNSS monitoring value on that date.

[0033] The fifth step: Perform grid processing on all discrete points to obtain the elevation monitoring value of the monitoring surface on the target date; Use the piecewise linear interpolation method to establish the ratio change function, and establish the combined measurement model accordingly.

[0034] Preferably, the calculation formula for solving the elevation value of each point in the monitoring area based on the GNSS monitoring value on that date in the fourth step is as follows:

[0035] E d = F(d)·G d (12)

[0036] Where, E d is an n×1 elevation vector representing the elevation of the discretized points on the required date, F(d) is an n×1 ratio function vector representing the ratio change function of each discretized point to GNSS, and G dIt is a scalar representing the elevation of the GNSS monitoring point at the target date to be determined, where d is the target date.

[0037] Preferably, the formula of the joint monitoring model is as follows:

[0038]

[0039] Where, E d is an n×1 elevation vector representing the elevations of discrete points at the date to be determined, S t is the elevation ratio vector of the previous period in the two flight dates, S t+1 is the elevation ratio vector of the latter period in the two flight dates, G d is a scalar representing the elevation of the GNSS monitoring point at the target date to be determined, t is the flight frequency of the unmanned aerial vehicle. D is the time interval between the two flight dates; d is the time interval between the target date to be determined and the previous period in the two flight dates.

[0040] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0041] The present invention uses the reference station in the GNSS monitoring system as a reference benchmark for the airborne receiver of the unmanned aerial vehicle, calculates the position coordinates of the airborne receiver of the unmanned aerial vehicle using the double-difference model, calculates the coordinates of the airborne receiver of the unmanned aerial vehicle at the exposure moment using the cubic spline interpolation function, combines the coordinate compensation value of the photo-taking record file of the unmanned aerial vehicle to reduce the coordinates to the camera center to obtain the image center coordinates, so that the image center coordinates are unified in the GNSS monitoring system, and then combines the image data with the image center coordinates and the collinearity equation to perform aerial triangulation to obtain the ground coordinates of the object in the image. Based on the point cloud derived from multi-period unmanned aerial vehicle image data and the GNSS monitoring sequence, a joint monitoring model is constructed, and the elevation deformation variable of the deformation area during the monitoring period is calculated using the joint monitoring model and displayed in combination with the DOM, intuitively showing the elevation deformation of the landslide surface monitoring area, realizing continuous and intuitive monitoring of landslide deformation. The method of unifying the unmanned aerial vehicle image coordinates and the GNSS monitoring station coordinate frame adopted by the present invention can effectively solve the problem of systematic errors existing between the two different monitoring means.

[0042] The calculation results of the joint measurement model constructed by the present invention have high reliability and practicability, and can be used as an effective supplementary means for GNSS monitoring to participate in landslide monitoring. Thus, not only can continuous high-precision point displacement information be obtained, but also it is convenient to obtain the continuous change trend of the surface shape, effectively making up for the problem of insufficient spatial sampling caused by the "point-like" monitoring mode of GNSS. It more intuitively shows the overall deformation characteristics of the landslide body, provides strong data support for researchers to master the landslide evolution process and formulate effective prevention and control means. Description of the Drawings

[0043] Figure 1Flow chart of the landslide deformation continuous monitoring method integrating UAV photogrammetry and GNSS in Embodiment 1;

[0044] Figure 2 File for UAV photo taking records;

[0045] Figure 3 Flow chart for UAV image coordinate calculation;

[0046] Figure 4 Schematic diagram of landslide monitoring area and vectorization;

[0047] Figure 5 Diagram of displacement sequence and elevation calculation results of GNSS monitoring points;

[0048] Figure 6 Schematic diagram of elevation deformation of the landslide at different times;

[0049] Figure 7 Principle diagram of the combined measurement model algorithm;

[0050] Figure 8 Diagram of calculation of the detection surface and elevation change results on January 31 and February 24, 2024;

[0051] Figure 9 Comparison diagram of the positions of GNSS monitoring stations in DOM and actual positions;

[0052] Figure 10 Schematic diagram of GNSS monitoring stations;

[0053] Figure 11 3D deformation field dataset of the landslide body;

[0054] Figure 12 Statistical chart of the calculation results of the combined measurement model. Specific implementation manners

[0055] To better illustrate the purpose, technical solution and advantages of the present invention, the present invention will be further described below in conjunction with specific embodiments. Those skilled in the art should understand that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0056] Embodiment 1

[0057] This embodiment provides a landslide deformation monitoring method integrating UAV photogrammetry and GNSS. As Figures 1 - 12 shown, the overall flow schematic diagram is as Figure 1 shown, and it includes the following steps:

[0058] 1. Obtain the elevation monitoring sequence data of GNSS monitoring stations based on the double-difference model

[0059] Deploy multiple GNSS monitoring points in the monitoring area, set up reference stations in stable areas, and install GNSS monitoring stations in landslide - unstable areas. Use the real - time kinematic carrier - phase differential (RTK) technology for continuous monitoring and record the three - dimensional coordinate changes of the GNSS monitoring stations. The data sampling interval is set according to the landslide activity level and can be hourly or calculated daily.

[0060] RTK technology generally uses more than two GNSS receivers for synchronous observation, receiving pseudorange and carrier - phase signals. By performing inter - station and inter - satellite differential processing on the observed data, a double - difference observation model is established to eliminate the main errors and obtain high - precision relative positioning results. At the same time, when the baseline distance is short, the correlation of the atmospheric delay errors suffered by the two stations at both ends of the baseline is strong, and the ionospheric delay error and tropospheric delay error during the propagation of satellite signals can be ignored.

[0061] At a certain moment t, the reference station b and the GNSS monitoring station r simultaneously observe satellites s and k. The double - difference model can be expressed as:

[0062]

[0063] In the formula: represents double - difference calculation; and are respectively the pseudorange and carrier - phase observation values of the GNSS receiver r for satellite s; and are respectively the pseudorange and carrier - phase observation values of the GNSS receiver r for satellite k; is the geometric distance between satellite s and GNSS receiver r, is the geometric distance between satellite k and GNSS receiver r; λ is the carrier wavelength; and are the carrier - phase ambiguities; and are respectively the pseudorange multipath error and carrier - phase multipath error; ε P and ε Φ are respectively the observation noises of the pseudorange and carrier - phase.

[0064] Deploy multiple GNSS receivers near the area to be monitored, and a regional GNSS monitoring system can be constructed. Usually, one or more GNSS reference stations are set up in stable areas, and multiple GNSS monitoring stations are arranged in areas with high landslide risks. Based on the position coordinates of the reference station b and simultaneously using and in formulas (1) and (2) to calculate the distance from the GNSS monitoring station to the reference station, and obtain the coordinates of the GNSS monitoring station through inverse coordinate calculation, and then obtain the displacement sequence information of the GNSS monitoring station (see Figure 5The left figure) and the elevation monitoring sequence (see Figure 5 The right figure).

[0065] Evaluating the stability state of a landslide through the displacement in the elevation direction of a GNSS monitoring station and providing early warning information for possible landslide disasters is the conventional processing method for landslide monitoring by existing GNSS monitoring stations. Its monitoring range is limited to point targets and it is difficult to reflect the dynamic changes of the overall landslide instability surface.

[0066] To address this problem, the present invention further involves UVA data collection and unifying the coordinate systems of the GNSS monitoring station coordinates and the UAV image coordinates. Furthermore, based on multi-period UAV image data and GNSS monitoring sequences, a joint monitoring model is constructed to enable continuous monitoring of the elevation changes in the landslide area, thereby accurately grasping the evolution process of the entire landslide surface, as follows.

[0067] II. Obtaining UAV on-board receiver data, exposure record file data, and image data through UAV photogrammetry

[0068] Use a UAV equipped with a high-precision camera to collect images of the landslide area. The present invention uses a DJI M300RTK UAV equipped with a Zenmuse P1 camera to collect UAV image data. The on-board receiver of the UAV is connected to a CORS station to perform fixed-frequency calculation of the coordinates of the UAV on-board receiver. At the moment when the UAV on-board camera takes a photo, an electrical pulse signal is generated, and the built-in sensor of the camera records this moment as GPS time and saves it as a photo record file. Taking the DJI UAV as an example, its exposure record file data is as Figure 2 shown, including data such as the photo point serial number, GPS seconds within a week (i.e., the photo-taking moment), GPS week, compensation value in the north direction, compensation value in the east direction, and compensation value in the elevation direction. Based on the coordinates of the on-board receiver and the exposure record file data, the coordinates of the image can be calculated.

[0069] III. Using the reference station in the GNSS monitoring system to provide a reference benchmark for the UAV on-board receiver, obtaining the coordinates of the UAV on-board receiver based on the double-difference model, and combining the exposure record file data to obtain the image center coordinates to achieve the unification of the coordinate frameworks of the GNSS monitoring station and the UAV

[0070] Since the data collected by GNSS forms an independent coordinate system by calculating the GNSS monitoring station coordinates through the reference station and the monitoring station using the double-difference model, the GNSS data adopts the global coordinate system (WGS84). The data collected by the UAV includes images and image coordinate data. The image POS data is calculated through the combination of the on-board receiver of the UAV and the CORS station using the double-difference model and the exposure record file, and then the ground objects captured in the image are restored to the ground coordinate system to produce the digital orthophoto map (DOM) and the digital terrain model (DSM). The generated DSM and DOM usually adopt the projected coordinate system (such as UTM or local coordinate system). It can be seen that the GNSS data and the UAV image data are two independent coordinate systems, and it is difficult to directly fuse the displacement data collected by GNSS with the image data of the UAV. Therefore, it is necessary to perform coordinate system conversion and unification on these two types of data, namely GNSS and UAV, to ensure the spatial matching of the data.

[0071] In order to unify the coordinate system, during the calculation process of converting the images and image coordinates collected by the UAV to the surface coordinate system through aerial triangulation and bundle adjustment (collinearity equation, Equation 3), the present invention uses the reference station in the GNSS monitoring system as a reference benchmark for the on-board receiver of the UAV, and uses the double-difference model to calculate the position of the on-board receiver of the UAV. Then, cubic spline interpolation and the exposure record file are used to calculate the image coordinates, and aerial triangulation is carried out in combination with the collinearity equation to restore the surface coordinates, ensuring that the coordinate system of the UAV data is consistent with that of the GNSS monitoring system. The specific principle is as follows:

[0072] There is a fixed geometric relationship between the aerial photographic film of UAV photogrammetry and the ground. The parameters determining the position relationship between the aerial camera and the film are the interior orientation elements, and the parameters describing the spatial position of the photographic light rays are the exterior orientation elements. Using these interior and exterior orientation elements, the geometric relationship between the image points on the aerial photographic film and the corresponding ground object points can be established, that is, the central projection collinearity equation, as shown in Equation 3 below:

[0073]

[0074] In the formula: x, y are the image point coordinates; f is the camera focal length; x0, y0 are the offsets of the principal point; a i , b i , c i (i = 1, 2, 3) are the direction cosines composed of the three angular elements of pitch roll angle ω, yaw angle κ; (X A , Y A , Z A ), (X S , Y S , Z S ) are the coordinates of the ground object point A and the photography center S in the ground photogrammetry coordinate system respectively.

[0075] The method for overall solving the photo orientation elements and the object space coordinates of the points to be determined using Equation (3) is the bundle adjustment of the regional network. However, during the process of bundle adjustment of the regional network, due to the strong correlation between the interior orientation elements f, x0, y0 of the photo and the exterior orientation elements X S , Y S , Z S , ω, κ, the normal equation obtained from Equation 3 by the least squares principle is ill-conditioned, resulting in unstable solutions for the unknowns sought.

[0076] By using GNSS technology, the corresponding coordinates in Equations (1) and (2) can be transformed into the initial position of the UAV aerial photo through interpolation and other methods. Taking the initial position of the photo as the weighted approximate value of (X S , Y S , Z S ) can greatly offset the correlation between the interior and exterior orientation elements, significantly improve the accuracy of object points, and restore accurate terrain features.

[0077] The mathematical models used for solving the coordinates of the monitoring stations in the GNSS monitoring system and the coordinates of the UAV-borne receiver are the same, both being the double-difference model shown in Equations (1) and (2). Therefore, the key to unifying the UAV photogrammetry coordinate system and the GNSS monitoring framework lies in using the information of the same reference station. Thus, in the present invention, the reference station in the GNSS monitoring system is used to provide a reference benchmark for the UAV-borne receiver, and the double-difference model is used to calculate the coordinates of the UAV-borne receiver.

[0078] Since the frequency of the UAV-borne receiver is fixed, and the exposure interval of the camera is not synchronized with this frequency, in order to obtain the coordinates of the image at the exposure moment, it is necessary to perform data interpolation processing by means of the photo-taking record file. Different interpolation methods will directly affect the accuracy of the coordinates. Cubic spline interpolation has good convergence, smooth curves, and reliable interpolation results. Therefore, in this embodiment, cubic spline interpolation is used to calculate the position of the UAV-borne receiver at the exposure point, and its basic principle is as follows:

[0079] The basic principle of this interpolation method is as follows:

[0080] It is known that n + 1 data points divide an interval [a, b] into n sub-intervals, and a cubic function S0(x) ～ S n-1 (x) is constructed on each sub-interval. Each cubic function is defined by the function in Equation (4):

[0081] S i (x) = a i + b i (x - xi ) + c i (x - x i ) 2 + d i (x - x i ) 3 , i = 0, 1, …, n - 1 (4)

[0082] Each cubic function has four coefficients a, b, c, and d. So, there are a total of 4n coefficients in equation (5). Therefore, 4n equations are needed to solve. Four necessary conditions are required for the solution:

[0083] First, all cubic functions must pass through the known nodes, that is:

[0084] S i (x i ) = y i , i = 0, 1, 2, …, n (5)

[0085] Second, except for the two end nodes, the function is continuous at all n - 1 nodes, that is:

[0086] S i (x i+1 ) = S i+1 (x i+1 ), i = 0, 1, 2, …, n - 2 (6)

[0087] Third, except for the two end nodes, the function is first - order differentiable at all n - 1 nodes, that is:

[0088] S′ i (x i+1 ) = S′ i+1 (x i+1 ), i = 0, 1, 2, …, n - 2 (7)

[0089] Fourth, except for the two end nodes, the function is second - order differentiable at all n - 1 nodes, that is:

[0090] S″ i (x i+1 ) = S″ i+1 (x i+1 ), i = 0, 1, 2, …, n - 2 (8)

[0091] From the above four conditions, 4n - 2 equations can be listed. The remaining two equations are composed of two boundary conditions, such as:

[0092] S″(x0) = 0 = S″(x n ) (9)

[0093] Define the second - order derivative M i = S″ i (xi ) A tridiagonal system of equations can be derived through the second - order continuity condition:

[0094]

[0095] where h i is the step size, and h i = x i+1 - x i . After solving for M i , the cubic spline interpolation function within the interval [x i , x i+1 is:

[0096]

[0097] Taking the DJI drone as an example, its photo - taking record file is as Figure 2 shown. Among them, columns 4, 5, and 6 are the coordinate compensation values between the on - board receiver of the drone and the exposure center at the exposure moment. These values are the angular element auxiliary solution results measured by the drone's IMU (Inertial Measurement Unit) system and can be used for reduction and correction when calculating the coordinates of the image at the exposure moment.

[0098] The solution process of the coordinates of the drone image at the exposure moment is as Figure 3 shown. As can be seen from Figure 3 , the specific steps for solving the coordinates of the sensor exposure center of the drone image at the exposure moment are as follows:

[0099] ① Construct a double - difference model through the reference station observation file, the drone on - board receiver observation file, and the satellite ephemeris to calculate the coordinates of the drone on - board receiver. That is, the base station is the reference station b in the GNSS monitoring system, and the drone on - board receiver is r; the on - board receiver can calculate the distance from the drone to the reference station according to and in equations (1) and (2), and calculate the coordinates of the drone on - board receiver by inverse coordinate calculation based on the reference station coordinates;

[0100] ② Interpolate and calculate the coordinates of the drone on - board receiver at that moment according to the exposure moment recorded in the photo - taking record file. As shown in equation (11), the known nodes are the coordinates of the on - board receiver after solution, h i is the distance between two adjacent points at the sampling interval of the on - board receiver, M i is the second - order derivative at the coordinate node of the on - board receiver, and S i (x) is the coordinate position of the receiver at the required exposure moment;

[0101] ③According to the exposure moment coordinate compensation values recorded in the photographing record file, the coordinates of the airborne receiver in the rectangular coordinate system are reduced to the camera exposure center to obtain the image center coordinates at the exposure moment. For example, the compensation value in the northeast direction of the 4th column is -552, so this value is subtracted during the calculation of the exposure moment coordinates for reduction and correction; for example, the compensation value in the northeast direction of the 5th column is -28, so this value is subtracted during the calculation of the exposure moment coordinates for reduction and correction; for example, the compensation value in the elevation direction of the 6th column is 250, so this value is added during the calculation of the exposure moment coordinates for reduction and correction.

[0102] Through the above steps ①②③, the image center coordinates at the exposure moment are obtained, and the initial position of the image to be found is used as the weighted approximate value of (X S , Y S , Z S ) in Equation 3. Together with the image data, aerial triangulation is carried out, and the ground point coordinates are restored using the collinearity equation of Equation 3, that is, the ground coordinates of the imaged object are obtained.

[0103] IV. Under the unified coordinate framework, the image data and the image center coordinates are used to obtain UAV-derived data through aerial triangulation. Then, based on the DEM difference of multiple periods of UAV-derived data, the deformation area is determined, and discrete points on different flight dates in the deformation area are obtained.

[0104] After the unification of the UAV and GNSS coordinate frameworks, the ground coordinates of the objects photographed by the UAV and the coordinate system of the GNSS monitoring station (the monitoring station is originally on the ground surface) have been unified, and further fusion of the monitoring results can be carried out.

[0105] First, the deformation area needs to be determined. The application of the Digital Elevation Model (DEM) is a key step in determining the elevation change area. By comparing the DEMs of different dates, the areas with elevation changes can be identified. In addition, means such as Digital Orthophoto Map (DOM) and on-site inspections are needed to further determine the research area. For example, in this embodiment, three main potential hazard areas are determined (as Figure 4 shown). Among them, the elevation change in the blue area is caused by snowmelt, and this area is located at the foot of the slope, so it is not within the scope of this monitoring. The elevation change in the green area is caused by human factors. The red area is the location of the GNSS monitoring station, which represents the landslide instability area to be monitored. In order to further analyze and monitor these areas, the scope of the monitoring area needs to be vectorized.

[0106] V. A joint measurement model is constructed with the GNSS elevation monitoring sequence data and the discrete points in the deformation area obtained by the UAV.

[0107] For the monitoring area, this study adopts the following ideas:

[0108] The first step: Calculate the average value of the second-level elevation monitoring values of GNSS by date as the elevation monitoring quantity for the day.

[0109] The second step: Discretize the monitoring area in each period of DEM at equal intervals.

[0110] The third step: Calculate the elevation ratio of each discretized point to the GNSS monitoring point on the flight date, and establish a ratio change function through interpolation.

[0111] The fourth step: Calculate the elevation ratio of each discretized point to GNSS on the target monitoring date, and use the GNSS monitoring value on this date to calculate the elevation value of each point in this area using Equation (12).

[0112] The fifth step: Perform grid processing on all discretized points to obtain the elevation monitoring value of the monitoring surface on the target date. The algorithm principle is as Figure 7 shown.

[0113] E d = F(d)·G d (12)

[0114] where E d is an n×1 elevation vector representing the elevation of the discretized points on the required date, expressed as Equation (13). F(d) is an n×1 ratio function vector representing the ratio change function of each discretized point to GNSS, expressed as Equation 14. G d is a scalar representing the elevation of the GNSS monitoring point on the required target date, n is the number of discretized points, and d is the target date.

[0115]

[0116] During the experiment, due to the long interval between flight dates and the relatively small number of flights, the results of multiple interpolation functions are unstable. Therefore, this paper uses the piecewise linear interpolation method to establish the ratio change function and establish a joint monitoring model as shown in Equation (15).

[0117]

[0118] where S t is the elevation ratio vector of the previous period in two flight dates, expressed by Equation (16). S t+1 is the elevation ratio vector of the latter period in two flight dates, expressed by Equation (17). is the elevation of each discretized point in the previous period of the two flight dates. is the elevation of each discretized point in the latter period of the two flight dates. G tis the elevation of the GNSS monitoring point in the previous period of the two flight dates. G t+1 is the elevation of the GNSS monitoring point in the latter period of the two flight dates. t is the flight frequency of the UAV. D is the time interval between the two flight dates; d is the time interval between the target date to be calculated and the previous period of the two flight dates.

[0119]

[0120] VI. Based on the combined measurement model, obtain the areal monitoring data for any date, construct the areal deformation dataset of the landslide, visualize the landslide monitoring results, and achieve continuous and intuitive monitoring of the landslide deformation.

[0121] Through the above method, using the multi-period UAV photogrammetry results and the GNSS time series monitoring values, the elevation monitoring quantity of the landslide instability area at any time can be calculated, and a relatively continuous areal monitoring result in time can be obtained through high-density calculation. It should be noted that the deformation quantity of the landslide is small within a short period of time, and the deformation quantity of the results obtained by high-density calculation for adjacent dates may be much smaller than the grid resolution. Therefore, the time interval can be adjusted according to the deformation quantity and the actual monitoring requirements.

[0122] Through the above method, using the multi-period UAV photogrammetry results and the GNSS time series monitoring values, the elevation monitoring quantity of the landslide instability area can be continuously obtained.

[0123] The present invention calculates the areal monitoring data on March 20, July 1, September 1, and November 1, 2023 according to the combined measurement model. The grid data is generated based on the calculated discrete points and is respectively differentiated from the DEM on April 1 to obtain the elevation deformation quantity within this time period. The areal monitoring results are combined with the DOM for display to intuitively show the elevation change and deformation of the landslide. The specific situation is as Figure 6 shown. A monitoring station B was installed in this area on December 28 to verify the accuracy and expansion ability of the combined measurement model (i.e., calculate the data outside the UAV flight time interval). The present invention calculates the areal monitoring data on January 3, 2024 and February 24, 2024. The differential results of the calculation results and the DEM on December 28, 2023 are shown in Figure 8 , the green area is the monitoring result on January 31, 2024, and the blue area is the monitoring result on February 24, 2024. It can be Figure 8 seen that the shape of the monitoring area calculated by this model is the same as the actual one, the deformation quantity is consistent with the monitoring quantity obtained by the monitoring station B, and the deformation process of the areal area is accurately restored.

[0124] Based on the coordinate unified framework and the combined measurement model, the present invention constructs a three-dimensional deformation field dataset of the landslide body using the combined monitoring model ( Figure 11) It more intuitively shows the overall deformation characteristics of the landslide body, effectively supplementing the problem of insufficient spatial sampling caused by the GNSS "point-like" monitoring mode.

[0125] Example 2

[0126] In this example, the experimental verification of the landslide monitoring effect is carried out based on the joint measurement model and the areal deformation data set constructed in Example 1. The specific test method is as follows:

[0127] ① To visually show the differences between the UAV and the GNSS monitoring framework in different situations and the effectiveness of the unified monitoring framework method, this paper uses the DJI M300RTK UAV equipped with the Zenmuse P1 camera to collect image data in three different flight modes.

[0128] Method 1 (single-point positioning): On February 14th, the single-point positioning flight mode was adopted, that is, relying on single-point positioning to obtain the initial position of the image without using any reference station information.

[0129] Method 2 (network RTK): In the UAV operation on April 1st, the network RTK mode was adopted, using the virtual reference station provided by the CORS station to solve the initial position of the UAV image in real time.

[0130] Method 3 (self-built RTK): In the flight operation on December 28th, the on-board receiver of the UAV was connected to the GNSS monitoring system cloud server through the 4G network to receive the regional reference station information (HF01), and the position of the UAV image was solved in real time.

[0131] By comparing the positions of the GNSS monitoring stations in the DOM generated from the UAV image coordinates before and after the solution of the above three methods with the actual position of the GNSS monitoring station A, the effectiveness of the coordinate framework unification method is verified. The results are as Figure 9 shown, and the specific values are shown in Table 1. From Figure 9 Combining the results in Table 1, it can be obtained that in the UAV operation on February 14th, the single-point positioning flight mode was used. The planar distance between the position of the GNSS monitoring station in the DOM before the solution and the actual monitoring station position was 4.69 meters, and this result is consistent with the error of single-point positioning. In the flight on December 28th, the self-built RTK flight mode was adopted, and the planar distance before the solution was reduced to 1.33 meters, reflecting higher positioning accuracy. However, for the network RTK adopted on April 1st, the difference before the solution reached 70.63 meters, which was mainly due to the difference between the CORS coordinate system and the GNSS monitoring framework.

[0132] Through the method described in this article, after the calculations on February 14th, April 1st, and December 28th, the positions of the GNSS monitoring stations in the DOM generated from the UAV coordinates were reduced to 0.2 meters, 0.07 meters, and 0.06 meters from their actual positions. This result indicates that the method of unifying the UAV image coordinates with the GNSS monitoring station coordinate framework in the present invention can greatly improve the positioning accuracy of GNSS monitoring stations.

[0133] In contrast, the accuracy on February 14th is poorer than that of the other two periods. The main reason is that when designing the UAV flight route, this area was at the edge of the route, and the monitoring station has a certain height, resulting in a certain deviation in the position of the GNSS monitoring station in the DOM. It should be noted that the GNSS monitoring station is a rigid whole, and the displacement at any position of this whole is the same. Therefore, as long as the actual position coordinates of the GNSS monitoring station fall within the range of the monitoring station in the DOM, its planar error has no impact on elevation monitoring. From the experimental results, within the allowable error range, this method can effectively combine the UAV photogrammetry coordinate system with the GNSS monitoring framework, and the error after calculation has little impact on joint monitoring.

[0134] Table 1 Distance between GNSS monitoring stations on DOM and actual GNSS monitoring stations (unit: meter)

[0135] Operation Date February 14th April 1st December 28th Flight Mode Single Point Positioning Network RTK Self - built RTK Before Resolution 4.69 70.63 1.33 After Resolution 0.2 0.07 0.06

[0136] ② To verify the reliability and accuracy of the calculation results of the joint measurement model, the areal monitoring data on March 20th, July 1st, September 1st, November 1st, January 31st, and February 24th, 2023 obtained on the calculation dates were statistically analyzed, and the elevation average values of different dates of the monitoring area were compared with the GNSS monitoring sequences during this time period to analyze their changing trends ( Figure 12 ), where the pole height of 1.5 meters has been subtracted from the monitoring sequence of monitoring station A ( Figure 10 ).

[0137] From Figure 12The results show that the elevation changes recorded by GNSS monitoring station A are consistent with the elevation change trend of the monitoring surface calculated by the combined measurement model. Taking April 1, 2023, the date of the second drone flight, as the reference date, during the monitoring period, the average elevation changes of the monitoring surface relative to the reference date at different time nodes are as follows: 0.10 m on March 20, -0.12 m on July 1, -0.15 m on September 1, and -0.18 m on November 1. Correspondingly, the GNSS monitoring values are: 0.088 m on March 20, -0.117 m on July 1, -0.132 m on September 1, and -0.182 m on November 1. This indicates that the elevation change trend of the monitoring surface is basically consistent with the change trend of the GNSS elevation monitoring values, with some minor differences. It should be noted that the described change amount of the monitoring surface is the average value of each point, which is different from the monitoring results of a single point of the GNSS monitoring station, and this is reasonable.

[0138] To verify the reliability and practicability of the calculation results of the combined measurement model, the research team installed a low-cost monitoring station B on the monitoring surface on December 28 for verification. The verification method is to obtain the calculated elevation values of the monitoring surface corresponding to the planar positions of monitoring station B on January 31, 2024, and February 24, 2024, and compare them with the actual elevation monitoring results of B. The results are as Figure 8 shown.

[0139] The deformation trend of the actual monitoring sequence of station B is also basically consistent with the calculation results of the combined measurement model. The specific numerical results are shown in Table 2.

[0140] Table 2 Calculation Results and Monitoring Results of Monitoring Station B (Unit: m)

[0141]

[0142] According to Table 2, the elevation change amounts of B calculated by the combined measurement model for the required dates are: 2 cm and 6 cm; while the actual elevation monitoring results of A01 are: 2 cm and 5.4 cm. The error between the two is less than 1 cm. This indicates that the calculation results of the combined measurement model of the present invention have high reliability and practicability, and can be used as an effective supplementary means for GNSS monitoring to participate in landslide monitoring. Thus, not only can continuous high-precision point displacement information be obtained, but also it is convenient to obtain the continuous change trend of the surface shape, providing strong data support for researchers to master the landslide evolution process and formulate effective prevention and control measures.

[0143] The above only discloses several specific embodiments of the present invention. Those skilled in the art can make various changes and modifications to the embodiments of the present invention without departing from the spirit and scope of the present invention. However, the embodiments of the present invention are not limited thereto, and any changes that can be thought of by those skilled in the art should fall within the protection scope of the present invention.

Claims

1. A continuous landslide deformation monitoring method integrating UAV photogrammetry and GNSS, characterized in that, It includes the following steps: (1) In the GNSS monitoring system, based on the position information of the reference station and the GNSS monitoring stations, the elevation monitoring sequence data of the GNSS monitoring stations is obtained based on the double-difference model; (2) UAV photogrammetry is used to obtain UAV airborne receiver data, exposure record file data, and image data; (3) The reference station in the GNSS monitoring system provides a reference benchmark for the UAV airborne receiver. Based on the double-difference model, the coordinates of the UAV airborne receiver are obtained. Combining the exposure record file data, the image center coordinates are obtained to achieve the unification of the coordinate frameworks of the GNSS monitoring stations and the UAV; (4) Under the unified coordinate framework, the image data and the image center coordinates are used to obtain UAV-derived data through aerial triangulation. Then, the deformation area is determined by DEM difference of multiple periods of UAV-derived data, and discrete points on different flight dates in the deformation area are obtained; (5) A joint measurement model is constructed with the GNSS elevation monitoring sequence data obtained in step (1) and the discrete points on different flight dates in the deformation area obtained in step (4); (6) Based on the joint measurement model, the areal monitoring data on any date is obtained, and an areal deformation dataset is constructed to visualize the landslide monitoring results and achieve continuous and intuitive monitoring of landslide deformation.

2. The continuous landslide deformation monitoring method integrating UAV photogrammetry and GNSS according to claim 1, characterized in that, The double-difference model in step (1) is: In the formula: represents double-difference calculation, r represents the GNSS monitoring station, and b represents the reference station; and are respectively the pseudorange and carrier phase observation values of the GNSS monitoring station r for the satellite s; and are respectively the pseudorange and carrier phase observation values of the GNSS monitoring station r for the satellite k; is the geometric distance between the satellite s and the GNSS monitoring station r, is the geometric distance between the satellite k and the GNSS monitoring station r; λ is the carrier wavelength; and are the carrier phase ambiguities; and are respectively the pseudorange multipath error and the carrier phase multipath error; ε P and ε Φ are the observation noises of the pseudorange and the carrier phase respectively; Based on the position coordinates of the reference station b and simultaneously based on the and in Equations (1) and (2), calculate the distance from the GNSS monitoring station to the reference station, obtain the coordinates of the GNSS monitoring station through inverse coordinate calculation, and then obtain the elevation monitoring sequence of the GNSS monitoring station.

3. The continuous landslide deformation monitoring method integrating UAV photogrammetry and GNSS according to claim 1, characterized in that The calculation formula of the double-difference model described in step (3) is the same as that of the double-difference model described in step (1), where b represents the reference station in the GNSS monitoring system, which is the same as the reference station in the double-difference model in step (1). The difference is that r represents the UAV-borne receiver; based on the coordinates of the reference station b, the distances from the UAV-borne receiver to the reference station are calculated from and in the double-difference model, and the coordinates of the UAV-borne receiver are obtained by inverse coordinate calculation.

4. The continuous landslide deformation monitoring method integrating UAV photogrammetry and GNSS according to claim 1, characterized in that, Since the frequency of the UAV airborne receiver is not synchronized with the camera exposure interval frequency, a cubic spline interpolation function is used to calculate the coordinates of the UAV airborne receiver at the exposure moment in step (3). Then, based on the coordinates of the UAV airborne receiver at the exposure moment and the exposure record file data, the image center coordinates are obtained; The cubic spline interpolation function is: This function divides an interval [a, b] into n sub-intervals through n + 1 data points, and constructs a cubic function S0(x)~S n-1 (x) on each sub-interval, and all functions pass through the known nodes within the interval; where the known nodes are the coordinates of the airborne receiver after solution, h i is the distance between two adjacent points at the sampling interval of the airborne receiver, and M i is the second derivative at the coordinate node of the airborne receiver, and S i (x) is the coordinate position of the receiver at the required exposure time.

5. The continuous landslide deformation monitoring method integrating UAV photogrammetry and GNSS according to claim 4, characterized in that, The method for obtaining the image center coordinates by combining the coordinates of the UAV airborne receiver at the exposure moment and the exposure record file data is as follows: Based on the exposure moment coordinate compensation value in the exposure record file data, the coordinates of the airborne receiver at the exposure moment in the rectangular coordinate system are reduced and corrected to obtain the image center coordinates at the exposure moment.

6. The continuous landslide deformation monitoring method integrating UAV photogrammetry and GNSS according to claim 4, characterized in that Based on the image center coordinates at the exposure moment, the corresponding coordinates in the double-difference model described in step (3) are interpolated and transformed into the initial position of the UAV aerial photograph; the initial position of the UAV aerial photograph is used as the weighted approximate value of (X , Y S , Z S , Z S ) in Equation 3, and together with the image data, aerial triangulation is carried out. The ground coordinates of the object photographed by the image are obtained by using the collinearity equation of Equation 3, that is, the UAV-derived data is obtained; the collinearity equation is as follows: Where: x and y are the image point coordinates; f is the camera focal length; x0 and y0 are the offsets of the principal point; a i , b i , c i (i = 1, 2, 3) are the three angular elements, the pitch angle , the roll angle ω, and the direction cosines formed by the yaw angle κ; (X A , Y A , Z A ), (X S , Y S , Z S ) are the coordinates of the ground object point A and the photography center S in the ground photogrammetry coordinate system, respectively.

7. The continuous landslide deformation monitoring method integrating UAV photogrammetry and GNSS according to claim 1, characterized in that, The step (5) of constructing a joint measurement model with the obtained GNSS elevation monitoring sequence data and the discrete points in the deformation area includes the following steps: The first step: The second-level elevation monitoring values of GNSS are averaged by date as the elevation monitoring quantity for that day; The second step: The monitoring area in each period of DEM is discretized at equal intervals; The third step: For each discretized point, the elevation ratio with the GNSS monitoring point on the flight date is obtained, and a ratio change function is established by interpolation; The fourth step: Calculate the elevation ratio of each discrete point to GNSS on the target monitoring date, and calculate the elevation value of each point in the area based on the GNSS monitoring value on that date; The fifth step: All discrete points are gridded to obtain the elevation monitoring value of the monitoring surface in that target date; A piecewise linear interpolation method is used to establish a ratio change function, and a joint monitoring model is established based on this.

8. The continuous landslide deformation monitoring method integrating UAV photogrammetry and GNSS according to claim 7, wherein, The calculation formula for calculating the elevation value of each point in the monitoring area based on the GNSS monitoring value on that date in the fourth step is as follows: E d = F(d)·G d (12) Among them, E d is an n×1 elevation vector representing the elevation of discrete points on the required date, F(d) is an n×1 proportional function vector representing the proportional change function of each discrete point with GNSS, and G d is a scalar representing the elevation of the GNSS monitoring point on the required target date, and d is the target date.

9. The continuous landslide deformation monitoring method integrating UAV photogrammetry and GNSS according to claim 7, characterized in that, The formula for the joint monitoring model is as follows: Among them, E d is an n×1 elevation vector representing the elevation of discrete points on the required date, S t is the elevation ratio vector of the previous period among the two flight dates, S t+1 is the elevation ratio vector of the latter period among the two flight dates, G d is a scalar representing the elevation of the GNSS monitoring point on the required target date, t is the flight frequency of the UAV, D is the time interval between the two flight dates; d is the time interval between the required target date and the previous period among the two flight dates.

Images