CN106969751A - A kind of method of the coal mining subsidence amount monitoring calculation based on unmanned aerial vehicle remote sensing - Google Patents

Abstract

The present invention relates to a kind of method of the coal mining subsidence amount monitoring calculation based on SUAV remote sensing technology, comprise the following steps：One, according to research area's topography and geomorphology, scope, position and ground resolution planning unmanned plane during flying course line, two, unmanned plane is manipulated by setting airline operation, and gather remote sensing images, three, a number of control point is laid according to research area's scope and geographical feature, and these points are measured one by one, four, point cloud generation processing is realized using professional remote sensing image processing software, and export research area's earth's surface three-dimensional coordinate, five, measurement result is controlled according to ground, the three-dimensional coordinate that aerial survey of unmanned aerial vehicle is obtained is corrected using 4 parametric surface fitting calibration models, six, it is compared using the terrain elevation values obtained after correction with topographic map height value before exploitation, calculate subsidence amount after exploitation.The method and step of the present invention is simple, reasonable in design and realizes that convenient, computational accuracy is high, and subsidence amount caused by coal mining effectively can be calculated.

Description

A kind of method of the coal mining subsidence amount monitoring calculation based on unmanned aerial vehicle remote sensing

Technical field

The invention belongs to mineral engineering coal mining subsidence monitoring technical field, put down more particularly, to a kind of vegetation is rare The monitoring calculation method of former area's coal mining subsidence amount.

Background technology

Coal mining subsidence is a kind of surface subsidence phenomenon caused by mining, big and the duration is long with influence area The characteristics of, its appearance is frequently accompanied by the destruction of surface movement and deformation, causes earth's surface collapse-fissure, Collapse Pit occur, collapse The geological disaster such as groove and association landslip, can cause different degrees of to facilities such as earth's surface traffic, building, communication and water conservancies Harm, seriously threatens mining area ecological environment and people life property safety.It can be ore deposit to carry out the monitoring of mining area coal mining land subsidence Area's land resource and building protection, damage appraisement and management provide foundation, and the also preventing and treating to coal mining surface collapse disaster has Significance.

At present, being used coal mining subsidence more set up surface strata movement observation station, utilize total powerstation etc. to carry out manual measurement Method is observed, but this conventional observation method inefficiency, cost are high and monitoring range is limited, are typically only used for surface strata movement Law study, and cannot be used for coal mining subsidence investigation and monitoring among.In recent years, land resources department has carried out profit Coal-mining subsidence investigation and monitoring are carried out with satellite remote sensing technology and aerial remote sens ing technique, but is existed because image resolution is low And the problem of None- identified smaller width surface cracks, there is also can not pass through image and carry out coal mining subsidence amount to calculate etc. Problem.The appearance and development of unmanned aerial vehicle remote sensing technology, a kind of new approaches are provided without the solution suspected of above mentioned problem.Open at present Application study work of the unmanned aerial vehicle remote sensing technology of exhibition in coal-mining subsidence monitoring shows that unmanned aerial vehicle remote sensing images can quickly have Effect ground identification subsidence crack, is a kind of good method of coal-mining subsidence monitoring.The parameter values such as coal mining subsidence amount are to evaluate The important indicator of Mining subsidence, but how to carry out the surface movement parameters such as coal mining subsidence amount using unmanned aerial vehicle remote sensing technology Calculate, there is not yet research report.

The content of the invention

The problem of it is an object of the invention to exist to prior art, is solved there is provided a kind of based on unmanned aerial vehicle remote sensing The method of coal mining subsidence amount monitoring calculation, this method step is simple, cost is low, speed is fast, result of calculation is reliable, can be rapid And the task of coal mining land subsidence monitoring is accurately finished, and then the wider region of no relief coal mining subsidence rare to vegetation Effectively monitored.

To realize that the technical scheme of above goal of the invention is as described below.

A kind of method of the coal mining subsidence amount monitoring calculation based on unmanned aerial vehicle remote sensing, it is characterised in that utilize small-sized nothing Man-machine remote sensing technology is calculated coal mining subsidence amount, and process is as follows：

Step 1: unmanned plane flight-line design：According to research area's topography and geomorphology, scope, position and ground resolution planning nothing The man-machine line of flight；

Unmanned plane course line mainly includes flying height, Duplication and course, and process is as follows：

Step 101, flying height are determined：According to formula

The unit that calculating obtains flying height H, H under definitely face resolution ratio correspondence is rice；In formula (1), GSD is ground Face resolution ratio, unit is rice；F is camera lens focal length, and unit is millimeter；A is camera pixel dimension, and unit is micron；

Step 102, Duplication are determined：Determined according to research area's topography and landform character, general setting region of no relief course and side It is respectively 60% and 40% to Duplication, mountain area course and sidelapping rate are respectively 80% and 60%；

Step 103, course are determined：Determined according to research area's shape along research area master is axial arranged；

Step 2: unmanned plane：The course line of design is imported into unmanned plane, between manipulation unmanned plane is by airline operation is set simultaneously Every collection research area's remote sensing image of taking photo by plane；Each sortie exports the remote sensing images collected and coordinate POS data after the completion of taking photo by plane, And carry out the inspection of degree of overlapping and the quality of image；

Step 3: remote sensing image processing：The remote sensing image and coordinate data that unmanned plane is obtained import Pix4D The professional remote sensing image processing softwares of mapper 2.0, successively complete image mosaic, the encryption of point cloud, DOM and DSM generations processing, and By each mesh point P ' in 1000cm grid density size export research area₁、P’₂、P’₃、…、P’_nPlane coordinates X, Y and big The high Z in ground_Greatly；

Step 4: ground control survey：A number of control point is laid according to research area's scope and geographical feature, and to this A little points are measured one by one, and specific practice is as follows：

Step 401, control measurement points are determined and laid：The each point P ' that aerial survey of unmanned aerial vehicle is obtained₁、P’₂、P’₃、…、 P’_nSpread selects m ground survey point P according to research area's area and topography and geomorphology in research Qu Tushang from these points₁、P₂、 P₃、…、P_m(m >=5), measurement point is evenly distributed in research area, by taking a fully-mechanized mining working as an example, the point of ground control survey Four corners and center of working face are distributed in, and quantity is no less than 5；

Step 402, control point survey：Respectively to ground survey point P₁、P₂、P₃、…、P_mSetting-out is carried out, and by the elevation fourth class The Control Surveying in High Accuracy requirement of the measurement of the level obtains each measurement point P to these ground survey point testings one by one₁、P₂、P₃、…、 P_mNormal high Z_Just；

Step 5: GPS height fittings are corrected：Measurement result is controlled according to ground, it is high using the GPS of 4 parametric surface fittings Journey fitting correction model is corrected to the three-dimensional coordinate that aerial survey of unmanned aerial vehicle is obtained, and process is as follows：

Step 501,4 parameters are calculated：According to formula

H_Just=H_Greatly-ζ_i(i=1,2,3 ..., m) (2)

Calculate each common point P₁、P₂、P₃、...、P_mCorresponding abnormal elevation value ζ_i, in formula (2), H_JustFor normal high level, H_GreatlyFor the earth high level, unit is rice；ζ_iFor each common point P₁、P₂、P₃、…、P_mCorresponding abnormal elevation value；And according to formula

ζ_i=a₀+a₁x_i+a₂y_i+a₃x_iy_i(i=1,2,3 ..., m) (3)

Calculating obtains a₀、a₁、a₂And a₃In the value of this 4 parameters, formula (3), a₀、a₁、a₂And a₃For unknown parameter, ζ_iFor Each common point P₁、P₂、P₃、...、P_mAbnormal elevation value, unit is rice, x_iAnd y_iFor each common point P₁、P₂、P₃、...、P_mIt is corresponding Plane coordinates, unit is rice；

Because the number of common point is more than 4, corresponding error equation is listed

v_i=a₀+a₁x_i+a₂y_i+a₃x_iy_i-ζ_i(i=1,2,3 ..., m) (4),

In formula (4), v_iFor error amount；Being write as matrix form is

V=XA- ζ (5),

In formula (5),

It can be tried to achieve according to principle of least square method

Step 502, accuracy assessment：The correction of GPS height fittings includes precision of inner coincidence and precision of exterior coincidence, and process is as follows：

Step 5021, precision of inner coincidence are determined：According to formula

V_i=ζ_i-ζ_i' (7)

Calculate the residual error V for obtaining model of fit_i, in formula (7), V_iFor regression criterion, ζ_iThe Fitting Calculation has been participated in for n Know actual value a little, ζ_i' it is the n match value for participating in the Fitting Calculation known point；And according to formula

The precision of inner coincidence M of digital simulation model_{It is interior}, in formula (8), M_{It is interior}For precision of inner coincidence, V_iFor regression criterion, n' is Participate in the known point quantity of the Fitting Calculation；

Step 5022, precision of exterior coincidence are determined：According to formula

The precision of exterior coincidence M of digital simulation model_Outside；In formula (9), M_OutsideFor precision of exterior coincidence, V_iFor regression criterion, m' is Check point quantity；

Step 503, elevation correction：Each parameter a that formula (6) is obtained₀、a₁、a₂And a₃Institute can be tried to achieve by substituting into formula (3) There is mesh point P '₁、P’₂、P’₃、...、P’_nHeight anomaly ζ_i(i=1,2,3 ..., n), by ζ_iIt is that can obtain to substitute into formula (2) All mesh point P '₁、P’₂、P’₃、...、P’_nNormal high H_Correction；

Step 6: sinkage is calculated：Carried out using the terrain elevation values obtained after correction with topographic map height value before exploitation Compare, calculate subsidence amount after exploitation；

The planar point P " obtained after GPS height fittings are corrected₁、P”₂、P”₃、...、P”_nCoordinate X, Y and normal high Z_Correction Machine plotting to research area's pre-existing topography figure on, and by 100 meters × 100 meters grid to the stacking chart carry out mesh generation, to score Analyse before each new grid points are adopted and adopt rear height value, according to formula

Δ H=H₀-H_Correction(10),

Can calculate before research paneling with subsidence value Δ H after exploitation, finally draw it is corresponding sink etc. It is worth line chart；In formula (10), Δ H is exploits subsidence value after preceding and exploitation, and unit is rice；H₀For former topographic map before exploitation On height value, unit is rice；H_CorrectionFor each planar point P " after being corrected through GPS height fittings₁、P”₂、P”₃、...、P”_nHeight Journey value, unit is rice.

Compared with prior art, calculating coal-mining subsidence measurer using the method for the invention has that method and step is simple, input Cost is low, monitoring cycle is short, monitoring accuracy high (coal-mining subsidence amount computational accuracy up to centimetre rank), monitoring range are big (can one Secondary property realizes the monitoring of the large area sinking land in coalmining areas), safe, good practical effect the advantages of, can be to earth's surface caused by mining Sinkage is effectively calculated, and is with a wide range of applications in measurement industry, mining industry, geology industry.Unmanned aerial vehicle remote sensing should be based on Coal mining subsidence amount monitoring calculation method not only reduce human and material resources and financial resources economic resources, shorten monitoring week Phase, improve operating efficiency, and can well be combined with social science and technology's progress, realized distant to coal mining surface collapse region Feel the quick obtaining of image, and combine field survey data, quickly and accurately complete the task of coal mining land subsidence monitoring, be The monitoring of wider coal mining land subsidence provides facility.

Brief description of the drawings

Fig. 1 is the FB(flow block) of the method for the invention.

Fig. 2 controls ground the schematic diagram that measurement point is laid for the present invention.

Fig. 3 is the schematic diagram calculated using the method for the invention sinkage.

Embodiment

The method of a kind of coal mining subsidence amount monitoring calculation based on unmanned aerial vehicle remote sensing as shown in Figure 1, with Ning Dongmei Exemplified by the fully-mechanized face of charcoal base ore deposit 805, the subsidence amount of the sinking land in coalmining areas is counted using SUAV Calculate, process is as described below.

Step 1: unmanned plane flight-line design：According to research area's topography and geomorphology, scope, position and ground resolution planning nothing The man-machine line of flight.

The embodiment of the present invention is taken photo by plane special small unmanned plane using the rotors of MD4-1000 tetra-, and sensor is public using Sony The SONY A7R cameras of production are taken charge of, the SUAV line of flight, which is related to, mainly includes flying height, Duplication and course design, tool Body is as follows：

Step 101, flying height are determined：Ground resolution GSD size is determined by coal mining surface subsidence fracture width, due to ground Table fracture width majority is more than 2 centimetres, therefore using ground resolution GSD as 2 centimetres of determination flying height, and SONY A7R cameras are inherently joined Number lens focus f is 35.0 millimeters, and pixel dimension a is 4.88 microns.According to formula

It is about 143 meters to calculate the flying height H obtained under definitely face resolution ratio correspondence.In formula (1), GSD is ground Resolution ratio, unit is rice；F is camera lens focal length, and unit is millimeter；A is camera pixel dimension, and unit is micron.

Step 102, Duplication are determined：Because the working face location of ore deposit 805 is the micro- mound topography and geomorphology in Plain, therefore unmanned plane Course and sidelapping rate are respectively 60% and 40%.

Step 103, course are determined：Because the ore deposit 805 work is rectangle, therefore the main shipping track of unmanned plane is moved towards along the working face Arrangement.

Step 2: unmanned plane：Unmanned plane is manipulated by setting airline operation, and gathers remote sensing images.

The course line of design is imported into MD4-1000 SUAVs, nothing is manipulated in noon sunny and breezy It is man-machine by setting airline operation, and be spaced take photo by plane collection the working face earth's surface remote sensing image.Each sortie is led after the completion of taking photo by plane Go out the remote sensing images collected and coordinate POS data, and carry out the inspection of degree of overlapping and the quality of image.

Step 3: remote sensing image processing：The generation processing of point cloud is realized using professional remote sensing image processing software, and export is ground Study carefully area's earth's surface three-dimensional coordinate.

All remote sensing images and coordinate POS data that unmanned plane is obtained import professional remote sensing image processing software Pix4D mapper 2.1, successively complete image mosaic, the encryption of point cloud, DOM and DSM generations, and by 1000cm grid density Size exports each mesh point P ' in the working face₁、P’₂、P’₃、...、P’_nPlane coordinates X, Y and geodetic height Z_Greatly(table 1).

The aerial survey of unmanned aerial vehicle height value of table 1

Step 4: ground control survey：A number of control point is laid according to research area's scope and geographical feature, and to this A little points are measured one by one.

Step 401, control measurement points are determined and laid：The each point P ' that aerial survey of unmanned aerial vehicle is obtained₁、P’₂、P’₃、...、 P’_nSpread is on the coal-face figure, and because the coal-face is rectangle and area is smaller, topography and geomorphology is simple, therefore choosing Select 5 ground survey point P₁、P₂、P₃、P₄And P₅(table 2), the point of ground control survey is distributed in four corners of working face with Heart position (as shown in Figure 2).

The ground control survey point coordinates of table 2

Step 402, control point survey：Respectively to ground survey point P₁、P₂、P₃、P₄And P₅Setting-out is carried out, and by the elevation fourth class The Control Surveying in High Accuracy requirement of the measurement of the level obtains each control measurement point P to these ground survey point testings one by one₁、P₂、P₃、 P₄And P₅Normal high Z_Just(table 3).

The ground control survey point coordinates of table 3

Step 5: GPS height fittings are corrected：Measurement result is controlled according to ground, using the parameter fitting calibration model of plane 4 The three-dimensional coordinate that aerial survey of unmanned aerial vehicle is obtained is corrected.

According to control measurement result, geodetic height is changed using the parameter fitting model of plane 4, the following of process

Step 501,4 parameters are calculated：According to formula

H_Just=H_Greatly-ζ_i(i=1,2,3,4,5) (2)

Each common point P can be calculated₁、P₂、P₃、P₄And P₅Corresponding abnormal elevation value ζ_iIn (table 4), formula (2), H_JustFor just Normal high level, H_GreatlyFor the earth high level, unit is rice；ζ_iFor each common point P₁、P₂、P₃、P₄And P₅Corresponding abnormal elevation value.

The common point abnormal elevation value of table 4

And according to formula

ζ_i=a₀+a₁x_i+a₂y_i+a₃x_iy_i(i=1,2,3,4,5) (3),

It can calculate and obtain a₀、a₁、a₂And a₃In the value of this 4 parameters, formula (3), a₀、a₁、a₂And a₃For unknown parameter, ζ_i For each common point P₁、P₂、P₃、...、P_mAbnormal elevation value, unit is rice, x_iAnd y_iFor each common point P₁、P₂、P₃、...、P_mCorrespondence Plane coordinates, unit is rice.Because the number of common point is more than 4, then corresponding error equation can be listed

v_i=a₀+a₁x_i+a₂y_i+a₃x_iy_i-ζ_i(i=1,2,3,4,5) (4),

In formula (4), v_iFor error amount.Being write as matrix form is

V=XA- ζ (5),

In formula (5),According to minimum Square law principle can be tried to achieve

Fit parameter values are shown in Table 5.

The GPS elevation plane fitting parameter computational charts of table 5

Step 502, accuracy assessment：The correction of GPS height fittings includes precision of inner coincidence and precision of exterior coincidence, and detailed process is such as Under：

Step 5021, precision of inner coincidence are determined：According to formula

Vi=Hi-H_i' (i=1,2,3,4,5) (7)

Calculate the residual error V for obtaining model of fit_iIn (table 6), formula (7), V_iFor regression criterion, unit is rice；Hi is n The actual value of the Fitting Calculation known point is participated in, unit is rice；H_i' it is the n match value for participating in the Fitting Calculation known point, unit is Rice.And according to formula

The precision of inner coincidence M of digital simulation model_{It is interior}For 0.2575 meter, in formula (8), M_{It is interior}For precision of inner coincidence, unit is Rice；V_iFor regression criterion, unit is rice；N' is the known point quantity for participating in the Fitting Calculation.

The residual error table of the GPS Height Fitting Models of table 6

Step 5022, precision of exterior coincidence are determined：This example is temporarily calculated without precision of exterior coincidence.

Step 503, elevation correction：Each parameter a that formula (6) is obtained₀、a₁、a₂And a₃Value substitute into formula (3) can ask Obtain all mesh point P '₁、P’₂、P’₃、...、P’_nHeight anomaly ζ_i(i=1,2,3 ..., n), by ζ_iSubstituting into formula (2) is All mesh point P ' can be obtained₁、P’₂、P’₃、...、P’_nNormal high H_Correction(table 7).

The aerial survey of unmanned aerial vehicle of table 7 corrects height value

Step 6: sinkage is calculated：Carried out using the terrain elevation values obtained after correction with topographic map height value before exploitation Compare, calculate subsidence amount after exploitation.

Each planar point P " obtained after GPS height fittings are corrected₁、P”₂、P”₃、...、P”_nCoordinate X, Y and normal high Z_CorrectionMachine plotting carries out grid to the stacking chart by 100 meters × 100 meters of grid and drawn on the working face pre-existing topography figure of ore deposit 805 Point, each new mesh point of comparative analysis is adopted preceding and adopts rear height value, according to formula

Δ H=H₀-H_Correction (10)

Can calculate before the working face mining with subsidence value Δ H (table 8) after exploitation, finally draw corresponding Study area's subsidence isogram (as shown in Figure 3).In formula (10), Δ H is to exploit subsidence value, unit after preceding and exploitation For rice；H₀For the height value before exploitation on former topographic map, unit is rice；H_CorrectionFor each planar point after being corrected through GPS height fittings P”₁、P”₂、P”₃、...、P”_nHeight value, unit is rice.

The coal mining subsidence value of table 8

The scope of the heavy basin of the fully-mechanized mining working surface collapse of ore deposit 805 is more than working face, and maximum sinking value appears in working face Centre, overlapping exploitation exacerbates destruction of the surface collapse to earth's surface.Obtained sinking is calculated using the method for the invention Value and the surface movement parameters such as subsidence factor with use probability integration process formula it is expected that it is basically identical, illustrate to utilize unmanned aerial vehicle remote sensing Technology can be advantageously applied to the monitoring of coal mining surface subsidence settling amount.

Described above is only the preferable case study on implementation of the present invention, and not the present invention is imposed any restrictions, every according to this hair Modification, change and equivalent the measure change for any simple properties that bright technical spirit above case study on implementation is made, still fall within In the protection domain of technical solution of the present invention.

Claims (1)

1. a kind of method of the coal mining subsidence amount monitoring calculation based on unmanned aerial vehicle remote sensing, it is characterised in that utilize miniature self-service Machine remote sensing technology is calculated coal mining subsidence amount, and process is as follows：

Step 1: unmanned plane flight-line design：According to research area's topography and geomorphology, scope, position and ground resolution planning unmanned plane The line of flight；

Unmanned plane course line mainly includes flying height, Duplication and course, and process is as follows：

Step 101, flying height are determined：According to formula

$H=\frac{f\×GSD}{a}---\left(1\right)$

The unit that calculating obtains flying height H, H under definitely face resolution ratio correspondence is rice；In formula (1), GSD is ground distributor Resolution, unit is rice；F is camera lens focal length, and unit is millimeter；A is camera pixel dimension, and unit is micron；

Step 102, Duplication are determined：Determined according to research area's topography and landform character, setting region of no relief course and sidelapping rate Respectively 60% and 40%, mountain area course and sidelapping rate are respectively 80% and 60%；

Step 103, course are determined：Determined according to research area's shape along research area master is axial arranged；

Step 2: unmanned plane：The course line of design is imported into unmanned plane, manipulation unmanned plane is by setting airline operation and is spaced boat Clap collection research area's remote sensing image；Each sortie exports the remote sensing images collected and coordinate data after the completion of taking photo by plane, and carries out The inspection of degree of overlapping and the quality of image；

Step 3: remote sensing image processing：The remote sensing image and coordinate data that unmanned plane is obtained import Pix4D mapper 2.0 professional remote sensing image processing softwares, successively complete image mosaic, the encryption of point cloud, DOM and DSM generations processing, and by 1000cm Grid density size export research area in each mesh point P '₁、P’₂、P’₃、…、P’_nPlane coordinates X, Y and geodetic height Z_Greatly；

Step 4: ground control survey：A number of control point is laid according to research area's scope and geographical feature, and to these points Measure one by one, specific practice is as follows：

Step 401, control measurement points are determined and laid：The each point P ' that aerial survey of unmanned aerial vehicle is obtained₁、P’₂、P’₃、…、P’_nSpread In research Qu Tushang, m ground survey point P is selected from these points according to research area's area and topography and geomorphology₁、P₂、P₃、…、P_m (m >=5), measurement point is evenly distributed in research area, and by taking a fully-mechanized mining working as an example, the point of ground control survey is distributed in work Make four corners and center in face, and quantity is no less than 5；

Step 402, control point survey：Respectively to ground survey point P₁、P₂、P₃、…、P_mSetting-out is carried out, and by elevation Leveling Survey The Control Surveying in High Accuracy requirement of measurement obtains each measurement point P to these ground survey point testings one by one₁、P₂、P₃、…、P_m's Normal high Z_Just；

Step 5: GPS height fittings are corrected：Measurement result is controlled according to ground, intended using the GPS elevations of 4 parametric surface fittings Close calibration model to be corrected the three-dimensional coordinate that aerial survey of unmanned aerial vehicle is obtained, process is as follows：

Step 501,4 parameters are calculated：According to formula

H_Just=H_Greatly-ζ_i(i=1,2,3 ..., m) (2)

Calculate each common point P₁、P₂、P₃、...、P_mCorresponding abnormal elevation value ζ_i, in formula (2), H_JustFor normal high level, H_GreatlyFor The earth high level, unit is rice；ζ_iFor each common point P₁、P₂、P₃、…、P_mCorresponding abnormal elevation value；And according to formula

ζ_i=a₀+a₁x_i+a₂y_i+a₃x_iy_i(i=1,2,3 ..., m) (3)

Calculating obtains a₀、a₁、a₂And a₃In the value of this 4 parameters, formula (3), a₀、a₁、a₂And a₃For unknown parameter, ζ_iFor each public affairs Concurrent P₁、P₂、P₃、...、P_mAbnormal elevation value, unit is rice, x_iAnd y_iFor each common point P₁、P₂、P₃、...、P_mCorresponding plane Coordinate, unit is rice；

Because the number of common point is more than 4, corresponding error equation is listed

v_i=a₀+a₁x_i+a₂y_i+a₃x_iy_i-ζ_i(i=1,2,3 ..., m) (4),

In formula (4), v_iFor error amount；Being write as matrix form is

V=XA- ζ (5),

In formula (5),

It can be tried to achieve according to principle of least square method

$A=\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\end{array}\right]={\left(X^{T}X\right)}^{-1}{X}^T\mathrm{\ζ}---\left(6\right);$

Step 502, accuracy assessment：The correction of GPS height fittings includes precision of inner coincidence and precision of exterior coincidence, and process is as follows：

Step 5021, precision of inner coincidence are determined：According to formula

V_i=ζ_i-ζ_i' (7)

Calculate the residual error V for obtaining model of fit_i, in formula (7), V_iFor regression criterion, ζ_iFor n participation the Fitting Calculation known point Actual value, ζ_i' it is the n match value for participating in the Fitting Calculation known point；And according to formula

The precision of inner coincidence M of digital simulation model_{It is interior}, in formula (8), M_{It is interior}For precision of inner coincidence, V_iFor regression criterion, n' is participation The known point quantity of the Fitting Calculation；

Step 5022, precision of exterior coincidence are determined：According to formula

The precision of exterior coincidence M of digital simulation model_Outside；In formula (9), M_OutsideFor precision of exterior coincidence, V_iFor regression criterion, m' is to check Point quantity；

Step 503, elevation correction：Each parameter a that formula (6) is obtained₀、a₁、a₂And a₃All grids can be tried to achieve by substituting into formula (3) Point P '₁、P’₂、P’₃、...、P’_nHeight anomaly ζ_i(i=1,2,3 ..., n), by ζ_iIt is that can obtain all nets to substitute into formula (2) Lattice point P '₁、P’₂、P’₃、...、P’_nNormal high H_Correction；

Step 6: sinkage is calculated：It is compared using the terrain elevation values obtained after correction with topographic map height value before exploitation, Calculate subsidence amount after exploitation；

The planar point P " obtained after GPS height fittings are corrected₁、P”₂、P”₃、...、P”_nCoordinate X, Y and normal high Z_CorrectionMachine plotting Mesh generation is carried out to the stacking chart onto research area's pre-existing topography figure, and by 100 meters × 100 meters of grid, comparative analysis is each New grid points are adopted preceding and adopt rear height value, according to formula

Δ H=H₀-H_Correction (10)

Can calculate research paneling before with subsidence value Δ H after exploitation, finally draw corresponding equivalent curve Figure；In formula (10), Δ H is exploits subsidence value after preceding and exploitation, and unit is rice；H₀For on the preceding former topographic map of exploitation Height value, unit is rice；H_CorrectionFor each planar point P " after being corrected through GPS height fittings₁、P”₂、P”₃、...、P”_nHeight value, Unit is rice.