Earth surface deformation forward and reverse modeling method and system based on mogi model

Abstract

The invention provides an earth surface deformation forward and reverse modeling method and system based on a mogi model. The method comprises the steps of obtaining an observation data file and a mogi point source parameter constraint file; determining deformation parameters of observation points from the observation data file, and obtaining regional optimal mogi point source parameters based on nonlinear least square inversion by combining the mogi point source parameter constraint file; position parameters of observation points are determined from the observation data file, mogi forward modeling is carried out based on regional optimal mogi point source parameters, and a theoretical earth surface deformation field is obtained; and calculating a residual index according to the theoretical earth surface deformation field and the deformation observation value of the observation point, and outputting the theoretical earth surface deformation field and the residual index which meet a preset evaluation requirement.

Description

Technical Field

[0001] This invention belongs to the field of Earth science, specifically relating to a method and system for forward and inverse modeling of surface deformation based on the Mogi model. Background Technology

[0002] Obtaining accurate surface deformation data is of great significance for studying volcanic activity and predicting volcanic eruptions. The Mogi model, proposed by Japanese scholar Kiyoo Mogi, is a model linking magma pressure sources and volcanic surface deformation. This model has important applications in simulating and interpreting deformation in numerous volcanic areas and has been widely used in deformation studies of volcanic regions worldwide. Furthermore, the model has been extended to studies of surface deformation caused by earthquake precursors, surface deformation caused by underground fluid injection and extraction, and deformation of mine collapses and goaf areas. Analyzing the surface deformation of the Mogi model can help understand the dynamic changes of underground Mogi point sources, which is helpful for quantitatively analyzing volcanic magma chamber activity, fluid migration, and mineral reservoir mining. Studying the surface deformation of the Mogi model can further predict volcanic activity, calculate fluid migration, and simulate mining-induced earthquakes. Currently, with the improvement of observation technology, especially the development of GNSS all-weather observation and InSAR and leveling observation technologies, data acquisition has become easier. In terms of data processing, nonlinear least squares is the main method for inverting the point source parameters of the Mogi model. For multi-source observation data (InSAR, GNSS, and leveling observations), few existing technologies provide surface deformation analysis using the Mogi model, and some lack inversion functionality. Summary of the Invention

[0003] To overcome the shortcomings of existing technologies in surface deformation analysis using the Mogi model without multi-source observation data and lacking inversion functions, this invention provides a forward and inverse surface deformation method and system based on the Mogi model. It relies on weighted nonlinear least squares inversion to determine point source parameters, sets upper and lower bounds and initial values ​​of parameters through the Mogi point source parameter constraint file, and improves the rationality of the evaluation process to verify the results. It can intuitively and quickly find the optimal solution or verify the reliability of the solution.

[0004] According to one aspect of this specification, a method for forward and inverse modeling of land surface deformation based on the Mogi model is provided, comprising:

[0005] Obtain the observation data file and the Mogi point source parameter constraint file;

[0006] The deformation parameters of the observation points are confirmed from the observation data file. Based on nonlinear least squares inversion, the optimal Mogi point source parameters for the region are obtained by combining the Mogi point source parameter constraint file.

[0007] The location parameters of the observation points are confirmed from the observation data file. Based on the optimal Mogi point source parameters in the region, Mogi forward modeling is performed to obtain the theoretical surface deformation field.

[0008] The residual index is calculated based on the theoretical surface deformation field and the deformation parameters at the observation points, and the theoretical surface deformation field and residual index that meet the preset evaluation requirements are output.

[0009] As a further technical solution, the inversion process, which combines the deformation parameter data of each observation point with the Mogi point source parameter constraint file, includes:

[0010] Read the deformation parameter data of each observation point in the observation data file, construct the observation data matrix, and construct the weight matrix according to the weight of each type of observation data file;

[0011] Determine the number of mogi point sources and determine the range of mogi parameter values ​​based on the mogi point source parameter constraint file;

[0012] Set a reference central meridian and a reference ellipsoid, and calculate the coordinates of each observation point and the coordinate range of each Mogi point source in the Gaussian plane rectangular coordinate system using Gaussian projection;

[0013] Based on the observation data matrix, weight matrix, and range of Mogi parameter values, the optimal Mogi point source parameters for the region are obtained using nonlinear least squares inversion.

[0014] As a further technical solution, the process based on nonlinear least squares inversion includes:

[0015] The sampling number for nonlinear least squares inversion is set, and the inversion is performed according to the minimization criterion of weighted least squares joint adjustment with weighted ratio factor to obtain the optimal mogi parameter matrix that satisfies the range of mogi parameter values ​​and the coordinate range of mogi point source.

[0016] The coordinates of the corresponding Mogi point sources in the Gaussian Cartesian coordinate system are converted into longitude and latitude in the matrix elements of the optimal Mogi parameter matrix. Then, all matrix elements are derived to obtain the optimal Mogi point source parameters for the region.

[0017] As a further technical solution, the mathematical expression of the minimization criterion for weighted least squares joint adjustment with weighted ratio factors is as follows:

[0018]

[0019] In the formula, min represents minimizing the objective function; W is the weight matrix; This is the deformation observation matrix; The design matrix includes the partial derivatives of the mogi point source parameters with respect to the deformation observations; This is the source parameter matrix for the Mogi points.

[0020] As a further technical solution, the parameters of the Mogi point source include longitude, latitude, the volume change of the Mogi point source, and the depth of the Mogi point source.

[0021] As a further technical solution, the steps for performing Mogi forward modeling to obtain the theoretical surface deformation field include:

[0022] Read several observation points and their corresponding location parameters from the observation data file;

[0023] Based on the optimal parameters of the regional mogi point sources, construct the parameter matrix of each mogi point source;

[0024] A reference central meridian and a reference ellipsoid are defined. Based on the position parameters of each observation point and the parameter matrix of each Mogi point source, the projected coordinates of each observation point and the Mogi point source in the Gaussian plane rectangular coordinate system are calculated by Gaussian projection.

[0025] Calculate the radial distance between each observation point and each Mogi point source on the Earth's surface. Based on the Mogi model, calculate the deformation field of each simulated Mogi point source. Superimpose the deformation fields simulated by all Mogi point sources to obtain the theoretical surface deformation field.

[0026] As a further technical solution, a method for forward and inverse modeling of surface deformation based on the Mogi model also includes a separate forward modeling process:

[0027] Confirm the location parameters of the observation points from the observation data file, obtain the actual Mogi point source parameters to replace the optimal Mogi point source parameters obtained by inversion in the region, perform Mogi forward modeling, and obtain the theoretical surface deformation field.

[0028] According to another aspect of this specification, a forward and inverse surface deformation modeling system based on the Mogi model is provided, comprising:

[0029] The data import module is used to acquire observation data files and Mogi point source parameter constraint files;

[0030] The inversion module is used to confirm the deformation parameters of the observation points from the observation data file, and obtain the optimal Mogi point source parameters in the region based on nonlinear least squares inversion in combination with the Mogi point source parameter constraint file.

[0031] The forward modeling module is used to confirm the location parameters of the observation points from the observation data file, perform Mogi forward modeling based on the optimal Mogi point source parameters in the region, and obtain the theoretical surface deformation field.

[0032] The evaluation and output module is used to calculate the residual index based on the theoretical surface deformation field and the deformation parameters of the observation points, and output the theoretical surface deformation field and residual index that meet the preset evaluation requirements.

[0033] According to another aspect of this specification, an electronic device is provided, including a memory and a processor, the memory storing program instructions executed by the processor, the processor invoking the program instructions to perform a forward and inverse method for land surface deformation based on the Mogi model.

[0034] According to another aspect of this specification, a non-transitory computer-readable storage medium is provided, the non-transitory computer-readable storage medium storing computer instructions that cause the computer to execute a forward and inverse method for surface deformation based on the Mogi model.

[0035] Compared with the prior art, the beneficial effects of the present invention are as follows: it relies on weighted nonlinear least squares inversion to determine point source parameters, sets the upper and lower bounds and initial values ​​of parameters through the mogi point source parameter constraint file, and provides the evaluation process to verify the rationality of the results. It can intuitively and quickly find the optimal solution or verify the reliability of the solution, overcome data noise and multi-source differences, and provide high-precision inversion solutions. It also provides a separate forward modeling process, filling the gap in the prior art of forward and inverse modeling of surface deformation based on the mogi model. Attached Figure Description

[0036] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the accompanying drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0037] Figure 1 A flowchart illustrating a forward and inverse method for land surface deformation based on the Mogi model, provided in an embodiment of the present invention;

[0038] Figure 2 This is a schematic diagram of the software interface for the forward modeling process in an embodiment of the present invention;

[0039] Figure 3 This is an example diagram of a portion of the code for the forward modeling software in an embodiment of the present invention;

[0040] Figure 4 This is a schematic diagram of the forward modeling results of InSAR, GNSS, and leveling observation data files in an embodiment of the present invention;

[0041] Figure 5 This is a schematic diagram of the software interface for the inversion process in an embodiment of the present invention;

[0042] Figure 6 This is an example diagram of part of the code for the inversion process software in an embodiment of the present invention;

[0043] Figure 7 This is a schematic diagram of the deformation rate results of the original line-of-sight direction of the Cerropreto geothermal field in an embodiment of the present invention;

[0044] Figure 8 This is a schematic diagram of the InSAR observed deformation values, simulated deformation values ​​of the optimal mogi point source inversion parameters, and corresponding fitting residuals after cropping and downsampling of the Cerropreto geothermal field in an embodiment of the present invention.

[0045] Figure 9 This is a schematic diagram of the structure of a surface deformation forward and inverse model based on the Mogi model in an embodiment of the present invention;

[0046] Figure 10 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention. Detailed Implementation

[0047] It should be noted that:

[0048] The terms “comprising” and “having”, and any variations thereof, in the specification, claims, and accompanying drawings of this invention are intended to cover a non-exclusive inclusion, such as a process, method, system, product, or apparatus that includes a series of steps or units, not necessarily limited to those explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0049] The block diagrams shown in the accompanying drawings are merely functional entities and do not necessarily correspond to physically independent entities. That is, these functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor devices and / or microcontroller devices. The flowcharts shown in the accompanying drawings are merely illustrative and do not necessarily include all content and operations / steps, nor do they necessarily have to be performed in the described order. For example, some operations / steps can be decomposed, while others can be combined or partially combined; therefore, the actual execution order may change depending on the specific circumstances.

[0050] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. In addition, the technical features of the various embodiments or individual embodiments provided by the present invention can be arbitrarily combined to form new technical solutions. Such combinations are not bound by the order of steps and / or structural composition patterns, but must be based on the ability of those skilled in the art to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed by the present invention.

[0051] like Figure 1 As shown, a method for forward and inverse modeling of surface deformation based on the Mogi model includes:

[0052] Step 1: Obtain the observation data file and the mogi point source parameter constraint file;

[0053] Step 2: Confirm the deformation parameters of the observation points from the observation data file, and obtain the optimal Mogi point source parameters for the region based on nonlinear least squares inversion using the Mogi point source parameter constraint file.

[0054] Step 3: Confirm the location parameters of the observation points from the observation data file, perform Mogi forward modeling based on the optimal Mogi point source parameters for the region, and obtain the theoretical surface deformation field;

[0055] Step 4: Calculate the residual index based on the theoretical surface deformation field and the deformation parameters of the observation points, and output the theoretical surface deformation field and residual index that meet the preset evaluation requirements.

[0056] Specifically, in step 1, the observation data file includes InSAR observation data file, GNSS observation data file or leveling observation data file, and the file includes, but is not limited to, the location parameters and deformation observation parameters of the observation points;

[0057] The mogi point source parameter constraint file includes lower / upper bounds for longitude, latitude, volume change, and depth.

[0058] Specifically, InSAR observation data files include longitude, latitude (geographic coordinate system coordinates), elevation, line-of-sight deformation value, azimuth, incident angle, and standard deviation.

[0059] Specifically, the GNSS observation data file includes the station name (observation point name), longitude, latitude, elevation, eastward deformation value, northward deformation value, vertical deformation value, eastward standard deviation, northward standard deviation, correlation coefficient, vertical error, and observation time.

[0060] Specifically, the leveling (LEV) observation data file includes the station name, longitude, latitude, elevation, vertical deformation value, standard deviation, and observation time.

[0061] Specifically, in the mogi point source parameter constraint file, the lower / upper bounds of longitude and latitude can be determined based on the longitude and latitude range of the high deformation value area; the lower / upper bounds of the mogi point source volume change can be determined based on the positive / negative value and average value of the deformation, with positive volume change corresponding to deformation uplift and negative volume change corresponding to deformation subsidence, and the volume change is directly proportional to the absolute value of the deformation; the lower / upper bounds of the mogi point source depth can be determined based on the average value of the deformation, and the point source depth is inversely proportional to the absolute value of the deformation.

[0062] Step 2, the inversion process combining the deformation parameter data of each observation point with the Mogi point source parameter constraint file, includes:

[0063] Step 2-1: Read the deformation parameter data of each observation point in the observation data file, construct the observation data matrix, and construct the weight matrix according to the weight of each type of observation data file;

[0064] Step 2-2: Determine the number of Mogi point sources and determine the range of Mogi parameter values ​​based on the Mogi point source parameter constraint file;

[0065] Steps 2-3: Set the reference central meridian and reference ellipsoid, and calculate the coordinates of each observation point and the coordinate range of each Mogi point source in the Gaussian plane rectangular coordinate system using Gaussian projection;

[0066] Steps 2-4: Based on the observation data matrix, weight matrix, and the range of values ​​for the mogi parameters, the optimal mogi point source parameters for the region are obtained using nonlinear least squares inversion.

[0067] Optionally, in steps 2-3, the reference ellipsoid can be selected from WGS-84, Krasovsky ellipsoid, or the 1975 reference ellipsoid.

[0068] Furthermore, in steps 2-4, the process based on nonlinear least squares inversion includes:

[0069] The sampling number for nonlinear least squares inversion is set, and the inversion is performed according to the minimization criterion of weighted least squares joint adjustment with weighted ratio factor to obtain the optimal mogi parameter matrix that satisfies the range of mogi parameter values ​​and the coordinate range of mogi point source.

[0070] The coordinates of the corresponding Mogi point sources in the Gaussian Cartesian coordinate system are converted into longitude and latitude in the matrix elements of the optimal Mogi parameter matrix. Then, all matrix elements are derived to obtain the optimal Mogi point source parameters for the region.

[0071] Furthermore, the mathematical expression for the minimization criterion of weighted least squares joint adjustment with weighted ratio factors is as follows:

[0072]

[0073] In the formula, min represents minimizing the objective function; This is the deformation observation matrix; The design matrix includes the partial derivatives of the mogi point source parameters with respect to the deformation observations; This is the source parameter matrix for the Mogi points.

[0074] The parameters of the Mogi point source include longitude, latitude, volume change of the Mogi point source, and depth of the Mogi point source.

[0075] Specifically, taking a fully contained observation data file (InSAR, GNSS, or leveling) as an example, the nonlinear least squares observation equation based on multi-source data is:

[0076] (1)

[0077] in, This is the observation matrix; The design matrix contains the partial derivatives of the model parameters (Mogi point source parameters) with respect to the deformation observations; The source parameter matrix for the Mogi points; Let be the error matrix. The observation matrix and error matrix are block matrices of the observation data matrix.

[0078] The expanded forms of each matrix are as follows:

[0079] (2)

[0080] in, For the first LOS deformation values ​​(line-of-sight deformation values) of each InSAR observation point. , and The first Vertical deformation values, eastward deformation values, and northward deformation values ​​at each GNSS observation point. For the first Vertical deformation values ​​at each leveling observation point; , and For the first The projected coordinates, depth, and volume changes of each mogi point source; , and These are observation data matrices corresponding to different observation data types. Optionally, in this invention, , and These correspond to InSAR, GNSS, and leveling observation data matrices, respectively. , and For the corresponding design matrix, , and For the corresponding error matrix, equation (1) can be rewritten as:

[0081] (3)

[0082] Introducing a weight matrix :

[0083] (4)

[0084] in, , and These represent the weights of observation data corresponding to different observation data types.

[0085] It is worth noting that the above matrices ( Taking an observation data file that fully includes InSAR, GNSS, or leveling data files as an example, in actual operation, the matrix elements are adjusted according to the actual observation data file type. For example, elements corresponding to observation data file types that do not have a practical basis in a specific environment are deleted.

[0086] Therefore, the mathematical expression for minimizing the weighted least squares joint adjustment following the weighted ratio factor is as follows:

[0087] (5)

[0088] Therefore, we can determine the location of the point source in the mogi point source parameter matrix. , ,depth and volume change Accurate determination of the parameters is crucial for inversion. Under the same deformation data, different initial parameter settings or inversion strategies may yield different solutions, leading to differences in the interpretation of subsurface dynamic processes (such as magmatic activity and fluid migration).

[0089] Preferably, under the same deformation parameter conditions, multiple sets of different initial parameters are set, and the inversion in step 2 and subsequent steps 3 and 4 are performed. Finally, the optimal solution of the inversion result that meets the preset evaluation requirements is selected as the optimal mogi point source parameter for the region.

[0090] Step 3, which involves performing Mogi forward modeling to obtain the theoretical surface deformation field, includes the following steps:

[0091] Step 3-1: Read several observation points and their corresponding location parameters from the observation data file;

[0092] Step 3-2: Confirm the number of Mogi point sources and construct the parameter matrix for each point source based on the optimal Mogi point source parameters for the region.

[0093] Step 3-3: Set the reference central meridian and reference ellipsoid. Based on the position parameters of each observation point and the parameter matrix of each Mogi point source, calculate the projected coordinates of each observation point and Mogi point source in the Gaussian plane rectangular coordinate system through Gaussian projection.

[0094] Steps 3-4: Calculate the radial distance from each observation point to each Mogi point source. Based on the Mogi model, the deformation field of each simulated Mogi point source is calculated separately, and the deformation fields simulated by all Mogi point sources are superimposed to obtain the theoretical surface deformation field.

[0095] Specifically, the application of the Mogi model is based on the assumption that the fluid pressure source is placed in an elastic half-space, and when the source radius is much smaller than the source depth, the fluid flux can be considered as an equivalent spherical source. In the case of a Poisson medium in the Earth's crust and considering the expansion of the Mogi point source volume with a radius of... In the case of an equivalent sphere, in steps 3-4, the deformation values ​​of the simulated Mogi point cloud are calculated based on the classic Mogi model, including the vertical deformation value of the ground surface, the deformation value in the east direction, and the deformation value in the north direction. The mathematical formula is as follows:

[0096] (6)

[0097] (7)

[0098] (8)

[0099] in, , and These are the vertical surface deformation value, the eastward deformation value, and the northward deformation value, respectively. It is the volume change of the Mogi point source; and These are the projected coordinates of the Mogi point source in the Gaussian Cartesian coordinate system; and These are the projected coordinates of the observation point in the Gaussian plane rectangular coordinate system; It is the radial distance from the Earth's surface.

[0100] Preferably, during forward modeling, the InSAR observation data file may also include InSAR observation incident angle and azimuth angle binary raster files, in degrees; and binary data description files, including the number of rows, the number of columns, the longitude of the top left pixel, the meridional step size, the latitude of the top left pixel, the latitudinal step size, and the pixel data format.

[0101] Optionally, according to the method of the present invention, given the known mogi parameters, a separate forward modeling process can be performed based on the observation data file. As a specific embodiment, the present invention implements a separate forward modeling process on a computer with a MATLAB environment that supports the various data file types involved in the method.

[0102] The software interface for the forward modeling process is as follows: Figure 2 As shown, the forward modeling operation steps are as follows:

[0103] Click the “InSAR Deformation” checkbox to start the forward modeling process of the InSAR observation data file. Import the binary raster file of the InSAR observation incident angle and azimuth angle and the binary data description file into the interface. Import the mogi point source parameters. Calculate the coordinates of the InSAR observation point and the mogi point source in the Gaussian plane rectangular coordinate system and the radial distance between each InSAR observation point and the mogi point source in the MATLAB environment. Click “Calculate” to calculate the simulated vertical deformation value, eastward deformation value and northward deformation value according to formulas (6), (7) and (8), and project them to the InSAR line of sight according to the incident angle and azimuth angle information. Click “Export” to export the InSAR observation forward modeling results.

[0104] Clicking the “GNSS Deformation” checkbox will start the forward modeling process for GNSS observation data files, and clicking the “Leveling Deformation” checkbox will start the forward modeling process for leveling observation data files. Following the same steps as the forward modeling process for InSAR observation data files, the GNSS observation forward modeling results and the leveling observation forward modeling results will be exported.

[0105] The software code implementation steps are as follows, with some code examples as shown below. Figure 3 As shown:

[0106] 1. Data Reading and Verification: The system receives various types of imported input data files, including:

[0107] The mogi point source parameters (.mog format) are read and stored in a list structure.

[0108] Leveling observation data files (.lev format) and GNSS observation data files (.gnss format) are read when the corresponding function group is enabled.

[0109] InSAR binary data description file (.rsc format), when the corresponding function group is enabled, its metadata (such as the number of rows and columns, latitude and longitude of the top left corner, cell size, etc.) is parsed into a raster information structure.

[0110] Azimuth raster data files (.img format) are parsed and their pixel values ​​are read into a one-dimensional array when the corresponding function group is enabled and the file is valid.

[0111] When the corresponding function group is enabled and the file is valid, the incident angle raster data file (.img format) is parsed, its metadata is analyzed, and its cell values ​​are read into a one-dimensional array.

[0112] 2. Parameter initialization and configuration, including: ellipsoid parameter settings: determine the semi-major and semi-minor axis values ​​of the reference ellipsoid based on the user interface selection; and central meridian: obtain the central meridian longitude value input from the user interface.

[0113] 3. Construction of the observation point coordinate matrix (Gaussian plane rectangular coordinate system) and the construction of the Mogi point source parameter matrix.

[0114] Preferably, the coordinate matrix is ​​constructed using the geographic coordinates (longitude and latitude) of the observation points (if they exist) in the leveling observation data file. If the leveling observation data file does not exist, the coordinate matrix is ​​constructed using the geographic coordinates (longitude and latitude) of the observation points in the GNSS observation data file (if it exists). If neither of the above exists, but the InSAR binary data description file has been successfully read, the longitude and latitude of the center point of each pixel are calculated based on its metadata (top-left corner coordinates, pixel size, number of rows and columns), and the coordinate matrix is ​​constructed accordingly.

[0115] 4. Coordinate transformation and forward modeling (using the MATLAB calculation engine):

[0116] Data preparation and transfer: Convert the data of the observation point coordinate matrix and mogi point source parameter matrix constructed in the above steps into an array format (mwArray) compatible with the MATLAB environment.

[0117] Calculation call: Calls a MATLAB function named forward_mogi to perform forward calculations based on array-formatted data.

[0118] Calculation results acquisition: After the MATLAB function completes the forward modeling calculation, it returns the deformation field of each observation point in the geographic coordinate system (including three deformation displacement components: eastward deformation value, northward deformation value, and vertical deformation value), which is a one-dimensional displacement result. These results are extracted from the MATLAB array format and stored in local memory.

[0119] 5. Results Reorganization and Output:

[0120] If the calculation is performed on InSAR data: the resulting one-dimensional displacement result array is reorganized into a two-dimensional matrix according to the grid size (number of columns × number of rows) of the InSAR data, representing the deformation values ​​in the east direction, north direction, and vertical direction, respectively.

[0121] If the calculation is for discrete point observation data (leveling or GNSS): directly associate the one-dimensional displacement result array with each discrete observation point, and organize it into an N×1 column vector (N is the number of observation points).

[0122] Finally, the system outputs a message indicating that the forward modeling calculation was successful.

[0123] 6. Export File Settings and Selection: Based on the observation data file type selected in the user interface (InSAR, GNSS, or leveling observation data file), the system will display the corresponding file save dialog box, allowing the user to select the storage path and filename for the exported file. The system supports exporting result files in multiple formats:

[0124] InSAR observation data files can be exported as a standard raster image format (.img) or a specific ROS format (.ros).

[0125] For GNSS observation data files, export them as specific GNSS format files (.gnss).

[0126] For leveling observation data files, export them as a specific leveling format file (.lev).

[0127] 7. File Existence Verification: After selecting the export path, the system checks whether the target file already exists. If it does, a warning is issued to the user and the export operation is aborted to prevent overwriting existing files.

[0128] 8. Associating and formatting the results data with metadata:

[0129] Exporting forward modeling results from leveling observation data files: Assign the vertical deformation value of each observation point obtained from the forward modeling to the corresponding observation point data object. Following a predefined file header format, write all leveling point data, including observation point name, latitude and longitude, elevation, three deformation displacement values ​​(eastward deformation value, northward deformation value, and vertical deformation value), error value, and time information, to the specified .lev file.

[0130] Exporting forward modeling results from GNSS observation data files: Assign the eastward, northward, and vertical deformation values ​​obtained from the forward modeling of each GNSS station to the corresponding GNSS observation point data object. Following a predefined file header format, write all GNSS station data, including observation point name, latitude and longitude, elevation, the three deformation displacement values ​​and their errors, correlation coefficient, and time information, to the specified .gnss file.

[0131] Exporting forward modeling results from InSAR observation data files: The entire three-dimensional deformation field obtained from the forward modeling calculation is converted from a matrix format in memory to a one-dimensional floating-point array arranged in row-major order to conform to the storage format of raster data. Optionally, depending on whether the user provides azimuth and incident angle data, there are two scenarios:

[0132] Scenario 1: If azimuth and incident angle data are not provided, convert the three one-dimensional arrays of the three deformation displacement values ​​into binary byte streams respectively. Use the raster data creation function to generate three independent .img format raster files for the three deformation displacement values.

[0133] Scenario 2: If both azimuth and incident angle data are provided, for each pixel in the deformation field, using its corresponding azimuth and incident angle values, the three deformation displacement values ​​are projected from the 3D ground surface onto the radar's line-of-sight (LOS) direction. The LOS-direction deformation value (line-of-sight direction deformation value) is calculated, and the calculated LOS-direction displacement value for each pixel is stored in a one-dimensional array. Simultaneously, a ROS data list containing the latitude, longitude, elevation, LOS displacement value, azimuth, incident angle, and error information for each pixel is constructed. (Based on the user-selected file extension:)

[0134] If the .ros format is selected, the above ROS data list will be written to a file in a specific text format.

[0135] If the .img format is selected, the LOS to deformed one-dimensional array is converted into a binary byte stream, and a raster image file containing the LOS deformation field is created.

[0136] 9. Operation Result Feedback

[0137] After any of the above export operations are completed, the system will display a corresponding success message or error warning to the user, depending on whether the operation was successful or not.

[0138] In the aforementioned software, taking a set of InSAR, GNSS, and leveling observations with a longitude range of 86.725°~87.075° and a latitude range of 32.900°~33.250° as an example, the given InSAR observation points are arranged in 35 rows * 35 columns, with a fixed azimuth angle of -10° and a fixed incident angle of 23°. The given GNSS observation points include a total of 42 stations, and the given leveling observation points include a total of 15 stations (observation points). The given Mogi point source parameters are 86.820°, 33.130°, and 4,500,000m. 3 6000m (Mogi point source 1) and 86.980°, 33.000°, 6000000m 3 6000m (mogi point source 2).

[0139] Use plotting tools to plot the forward modeling results of InSAR, GNSS, and leveling observations together. Figure 4 . Figure 4 The results also include forward modeling data from InSAR, GNSS, and leveling observations. InSAR data is represented by colored rectangular blocks, GNSS data by black vector arrows with blue error ellipses, and leveling data by pink arrows with blue error bands. It can be seen that the surface deformation fields generated by the two established Mogi point sources are clearly identifiable, and the forward modeling results are intuitive and straightforward.

[0140] As a specific embodiment, the present invention implements a method for forward and inverse modeling of land surface deformation based on the Mogi model in the computer environment mentioned in the above embodiments.

[0141] The software interface for the inversion process is as follows: Figure 5 As shown, the code implementation steps are as follows, with some code as shown below. Figure 6 As shown:

[0142] 1. Multi-source observation data reading and analysis:

[0143] The system reads various types of input data files, including:

[0144] The mogi model constraint parameter file (.moc format) is parsed and stored as a parameter list.

[0145] The InSAR observation data file [specifically, the LOS (line of sight) deformation data file (.los or .ros format)] contains the latitude and longitude, LOS deformation value, azimuth, angle of incidence, and error information for each observation point.

[0146] The leveling observation data file (.lev format) contains the latitude, longitude, elevation, vertical deformation value, and error information for each observation point.

[0147] The GNSS observation data file (.gnss format) contains information such as latitude and longitude, elevation, eastward, northward, and vertical deformation components, errors, and correlation coefficients for each observation point.

[0148] 2. Inversion parameter initialization and configuration:

[0149] Ellipsoid parameter settings: Select the values ​​of the semi-major and semi-minor axes of the reference ellipsoid according to the user interface selection.

[0150] Central Meridian: Obtains the longitude value of the central meridian input from the user interface.

[0151] Sampling parameter settings: Obtain the number of samples set by the user.

[0152] Weight matrix construction: Based on the weight values ​​set by the user for InSAR, GNSS and leveling data, a 1x3 weight matrix is ​​constructed to balance the contributions of different data sources in joint inversion.

[0153] Construction of the observation data matrix:

[0154] Leveling observation data matrix: Construct an Nx7 matrix from the latitude, longitude, observed values, errors, and other information of each leveling observation point.

[0155] GNSS observation data matrix: This matrix is ​​constructed by combining the latitude and longitude of each GNSS observation point, deformation components in three directions, errors, correlation coefficients, and other information into an Nx12 matrix. .

[0156] InSAR observation data matrix: For each InSAR observation point, the projection coefficients of its LOS direction with respect to east, north, and vertical deformation are calculated based on its azimuth and incident angle. The latitude and longitude, LOS deformation value, error, and three projection coefficients of each point are then constructed into an Nx7 matrix. .

[0157] The initial parameter matrix for Mogi point sources is constructed from the longitude, latitude, depth, and volume variation parameters of each point source read from the constraint file, forming an Nx4 matrix. .

[0158] in, , , and The expanded form is as follows: These represent the east-west, north-south, and vertical projection coefficients of the InSAR observation point, respectively. The meanings of the other variables are explained in the section above on nonlinear least squares inversion.

[0159]

[0160]

[0161] 3. Call the inversion algorithm engine to perform calculations, including:

[0162] Data preparation and transmission: The constructed observation data matrix, initial parameter matrix of the mogi point source, weight matrix, central meridian, number of samples, and other parameters are converted into an array format (mwArray) compatible with the MATLAB environment.

[0163] Perform the inversion calculation: Call the inversion algorithm function named Inversion_mogi. This function takes all the above parameters as input and performs a nonlinear inversion calculation based on the mogi model according to Equation 5.

[0164] Calculation results acquisition: After the inversion algorithm is completed, the result returned is the optimized mogi point source parameters: including the optimized longitude, latitude, depth, and volume change parameters of each inverted point source.

[0165] Model predictions and residuals for various data types: Returns the residuals between model predictions and observed values ​​for InSAR, GNSS, and leveling data.

[0166] Goodness-of-fit index: Returns a value (root mean square error) that represents the overall goodness of fit.

[0167] 4. Inversion Results Output and Residual Analysis

[0168] Optimize mogi parameter output: Write the optimized mogi point source parameters (longitude, latitude, depth, volume change) to the specified output file (.mog format).

[0169] GNSS Results Output: The model prediction values ​​(East, North, Vertical) and residual values ​​of the GNSS deformation components are assigned to the corresponding GNSS data objects. The results, including the model prediction values ​​and residuals, are output to two separate GNSS format files (.gnss).

[0170] Leveling results output: The model prediction values ​​and residual values ​​of leveling deformation are assigned to the corresponding leveling data objects.

[0171] The results, including model predictions and residuals, are output to two separate level format files (.lev).

[0172] InSAR output: The model predictions and residuals of the InSAR LOS deformation are assigned to the corresponding InSAR data objects. The results, including the model predictions and residuals, are output to two separate LOS format files (.los).

[0173] Root Mean Square Error (RMSE) Calculation: Based on the residuals of all observation data (three components of InSAR points, leveling points, and GNSS stations), a comprehensive root mean square error is calculated and displayed in the user interface as an overall evaluation of the inversion fitting effect.

[0174] 5. Feedback upon completion of calculation:

[0175] After all output operations are completed, the system will display a message to the user indicating that the calculation is complete.

[0176] This embodiment uses InSAR observations of the Cerro Preto geothermal field at the southern end of the San Andreas Fault System as an example. This geothermal field is the world's largest, primarily hydrothermal, geothermal field. The acquired satellite imagery comes from the Sentinel-1A T173 descending orbit, with observations from October 29, 2014 to May 4, 2025. The longitude range is -115.611° to -114.678°, the latitude range is 32.258° to 32.583°, the azimuth range is 190.469° to 191.019°, the incident angle range is 30.558° to 36.692°, the number of rows is 1749, and the number of columns is 609. The original line-of-sight deformation rate results are as follows: Figure 7 As shown, the deformation rate ranges from -14.769 to 59.736 mm / yr (negative values ​​represent deformation closer to the satellite, and positive values ​​represent deformation farther from the satellite). The given Mogi point source parameters are constrained to -115.270 to -115.080°, 32.340 to 32.490°, -100000000 to -10000 m³, and 100 to 10000 m (point source 1 and point source 2).

[0177] First, import the cropped and downsampled InSAR observation data file into the interface, and import the mogi point source parameter constraint file; set the reference ellipsoid "Reference Ellipsoid" to "WGS-84", set the reference central meridian "L0" to "86.725°", and calculate the coordinates of the InSAR observation points and the upper and lower bounds of the mogi point source in the Gaussian plane rectangular coordinate system. and , The radial distances on the ground between each InSAR observation point and the upper and lower bounds of the Mogi point source coordinates. ;

[0178] Set the inversion parameter "Sample Size" to "100";

[0179] The data weight scaling factor can be omitted because there is no GNSS or leveling observation data in this example. The weight matrix... It is a diagonal matrix of [1,0,0].

[0180] The observation equation is constructed according to formula (3), and the optimal solution of the mogi point source parameters is obtained by minimizing the formula (5). Table 1 gives the optimal mogi point source inversion parameters for this example.

[0181] Table 1 Optimal Mogi point source inversion parameters

[0182]

[0183] Figure 8 The InSAR observed deformation values ​​(a), simulated deformation values ​​(b), and corresponding fitting residuals (c) of the Cerro Preto geothermal field after trimming and downsampling, generated using a surface deformation forward and inverse retrieval method based on the Mogi model provided by this invention, are shown. It can be seen that the location of the optimal Mogi point source precisely corresponds to the high-value center region of the surface deformation field (Table 1). Figure 8 a) fully conforms to the Mogi model theory. Further model validation shows that the fitting residuals are uniformly distributed and no significant systematic bias was detected. Figure 8 c) The inversion results are accurate. This strongly demonstrates the effectiveness and reliability of this method (especially its multi-source data joint inversion framework based on weighted nonlinear least squares) in accurately inverting the location and parameters of underground point sources. Its core advantage lies in its ability to overcome data noise and multi-source differences, providing high-precision inversion solutions.

[0184] After being integrated into a highly efficient and intuitive tool platform, the method proposed in this invention allows users to intuitively set or adjust the position of the Mogi point source for forward modeling. , ,depth and volume change The value of is calculated in real time, and the corresponding theoretical deformation field is calculated ( , and The system instantly displays the results on a graphical interface. For inversion, it can instantly display the fitting relationship (residual plot) between the theoretical deformation field and the observed data. Users can easily modify the upper and lower bounds of parameters, set weights (for multi-source data), and dynamically view changes in the fitting effect. This greatly optimizes the inversion process, enabling intuitive exploration and efficient determination of Mogi point source parameters, and solving the shortcomings of existing technologies in terms of interactivity, visual feedback, and efficiency.

[0185] The implementation of the various embodiments of this invention is based on programmed processing through a system with processor functionality. Therefore, in practical engineering, the technical solutions and functions of the various embodiments of this invention are encapsulated into various modules. Based on this reality, and building upon the above embodiments, the embodiments of this invention provide a surface deformation forward and inverse modeling system based on the Mogi model. This system is used to execute a surface deformation forward and inverse modeling method based on the Mogi model from the above method embodiments.

[0186] See Figure 9 The system includes:

[0187] The data import module is used to acquire observation data files and Mogi point source parameter constraint files; the inversion module is used to confirm the deformation parameters of the observation points from the observation data files, and obtain the optimal Mogi point source parameters for the region based on nonlinear least squares inversion using the Mogi point source parameter constraint files; the forward modeling module is used to confirm the location parameters of the observation points from the observation data files, and perform Mogi forward modeling based on the optimal Mogi point source parameters for the region to obtain the theoretical surface deformation field.

[0188] The evaluation and output module is used to calculate the residual index based on the theoretical surface deformation field and the deformation parameters of the observation points, and output the theoretical surface deformation field and residual index that meet the preset evaluation requirements.

[0189] It should be noted that the system embodiments provided by the present invention are used not only to implement the methods in the above method embodiments, but also to implement the methods in other method embodiments provided by the present invention. The only difference is that corresponding functional modules are set. The principle is basically the same as that of the above system embodiments provided by the present invention. As long as those skilled in the art can improve the modules in the above system embodiments by referring to the specific technical solutions in other method embodiments and combining technical features to obtain corresponding technical means and technical solutions composed of these technical means, on the basis of the above system embodiments, and on the premise of ensuring the practicality of the technical solutions, they can obtain corresponding system-like embodiments for implementing the methods in other method-like embodiments.

[0190] Preferably, the data import module is also used to obtain the mogi point source parameter file, providing a data foundation for a separate forward modeling process.

[0191] The method in this embodiment of the invention is implemented using an electronic device; therefore, it is necessary to introduce the relevant electronic device. For this purpose, embodiments of the present invention provide an electronic device, such as... Figure 10As shown, the electronic device includes: at least one processor, a communication interface, at least one memory, and a communication bus, wherein the at least one processor, the communication interface, and the at least one memory communicate with each other via the communication bus. The at least one processor invokes logical instructions stored in the at least one memory to execute all or part of the steps of the methods provided in the foregoing method embodiments.

[0192] Furthermore, when the logical instructions in at least one of the aforementioned memories are implemented as software functional units and sold or used as independent products, they are stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, or a part of the technical solution, is embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (a personal computer, server, or network device) to execute all or part of the steps of the methods described in the various method embodiments of the present invention. The aforementioned storage medium includes: USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks—various media for storing program code.

[0193] The system embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units, located in one place, or distributed across multiple network units. The purpose of this embodiment is achieved by selecting some or all of the modules according to actual needs. Those skilled in the art will understand and implement this without any inventive effort.

[0194] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

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

[0196] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0197] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0198] Based on the same technical concept as the foregoing embodiments, the present invention provides a non-transitory computer-readable storage medium that stores computer instructions that cause the computer to execute a forward and inverse method for surface deformation based on the Mogi model.

[0199] In summary, this invention relates to a method and system for forward and inverse modeling of surface deformation based on the Mogi model. At a given InSAR, GNSS, or leveling observation point, the surface deformation field is forward modeled according to given Mogi point source parameters. Based on given InSAR, GNSS, or leveling deformation observations and constrained by the Mogi point source parameters, the Mogi point source parameters for the study area are inverted, and the corresponding simulated deformation and fitting residuals are calculated. This method and system fill the gap in existing technologies for forward and inverse modeling of surface deformation based on the Mogi model, providing an efficient and intuitive tool platform integrating data visualization, forward simulation, inverse analysis, and result output, significantly improving the model's constraint accuracy and application efficiency.

[0200] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the technical solutions of the embodiments of the present invention.

Claims

1. A surface deformation forward and inverse method based on a mogi model, characterized in that, The method comprises the following steps: Obtaining an observation data file and a mogi point source parameter constraint file; Confirming deformation parameters of observation points from the observation data file, and combining the mogi point source parameter constraint file to obtain regional optimal mogi point source parameters based on nonlinear least squares inversion; Confirming position parameters of the observation points from the observation data file, and obtaining a theoretical ground surface deformation field based on mogi forward calculation of the regional optimal mogi point source parameters; Calculating a residual index according to the theoretical ground surface deformation field and the deformation parameters of the observation points, and outputting the theoretical ground surface deformation field and the residual index that meet preset evaluation requirements.

2. The Mogi model-based surface deformation forward and inversion method according to claim 1, characterized in that, The process of combining the deformation parameter data of each observation point with the mogi point source parameter constraint file for inversion comprises the following steps: Reading the deformation parameter data of each observation point in the observation data file, constructing an observation data matrix, and constructing a weight matrix according to the weights of each type of observation data file; Determining the number of mogi point sources, and determining the mogi parameter value range according to the mogi point source parameter constraint file; Setting a reference central meridian and a reference ellipsoid, and calculating the coordinates of each observation point and the coordinate range of each mogi point source in the Gauss plane rectangular coordinate system through Gauss projection; Based on the observation data matrix, the weight matrix and the mogi parameter value range, regional optimal mogi point source parameters are obtained based on nonlinear least squares inversion.

3. The Mogi model-based surface deformation forward and inversion method of claim 2, wherein, The process of nonlinear least squares inversion comprises the following steps: Setting the sampling number of nonlinear least squares inversion, and performing inversion according to the minimization criterion of weighted least squares combined adjustment with weight ratio factors to obtain an optimal mogi parameter matrix that meets the mogi parameter value range and the coordinate range of the mogi point source; Converting the coordinates of the mogi point source in the Gauss plane rectangular coordinate system in the matrix elements of the optimal mogi parameter matrix into longitude and latitude, and then deriving all the matrix elements to obtain regional optimal mogi point source parameters.

4. The Mogi model-based surface deformation forward and inversion method of claim 3, wherein, The mathematical expression of the minimization criterion of the weighted least squares combined adjustment with weight ratio factors is as follows: ； In the formula, min represents the minimization objective function. W is a weight matrix; is a deformation observation matrix; is a design matrix, containing the partial derivatives of the mogi point source parameters with respect to the deformation observations; is a mogi point source parameter matrix.

5. The Mogi model-based surface deformation forward and inversion method of claim 1, wherein, The mogi point source parameters include longitude, latitude, mogi point source volume change and mogi point source depth.

6. The Mogi model-based surface deformation forward and inversion method of claim 1, wherein, The process of mogi forward calculation to obtain the theoretical ground surface deformation field comprises the following steps: Reading a plurality of observation points and their corresponding position parameters from the observation data file; Based on the regional optimal mogi point source parameters, constructing a parameter matrix of each mogi point source; Setting a reference central meridian and a reference ellipsoid, and calculating the projection coordinates of each observation point and mogi point source in the Gauss plane rectangular coordinate system based on the position parameters of each observation point and the parameter matrix of each mogi point source; Calculating the ground surface radial distance between each observation point and each mogi point source, calculating the deformation field of each mogi point source based on the mogi model, superimposing the deformation fields of all mogi point sources to obtain the theoretical ground surface deformation field.

7. The Mogi model-based surface deformation forward and inversion method according to claim 6, characterized in that, The method further comprises a separate forward calculation process: The position parameters of the observation points are confirmed from the observation data file, actual mogi point source parameters are obtained to replace the regional optimal mogi point source parameters obtained through inversion, mogi forward calculation is performed, and a theoretical ground surface deformation field is obtained.

8. A surface deformation forward and inversion system based on mogi model, characterized in that, The method comprises the following steps: a data import module is configured to obtain an observation data file and a mogi point source parameter constraint file; an inversion module is configured to confirm deformation parameters of observation points from the observation data file, and obtain regional optimal mogi point source parameters through inversion based on a non-linear least square method in combination with the mogi point source parameter constraint file; a forward calculation module is configured to confirm position parameters of the observation points from the observation data file, and perform mogi forward calculation based on the regional optimal mogi point source parameters to obtain a theoretical ground surface deformation field; an evaluation and output module is configured to calculate a residual index according to the theoretical ground surface deformation field and the deformation parameters of the observation points, and output the theoretical ground surface deformation field and the residual index that meet preset evaluation requirements.

9. An electronic device, comprising: The computer readable storage medium stores computer instructions, and the computer instructions enable the computer to execute the method in any one of claims 1 to 7.

10. A non-transitory computer-readable storage medium, comprising: The non-transitory computer readable storage medium stores computer instructions, and the computer instructions enable the computer to execute the method in any one of claims 1 to 7.

Images