CN115143942A - Satellite photogrammetry earth positioning method based on photon point cloud assistance - Google Patents

Satellite photogrammetry earth positioning method based on photon point cloud assistance Download PDF

Info

Publication number
CN115143942A
CN115143942A CN202210843831.XA CN202210843831A CN115143942A CN 115143942 A CN115143942 A CN 115143942A CN 202210843831 A CN202210843831 A CN 202210843831A CN 115143942 A CN115143942 A CN 115143942A
Authority
CN
China
Prior art keywords
elevation
grid
point cloud
point
points
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.)
Granted
Application number
CN202210843831.XA
Other languages
Chinese (zh)
Other versions
CN115143942B (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.)
Guangdong University of Technology
Original Assignee
Guangdong University of Technology
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 Guangdong University of Technology filed Critical Guangdong University of Technology
Priority to CN202210843831.XA priority Critical patent/CN115143942B/en
Publication of CN115143942A publication Critical patent/CN115143942A/en
Application granted granted Critical
Publication of CN115143942B publication Critical patent/CN115143942B/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
    • G01C11/00—Photogrammetry or videogrammetry, e.g. stereogrammetry; Photographic surveying
    • G01C11/02—Picture taking arrangements specially adapted for photogrammetry or photographic surveying, e.g. controlling overlapping of pictures
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B11/00—Measuring arrangements characterised by the use of optical techniques
    • G01B11/02—Measuring arrangements characterised by the use of optical techniques for measuring length, width or thickness
    • G01B11/06—Measuring arrangements characterised by the use of optical techniques for measuring length, width or thickness for measuring thickness ; e.g. of sheet material
    • G01B11/0608—Height gauges
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88—Lidar systems specially adapted for specific applications
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G06T7/337—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving reference images or patches
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10028—Range image; Depth image; 3D point clouds
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Electromagnetism (AREA)
  • Multimedia (AREA)
  • Image Processing (AREA)
  • Optical Radar Systems And Details Thereof (AREA)

Abstract

A satellite photogrammetry is a positioning method to the ground based on photon point cloud assistance, which belongs to the field of satellite photogrammetry and remote sensing. The invention solves the problem of low positioning precision of the existing combined processing method. The method mainly comprises the following steps: acquiring stereo image data of a target area and generating a photogrammetric point cloud of a ground object; acquiring and filtering photon point cloud, and generating an elevation sampling point according to the filtered photon point cloud; registering the elevation sampling point cloud and the photogrammetric point cloud to obtain a transformation matrix; obtaining control points according to the transformation matrix; and correcting the system error of the photogrammetric point cloud according to the obtained control point. The method can be applied to the fields of satellite photogrammetry and remote sensing.

Description

Satellite photogrammetry earth positioning method based on photon point cloud assistance
Technical Field
The invention belongs to the field of satellite photogrammetry and remote sensing, and particularly relates to a satellite photogrammetry and ground positioning method based on photon point cloud assistance.
Background
With the development of optical remote sensing satellites and digital photogrammetry technologies, satellite photogrammetry has become a remote sensing mapping technology widely applied to the field of national defense economic construction. The ground resolution of the existing optical remote sensing satellites such as the third resource, the seventh high-resolution, the wordview3/4, the Pleiades and the Cartosat-3 are all fully stepped into the level of decimeters. Various high-resolution remote sensing satellites are transmitted and operated sequentially, and image data with abundant quantity and excellent quality are provided for various industry fields. Meanwhile, scholars at home and abroad carry out a great deal of research on a theoretical method for topographic mapping by using remote sensing images, and a mature satellite photogrammetry technical system is formed.
In recent years, the emerging satellite laser height measurement technology has made great progress, the earth observation positioning capability is greatly enhanced, and strong data support is provided for the fields of topographic mapping and the like. The ICESat-2 satellite adopts an Advanced Terrain Laser Altimeter System (ATLAS), has strong Laser detection capability, and can acquire ground surface Laser point cloud data with higher observation frequency and higher positioning accuracy.
Because the optical satellite is influenced by various factors such as attitude and orbit observation errors, platform flutter, atmospheric refraction, internal splicing and the like in the imaging process, a positioning result contains complex system errors, a certain number of uniformly distributed ground control points need to be acquired, and the positioning accuracy is improved through block adjustment. However, it is very difficult to obtain enough high-precision ground control points for areas such as mountain forests, deserts and plateaus, and the operation cost is high. The satellite laser height measurement technology can automatically acquire high-precision height measurement data in orbit without being influenced by landform. Therefore, the ICESat-2 photon point cloud and the optical remote sensing image are processed in a combined mode, the traditional manual field measurement point is replaced by the laser height measurement point, and the efficiency of extracting the earth surface point cloud by the satellite photogrammetry technology can be greatly improved. For the joint processing of satellite-borne laser radar data and satellite remote sensing images, the existing research shows that the positioning precision of the satellite photogrammetry technology can be obviously improved by fusing laser height measurement data. The current commonly used joint processing scheme is to screen out laser points in a flat area from point clouds to be used as control points, and then in an image space coordinate, the control points are used as weighted observation values to carry out block adjustment. The scheme can effectively improve the positioning precision of satellite photogrammetry without field control points, but faces the following problems to be further solved:
(1) Different from the traditional pulse type laser radar, the current laser height measurement satellite adopts a novel photon counting laser radar, and the observation data of the laser height measurement satellite is in a point cloud form. Because the laser radar has certain penetrability to the vegetation, the photon point cloud in the mountain area contains radiation information from crown reflection and surface reflection, and the optical remote sensing image mainly shows a crown target.
(2) Establishing an accurate corresponding relation between a laser point and an image is a precondition for carrying out combined processing on two kinds of data, and for the precondition, an image space registration scheme is generally adopted at present, namely, the laser point in a flat area is selected firstly, an original image positioning model or an image positioning model after adjustment processing of a free net is utilized, and an image space projection point of the laser point is calculated to serve as an image corresponding point.
In summary, the photon point cloud contains radiation information reflected from the earth surface and an image positioning model in an image registration scheme has errors, so that the positioning accuracy of the conventional joint processing method is low. Therefore, when the satellite image data and the photon radar data are jointly processed, the crown elevation needs to be extracted from the photon point cloud in a classified mode, and consistency of the target elevation is guaranteed. Then, researching the photon point cloud filtering classification algorithm to extract accurate surface elevation points has important significance for the combined processing scheme. In addition, in view of the defects of the existing combined processing scheme, it is necessary to deeply research an object-side registration scheme of a laser point and a remote sensing image, and necessary conditions are provided for the combined processing of satellite image data and photon radar data.
Disclosure of Invention
The invention aims to solve the problem of low positioning accuracy of the existing joint processing method, and provides a satellite photogrammetry earth positioning method based on photon point cloud assistance.
The technical scheme adopted by the invention for solving the technical problems is as follows:
a satellite photogrammetry earth positioning method based on photon point cloud assistance specifically comprises the following steps:
step one, acquiring stereoscopic image data of a target area;
secondly, generating a photogrammetric point cloud of the ground object by utilizing the three-dimensional image data;
acquiring an ICESat-2 photon point cloud in a stereopair region;
step four, performing grid filtering on the photon point cloud obtained in the step three, and filtering noise points to obtain a filtered photon point cloud;
generating elevation sampling points at equal intervals according to the filtered photon point cloud, and respectively calculating the ground elevation value of each elevation sampling point;
registering the elevation sampling point cloud and the photogrammetric point cloud generated in the second step according to the ground elevation value of each elevation sampling point to obtain a transformation matrix; then obtaining elevation sampling points serving as control points according to the transformation matrix;
and step seven, correcting the system error of the photogrammetric point cloud according to the obtained control point.
The beneficial effects of the invention are:
1. the invention provides a photon point cloud filtering algorithm based on grid continuity constraint, and the photon point cloud filtering and elevation extracting effects are excellent. The ICESat-2/ATL03 data test result shows that the algorithm can efficiently filter noise points and extract earth surface elevation points, the filtering success rate is higher than 94%, the elevation extraction precision is better than 5.4m, and the stable processing effect can be achieved on photon point cloud data of different terrain changes, terrain types and laser energy.
2. The joint processing algorithm based on the linear and planar point cloud registration provided by the invention obviously improves the positioning accuracy of the photogrammetric point cloud. Test results show that after ICESat-2 photon point cloud assistance is introduced, the X direction precision, the Y direction precision and the Z direction precision of the point cloud measured by the resource No. three camera respectively reach 0.46m, 2.43m and 2.13m, the X direction precision, the Y direction precision and the Z direction precision of the point cloud measured by the high-resolution No. seven camera respectively reach 0.36m, 1.34m and 0.59m, and the positioning precision is obviously improved.
Drawings
FIG. 1 is a flow chart of the method of the present invention;
FIG. 2 is a flow chart of ICESat-2 photon point cloud filtering and elevation extraction;
FIG. 3 is a schematic diagram illustrating grid continuity determination;
FIG. 4 is a schematic view of a point cloud trending process;
FIG. 5 is a schematic view of elevation sampling points;
FIG. 6a is a schematic diagram of adjacent points of an elevation sampling point;
FIG. 6b is a schematic diagram of the statistics of the elevation frequency of the neighborhood point cloud;
FIG. 7 is a flow chart of registration of a photon point cloud and a photogrammetric point cloud;
FIG. 8 is a schematic view of elevation similarity calculation;
FIG. 9a is a schematic view of a resource number three photogrammetry point cloud and an ICESat-2 point cloud;
FIG. 9b is a schematic diagram of the resource No. three photogrammetry point cloud system before error correction;
FIG. 9c is a schematic view of the resource No. three photogrammetric point cloud system after error correction;
FIG. 10a is a schematic view of a top seven photogrammetric point cloud and an ICESat-2 point cloud;
FIG. 10b is a schematic diagram of a high resolution seven-camera measurement point cloud system before error correction;
fig. 10c is a schematic diagram of the high-resolution seven-camera measurement point cloud system after error correction.
Detailed Description
First embodiment this embodiment will be described with reference to fig. 1. The satellite photogrammetry ground positioning method based on photon point cloud assistance in the embodiment specifically comprises the following steps:
step one, acquiring stereo image data of a target area;
secondly, generating a photogrammetric point cloud of the ground object by utilizing the three-dimensional image data;
acquiring ICESat-2 photon point cloud in a stereopair area (namely a target area);
step four, performing grid filtering on the photon point cloud obtained in the step three, and filtering noise points to obtain a filtered photon point cloud;
step five, generating elevation sampling points at equal intervals according to the filtered photon point cloud, and then respectively calculating the ground elevation value of each elevation sampling point;
registering the elevation sampling point cloud and the photogrammetric point cloud generated in the second step according to the ground elevation value of each elevation sampling point to obtain a transformation matrix; then obtaining elevation sampling points as control points according to the transformation matrix;
and step seven, correcting the system error of the photogrammetric point cloud according to the obtained control point.
The second embodiment is as follows: this embodiment will be described with reference to fig. 2. The difference between this embodiment and the first embodiment is that the specific process of the fourth step is:
step four, grid division
Carrying out grid division on the whole photon point cloud profile range to form a grid with the size of w multiplied by h;
step two, counting the number of points contained in each grid in the grid respectively, and for any grid in a column, selecting the first t grids containing the most points from all grids in the column as candidate grids of the column:
Figure BDA0003751417050000041
wherein G is j Representing the extracted candidate grid in the jth column,
Figure BDA0003751417050000042
is the number of dots contained in the mth grid in the jth column, m is the total number of grids in that column,
Figure BDA0003751417050000043
representing that t grids containing the most number of points are selected from the jth column;
in the same way, respectively obtaining candidate grids of each column;
step four and three, detecting the continuity of the grid
Step four, taking a certain candidate grid in the j-th row as a grid to be detected, and taking the candidate grid in the j-1-th row, the j-2-th row, \8230, the candidate grid in the j-k 'th row, the j + 1-th row, the j + 2-th row, \8230, the candidate grid in the j + k' th row as a neighborhood range of the grid to be detected (namely, selecting the candidate grid in the grid of each adjacent k 'row at the left and right sides of the grid to be detected as the neighborhood range, and the selected left side row number and the right side row number can be unequal), wherein the candidate grid in the j-k' th row, the j-1-th row, the j-2-th row, \8230, the candidate grid in the j-k 'th row forms a left connected domain, and the candidate grid in the j + 1-th row, the j + 2-th row, \8230, the j + k' th row forms a right connected domain; as shown in fig. 3.
Figure BDA0003751417050000044
Figure BDA0003751417050000045
In the formula (I), the compound is shown in the specification,
Figure BDA0003751417050000046
is the set of grids in the left connected domain,
Figure BDA0003751417050000047
is a grid set in a right connected domain, (i, j) is a row and column serial number of a grid to be detected in the grid,
Figure BDA0003751417050000048
for the m1 st grid in the set of grids in the left connected domain,
Figure BDA0003751417050000049
the grid is the (n 1) th grid in the grid set in the right connected domain;
step four, step two, calculating the continuity score Continue _ score of the grid to be detected:
Figure BDA00037514170500000410
in the formula (I), the compound is shown in the specification,
Figure BDA0003751417050000051
is the number of grids in the left connected domain,
Figure BDA0003751417050000052
the number of grids in the right connected domain;
step four, calculating the continuity score of each candidate grid in the jth column of the grid in a traversing manner, and taking the candidate grid with the highest continuity score and t' grids adjacent to the upper part and the lower part of the candidate grid with the highest continuity score as signal grids in the jth column;
step four, step three, the process from the step four, step three, to the step four, step three is repeated, obtain the signal grid of each row in the grid separately;
reducing the width and height of the unit grid, and performing grid division on the acquired signal grid point cloud section range based on the reduced width and height of the unit grid;
step four, according to the grid division result obtained in step four, step three, carry out the course of step two and step four repeatedly; and stopping until the width and the height of the reduced unit grid are smaller than the given minimum value, and taking the points contained in the signal grid obtained by the last iteration as the filtered point cloud.
Other steps and parameters are the same as those in the first embodiment.
The third concrete implementation mode: the present embodiment is different from the first or second embodiment in that the value of t is 3.
In this embodiment, the value of t is set to 3 in consideration of the requirement of continuity and fault tolerance of the signal points distributed in the elevation direction.
Other steps and parameters are the same as those in the first or second embodiment.
The fourth concrete implementation mode: the present embodiment is different from the first to third embodiments in that t' has a value of 1.
In the present embodiment, 1 grid adjacent to each other above and below the candidate grid having the highest continuity score can be retained in consideration of the steep slope topography and the truncation effect of the grid edge.
Other steps and parameters are the same as those in one of the first to third embodiments.
The fifth concrete implementation mode: the present embodiment is different from one of the first to fourth embodiments in that a specific method of reducing the width and height of the unit cell is:
width i’+1 ＝width i’ /r w (5)
height i’+1 ＝height i’ /r h (6)
in the formula, width i’ Unit grid width, height for the ith' time filtering i’ Unit grid height, r, for the ith' time filtering w Is a scaling factor of unit grid width, r h Is a scaling factor per unit grid height.
r w And r h Generally, 2 is adopted.
Other steps and parameters are the same as in one of the first to fourth embodiments.
The sixth specific implementation mode: the difference between this embodiment and one of the first to fifth embodiments is that the specific process of the fifth step is:
step five, point cloud trend reduction processing
Performing grid division on the filtered photon point cloud (the value range of the grid width and the grid height is generally within 15 meters) to form a grid with the size of w0 multiplied by h 0;
calculating median of elevation values of all points in a jth grid in the grid, and taking the median as a trend constant term of the grid; then, the trend constant term of the column is subtracted from the elevation value of each point in the column, so as to obtain the elevation value after trend reduction corresponding to each point in the column, as shown in fig. 4:
Figure BDA0003751417050000061
Figure BDA0003751417050000062
wherein, C _ detrend j The trend constant term for column j, grid j The set of elevation values for the points contained in the jth column in the w0 x h0 grid,
Figure BDA0003751417050000063
representing the elevation value of the ith "point in the jth column of the grid in the grid, I" =1,2, \ 8230where I, I is the total number of points in the jth column of the grid,
Figure BDA0003751417050000064
the elevation value after the ith' point falling trend in the jth column of grids is obtained;
in the same way, the elevation value of each point in the grid after the trend is reduced is respectively obtained;
fifthly, generating equal-interval elevation sampling points
In the direction of the distance along the track of the point cloud after the trend is reduced, performing equal-interval interpolation to form elevation sampling points according to the minimum value and the maximum value of the distance along the track of the point cloud after the trend is reduced, as shown in fig. 5;
Figure BDA0003751417050000065
Figure BDA0003751417050000066
in the formula (I), the compound is shown in the specification,
Figure BDA0003751417050000067
to maximize the distance of the point cloud along the track after trending down,
Figure BDA0003751417050000068
minimum distance of point cloud along the track after trend reduction, d sample Is the distance between two adjacent elevation sampling points, N sample For the total number of elevation sampling points,
Figure BDA0003751417050000069
distance along the track for the kth elevation sampling point;
in order to make the extracted elevation curve conform as closely as possible to the actual terrain, d sample The value of (A) should not exceed 10 m.
Step five and three, acquiring adjacent points
As shown in fig. 6a, for the kth elevation sampling point, according to a given neighborhood distance, a set of neighboring points of the kth elevation sampling point is selected:
Figure BDA00037514170500000610
wherein the content of the first and second substances,
Figure BDA0003751417050000071
a set of adjacent points for the kth elevation sampling point,
Figure BDA0003751417050000072
representation collection
Figure BDA0003751417050000073
The nth point in the inner part of the cell,
Figure BDA0003751417050000074
being points within a set of adjacent points
Figure BDA0003751417050000075
Along the rail distance, d ε Is the neighborhood distance;
step five four, counting the frequency distribution of the adjacent point set
According to collections
Figure BDA0003751417050000076
The maximum value and the minimum value of the point cloud elevation after the internal descending trend are set, the group distance is set, and a set of adjacent points is counted
Figure BDA0003751417050000077
A frequency distribution histogram of (a);
fifthly, calculating a point cloud classification threshold
Figure BDA0003751417050000078
Wherein bin j′ The frequency number of j' th interval in the frequency distribution histogram, N bin Is the number of intervals in the frequency distribution histogram, p is the threshold coefficient,
Figure BDA0003751417050000079
a point cloud classification threshold representing a kth elevation sampling point;
step five and six, point cloud classification
In the frequency distribution histogram of the adjacent point set, according to the direction of change from large to small of the elevation value, taking the interval of which the first frequency value is greater than the classification threshold value as a ground object elevation interval, and taking the point cloud falling in the ground object elevation interval as a ground object point cloud, as shown in fig. 6 b;
fifthly, calculating the ground elevation of the elevation sampling point
Calculating the average value of the elevation values of the point clouds falling in the elevation interval of the ground object after the descending trend, and taking the calculated average value as the ground object elevation value of the elevation sampling point;
Figure BDA00037514170500000710
Figure BDA00037514170500000711
in the formula (I), the compound is shown in the specification,
Figure BDA00037514170500000712
is the ground object elevation value, bin, of the kth elevation sampling point surface Is a ground object elevation interval, N surface Is a ground object elevation areaThe number of point clouds in the space between,
Figure BDA00037514170500000713
is the height value after the descending trend of the first point cloud in the height interval of the ground object, C _ detrend l A trend constant item corresponding to the first point cloud;
and fifthly, extracting the ground feature elevation value of each elevation sampling point by adopting the method from the fifth step three to the fifth step seven.
Other steps and parameters are the same as those in one of the first to fifth embodiments.
The seventh concrete implementation mode: this embodiment will be described with reference to fig. 7. The difference between this embodiment and one of the first to sixth embodiments is that the specific process of the step six is as follows:
sixthly, dividing the elevation sampling point cloud by a fixed distance length along the track to obtain each section of elevation sampling point cloud;
sixthly, setting a parameter dx representing the coordinate offset in the x direction and a parameter dy representing the coordinate offset in the y direction on a point cloud xoy plane (a plane rectangular coordinate system) for any section of elevation sampling point cloud, and searching a range L according to the x direction dx Y-direction search range L dy And x-direction search step size alpha dx Y-direction search step length α dy Generating all combinations Θ { (dx) of all values of the parameters dx and dy 1 ,dy 1 ),9dx 2 ,dy 2 ),…,(dx j″ ,dy j″ ),…}；
Figure BDA0003751417050000081
Figure BDA0003751417050000082
Sixthly, interpolation of elevation values
According to coordinate offset (dx) 1 ,dy 1 ) After the elevation sampling point cloud is integrally translated, the elevation sampling point cloud is translatedLater elevation sampling point p a Plane coordinates of
Figure BDA0003751417050000083
As a circle center, a photogrammetric point (q) within the radius r of the circle center is to be positioned 1 ,q 2 ,…,q m′ ) Reusing neighborhood points (q) as neighborhood points 1 ,q 2 ,…,q m′ ) Performing inverse distance interpolation to obtain interpolation point
Figure BDA0003751417050000084
The formula for interpolation is as follows:
Figure BDA0003751417050000085
Figure BDA0003751417050000086
Figure BDA0003751417050000087
in the formula (I), the compound is shown in the specification,
Figure BDA0003751417050000088
for the translated elevation sampling point p a The plane coordinates of the optical axis of the optical fiber,
Figure BDA0003751417050000089
for photogrammetric measuring points q b Is determined by the three-dimensional coordinates of (a),
Figure BDA00037514170500000810
is p a 、q b The planar distance between the two points is,
Figure BDA00037514170500000811
is a set of neighborhood points, w b Is the weight of the inverse distance interpolation,
Figure BDA00037514170500000812
interpolated values for elevation;
sixthly, calculating the elevation similarity, as shown in fig. 8:
Figure BDA00037514170500000813
Figure BDA00037514170500000814
Figure BDA00037514170500000815
in the formula (I), the compound is shown in the specification,
Figure BDA00037514170500000816
for the translated elevation sampling point p a N' is the number of sampling points in the section of elevation sampling point cloud, and rho is elevation similarity;
sixthly, repeating the process from the sixth step to the sixth step, traversing the coordinate offset to each group of coordinate offsets in the combined theta, and taking the maximum elevation similarity rho max Corresponding coordinate offset parameter combination as the optimal estimation value under the current search range and search step length
Figure BDA0003751417050000091
Sixthly, continuously adjusting the search range and the search step length of the elevation sampling point cloud according to the scale factor, repeating the processes from the sixth step to the sixth step after each adjustment, and obtaining the corresponding optimal estimation value and the corresponding maximum elevation similarity under each group of search range and search step length;
up to
Figure BDA0003751417050000092
Terminating the iteration when or up to a set maximum number of iterations, and endingThe optimal estimation value obtained by one iteration is used as the final x-direction coordinate offset and the y-direction coordinate offset;
wherein the content of the first and second substances,
Figure BDA0003751417050000093
for the maximum elevation similarity obtained for the ith iteration,
Figure BDA0003751417050000094
maximum elevation similarity obtained for the (i + 1) th iteration;
sixthly, translating the elevation sampling point cloud according to the final x-direction coordinate offset and the y-direction coordinate offset;
sixthly, in a similar way, translating each section of elevation sampling point cloud by adopting the method from the sixth step to the sixth step;
sixthly, randomly sampling a plurality of sections of elevation sampling point clouds from the elevation sampling point clouds after each section of translation, and registering elevation sampling points in the plurality of sections of sampled elevation sampling point clouds and photogrammetric point clouds by utilizing an ICP (inductively coupled plasma) algorithm to obtain a transformation matrix with the minimum error distance between the elevation sampling points and the photogrammetric point clouds;
and after rigid transformation is carried out on all elevation sampling points by using the transformation matrix, the elevation sampling points with error distances not meeting the threshold requirement are removed, and the remaining elevation sampling points are used as control points.
Other steps and parameters are the same as those in one of the first to sixth embodiments.
The specific implementation mode is eight: the difference between this embodiment and one of the first to seventh embodiments is that the search range and the search step length are adjusted according to the scaling factor, and the specific process is as follows:
Figure BDA0003751417050000095
Figure BDA0003751417050000096
Figure BDA0003751417050000097
Figure BDA0003751417050000098
in the formula (I), the compound is shown in the specification,
Figure BDA0003751417050000099
the values of the scale factors of the search range are all 2,
Figure BDA00037514170500000910
the scale factor for the search step length is generally 0.5,
Figure BDA00037514170500000911
for the x-direction search range of the ith iteration,
Figure BDA00037514170500000912
the range is searched for in the y-direction for the ith iteration,
Figure BDA00037514170500000913
for the x-direction search step of the ith iteration,
Figure BDA00037514170500000914
the step size is searched for the y-direction of the ith iteration,
Figure BDA00037514170500000915
the x-direction search range for the (i + 1) th iteration,
Figure BDA0003751417050000101
the y-direction search range for the (i + 1) th iteration,
Figure BDA0003751417050000102
the step size is searched for in the x direction for the (i + 1) th iteration,
Figure BDA0003751417050000103
the step size is searched for the y direction of the (i + 1) th iteration.
Other steps and parameters are the same as those in one of the first to seventh embodiments.
The specific implementation method nine: the difference between this embodiment and the first to eighth embodiments is that the specific process of the seventh step is:
step seven, establishing a fitting model of the elevation system error
Set of observations for control points
Figure BDA0003751417050000104
Carrying out barycenter transformation on the control points:
Figure BDA0003751417050000105
Figure BDA0003751417050000106
Figure BDA0003751417050000107
in the formula (I), the compound is shown in the specification,
Figure BDA0003751417050000108
for the observed values of the set of centrolized control points,
Figure BDA0003751417050000109
for known values of the centrobaric control point set,
Figure BDA00037514170500001010
is a known value of the set of control points, (x) q ,y q ,z q ) For photogrammetry the coordinates of the point cloud,
Figure BDA00037514170500001011
the coordinate of the C-th photogrammetric point in the photogrammetric point cloud is shown, and C is the number of photogrammetric points contained in the photogrammetric point cloud;
taking the z direction as an example, a second order polynomial fit model for the elevation system error:
Figure BDA00037514170500001012
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA00037514170500001013
for the observed value of the coordinates of the d-th control point in the barycentric control point set,
Figure BDA00037514170500001014
Figure BDA00037514170500001015
and
Figure BDA00037514170500001016
is the coefficient of a quadratic polynomial,
Figure BDA00037514170500001017
for a known value of the coordinates of the d-th control point in the set of centrolized control points, v d Is an intermediate variable;
let equation (29) be written in matrix form:
Figure BDA00037514170500001018
Figure BDA00037514170500001019
Figure BDA00037514170500001020
wherein D is the number of control points contained in the control point set;
according to the least square principle, the adjustment solution result of equation (30) is:
Figure BDA0003751417050000111
seventhly, similarly, obtaining quadratic polynomial coefficients in the x direction and the y direction by adopting the method in the step seventeenth;
seventhly, substituting the coordinates of each photogrammetric point into a quadratic polynomial fitting model of the elevation system error in the x direction, the y direction and the z direction and a planar quadratic polynomial fitting model to obtain a system error estimation value of each photogrammetric point in the x direction, the y direction and the z direction;
and correcting the coordinates of each photogrammetric point according to the system error estimation value.
Zeroth, first and second order polynomials are commonly used fitting models, and their specific form is as follows:
Figure BDA0003751417050000112
Figure BDA0003751417050000113
Figure BDA0003751417050000114
in the formula, a i Is a coefficient of a polynomial;
observed value for control point
Figure BDA0003751417050000115
Assuming its plane coordinates
Figure BDA0003751417050000116
Is related to the systematic error (Δ x, Δ y, Δ z) in generalThe functional relationship between the two is as follows:
Figure BDA0003751417050000117
Figure BDA0003751417050000118
Figure BDA0003751417050000119
in the formula, epsilon is a random error;
correcting the coordinate (x) of photogrammetric point cloud according to the estimated value q +Δx,y q +Δy,z q +Δz)。
Other steps and parameters are the same as those in one to eight of the embodiments.
Examples
In order to verify the effectiveness of the method, the invention adopts the resource No. three image and the high-resolution No. seven image to carry out experiments, and the experimental data parameters of the satellite remote sensing image are shown as the following table 1:
TABLE 1 satellite remote sensing image experiment data parameters
Figure BDA00037514170500001110
Figure BDA0003751417050000121
As shown in table 2, the data is obtained by performing ICESat-2 point cloud filtering and elevation extraction on the resource number three image:
TABLE 2
Figure BDA0003751417050000122
As shown in table 3, the data of the precision of the ICESat-2 point cloud filtering and elevation extraction for the top-scoring seven images are:
TABLE 3
Figure BDA0003751417050000123
Figure BDA0003751417050000131
As can be seen from tables 2 and 3, the filtering success rate of the ICESat-2 point cloud of the resource No. three image is 98.7% at most and 94.8% at least, and the filtering success rate of the ICESat-2 point cloud of the high-resolution No. seven image is 99.11% at most and 95.73% at most, which indicates that the point cloud filtering algorithm can obtain good filtering effect on photon data of different terrain variations, terrain types and laser intensities. The RMSE value of the ground object elevation extracted from the ICESat-2 point cloud of the resource No. three image is 5.36m at most, 3.04m at least, the MAD value is 4.20m at most, 2.23m at least, the RMSE value of the ground object elevation extracted from the ICESat-2 point cloud of the high-resolution No. seven image is 4.06m at most, 2.58m at most, the MAD value is 3.06m at most, and 1.64m at least, which proves that the photon point cloud elevation extraction algorithm has effectiveness, and the extracted ground object point elevation can be used as an elevation control point for satellite photogrammetry.
As shown in table 4, the positioning accuracy data of the resource No. three photogrammetric point clouds is:
TABLE 4
Figure BDA0003751417050000132
As shown in table 5, the positioning accuracy data of the high-resolution seven- # photogrammetry point cloud is:
TABLE 5
Figure BDA0003751417050000133
Figure BDA0003751417050000141
As can be seen from tables 4 and 5, the precision of the resource No. three photogrammetric point cloud in the X direction is improved from 31.49m to 0.46m, the precision of the resource No. three photogrammetric point cloud in the Y direction is improved from 10.83m to 2.43m, the precision of the height Z direction is improved from 7.28m to 2.11m, the precision of the height No. seven photogrammetric point cloud in the X direction is improved from 2.95m to 1.34m, the precision of the Y direction is not obviously improved, and the precision of the height Z direction is improved from 12.34m to 0.59m. Fig. 9a shows the resource point cloud obtained by photogrammetry No. three and the ICESat-2 point cloud, fig. 9b shows the resource point cloud obtained by photogrammetry No. three before the system error correction, and fig. 9c shows the resource point cloud obtained by photogrammetry No. three after the system error correction. Fig. 10a shows a schematic diagram of the seven-high point cloud and the ICESat-2 point cloud, fig. 10b shows a schematic diagram of the seven-high point cloud before system error correction, and fig. 10c shows a schematic diagram of the seven-high point cloud after system error correction. Therefore, the photon point cloud assisted satellite photogrammetry ground positioning technology has a remarkable effect on improving the point cloud positioning accuracy of photogrammetry in mountain areas.
The above-described calculation examples of the present invention are merely to explain the calculation model and the calculation flow of the present invention in detail, and are not intended to limit the embodiments of the present invention. It will be apparent to those skilled in the art that other variations and modifications of the present invention can be made based on the above description, and it is not intended to be exhaustive or to limit the invention to the precise form disclosed, and all such modifications and variations are possible and contemplated as falling within the scope of the invention.

Claims (9)

1. A satellite photogrammetry geostational positioning method based on photon point cloud assistance is characterized by comprising the following steps:
step one, acquiring stereoscopic image data of a target area;
step two, generating a photogrammetric point cloud of a ground object by utilizing the stereo image data;
acquiring an ICESat-2 photon point cloud in a stereopair region;
step four, performing grid filtering on the photon point cloud obtained in the step three, and filtering noise points to obtain a filtered photon point cloud;
generating elevation sampling points at equal intervals according to the filtered photon point cloud, and respectively calculating the ground elevation value of each elevation sampling point;
registering the elevation sampling point cloud and the photogrammetric point cloud generated in the second step according to the ground elevation value of each elevation sampling point to obtain a transformation matrix; then obtaining elevation sampling points serving as control points according to the transformation matrix;
and step seven, correcting the system error of the photogrammetric point cloud according to the obtained control point.
2. The method for positioning the earth through satellite photogrammetry based on photon point cloud assistance as claimed in claim 1, wherein the specific process of the fourth step is as follows:
step four, grid division
Carrying out grid division on the whole photon point cloud profile range to form a grid with the size of w multiplied by h;
step two, counting the number of points contained in each grid in the grid respectively, and for any grid in a column, selecting the first t grids containing the most points from all grids in the column as candidate grids of the column:
Figure FDA0003751417040000011
wherein, G j Representing the extracted candidate grid in column j,
Figure FDA0003751417040000012
is the number of points contained within the mth grid in column j, m is the total number of grids in that column,
Figure FDA0003751417040000013
representing that t grids containing the most number of points are selected from the jth column;
in the same way, respectively obtaining candidate grids of each column;
step four and three, detecting the continuity of the grid
Step three, taking a certain candidate grid in the j column as a grid to be detected, taking the candidate grid in the j-1 column, the j-2 column, \8230, the candidate grid in the j-k 'column, the j +1 column, the j +2 column, \8230, the candidate grid in the j + k' column as the neighborhood range of the grid to be detected, forming a left connected domain by the grid to be detected and the candidate grid in the j-1 column, the j-2 column, \8230, the j-k 'column, forming a right connected domain by the grid to be detected and the candidate grid in the j +1 column, the j +2 column, \8230andthe j + k' column;
Figure FDA0003751417040000014
Figure FDA0003751417040000021
in the formula (I), the compound is shown in the specification,
Figure FDA0003751417040000022
is the set of grids within the left connected domain,
Figure FDA0003751417040000023
is a grid set in a right connected domain, (i, j) is a row and column serial number of a grid to be detected in the grid,
Figure FDA0003751417040000024
for the m1 st grid in the set of grids in the left connected domain,
Figure FDA0003751417040000025
the n1 grid in the grid set in the right connected domain is obtained;
step four, step three, step two, calculate continuity score Continue _ score of the grid to be detected:
Figure FDA0003751417040000026
in the formula (I), the compound is shown in the specification,
Figure FDA0003751417040000027
is the number of grids in the left connected domain,
Figure FDA0003751417040000028
the number of grids in the right connected domain;
step four, calculating the continuity score of each candidate grid in the jth column of the grid in a traversing manner, and taking the candidate grid with the highest continuity score and t' grids adjacent to the upper part and the lower part of the candidate grid with the highest continuity score as signal grids in the jth column;
step four, step three, repeat the process from step four three one to step four three, obtain the signal grid of each row in the grid separately;
reducing the width and height of the unit grid, and performing grid division on the acquired signal grid point cloud section range based on the reduced width and height of the unit grid;
step four, according to the grid division result obtained in step four, step three, carry out the course of step two and step four repeatedly; and stopping until the width and the height of the reduced unit grid are smaller than the given minimum value, and taking the point contained in the signal grid obtained by the last iteration as the filtered point cloud.
3. The method for satellite photogrammetry-based satellite positioning based on photon point cloud assistance as claimed in claim 2, wherein the value of t is 3.
4. The method for positioning the earth through satellite photogrammetry based on photon point cloud assistance as recited in claim 3, wherein the value of t' is 1.
5. The method for positioning the earth through satellite photogrammetry based on photon point cloud assistance as claimed in claim 4, wherein the specific method for reducing the width and height of the unit grid is as follows:
width i’+1 ＝width i’ /r w (5)
height i’+1 ＝height i’ /r h (6)
in the formula, width i’ Unit grid width, height for the ith' time filtering i’ Unit grid height, r, for the ith' time filtering w Is a scaling factor of unit grid width, r h Is a scaling factor per unit grid height.
6. The method for positioning the earth through satellite photogrammetry based on photon point cloud assistance as claimed in claim 5, wherein the concrete process of the fifth step is as follows:
step five, point cloud trend reduction processing
Performing grid division on the filtered photon point cloud to form a grid with the size of w0 multiplied by h 0;
calculating median of elevation values of all points in a jth grid in the grid, and taking the median as a trend constant term of the grid; and then respectively subtracting the trend constant term of the column from the elevation value of each point in the column to obtain the elevation value after trend reduction corresponding to each point in the column:
Figure FDA0003751417040000031
Figure FDA0003751417040000032
wherein, C _ detrend j In column jTrend constant term, grid j The set of elevation values for the points contained in the jth column in the w0 x h0 grid,
Figure FDA0003751417040000033
representing the elevation value of the ith "point in the jth column of the grid in the grid, I" =1,2, \ 8230where I, I is the total number of points in the jth column of the grid,
Figure FDA0003751417040000034
the elevation value after the ith' point falling trend in the jth grid is taken as the elevation value;
in the same way, the elevation value of each point in the grid after the trend is reduced is respectively obtained;
fifthly, generating equal-interval elevation sampling points
In the distance direction of the point cloud along the track after the trend is reduced, performing equal-interval interpolation to form elevation sampling points according to the minimum value and the maximum value of the distance of the point cloud along the track after the trend is reduced;
Figure FDA0003751417040000035
Figure FDA0003751417040000036
in the formula (I), the compound is shown in the specification,
Figure FDA0003751417040000037
to maximize the distance of the point cloud along the track after trending down,
Figure FDA0003751417040000038
minimum distance of point cloud along the track after trend reduction, d sample Is the distance between two adjacent elevation sampling points, N sample For the total number of elevation sampling points,
Figure FDA0003751417040000039
distance along the track for the kth elevation sampling point;
step five and three, acquiring adjacent points
For the kth elevation sampling point, selecting a set of adjacent points of the kth elevation sampling point according to a given neighborhood distance:
Figure FDA00037514170400000310
wherein the content of the first and second substances,
Figure FDA00037514170400000311
a set of adjacent points for the kth elevation sampling point,
Figure FDA00037514170400000312
representation collection
Figure FDA00037514170400000313
The nth point in the inner part of the cell,
Figure FDA00037514170400000314
for points within a set of adjacent points
Figure FDA00037514170400000315
Along the rail distance, d ε Is the neighborhood distance;
step five and four, counting the frequency distribution of the adjacent point set
According to collections
Figure FDA0003751417040000041
The maximum value and the minimum value of the point cloud elevation after the internal descending trend, the group distance is set, and the adjacent point set is counted
Figure FDA0003751417040000042
A frequency distribution histogram of (a);
fifthly, calculating a point cloud classification threshold
Figure FDA0003751417040000043
Wherein bin j′ The frequency of j' th interval in the frequency distribution histogram, N bin Is the number of intervals in the frequency distribution histogram, p is the threshold coefficient,
Figure FDA0003751417040000044
a point cloud classification threshold representing a kth elevation sampling point;
fifthly, point cloud classification
In the frequency distribution histogram of the adjacent point set, according to the change direction of the elevation value from large to small, taking the interval of which the first frequency value is greater than the classification threshold value as a ground object elevation interval, and taking the point cloud falling in the ground object elevation interval as a ground object point cloud;
fifthly, calculating the ground elevation of the elevation sampling point
Calculating the average value of the elevation values of the point clouds falling in the elevation interval of the ground object after the descending trend, and taking the calculated average value as the ground object elevation value of the elevation sampling point;
Figure FDA0003751417040000045
in the formula (I), the compound is shown in the specification,
Figure FDA0003751417040000046
is the ground object elevation value, bin, of the kth elevation sampling point surface Is a ground object elevation interval, N surface The number of point clouds in the elevation interval of the ground object,
Figure FDA0003751417040000047
is the height value after the descending trend of the first point cloud in the height interval of the ground object, C _ detrend l A trend constant item corresponding to the first point cloud;
and fifthly, extracting the ground feature elevation value of each elevation sampling point by adopting the method from the fifth step three to the fifth step seven.
7. The method for satellite photogrammetry to ground positioning based on photon point cloud assistance as claimed in claim 6, wherein the specific process of the sixth step is as follows:
sixthly, dividing the elevation sampling point cloud by using the fixed distance length along the track to obtain each section of elevation sampling point cloud;
sixthly, for any section of elevation sampling point cloud, respectively setting a parameter dx representing the coordinate offset in the x direction and a parameter dy representing the coordinate offset in the y direction on a point cloud xoy plane, and searching a range L according to the x direction dx Y-direction search range L dy And x-direction search step size alpha dx Y-direction search step length alpha dy Generating all combinations Θ { (dx) of all values of the parameters dx and dy 1 ,dy 1 ),(dx 2 ,dy 2 ),…,(dx j″ ,dy j″ ),…}；
Figure FDA0003751417040000051
Figure FDA0003751417040000052
Sixthly, interpolation of elevation values
According to coordinate offset (dx) 1 ,du 1 ) After the elevation sampling point cloud of the section is integrally translated, the translated elevation sampling point p is used a Plane coordinates of
Figure FDA0003751417040000053
As the center of circle, the photogrammetric point (q) within the radius r of the center of circle 1 ,q 2 ,…,q m′ ) Reusing neighborhood points (q) as neighborhood points 1 ,q 2 ,…,q m′ ) Performing inverse distance interpolation to obtain interpolation point
Figure FDA0003751417040000054
The formula for interpolation is as follows:
Figure FDA0003751417040000055
Figure FDA0003751417040000056
Figure FDA0003751417040000057
in the formula (I), the compound is shown in the specification,
Figure FDA0003751417040000058
for the translated elevation sampling point p a The plane coordinates of (a) are determined,
Figure FDA0003751417040000059
for photogrammetric measuring points q b Is determined by the three-dimensional coordinates of (a),
Figure FDA00037514170400000510
is p a 、q b The planar distance between the two points is,
Figure FDA00037514170400000511
is a set of neighborhood points, w b Is the weight of the inverse distance interpolation,
Figure FDA00037514170400000512
interpolated values for elevation;
sixthly, calculating the elevation similarity
Figure FDA00037514170400000513
Figure FDA00037514170400000514
Figure FDA00037514170400000515
In the formula (I), the compound is shown in the specification,
Figure FDA00037514170400000516
for the translated elevation sampling point p a N' is the number of sampling points in the section of elevation sampling point cloud, and rho is elevation similarity;
sixthly, repeating the process from the sixth step to the sixth step, traversing the coordinate offset to each group of coordinate offsets in the combined theta, and taking the maximum elevation similarity rho max Corresponding coordinate offset parameter combination as the optimal estimation value under the current search range and search step length
Figure FDA00037514170400000517
Sixthly, continuously adjusting the search range and the search step length of the elevation sampling point cloud according to the scale factor, repeating the processes from the sixth step to the sixth step after each adjustment, and obtaining the corresponding optimal estimation value and the corresponding maximum elevation similarity under each group of search range and search step length;
up to
Figure FDA00037514170400000518
Stopping iteration when the set maximum iteration times are reached, and taking the optimal estimation value obtained by the last iteration as the final x-direction coordinate offset and the final y-direction coordinate offset;
wherein the content of the first and second substances,
Figure FDA0003751417040000061
for the maximum elevation similarity obtained for the ith iteration,
Figure FDA0003751417040000062
maximum elevation similarity obtained for the (i + 1) th iteration;
sixthly, translating the elevation sampling point cloud according to the final x-direction coordinate offset and the y-direction coordinate offset;
sixthly, in a similar way, translating each section of elevation sampling point cloud by adopting the method from the sixth step to the sixth step;
sixthly, randomly sampling a plurality of sections of elevation sampling point clouds from the translated elevation sampling point clouds, and registering the elevation sampling points in the plurality of sections of the sampled elevation sampling point clouds and the photogrammetric point clouds by utilizing an ICP (inductively coupled plasma) algorithm to obtain a transformation matrix with the minimum error distance between the elevation sampling points and the photogrammetric point clouds;
and after rigid transformation is carried out on all elevation sampling points by using the transformation matrix, eliminating the elevation sampling points of which the error distances do not meet the threshold requirement, and taking the remaining elevation sampling points as control points.
8. The method for positioning the earth through satellite photogrammetry based on photon point cloud assistance as claimed in claim 7, wherein the scaling factor adjusts the search range and the search step length by the following specific processes:
Figure FDA0003751417040000063
Figure FDA0003751417040000064
Figure FDA0003751417040000065
Figure FDA0003751417040000066
in the formula (I), the compound is shown in the specification,
Figure FDA0003751417040000067
for the scale factor of the search range,
Figure FDA0003751417040000068
to search for the scaling factor of the step size,
Figure FDA0003751417040000069
the x-direction search range for the ith iteration,
Figure FDA00037514170400000610
the range is searched for in the y-direction for the ith iteration,
Figure FDA00037514170400000611
for the x-direction search step of the ith iteration,
Figure FDA00037514170400000612
the step size is searched for the y-direction of the ith iteration,
Figure FDA00037514170400000613
the x-direction search range for the (i + 1) th iteration,
Figure FDA00037514170400000614
the range is searched for in the y direction for the (i + 1) th iteration,
Figure FDA00037514170400000615
the step size is searched for in the x direction for the (i + 1) th iteration,
Figure FDA00037514170400000616
the step size is searched for the y direction of the (i + 1) th iteration.
9. The method for positioning the earth through satellite photogrammetry based on photon point cloud assistance as claimed in claim 8, wherein the concrete process of the seventh step is as follows:
step seven, establishing a fitting model of the elevation system error
Set of observations for control points
Figure FDA00037514170400000617
And (3) carrying out center of gravity on the control points:
Figure FDA0003751417040000071
Figure FDA0003751417040000072
Figure FDA0003751417040000073
in the formula (I), the compound is shown in the specification,
Figure FDA0003751417040000074
as an observation of the set of centrobaric control points,
Figure FDA0003751417040000075
for known values of the centrobaric control point set,
Figure FDA0003751417040000076
is a known value of the set of control points, (x) q ,y q ,z q ) For photogrammetry the coordinates of the point cloud,
Figure FDA0003751417040000077
the coordinate of the C-th photogrammetric point in the photogrammetric point cloud is shown, and C is the number of photogrammetric points contained in the photogrammetric point cloud;
taking the z direction as an example, fitting a model to a quadratic polynomial of the elevation system error:
Figure FDA0003751417040000078
wherein, the first and the second end of the pipe are connected with each other,
Figure FDA0003751417040000079
for the observed value of the coordinates of the d-th control point in the set of barycentric control points,
Figure FDA00037514170400000710
and
Figure FDA00037514170400000711
is a coefficient of a quadratic polynomial and is,
Figure FDA00037514170400000712
for a known value of the coordinates of the d-th control point in the set of centrolized control points, v d Is an intermediate variable;
let equation (29) be written in matrix form:
Figure FDA00037514170400000713
Figure FDA00037514170400000714
Figure FDA00037514170400000715
wherein D is the number of control points contained in the control point set;
according to the least square principle, the adjustment solution result of equation (30) is:
Figure FDA00037514170400000716
seventhly, similarly, obtaining quadratic polynomial coefficients in the x direction and the y direction by adopting the method in the step seventeenth;
seventhly, substituting the coordinates of each photogrammetric point into a quadratic polynomial fitting model of the elevation system error in the x direction, the y direction and the z direction to obtain a system error estimated value of each photogrammetric point in the x direction, the y direction and the z direction;
and correcting the coordinates of each photogrammetric point according to the system error estimation value.
CN202210843831.XA 2022-07-18 2022-07-18 Satellite photogrammetry earth positioning method based on photon point cloud assistance Active CN115143942B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210843831.XA CN115143942B (en) 2022-07-18 2022-07-18 Satellite photogrammetry earth positioning method based on photon point cloud assistance

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210843831.XA CN115143942B (en) 2022-07-18 2022-07-18 Satellite photogrammetry earth positioning method based on photon point cloud assistance

Publications (2)

Publication Number Publication Date
CN115143942A true CN115143942A (en) 2022-10-04
CN115143942B CN115143942B (en) 2023-07-28

Family

ID=83411933

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210843831.XA Active CN115143942B (en) 2022-07-18 2022-07-18 Satellite photogrammetry earth positioning method based on photon point cloud assistance

Country Status (1)

Country Link
CN (1) CN115143942B (en)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102645209A (en) * 2012-04-24 2012-08-22 长江勘测规划设计研究有限责任公司 Joint positioning method for spatial points by means of onboard LiDAR point cloud and high resolution images
CN103017739A (en) * 2012-11-20 2013-04-03 武汉大学 Manufacturing method of true digital ortho map (TDOM) based on light detection and ranging (LiDAR) point cloud and aerial image
CN104123730A (en) * 2014-07-31 2014-10-29 武汉大学 Method and system for remote-sensing image and laser point cloud registration based on road features
CN112529946A (en) * 2020-12-04 2021-03-19 中南大学 High discrete body model optimization method and system based on elevation data, electronic equipment and readable storage medium

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102645209A (en) * 2012-04-24 2012-08-22 长江勘测规划设计研究有限责任公司 Joint positioning method for spatial points by means of onboard LiDAR point cloud and high resolution images
CN103017739A (en) * 2012-11-20 2013-04-03 武汉大学 Manufacturing method of true digital ortho map (TDOM) based on light detection and ranging (LiDAR) point cloud and aerial image
CN104123730A (en) * 2014-07-31 2014-10-29 武汉大学 Method and system for remote-sensing image and laser point cloud registration based on road features
CN112529946A (en) * 2020-12-04 2021-03-19 中南大学 High discrete body model optimization method and system based on elevation data, electronic equipment and readable storage medium

Also Published As

Publication number Publication date
CN115143942B (en) 2023-07-28

Similar Documents

Publication Publication Date Title
US20200103530A1 (en) Method for extracting elevation control point with assistance of satellite laser altimetry data
CN105866762B (en) Laser radar automatic calibrating method and device
CN104111456B (en) A kind of line of high-speed railway Ground Deformation high-resolution InSAR monitoring methods
CN103673995A (en) Calibration method of on-orbit optical distortion parameters of linear array push-broom camera
CN110889899B (en) Digital earth surface model generation method and device
CN109100719A (en) Combine plotting method with the topographic map of optical image based on satellite-borne SAR image
CN112068136A (en) Azimuth deformation monitoring method based on amplitude offset
CN113393571B (en) Cloud-free satellite image generation method and device
CN113358091A (en) Method for producing digital elevation model by using three-linear array three-dimensional satellite image
CN107393004A (en) A kind of method and device for obtaining building amount of demolition in power transmission line corridor
CN108562900B (en) SAR image geometric registration method based on elevation correction
CN108444451B (en) Planet surface image matching method and device
CN113238228B (en) Three-dimensional earth surface deformation obtaining method, system and device based on level constraint
CN103065295A (en) Aviation and ground lidar data high-precision automatic registering method based on building angular point self-correction
CN114821522A (en) Urban road cross slope and super height value calculation method based on vehicle-mounted laser point cloud data
Huang et al. Multi-view large-scale bundle adjustment method for high-resolution satellite images
CN115143942A (en) Satellite photogrammetry earth positioning method based on photon point cloud assistance
Gao et al. Automatic geo-referencing mobile laser scanning data to UAV images
CN117092621A (en) Hyperspectral image-point cloud three-dimensional registration method based on ray tracing correction
CN113344866B (en) Point cloud comprehensive precision evaluation method
CN114964176A (en) Landform surveying and mapping method for permanent shadow area of moon
CN115830476A (en) Terrain factor space downscaling method
CN115546264A (en) Satellite-borne InSAR image fine registration and stereo measurement method
CN113160374A (en) Three-dimensional calculation method for volume change of gully based on terrain point cloud
CN111308469B (en) Building elevation measurement method based on PSInSAR technology

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