CN110440753B - High-precision DEM aviation gravity remote zone terrain correction method considering earth curvature - Google Patents

High-precision DEM aviation gravity remote zone terrain correction method considering earth curvature Download PDF

Info

Publication number
CN110440753B
CN110440753B CN201910746283.7A CN201910746283A CN110440753B CN 110440753 B CN110440753 B CN 110440753B CN 201910746283 A CN201910746283 A CN 201910746283A CN 110440753 B CN110440753 B CN 110440753B
Authority
CN
China
Prior art keywords
measuring point
dem
elevation
earth
coordinate system
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910746283.7A
Other languages
Chinese (zh)
Other versions
CN110440753A (en
Inventor
冯治汉
郭培虹
熊盛青
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
XI'AN CENTER OF GEOLOGICAL SURVEY CGS
Original Assignee
XI'AN CENTER OF GEOLOGICAL SURVEY CGS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by XI'AN CENTER OF GEOLOGICAL SURVEY CGS filed Critical XI'AN CENTER OF GEOLOGICAL SURVEY CGS
Priority to CN201910746283.7A priority Critical patent/CN110440753B/en
Publication of CN110440753A publication Critical patent/CN110440753A/en
Application granted granted Critical
Publication of CN110440753B publication Critical patent/CN110440753B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C5/00—Measuring height; Measuring distances transverse to line of sight; Levelling between separated points; Surveyors' levels
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V9/00—Prospecting or detecting by methods not provided for in groups G01V1/00 - G01V8/00
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/29—Geographical information databases

Landscapes

  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Remote Sensing (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

A high-precision DEM aviation gravity remote zone terrain correction method considering the curvature of the earth comprises the following steps: and S1, acquiring the position information of the measuring point Q. S2, extracting the average elevation value of all DEM units in the ground improvement range of the measuring point Q from a pre-constructed ground improvement elevation database after coordinate conversion and rotation until the measuring point Q is positioned at the earth north pole according to the position information of the measuring point Q
Figure DDA0002165681450000011
S3, according to the height average value of all DEM units in the measuring point Q ground change range
Figure DDA0002165681450000012
And combining a preset topographic correction value calculation formula of the measuring point Q to obtain the topographic correction value of the measuring point Q. And S4, accurately correcting the gravity measurement of the measuring point Q according to the terrain correction value of the measuring point Q. The accuracy of terrain correction in a gravity remote area is greatly improved, and meanwhile, the calculated amount is small and the calculation speed is high; the long-term work can splice the ground change elevation database finished in different periods so as to update and maintain the elevation database.

Description

High-precision DEM aviation gravity remote zone terrain correction method considering earth curvature
Technical Field
The invention relates to the technical field of geophysical exploration, in particular to a high-precision DEM aeronautical gravity remote terrain correction method considering the curvature of the earth.
Background
The terrain correction in the gravity survey refers to correcting the mass from the ground level to the natural terrain, and the correction is respectively carried out by distinguishing a near middle area from a far area. The near-middle region terrain improvement is to approximate the earth to be an infinite flat plate near a measuring point (within 2km of radius), eliminate the influence of surplus or loss of substances higher or lower than the plane of the measuring point and calculate by using a plane vertical square column formula; the remote terrain improvement is carried out by a regional gravity database in a region which is far away from a measuring point (within 166.7km beyond the radius of 2 km), and each unit can correct the regional gravity database before the gravity database is built.
The method has the advantages that the method is widely distributed in mountains and rich in mineral products in China, a large amount of rapid and efficient aerial gravity measurement needs to be carried out, but the influence of terrain change on gravity values is very large, and a high-precision aerial gravity remote terrain correction method is needed. In addition, the terrain correction range specified in regional gravity survey Specification (DZ/T0082-2006) of China is unified to 166.7km, and for comparison and joint explanation with the ground gravity, the terrain correction range of the aviation gravity is also unified to 166.7 km. The gravity anomaly precision calculated in the calculation range without considering the influence of the earth curvature cannot meet the requirement of gravity interpretation, and the difference between the gravity anomaly precision and the gravity value observed on the ground is large. Therefore, a high-precision DEM aviation gravity remote zone terrain correction method considering the curvature of the earth is needed.
In the current domestic gravity measurement, RGIS is adopted to calculate the remote terrain correction value, 1km multiplied by 1km node elevation is used to calculate the 5 'multiplied by 5' sector average elevation, and the terrain correction precision is about +/-0.2 multiplied by 10-5m/s2. However, there are two drawbacks to the RGIS software: firstly, high-precision DEM data cannot be used, and the requirement of current high-precision gravity measurement cannot be met; secondly, the measuring point is required to be positioned on the ground, and the requirement of aviation gravity measurement cannot be met. The domestic aviation gravity measurement is only in the test and preliminary popularization stage. At present, foreign aviation gravity measurement is mainly implemented in a plain area, although high-precision DEM data is adopted to calculate a gravity ground correction value, the terrain correction range is small, and the influence of the curvature of the earth is not considered during calculation.
And because mass data are obtained by aviation gravity measurement, the calculation speed of the currently adopted calculation method is low, and the production requirement is difficult to meet. Therefore, a method for quickly and efficiently calculating the terrain correction value of the aviation gravity remote area based on high-precision DEM data and by considering the influence of the earth curvature is needed.
Disclosure of Invention
Technical problem to be solved
In order to solve the problems in the prior art, the invention provides a high-precision DEM aviation gravity remote zone terrain correction method considering the curvature of the earth. The accuracy of terrain correction in a gravity remote area is greatly improved, and meanwhile, the calculated amount is small and the calculation speed is high; the long-term work can splice the ground change elevation database finished in different periods so as to update and maintain the elevation database.
(II) technical scheme
In order to achieve the purpose, the invention adopts the main technical scheme that:
a high-precision DEM aviation gravity remote zone terrain correction method considering the curvature of the earth comprises the following steps:
and step S1, acquiring the position information of the measuring point Q.
Step S2, according to the position information of the measuring point Q, extracting the average elevation value of all DEM units in the measuring point Q changing range from the pre-constructed ground changing elevation database after coordinate conversion and rotation until the measuring point Q is positioned at the earth north pole
Figure BDA0002165681430000022
Step S3, according to the height average value of all DEM units in the measuring point Q ground change range
Figure BDA0002165681430000023
And combining a preset topographic correction value calculation formula of the measuring point Q to obtain the topographic correction value of the measuring point Q.
And step S4, accurately correcting the gravity measurement of the measuring point Q according to the terrain correction value of the measuring point Q.
As an improvement of the method, a preset measuring point Q terrain correction value calculation formula comprises
Figure BDA0002165681430000021
Wherein σ is the ground density; rqR + Z, wherein R is the average radius of the earth, and Z is the elevation of a measuring point Q; e (theta, r) ═ 2-r2-r cosθ-3cos2θ)+3sin2θcosθln(r-cosθ+l)，
Figure BDA0002165681430000031
r is the distance from the geocentric, and theta is the latitude; theta1、θ2Representing the latitude range of the change range of the measuring point Q.
As a modification of the method of the present invention, after step S1 and before step S2, the method further comprises:
and S11, acquiring a digital elevation model DEM.
And S12, performing coordinate conversion and coordinate rotation to enable the measuring point Q to be positioned at the north pole of the earth and obtain a new space geodetic coordinate of the middle point P on the top surface of the DEM unit.
S13, calculating the average elevation value of all DEM units in the ground change range of the measuring point Q according to the new space geodetic coordinates of the middle point P on the top surface of the DEM unit
Figure BDA0002165681430000032
S14, establishing measuring point Q position information and elevation average value
Figure BDA0002165681430000033
An associated change of elevation database.
As a modification of the method of the present invention, step S12 includes:
and S121, converting the measuring point Q described by the projection coordinate system and the top surface midpoint P of the DEM unit into coordinates described by a space geodetic coordinate system.
And S122, converting the measuring point Q described by the space geodetic coordinate system and the top surface midpoint P of the DEM unit into coordinates described by a space rectangular coordinate system.
And S123, rotating the space rectangular coordinate system until the measuring point Q is positioned at the north pole of the earth, obtaining a new space rectangular coordinate system, and calculating the rotation angle of the space rectangular coordinate system.
And S124, calculating a new space rectangular coordinate of the central point P on the top surface of the DEM unit according to the space rectangular coordinate of the central point P on the top surface of the DEM unit and the rotation angle of the space rectangular coordinate system.
And S125, calculating a new space geodetic coordinate of the middle point P on the top surface of the DEM unit according to the new space rectangular coordinate of the middle point P on the top surface of the DEM unit.
As an improvement of the method, the origin of the space rectangular coordinate system is the geocentric, and the Z axis passes through the North Pole of the earth.
As a modification of the method of the invention, in step S13, the average elevation value of all DEM units in the area of the change of the measuring point Q is calculated
Figure BDA0002165681430000034
The method comprises the following steps: calculating the average elevation value of all DEM units in the range of measuring point Q longitude of 0-2 pi and latitude of 20-166.7km away from the arctic arc length
Figure BDA0002165681430000035
(III) advantageous effects
The invention has the beneficial effects that:
1. based on the high-precision DEM, the influence of the curvature of the earth on a calculation result is considered, so that the correction precision of the terrain in the ground gravity remote area is greatly improved.
2. The ground change elevation data only correspond to the coordinates of the measuring points and are irrelevant to the elevations of the measuring points, so that the ground change elevation data are suitable for gravity measurement work of the measuring points at different elevations.
3. The method avoids huge calculation workload caused by direct calculation of land improvement, adopts a land improvement elevation database of a land improvement range to be formed or used firstly, and uses a spherical shell sector to calculate the land improvement of a gravity remote area by reading the average elevation in the database so as to achieve the aim of rapid calculation.
4. The long-term work can splice the ground change elevation database finished in different periods so as to update and maintain the elevation database.
Drawings
The invention is described with the aid of the following figures:
FIG. 1 is a schematic diagram of the earth as approximately a sphere in an embodiment of the present invention;
FIG. 2 is a schematic diagram of ellipsoid integration after transformation of spherical coordinates according to an embodiment of the present invention;
FIG. 3 is a flow chart of a high-precision DEM aviation gravity remote zone terrain correction method considering the curvature of the earth in the embodiment of the invention.
Detailed Description
For the purpose of better explaining the present invention and to facilitate understanding, the present invention will be described in detail by way of specific embodiments with reference to the accompanying drawings.
As shown in fig. 1, the earth is approximated as a sphere, and the center of the earth is the center of the sphere. As shown in fig. 2, to account for earth curvature and based on high precision DEM data,calculating the terrain correction value of an aviation gravity far zone, and calculating by the applicant through coordinate conversion and ellipsoid integral
Figure BDA0002165681430000041
Point fan-shaped column pair aviation gravity measuring point
Figure BDA0002165681430000042
The gravity influence value of (a).
The sphere center is taken as the origin of coordinates, the Z axis passes through the top end (north pole) of the sphere, and the aviation gravity measuring point Q is on the Z axis. The area from the earth level to the free surface of the earth and around 20-166.7km of an aviation gravity measuring point is divided into theta along the meridian direction0,θ1,θ2,…,θnWeft loops and
Figure BDA0002165681430000043
divided into m × n crust segments. The latitude, longitude and elevation ranges of sector (i, j) with the top surface midpoint P are
Figure BDA0002165681430000051
The terrain influence value of the aviation gravity measuring point Q is as follows:
Figure BDA0002165681430000052
wherein R is the average radius of the earth, and 6371.025km is taken; z is the elevation of an aviation gravity measuring point Q; h is1And h2The elevation of the bottom surface and the elevation of the top surface of the sector (i, j) respectively; g is the gravitational constant of 6.67X 10-8cm3/(g·s2)；σijThe average density of the segments (i, j) is generally 2.67g/cm3(ii) a r is the distance from the geocentric.
The integration result of equation (1) is:
Figure BDA0002165681430000053
wherein E (theta),r)＝(2-r2-r cosθ-3cos2θ)+3sin2θcosθln(r-cosθ+l)，
Figure BDA0002165681430000054
As can be seen from equation (2), if the aviation gravity measures the point
Figure BDA0002165681430000055
And when the DEM is positioned at the north pole of the earth, the gravity influence value of each DEM on the measuring point can be calculated as long as the longitude, the latitude and the elevation of each DEM in the new coordinate system are solved.
And calculating the gravity influence value of the aerial gravity measuring point Q in actual work based on the direct calculation method of the high-precision DEM aerial gravity far-zone terrain correction value considering the curvature of the earth. The method specifically comprises the following steps:
and step S1, acquiring position information of the measuring point Q and a digital elevation model DEM, and performing coordinate conversion and coordinate rotation to enable the measuring point Q to be positioned at the north pole of the earth and obtain new space geodetic coordinates of a middle point P on the top surface of the DEM unit.
The step S1 includes:
step S11, converting the plane projection coordinate system into a spatial geodetic coordinate system: converting a measuring point Q described by a projection coordinate system and a central point P on the top surface of the DEM unit into a space geodetic coordinate system
Figure BDA0002165681430000056
The coordinates of the description. Such that the coordinates of both Q and P are in longitude
Figure BDA0002165681430000061
Latitude θ and altitude elevation h.
In actual operation, the coordinate of the measuring point Q is generally described by using a plane projection coordinate, such as a Gaussian 6-degree projection coordinate, and the like, and some coordinate data of the high-precision DEM unit are described by using a plane projection coordinate system and some coordinate data of the high-precision DEM unit are described by using a space geodetic coordinate system. When the gravity remote terrain correction value is calculated, firstly, coordinate conversion is carried out, a coordinate system is unified, and then the gravity terrain correction calculation can be carried out. The invention uniformly converts the projection coordinates of the Gaussian 6-degree band into a space geodetic coordinate system. The inverse gaussian projection formula is:
Figure BDA0002165681430000062
Figure BDA0002165681430000063
wherein L is0Is the central meridian longitude;
Figure BDA0002165681430000064
Mf＝a(1-e2)(1-e2sin2Bf)-3/2；tf＝tanBf；
Figure BDA0002165681430000065
Bfcalculating the latitude of the bottom point, namely the latitude corresponding to the meridian arc length when X is equal to X, according to the meridian arc length formula in an iterative way:
Figure BDA0002165681430000066
step S12, converting the spatial geodetic coordinate system into a spatial rectangular coordinate system: will be in longitude
Figure BDA0002165681430000067
And converting the measuring point Q represented by the latitude theta and the altitude elevation h and the central point P on the top surface of the DEM into coordinates (x, y and z) of a space rectangular coordinate system. Such that the coordinates of Q and P are both represented by x, y and z. Preferably, the coordinate origin of the rectangular space coordinate system is the earth center, and the Z axis passes through the north pole of the earth.
Figure BDA0002165681430000068
Because the three-dimensional coordinate rotation by using the geodetic coordinate system is difficult, the space geodetic coordinate system is converted into a space rectangular coordinate system to prepare for coordinate rotation.
And step S13, rotating the space rectangular coordinate system until the measuring point Q is positioned at the north pole of the earth, obtaining a new space rectangular coordinate system, and calculating the rotation angle of the space rectangular coordinate system. The gravity calculation is simplified. The method comprises the following steps:
and respectively rotating the X, Y and Z axes of the space rectangular coordinate system by alpha, beta and gamma to obtain a new space rectangular coordinate system X ', Y' and Z ', wherein the coordinates of the original measuring point Q (X, Y and Z) in the new coordinate system are changed into Q (X', Y 'and Z'). When the coordinate system is rotated, first, γ is rotated around the z-axis, and then the coordinate of Q (x, y, z) is (x)1,y1,z)：
Figure BDA0002165681430000071
Then, after the coordinate system rotates around the x axis by alpha, the coordinate of the point Q is (x)1,y′,z2)：
Figure BDA0002165681430000072
After the coordinate system is rotated by β around the y-axis, the coordinate of the point Q is (x ', y ', z '):
Figure BDA0002165681430000073
let x '═ y' ═ 0,
Figure BDA0002165681430000074
solving the above equations (7), (8), (9) to obtain the rotation angles α, β, γ respectively as:
Figure BDA0002165681430000075
and S14, calculating a new space rectangular coordinate of the top surface midpoint P of the DEM unit according to the space rectangular coordinate of the top surface midpoint P of the DEM unit and the rotation angle of the space rectangular coordinate system.
The spatial rectangular coordinate P (x, y, z) of the middle point of the top surface of the DEM is represented as P (x ', y ', z ') in a new spatial rectangular coordinate system:
Figure BDA0002165681430000076
and S15, calculating a new space geodetic coordinate of the central point P on the top surface of the DEM unit according to the new space rectangular coordinate of the central point P on the top surface of the DEM unit.
Figure BDA0002165681430000077
Wherein the content of the first and second substances,
Figure BDA0002165681430000078
and the longitude of the central point P of the top surface of the DEM unit in the new space geodetic coordinate system is shown, theta 'is the latitude, and r' is the distance from the geocenter.
And step S2, calculating the gravity influence value of the measuring point Q according to the new space geodetic coordinates of the top surface middle points P of all DEM units in the ground changing range of the measuring point Q by combining a formula (2).
Although the calculation of the gravity influence value at the measuring point is based on the high-precision DEM data and takes the influence of the curvature of the earth into consideration, the ground change value of the measuring point is directly calculated, so that huge calculation workload is caused, the calculation speed is low, and the production requirement is difficult to meet.
For this purpose, the applicant optimizes step S2 in the above direct calculation method of high-precision DEM airborne gravity far-field terrain correction value considering the curvature of the earth, and the step S2 after optimization includes:
step S21, calculating the average elevation value of all DEM units in the ground change range of the measuring point Q according to the new space geodetic coordinates of the middle point P on the top surface of the DEM unit
Figure BDA0002165681430000081
Average elevation value
Figure BDA0002165681430000082
The remote terrain correction value sigma delta g of a measuring point is obtained by accumulating the gravity values calculated by each fan-shaped blocki,j. Then, the floor-modifying value sigma delta g is back-calculated in the floor-modifying range by the formula (1)i,jCorresponding elevation h2(taking h in general)10). Here h is2It is the average elevation over the range of surveyed point changes. Due to elevation mean
Figure BDA0002165681430000083
The elevation of the sector is not directly calculated, but the elevation is obtained by solving a nonlinear equation through the calculated ground change value, errors mainly come from coordinate conversion, integral calculation and the nonlinear equation, and the method is high in calculation speed, controllable in errors and high in precision.
Step S22, according to the height average value of the midpoint P of the top surfaces of all DEM units in the measuring point Q ground changing range
Figure BDA0002165681430000084
The gravity influence value of the measurement point Q is calculated in combination with the modification (13) of the formula (2).
Figure BDA0002165681430000085
Wherein, the density is changed by sigma, and is generally 2.67g/cm3；Rq＝R+Z；θ1、θ2Representing the latitude range of the change range of the measuring point Q.
Above average value according to elevation
Figure BDA0002165681430000086
And the gravity influence value of the measuring point Q is calculated, so that the calculation workload is relatively small, and the calculation speed is high.
The applicant further optimizes the method for correcting the terrain in the high-precision DEM aeronautical gravity far zone considering the curvature of the earth, and the proposed method for correcting the terrain in the aeronautical gravity far zone, as shown in FIG. 3, comprises the following steps:
and step S1, constructing and modifying the elevation database. The method comprises the following steps:
and step S11, acquiring position information of the measuring point Q and a digital elevation model DEM, and performing coordinate conversion and coordinate rotation to enable the measuring point Q to be positioned at the north pole of the earth and obtain new space geodetic coordinates of a middle point P on the top surface of the DEM unit.
Step S12, calculating the average elevation value of all DEM units in the ground change range of the measuring point Q according to the new space geodetic coordinates of the middle point P on the top surface of the DEM unit
Figure BDA0002165681430000091
Step S13, establishing measuring point Q position information and elevation average value
Figure BDA0002165681430000092
An associated change of elevation database.
Step S2, according to the position information of the measuring point Q, extracting the average elevation value of all DEM units in the measuring point Q ground change range from the constructed ground change elevation database after coordinate conversion and rotation until the measuring point Q is positioned at the earth north pole
Figure BDA0002165681430000093
Step S3, according to the height average value of all DEM units in the measuring point Q ground change range
Figure BDA0002165681430000094
And (13) calculating the gravity influence value of the measuring point Q.
And step S4, accurately correcting the gravity measurement of the measuring point Q according to the gravity influence value of the measuring point Q.
Preferably, in step S12, the average elevation value of all DEM units in the variation range of the measuring point Q is calculated
Figure BDA0002165681430000095
Includes calculating the longitude range of the measuring point Q to be 0-2 pi, and the latitude range to be 20-166.7km from the arc length of the north poleAverage elevation value of all DEM units in enclosure
Figure BDA0002165681430000096
The optimized method is based on the high-precision DEM, and simultaneously considers the influence of the curvature of the earth on the calculation result; the ground change elevation data only correspond to the coordinates of the measuring points and are irrelevant to the elevations of the measuring points, so that the ground change elevation data are suitable for gravity measurement work of the measuring points at different elevations; the method avoids huge calculation workload caused by direct calculation of the land improvement, adopts a land improvement elevation database of a land improvement range to be formed or used firstly, and uses spherical shell segments to calculate the land improvement value of a gravity remote area by reading the average elevation in the database so as to achieve the aim of quick calculation; the long-term work can splice the ground change elevation database finished in different periods so as to update and maintain the elevation database.
In conclusion, the method provided by the invention fills the blank of correcting the terrain of the aviation gravity remote area in consideration of the curvature of the earth; the defect that the high-precision DEM cannot be used by the RGIS is overcome, and the correction precision of the terrain in the ground gravity far zone is greatly improved; by pre-calculating and changing the database, the speed of terrain correction calculation can be greatly improved on the premise of ensuring the calculation precision.
The method is realized by a computer program, and the operation body is essentially a control terminal.
It should be understood that the above description of specific embodiments of the present invention is only for the purpose of illustrating the technical lines and features of the present invention, and is intended to enable those skilled in the art to understand the contents of the present invention and to implement the present invention, but the present invention is not limited to the above specific embodiments. It is intended that all such changes and modifications as fall within the scope of the appended claims be embraced therein.

Claims (5)

1. A high-precision DEM aviation gravity remote zone terrain correction method considering the curvature of the earth is characterized by comprising the following steps:
step S1, acquiring the position information of the measuring point Q;
step S2,According to the position information of the measuring point Q, extracting the average elevation value of all DEM units in the ground improvement range of the measuring point Q from a pre-constructed ground improvement elevation database after coordinate conversion and rotation until the measuring point Q is positioned at the north pole of the earth
Figure FDA0002920027560000011
Step S3, according to the height average value of all DEM units in the measuring point Q ground change range
Figure FDA0002920027560000012
Combining a preset topographic correction value calculation formula of the measuring point Q to obtain a topographic correction value of the measuring point Q;
step S4, according to the terrain correction value of the measuring point Q, the gravity measurement of the measuring point Q is corrected accurately;
after step S1, before step S2, the method further includes:
s11, acquiring a digital elevation model DEM;
s12, performing coordinate conversion and coordinate rotation to enable the measuring point Q to be located at the north pole of the earth and obtain a new space geodetic coordinate of the middle point P on the top surface of the DEM unit;
s13, calculating the average elevation value of all DEM units in the ground change range of the measuring point Q according to the new space geodetic coordinates of the middle point P on the top surface of the DEM unit
Figure FDA0002920027560000013
S14, establishing measuring point Q position information and elevation average value
Figure FDA0002920027560000014
An associated change of elevation database.
2. The method of claim 1, wherein the predetermined survey point Q terrain correction value calculation formula comprises:
Figure FDA0002920027560000015
wherein σ is the ground density; rqR + Z, wherein R is the average radius of the earth, and Z is the elevation of a measuring point Q; e (theta, r) ═ 2-r2-r cosθ-3cos2θ)+3sin2θcosθln(r-cosθ+l)，
Figure FDA0002920027560000016
r is the distance from the geocentric, and theta is the latitude; theta1、θ2Representing the latitude range of the change range of the measuring point Q.
3. The method according to claim 1, wherein the step S12 includes:
s121, converting a measuring point Q described by a projection coordinate system and a central point P on the top surface of the DEM unit into coordinates described by a space-earth coordinate system;
s122, converting a measuring point Q described by the space geodetic coordinate system and a central point P on the top surface of the DEM unit into coordinates described by a space rectangular coordinate system;
s123, rotating the space rectangular coordinate system until the measuring point Q is positioned at the north pole of the earth to obtain a new space rectangular coordinate system, and calculating the rotation angle of the space rectangular coordinate system;
s124, calculating a new space rectangular coordinate of the midpoint P on the top surface of the DEM unit according to the space rectangular coordinate of the midpoint P on the top surface of the DEM unit and the rotation angle of the space rectangular coordinate system;
and S125, calculating a new space geodetic coordinate of the middle point P on the top surface of the DEM unit according to the new space rectangular coordinate of the middle point P on the top surface of the DEM unit.
4. The method of claim 3, wherein the rectangular spatial coordinate system has an origin of earth centered and a Z-axis passing through Earth north.
5. The method of claim 1, wherein in step S13, the average elevation value of all DEM units within the range of changes of the measurement point Q is calculated
Figure FDA0002920027560000021
Comprises that
Calculating the average elevation value of all DEM units in the range of measuring point Q longitude of 0-2 pi and latitude of 20-166.7km away from the arctic arc length
Figure FDA0002920027560000022
CN201910746283.7A 2019-08-13 2019-08-13 High-precision DEM aviation gravity remote zone terrain correction method considering earth curvature Active CN110440753B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910746283.7A CN110440753B (en) 2019-08-13 2019-08-13 High-precision DEM aviation gravity remote zone terrain correction method considering earth curvature

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910746283.7A CN110440753B (en) 2019-08-13 2019-08-13 High-precision DEM aviation gravity remote zone terrain correction method considering earth curvature

Publications (2)

Publication Number Publication Date
CN110440753A CN110440753A (en) 2019-11-12
CN110440753B true CN110440753B (en) 2021-04-09

Family

ID=68435178

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910746283.7A Active CN110440753B (en) 2019-08-13 2019-08-13 High-precision DEM aviation gravity remote zone terrain correction method considering earth curvature

Country Status (1)

Country Link
CN (1) CN110440753B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112462443B (en) * 2020-11-13 2021-07-20 中国自然资源航空物探遥感中心 Synchronous terrain correction method and device for aerial gravity measurement
CN114140397A (en) * 2021-11-13 2022-03-04 北京中勘迈普科技有限公司 Method and system for correcting gravity near-zone terrain by full-digital ground imaging method

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1647138A (en) * 2002-04-22 2005-07-27 Dgs计算机株式会社 Digital altimetric map drawing method and device

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1007293B (en) * 1985-09-20 1990-03-21 乔斯塔公司 Utilize the topomap location of satellite and storage and the system of message transfer
CN202748482U (en) * 2012-06-06 2013-02-20 中国地质调查局西安地质调查中心 Gravity measurement short-range terrain correction device and gravity measurement short-range terrain correction system
CN104155699B (en) * 2014-08-13 2017-02-22 昆明理工大学 Method of positioning and detecting high density concealed ore body in full spatial domain by tunnel gravity
JP6468606B2 (en) * 2016-02-01 2019-02-13 滋樹 水谷 Method for estimating surface density of gravity deviation data
CN207992458U (en) * 2017-12-24 2018-10-19 航天恒星科技有限公司 Carbon global position system under complicated observation mode

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1647138A (en) * 2002-04-22 2005-07-27 Dgs计算机株式会社 Digital altimetric map drawing method and device

Also Published As

Publication number Publication date
CN110440753A (en) 2019-11-12

Similar Documents

Publication Publication Date Title
CN102506824B (en) Method for generating digital orthophoto map (DOM) by urban low altitude unmanned aerial vehicle
CN103679711B (en) A kind of remote sensing satellite linear array push sweeps optics camera outer orientation parameter calibration method in-orbit
CN105158760B (en) Method for inverting underground fluid volume change and three dimension surface deformation using InSAR
CN106871927A (en) A kind of UAV electro-optical pod's alignment error Calibration Method
CN108226982B (en) Single linear array satellite laser combined high-precision positioning processing method
CN110440753B (en) High-precision DEM aviation gravity remote zone terrain correction method considering earth curvature
CN103673995A (en) Calibration method of on-orbit optical distortion parameters of linear array push-broom camera
CN104296751A (en) Layout design method of multi-star sensor configuration layout
CN113538595B (en) Method for improving geometric precision of remote sensing stereo image by using laser height measurement data in auxiliary manner
CN103207419A (en) Three-dimensional measurement method for tunnel rock formation attitude
CN111256730A (en) Earth mass balance correction calculation method for low-altitude oblique photogrammetry technology
CN110146052B (en) Plane normal astronomical directional measurement method and system based on total station
CN106875330B (en) Method for rotating plane model into spherical model
CN114608531A (en) GNSS continuous operation reference station pier mark inclination measuring method
CN110030968A (en) A kind of ground shelter measurement of elevation method based on spaceborne stereoptics image
CN113252009A (en) Earth and stone calculation method based on unmanned aerial vehicle aerial survey technology
CN113255740A (en) Multisource remote sensing image adjustment positioning precision analysis method
CN111307125B (en) Inclined-axis cylindrical projection method based on GNSS and ground ranging combined adjustment
CN110310370B (en) Method for point-plane fusion of GPS (Global positioning System) and SRTM (short Range TM)
CN115236759B (en) Hexagonal grid subdivision method for determining earth gravity field
CN111060140A (en) Polar region inertial navigation error obtaining method under earth ellipsoid model
CN105928513B (en) A kind of airborne synthetic aperture radar movement parameter measurement method based on position and attitude measuring system
CN106767824B (en) Method for calculating relative position of double detectors on surface of extraterrestrial celestial body
CN113532377B (en) Method for assisting adjustment of area network by using high-resolution seven-grade laser height measurement data
CN111914396B (en) Sub-grid terrain three-dimensional earth surface solar radiation forced effect rapid parameterization method based on high-resolution DEM data

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant