CN104913780A - GNSS-CCD-integrated zenith telescope high-precision vertical deflection fast measurement method - Google Patents

Abstract

The invention discloses a GNSS-CCD-integrated zenith telescope high-precision vertical deflection fast measurement method. At a measurement station, the optical axis of a zenith telescope faces to the zenith and takes a picture of a zenith zone so that a CCD fixed star image is obtained, fixed star image coordinates are calculated, exposure epoch is controlled by geodetic coordinates and a GNSS time signal measured by GNSS, fixed star information in the zenith zone is acquired in an appropriate star catalogue according to the measurement station vertical deflection calculated by an EGM2008 geopotential model, coupling identification of fixed stars in a zenith tangent plane zone of a celestial sphere in the star catalogue and in the CCD fixed star image is realized, zenith astronomical coordinates are calculated by least spuare iterative computation, and vertical deflection of the observed point is calculated according to geodetic coordinates. Compared with the prior art, the method provided by the invention has the characteristics of simple operation, time and labor saving and high measurement precision and is suitable for high-precision vertical deflection fast measurement.

Description

The high precision plumb line deviation method for fast measuring of integrated GNSS and CCD zenith telescope

Technical field

The present invention relates to a kind of astronomical geodesy method of plumb line deviation, particularly relate to a kind of high precision plumb line deviation method for fast measuring of integrated GNSS and CCD zenith telescope.

Background technology

Angle between normal line vector n on ground on any gravity vector g and corresponding ellipsoid is defined as the plumb line deviation of this point.Plumb line deviation characterizes the direction of gravity, is that geodetic surveying reduction, earth gravity field and geoid's model are refined, study the necessary earth physical quantity of the figure of the earth, in geodetic surveying and space technology etc., has important use.

Plumb line deviation is the performance of earth interior mass distribution, and can be used for inverting air abnormal refraction research, monitoring earth interior mass transfer and energy accumulation, the disaster such as monitor earthquake, volcanic explosion, also can be used for Underground abnormal quality, carries out resource exploration.In military field, plumb line deviation is the required basic data such as guided missile, rocket launching.

In prior art, conventional plumb line deviation observation procedure, mainly comprises traditional astronomical geodesy method, Gravimetric Method, astrogravimetric method and GPS leveling measuring method.But these method ubiquity basic data acquisition processes above-mentioned are complicated, or just cannot effectively obtain at all, late time data processing procedure very complicated, shortcoming or the deficiency such as waste time and energy:

Tradition astronomical geodesy method, adopts the astronomical theodolites such as T4, J05 usually, is equipped with a set of punctual, time dissemination system.Telgte method is commonly used in the observation of its astronomical latitude, and the observation of astronomical longitude is commonly used double star equal altitude method and includes radio time signal simultaneously.The plumb line deviation precision that this method obtains can reach 0.3 ".

But for this traditional astronomical geodesy method, usually, it is in order to make accuracy of observation enough high, the volume of the astronomical theodolite used all very greatly, very heavy; And the operation of instrument is too high to observer's skill requirement, observation time is very long, efficiency is extremely low.And its conventional punctual, time dissemination system, not only complicated operation, precision is lower, and easily makes mistakes.

Gravimetric Method, mainly by means of the gravity anomaly on geoid surface and earth ellipsoid face, utilizes Stokes or VeningMeinesz formulae discovery plumb line deviation.

But the prerequisite of this method does not have disturbance material outside supposition geoid surface, and its prerequisite is, needs global gravity anomaly data known, and this is very difficult and not real.Therefore, the method is not also independently applied, and the measuring accuracy of the plumb line deviation of the method is also on the low side.

Astrogravimetric method, its essence is in conjunction with astronomical geodesy method and Gravimetric Method to determine plumb line deviation.First the method needs the plumb line deviation value of the astrogeodetic point of some, then utilizes Gravimetric Method to encrypt gravity anomaly data, adopts least square collocation or Multi-surface fitting method to determine the plumb line deviation value of other points.

The major defect of this method is, the plumb line deviation precision obtained is on the low side.

GPS leveling measuring method, needs to carry out GPS measurement and measurement of the level simultaneously, obtains the difference of base length, geodetic azimuth and height anomaly, just can obtain plumb line deviation.

But the use of this method is had ready conditions restriction, such as, the difficulty of the region that some topographic relief is larger being carried out to precise leveling is too large, baseline can not be long, and due to the impact of base length and geodetic azimuth, the plumb line deviation precision causing it to obtain is on the low side.

Summary of the invention

The object of the invention is, provide one to have pervasive characteristic, and the data acquisition of high precision plumb line deviation is time saving and energy saving, the high precision plumb line deviation method for fast measuring of relatively simple a kind of integrated GNSS and the CCD zenith telescope of the disposal route of data.

The technical scheme that the present invention is adopted for achieving the above object is, a kind of high precision plumb line deviation method for fast measuring of integrated GNSS and CCD zenith telescope, is characterized in that, comprise the following steps:

The first step, survey station point sets up CCD zenith telescope, its optical axis is made to point to zenith direction, zenith region is taken at 0 ° and 180 ° of directions continuously, obtain CCD star chart picture, and the tilt data of zenith telescope is obtained by high-precision electronic level meter, the terrestrial coordinate of observation time information and survey station point is obtained by GNSS;

Second step, reads the CCD star chart picture of the FITS form of CCD zenith telescope shooting, then carries out ground unrest elimination;

3rd step, the CCD star chart picture after eliminating noise carries out the automatic search in fixed star astrology region;

4th step, uses the correction Two-Dimensional Moment method improved, calculates the image coordinate of the center of energy in the fixed star astrology region that above-mentioned automatic search goes out;

5th step, processes applicable star catalogue, the terrestrial coordinate of the temporal information obtained by GNSS and survey station point, calculates the rough astronomic coordinates of survey station point in conjunction with EGM2008.To intercept as corresponding region with CCD star chart in star catalogue according to result, and the independent equatorial coordinate system in calculation exposure moment zenith region and apparent magnitude information, and independent equatorial coordinate system is projected on section be converted to planimetric coordinates;

6th step, the fixed star projected in plane the astrology in CCD image and star catalogue identifies, mates;

7th step, according to the star pair that the match is successful, utilizes the coordinate conversion coefficient in least square method computational photogrammetry projective transformation formula, and the initial equatorial coordinate of iterative computation survey station point, then be converted to the initial astronomic coordinates of survey station point;

8th step, the zenith telescope inclination correction data obtain high-precision electronic level meter and Ghandler motion correct data and are converted on meridian direction and direction at the tenth of the twelve Earthly Branches, the fourth of the twelve Earthly Branches, determine the accurate astronomic coordinates of survey station point;

9th step, in conjunction with the survey station point terrestrial coordinate that GNSS measurement data calculates, utilizes Helmert formula, can calculate the plumb line deviation of survey station point.

The technique effect that technique scheme is directly brought is, plumb line deviation method for fast measuring is based on GNSS and CCD zenith telescope, on the basis meeting measuring accuracy and resolution needs, utilize the data acquisition and disposal route that improve, reduce the time that plumb line deviation acquisition process consumes, improve measurement efficiency, high measuring accuracy can be met and without requirements such as region restrictions.

Be preferably, it is 3 × 3 median filtering methods that above-mentioned background noise eliminates the method adopted, and concrete steps are as follows:

Step one, determines 3 × 3 neighborhoods put centered by certain pixel;

Step 2, by the gray-scale value (x of each pixel in this neighborhood ₁, x ₂, x ₃..., x _n) by sorting from small to large:

x ₁≤x ₂≤x ₃...≤x _n；

Step 3, goes out the intermediate value Y of the gray-scale value of each pixel by formulae discovery:

Step 4, replaces the gray-scale value of original central point with intermediate value Y, carry out medium filtering;

Step 5, determines 3 × 3 neighborhoods put centered by one other pixel, repeats above-mentioned steps two to step 4, until all pixel filter completes.

The technique effect that this optimal technical scheme is directly brought is, before and after denoising, fixed star discrimination improves about 60%, adds observed quantity, largely reducing the impact of random noise for CCD fixed star image procossing.The standard deviation of plumb line deviation meridional component and fourth of the twelve Earthly Branches component at the tenth of the twelve Earthly Branches on average reduces 0.016 " and 0.021 " respectively before denoising, and inner symbol precision increases.

3 × 3 medium filterings are a kind of non-linear processing methods removing noise based on sequencing statistical theory, and ultimate principle is that the Mesophyticum of each point gray-scale value in 3 × 3 neighborhoods of this point of the gray-scale value in digital picture is replaced.

Here neighborhood is commonly called window, after window moves in the picture up and down, utilizes medium filtering can carry out denoising to image well.

Further preferably, the automatic search in above-mentioned fixed star astrology region adopts algorithm of region growing to carry out.

The technique effect that this optimal technical scheme is directly brought is, employing algorithm of region growing carries out the astrology automatic search in CCD star chart picture, and its method is simple, and the time of required cost is shorter.

Reason is, the ultimate principle of algorithm of region growing is setting minimal gray threshold value, if certain region is made up of some neighbors that gray-scale value is close in ccd image, and in this region, each grey scale pixel value is all greater than minimal gray threshold value, then judge that this region is as a fixed star astrology region.

Further preferably, the correction Two-Dimensional Moment method of above-mentioned improvement is in correction Two-Dimensional Moment method in the past, and use iterative calculation method to improve threshold value, its step is as follows:

1st step, extracts maximum gradation value g in CCD star chart picture _maxwith minimum gradation value g _min, order

$T^0=\frac{g_{\min }+g_{\max }}{2}---\left(2\right);$

2nd step, according to threshold value T ^kcCD star chart picture is divided into astrology region and background area, wherein, astrology regional location (x _i, y _j) gray-scale value g _k+1(x _i, y _j) > 0, and the gray-scale value of background area is 0, as shown by the equation:

Then, the average gray value in astrology region is tried to achieve by formula (4):

$g_{k+1}=\frac{\underset{g_k(x_i,y_j)>T^k}{\mathrm{\Σ}}(g_k(x_i,y_j)-T^k)}{\underset{g_k(x_i,y_j)>T^k}{\mathrm{\Σ}1}}---\left(4\right)$

Order formula (4) is substituted into, can obtain:

$T^{k+1}=\frac{\underset{g_k(x_i,y_j)>T^k}{\mathrm{\Σ}}(g_k(x_i,y_j)-T^k)}{\underset{g_k(x_i,y_j)>T^k}{\mathrm{\Σ}2}}---\left(5\right)$

3rd step, arranges a positive number μ, if | T ^k-T ^k+1| < μ, then termination of iterations computation process, try to achieve optimal threshold by formula (6):

T＝T ^k+B (6)

If | T ^k-T ^k+1|>=μ, then make T ^k=T ^k+1, proceed the iterative computation of the 2nd step.

4th step, by formula (7) and formula (8), obtain astrology regional center coordinate in CCD star chart picture:

$x_0=\frac{\mathrm{\&Sigma;\&Sigma;}{x}_i(g(x_i,y_j)-T)}{\mathrm{\&Sigma;\&Sigma;}(g(x_i,y_j)-T)}---\left(7\right)$

$y_0=\frac{\mathrm{\&Sigma;\&Sigma;}{y}_i(g(x_i,y_j)-T)}{\mathrm{\&Sigma;\&Sigma;}(g(x_i,y_j)-T)}---\left(8\right).$

The technique effect that this optimal technical scheme is directly brought is, on the basis of the square of two dimension correction in the past algorithm, threshold value T is improved by process of iteration, astrology centralized positioning precision is greatly improved, positioning precision standard deviation controls, in 0.1-0.15 pixel, to control 0.1 the impact of plumb line deviation measuring accuracy " in.Further preferably, the above-mentioned process being suitable for star catalogue is the NOVAS-F program by calling in AURIGA, and in conjunction with the temporal information that GNSS obtains, to the star place observed epoch of observation and be suitable for each influent factor that star catalogue fixed star epoch mean place there are differences and correct, to obtain the apparent place of fixed star in star catalogue;

Then, intercept region corresponding with CCD fixed star star image in star catalogue, intercepting process is divided into two steps:

The first step, before integrated CCD and GPS numeral zenith telescope is observed, first can intercept by the fixed star information in survey station point zenith region this observation period in star catalogue, star catalogue intercepts scope so that the fixed star that may appear in visual field, zenith region during observation is all covered as principle, its survey station point equatorial coordinate (α ₂, δ ₂) tried to achieve by following formula (9):

δ ₂＝Φ ₂

α ₂＝Λ ₂+GAST (9)

In above formula (9), (Φ ₂, Λ ₂) be the initial astronomic coordinates of survey station point, this initial astronomic coordinates is the plumb line deviation adopting GNSS geodetic surveying result and EGM2008 gravity field model to calculate survey station point, then pushes away to obtain the astronomic coordinates of survey station point, gets that approximate value obtains;

Described star catalogue intercepts scope and is greater than digital zenith telescope visual field;

After intercepting scope is determined, use apparent position calculation of star model to intercept the fixed star information that may appear within the scope of this in the observation period, and these information are stored;

Second step, before carrying out the work of fixed star coupling, need intercept above star catalogue region further, and enable this region completely corresponding with CCD fixed star image-region, it is as follows that it intercepts condition:

α _min＜α ₁＜α _max

δ _min＜δ ₁＜δ _max(10)

In above formula (10):

(α ₁, δ ₁) be the equatorial coordinate of fixed star in visual field, zenith region;

α _min, α _max, δ _min, δ _maxdepend on the equatorial coordinate (α of survey station point ₂, δ ₂) and the visual field Fi of digital zenith telescope, its computing formula is as follows:

α _min＝α ₂-Fi/2

α _max＝α ₂+Fi/2 (11)

δ _min＝δ ₂-Fi/2

δ _max＝δ ₂+Fi/2

Described being projected to by independent equatorial coordinate system plane is converted to planimetric coordinates calculates by projective transformation formula (12) and (13):

$l=\frac{\tan ({\mathrm{\&alpha;}}_1-{\mathrm{\&alpha;}}_2)}{\cos (q-{\mathrm{\&delta;}}_2)}---\left(12\right)$

m＝tan(q-δ ₂) (13)

In above formula (12) and (13):

cotq＝cotδ ₁cos(α ₁-α ₂)；

(l, m) is fixed star planimetric coordinates;

(α ₁, δ ₁) be the equatorial coordinate of fixed star in visual field, zenith region _.

Illustrate: above-mentioned " star place that epoch of observation observes and be suitable for each influent factor that star catalogue fixed star epoch mean place there are differences " refers in the processing procedure being suitable for star catalogue, based on the fixed star in the zenith region that the zenith telescope be positioned on earth surface survey station point observes, be subject to stellar proper motion, the precession of the equinoxes, nutating, annual parallax, annual aberration, diurnal parallax, diurnal aberration, the impact of the factors such as astronomical refraction, the star place observed epoch of observation and applicable star catalogue fixed star epoch mean place there are differences, the present invention calls the NOVAS-F program in AURIGA, the temporal information obtained in conjunction with GNSS and the terrestrial coordinate of survey station point, above influent factor is corrected, obtain independent equatorial coordinate system and the apparent magnitude information in zenith time of exposure region.

The technique effect that this optimal technical scheme is directly brought is, 1, new star catalogue intercept method is divided into two steps, before observation and in observation, twice pair of star catalogue intercepts respectively, the technique effect that this optimal technical scheme is directly brought is, namely realize intercepting to the fixed star that may appear at zenith region before observation, only need intercept out in star catalogue in the first step when then processing in real time observation data and intercept, make fixed star identification be matched to power and significantly improve, the whole identification matching process used time is less than 0.3 second;

2, star catalogue intercepting scope must be greater than digital zenith telescope visual field, can guarantee that the fixed star in digital zenith telescope visual field is unlikely to occur any omission; Usually, above-mentioned star catalogue can be intercepted 2 times that scope is chosen as digital zenith telescope visual field, like this, be convenient to practical operation, can increase work efficiency.

Further preferably, the identification of above-mentioned fixed star, coupling adopt quadrilateral algorithm to carry out, its matching principle is that the quadrilateral of four bright stars composition in CCD star chart picture is identical with the four edges ratio of the quadrilateral that four bright stars in star catalogue form, i.e. formula (14):

$\frac{S_{\mathrm{star}1}}{S_{\mathrm{cat}1}}=\frac{S_{\mathrm{star}2}}{S_{\mathrm{cat}2}}=\frac{S_{\mathrm{star}3}}{S_{\mathrm{cat}3}}=\frac{S_{\mathrm{star}4}}{S_{\mathrm{cat}4}}---\left(14\right)$

Order (i=1,2,3,4), then side ratio SC _istandard deviation calculate by formula (15):

$\mathrm{\&sigma;}=\sqrt{\frac{\mathrm{\&Sigma;}{({\mathrm{SC}}_i-\mathrm{\&nu;})}^2}{3}}---\left(15\right)$

In above formula (15), $\mathrm{\&nu;}=\mathrm{\&Sigma;}\left({\mathrm{SC}}_i\right)/4;$

Then think that four bright stars that σ is minimum in CCD star chart picture mate with the bright star of four in star catalogue;

According to four bright stars that the match is successful, calculate initial coordinate transformation model by formula (16), (17):

$X=\frac{a_{11}+a_{12}x+a_{13}y}{1+a_{31}x+a_{32}y}---\left(16\right)$

$Y=\frac{a_{21}+a_{22}x+a_{23}y}{1+a_{31}x+a_{32}y}---\left(17\right)$

In above formula (16), (17):

(X, Y) is the planimetric coordinates after bright star projection in star catalogue;

(x, y) is the image coordinate of the astrology in CCD star chart picture;

Then, utilize least square method, calculate initial coordinate transformation model by formula (18):

According to other astrology coordinates calculated in the initial coordinate transformation model of gained, CCD star chart picture except four bright stars, and in star catalogue, the planimetric coordinates of other stars completes the coupling of institute's any stars;

Recalculate Coordinate Transformation Models according to the matching result of all stars again, after rejecting the poor fixed star of matching result, calculate final Coordinate Transformation Models.

The technique effect that this optimal technical scheme is directly brought is, Quadrilateral Method implements easier and identifies that matching precision is higher.

Further preferably, the calculating of above-mentioned astronomic coordinates refers to, obtain fixed star in the image coordinate of the astrology in CCD star chart picture and star catalogue project after planimetric coordinates Coordinate Transformation Models after, utilize 0 ° and 180 ° of directions to photograph CCD star chart picture, carry out that survey station point astronomic coordinates calculates as follows:

(1) step, supposes that survey station point is positioned at the center of ccd image, i.e. (x, y) _z=(0,0);

(2) step, according to the Coordinate Transformation Models calculated, obtains the initial plane coordinate (X, Y) of survey station point _z, the more initial equatorial coordinate of survey station point is gone out according to the inverse operation formulae discovery of photogrammetric projective transformation formula;

(3) step, respectively according to the initial equatorial coordinate of survey station point of the CCD fixed star Image Acquisition of 0 ° and 180 ° direction shooting, calculates the initial astronomic coordinates (Φ of survey station point by following formula (19), (20) ₁, Λ ₁), (Φ ₂, Λ ₂):

Φ＝δ (19)

Λ＝α-GAST (20)

In above-mentioned formula (19), (20):

(Φ, Λ) is survey station point astronomic coordinates;

(δ, α) is survey station point equatorial coordinate;

GAST is the Greenwich hour angle calculated according to GNSS temporal information;

Then, (Φ is calculated by formula (21), (22) ₁, Λ ₁), (Φ ₂, Λ ₂) mean value:

$\mathrm{\&Phi;}=\frac{{\mathrm{\&Phi;}}_1+{\mathrm{\&Phi;}}_2}{2}---\left(21\right)$

$\mathrm{\&Lambda;}=\frac{{\mathrm{\&Lambda;}}_1+{\mathrm{\&Lambda;}}_2}{2}---\left(22\right);$

(4) step, by (Φ, Λ) inverse to equatorial coordinate, recalculates survey station point planimetric coordinates (X, Y) according to photogrammetric projective transformation formula _z1, (X, Y) _z2;

The CCD image coordinate (x, y) of survey station point is calculated again according to Formula of Coordinate System Transformation _z1, (x, y) _z2;

(5) step, repeats above iterative process, until the difference of the astronomic coordinates of twice adjacent calculation is less than a certain given numerical value, then computation process terminates, and gets the astronomic coordinates of average as observation station of the result of calculation of last twice astronomic coordinates;

(6) step, according to the instrument horizontal data that high-precision electronic level meter obtains, calculates the horizontal corrected value of instrument by formula (23), (24):

δΦ _l＝cos(θ+β)n ₁-sin(θ+β)n ₂(23)

δΛ _l＝sin(θ+β)n ₁+cos(θ+β)n ₂/cosΦ (24)

In above-mentioned formula (23), (24):

θ is the coordinate axis of CCD coordinate systems in image and the angle of east-west direction;

β is the coordinate axis of CCD coordinate systems in image and the angle of high-precision electronic level meter axle;

N ₁, n ₂for the instrument tilting value of high-precision electronic level meter record;

Calculate Ghandler motion corrected value by formula (25), (26) again, be subject to Ghandler motion impact with the CCD image data obtained instrument and carry out Ghandler motion correction:

δΦ _p＝-(x _pcosΛ-y _psinΛ) (25)

δΛ _p＝-(x _psinΛ+y _pcosΛ)tanΦ (26)；

In above-mentioned formula (25), (26):

(x _p, y _p) be the coordinate of instantaneous pole;

(Φ, Λ) is the astronomic coordinates of survey station point.

Above-mentioned Helmert formula is formula (27) and formula (28):

ξ＝Φ-φ (27)

η＝(Λ-λ)cosφ (28)；

In above-mentioned formula (27), (28):

ξ is the meridional component of plumb line deviation;

η is the fourth of the twelve Earthly Branches component at the tenth of the twelve Earthly Branches of plumb line deviation;

(Φ, Λ) is survey station point astronomic coordinates;

(φ, λ) is survey station point terrestrial coordinate.

The technique effect that this optimal technical scheme is directly brought is, use the plumb line deviation of the CCD fixed star processing result image iterative computation survey station point of 0 ° and 180 ° direction shooting, the computational accuracy of survey station point coordinate can be improved, eliminate the impact that digital zenith telescope optical axis not exclusively overlaps with mechanical axis.

In sum, the present invention, relative to prior art, has following beneficial effect:

1, method of the present invention easy and simple to handle, take a short time, to consume financial resource and material resource few, the precision of the plumb line deviation obtained is high; And there is pervasive characteristic, be suitable for the Quick Measurement of high precision plumb line deviation.

2, star catalogue of the present invention intercepts, and before observation and in observation, twice pair of star catalogue intercepts respectively, and make fixed star image recognition be matched to power and significantly improve, the whole identification matching process used time is less than 0.3 second.

3, the disposal route of method data of the present invention is relatively simple: single observation data processing procedure controls within 30s, and single observation computational accuracy is ± 0.3 " within.

Accompanying drawing explanation

Fig. 1 is FB(flow block) of the present invention;

Fig. 2 be on the astronomical longitude of the continuous 75 groups of observed results of embodiment 1 and astronomical latitude direction or the difference v of value and every group observations _icurve map;

Fig. 3 be the embodiment of the present invention 1 astronomical latitude group in the group internal standard difference distribution situation histogram of observed reading;

Fig. 4 be the embodiment of the present invention 1 astronomical longitude group in the group internal standard difference distribution situation histogram of observed reading.

Embodiment

Below in conjunction with drawings and Examples, the present invention is described in detail.

As shown in Figure 1, the high precision plumb line deviation method for fast measuring of integrated GNSS and CCD zenith telescope of the present invention, is characterized in that, comprise the following steps:

The first step, survey station point sets up CCD zenith telescope, its optical axis is made to point to zenith direction, zenith region is taken at 0 ° and 180 ° of directions continuously, obtain CCD star chart picture, and the tilt data of zenith telescope is obtained by high-precision electronic level meter, the terrestrial coordinate of observation time information and survey station point is obtained by GNSS;

Second step, reads the CCD star chart picture of the FITS form of CCD zenith telescope shooting, then carries out ground unrest elimination;

3rd step, the CCD star chart picture after eliminating noise carries out the automatic search in fixed star astrology region;

4th step, uses the correction Two-Dimensional Moment method improved, calculates the image coordinate of the center of energy in the fixed star astrology region that above-mentioned automatic search goes out;

5th step, processes applicable star catalogue, the terrestrial coordinate of the temporal information obtained by GNSS and survey station point, calculates the rough astronomic coordinates of survey station point in conjunction with EGM2008.According to survey station astronomic coordinates result of calculation, star catalogue is intercepted, the independent equatorial coordinate system in calculation exposure moment zenith region and apparent magnitude information, and independent equatorial coordinate system is projected in plane be converted to planimetric coordinates;

6th step, the fixed star projected in plane the astrology in CCD image and star catalogue identifies, mates;

7th step, according to the star pair that the match is successful, utilizes least square method, the coordinate conversion coefficient in iterative computation photogrammetric projective transformation formula and the initial equatorial coordinate of survey station point, then is converted to the initial astronomic coordinates of survey station point;

8th step, the zenith telescope inclination correction data obtain high-precision electronic level meter and Ghandler motion correct data and are converted on meridian direction and direction at the tenth of the twelve Earthly Branches, the fourth of the twelve Earthly Branches, determine the accurate astronomic coordinates of survey station point;

9th step, in conjunction with the survey station point terrestrial coordinate that GNSS measurement data calculates, according to Helmert formula, can calculate the plumb line deviation of survey station point.

It is 3 × 3 median filtering methods that above-mentioned background noise eliminates the method adopted, and concrete steps are as follows:

Step one, determines 3 × 3 neighborhoods put centered by certain pixel;

Step 2, by the gray-scale value (x of each pixel in this neighborhood ₁, x ₂, x ₃..., x _n) by sorting from small to large:

x ₁≤x ₂≤x ₃...≤x _n；

Step 3, goes out the intermediate value Y of the gray-scale value of each pixel by formulae discovery:

Step 4, replaces the gray-scale value of original central point with intermediate value Y, carry out medium filtering;

Step 5, determines 3 × 3 neighborhoods put centered by one other pixel, repeats above-mentioned steps two to step 4, until all pixel filter completes.

The automatic search in above-mentioned fixed star astrology region adopts algorithm of region growing to carry out.

The correction Two-Dimensional Moment method of above-mentioned improvement is in correction Two-Dimensional Moment method in the past, and use iterative calculation method to improve threshold value, its step is as follows:

1st step, extracts maximum gradation value g in CCD star chart picture _maxwith minimum gradation value g _min, order

$T^0=\frac{g_{\min }+g_{\max }}{2}---\left(2\right);$

2nd step, according to threshold value T ^kcCD star chart picture is divided into astrology region and background area, wherein, astrology regional location (x _i, y _j) gray-scale value g _k+1(x _i, y _j) > 0, and the gray-scale value of background area is 0, as shown by the equation:

Then, the average gray value in astrology region is tried to achieve by formula (4):

$g_{k+1}=\frac{\underset{g_k(x_i,y_j)>T^k}{\mathrm{\&Sigma;}}(g_k(x_i,y_j)-T^k)}{\underset{g_k(x_i,y_j)>T^k}{\mathrm{\&Sigma;}1}}---\left(4\right)$

Order formula (4) is substituted into, can obtain:

$T^{k+1}=\frac{\underset{g_k(x_i,y_j)>T^k}{\mathrm{\&Sigma;}}(g_k(x_i,y_j)-T^k)}{\underset{g_k(x_i,y_j)>T^k}{\mathrm{\&Sigma;}2}}---\left(5\right)$

3rd step, arranges a positive number μ, if | T ^k-T ^k+1| < μ, then termination of iterations computation process, try to achieve optimal threshold by formula (6):

T＝T ^k+B (6)

If | T ^k-T ^k+1|>=μ, then make T ^k=T ^k+1, proceed the iterative computation of second step.

4th step, by formula (7) and formula (8), obtain astrology regional center coordinate in CCD star chart picture:

$x_0=\frac{\mathrm{\&Sigma;\&Sigma;}{x}_i(g(x_i,y_j)-T)}{\mathrm{\&Sigma;\&Sigma;}(g(x_i,y_j)-T)}---\left(7\right)$

$y_0=\frac{\mathrm{\&Sigma;\&Sigma;}{y}_i(g(x_i,y_j)-T)}{\mathrm{\&Sigma;\&Sigma;}(g(x_i,y_j)-T)}---\left(8\right)$

The above-mentioned process being suitable for star catalogue is the NOVAS-F program by calling in AURIGA, the temporal information obtained in conjunction with GNSS and the terrestrial coordinate of survey station point, to the star place observed epoch of observation and be suitable for each influent factor that star catalogue fixed star epoch mean place there are differences and correct, to obtain independent equatorial coordinate system and TV star's information in zenith time of exposure region;

Above-mentioned equatorial coordinate is the star place coordinate obtained after star catalogue process;

Above-mentioned image plane coordinate is the astrology position coordinates that CCD fixed star image procossing obtains;

Above-mentioned being projected to by independent equatorial coordinate system plane is converted to planimetric coordinates calculates by projective transformation formula (9) and (10):

$l=\frac{\tan (\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_0)\cos q}{\cos (q-{\mathrm{\&delta;}}_0)}---\left(9\right)$

m＝tan(q-δ ₀) (10)

In above formula (9) and (10):

cotq＝cotδcos(α-α ₀)；

(l, m) is fixed star planimetric coordinates;

(δ ₀, α ₀) be initial astronomic coordinates;

(δ, α) is independent equatorial coordinate system.

The identification of above-mentioned fixed star, coupling adopt quadrilateral algorithm to carry out, and its matching principle is that the quadrilateral of four bright star compositions in CCD star chart picture is identical with the four edges ratio of the quadrilateral that four bright stars in star catalogue form, i.e. formula (11):

$\frac{S_{\mathrm{star}1}}{S_{\mathrm{cat}1}}=\frac{S_{\mathrm{star}2}}{S_{\mathrm{cat}2}}=\frac{S_{\mathrm{star}3}}{S_{\mathrm{cat}3}}=\frac{S_{\mathrm{star}4}}{S_{\mathrm{cat}4}}---\left(11\right)$

Order (i=1,2,3,4), then side ratio SC _istandard deviation calculate by formula (12):

$\mathrm{\&sigma;}=\sqrt{\frac{\mathrm{\&Sigma;}{({\mathrm{SC}}_i-\mathrm{\&nu;})}^2}{3}}---\left(12\right)$

In above formula (12), ν=Σ (SC _i)/4;

Then think that four bright stars that σ is minimum in CCD star chart picture mate with the bright star of four in star catalogue;

According to four bright stars that the match is successful, calculate initial coordinate transformation model by formula (13), (14):

$X=\frac{a_{11}+a_{12}x+a_{13}y}{1+a_{31}x+a_{32}y}---\left(13\right)$

$Y=\frac{a_{21}+a_{22}x+a_{23}y}{1+a_{31}x+a_{32}y}---\left(14\right)$

In above formula (13), (14):

(X, Y) is the planimetric coordinates after bright star projection in star catalogue;

(x, y) is the image coordinate of the astrology in CCD star chart picture;

Then, utilize least square method, calculate initial coordinate transformation model by formula (15):

According to other astrology coordinates calculated in the initial coordinate transformation model of gained, CCD star chart picture except four bright stars, and in star catalogue, the planimetric coordinates of other stars completes the coupling of institute's any stars;

Recalculate Coordinate Transformation Models according to the matching result of all stars again, after rejecting the poor fixed star of matching result, calculate final Coordinate Transformation Models.

The calculating of above-mentioned astronomic coordinates refers to, obtain fixed star in the image coordinate of the astrology in CCD star chart picture and star catalogue project after planimetric coordinates Coordinate Transformation Models after, utilize 0 ° and 180 ° of directions to photograph CCD star chart picture, carry out that survey station point astronomic coordinates calculates as follows:

(1) step, supposes that survey station point is positioned at the center of CCD star chart picture, i.e. (x, y) _z=(0,0);

(2) step, according to the Coordinate Transformation Models calculated, obtains the initial plane coordinate (X, Y) of survey station point _z, the more initial equatorial coordinate of survey station point is gone out according to the inverse operation formulae discovery of photogrammetric projective transformation formula;

(3) step, respectively according to the initial equatorial coordinate of survey station point of the CCD fixed star Image Acquisition of 0 ° and 180 ° direction shooting, calculates the initial astronomic coordinates (Φ of survey station point ₁, Λ ₁), (Φ ₂, Λ ₂), as shown in formula (16), (17):

Φ＝δ (16)

Λ＝α-GAST (17)

In above-mentioned formula (16), (17):

(Φ, Λ) is survey station point astronomic coordinates;

(δ, α) is survey station point equatorial coordinate;

GAST is the Greenwich hour angle calculated according to GNSS temporal information;

Then, (Φ is calculated by formula (18), (19) ₁, Λ ₁), (Φ ₂, Λ ₂) mean value:

$\mathrm{\&Phi;}=\frac{{\mathrm{\&Phi;}}_1+{\mathrm{\&Phi;}}_2}{2}---\left(18\right)$

$\mathrm{\&Lambda;}=\frac{{\mathrm{\&Lambda;}}_1+{\mathrm{\&Lambda;}}_2}{2}---\left(19\right);$

(4) step, by (Φ, Λ) inverse to equatorial coordinate, recalculates survey station point planimetric coordinates (X, Y) according to photogrammetric projective transformation formula _z1, (X, Y) _z2;

The CCD image coordinate (x, y) of survey station point is calculated again according to Formula of Coordinate System Transformation _z1, (x, y) _z2;

(5) step, repeats above iterative process, until the difference of the astronomic coordinates of twice adjacent calculation is less than a certain given numerical value, then computation process terminates, and gets the astronomic coordinates of average as observation station of the result of calculation of last twice astronomic coordinates;

(6) step, according to the instrument horizontal data that high-precision electronic level meter obtains, calculates the horizontal corrected value of instrument by formula (20), (21):

δΦ _l＝cos(θ+β)n ₁-sin(θ+β)n ₂(20)

δΛ _l＝sin(θ+β)n ₁+cos(θ+β)n ₂/cosΦ (21)

In above-mentioned formula (20), (21):

θ is the coordinate axis of CCD coordinate systems in image and the angle of east-west direction;

β is the coordinate axis of CCD coordinate systems in image and the angle of high-precision electronic level meter axle;

N ₁, n ₂for the instrument tilting value of high-precision electronic level meter record;

Calculate Ghandler motion corrected value by formula (22), (23) again, be subject to Ghandler motion impact with the CCD image data obtained instrument and carry out Ghandler motion correction:

δΦ _p＝-(x _pcosΛ-y _psinΛ) (22)

δΛ _p＝-(x _psinΛ+y _pcosΛ)tanΦ (23)；

In above-mentioned formula (22), (23):

(x _p, y _p) be the coordinate of instantaneous pole;

(Φ, Λ) is the astronomic coordinates of survey station point.

Above-mentioned Helmert formula is formula (24) and formula (25):

ξ＝Φ-φ (24)

η＝(Λ-λ)cosφ (25)；

In above-mentioned formula (24), (25):

ξ is the meridional component of plumb line deviation;

η is the fourth of the twelve Earthly Branches component at the tenth of the twelve Earthly Branches of plumb line deviation;

(Φ, Λ) is survey station point astronomic coordinates;

(φ, λ) is survey station point terrestrial coordinate.

Embodiment 1

For certain university campus, adopt the high precision plumb line deviation method for fast measuring of above-mentioned integrated GNSS and CCD zenith telescope, plumb line deviation Quick Measurement is carried out to this university campus.

Concrete measuring process is as follows:

From on Dec 19,6 days to 2013 June in 2013, in experimentation, choose weather conditions therebetween better, 34 day night that observing environment is comparatively suitable, carry out the data acquisition and procession of integrated GNSS and CCD zenith telescope, total observation group number is 275 groups, often organizing average observed duration is 15min, often organizing average observed number of times is 26 times, often organizing average star number used is 724, in group, average observed precision is 0.246 " (astronomical latitude), 0.269 " (astronomical longitude), average observed group number is 11 groups the whole night, average observed precision is 0.10 " (astronomical latitude) the whole night, 0.12 " (astronomical longitude).

The astronomy group internal standard difference distribution situation of above-mentioned whole process 275 group observations is as follows:

In 275 group observationses, in the group of astronomical latitude direction, observed reading standard deviation is 0 ~ 0.05 " between have 16 groups, 0.05 ~ 0.1 " between have 9 groups, 0.1 ~ 0.15 " between have 42 groups, 0.15 ~ 0.2 " between have 51 groups, 0.2 ~ 0.25 " between have 74 groups, 0.25 ~ 0.3 " between have 53 groups, 0.3 ~ 0.35 " between have 15 groups, 0.35 ~ 0.4 " between have 6 groups, 0.4 ~ 0.45 " between have 4 groups, 0.45 ~ 0.5 " between have 2 groups, 0.5 ~ 0.55 " between have 2 groups, 0.55 ~ 0.6 " between have 1 group, it can thus be appreciated that have 245 groups in 275 groups of astronomical latitude observed reading standard deviations 0.3 " within, account for 89% of sum, have 30 groups 0.3 " more than, account for 11% of sum.

In the group of astronomical longitude direction, observed reading standard deviation is 0 ~ 0.05 " between have 12 groups, 0.05 ~ 0.1 " between have 5 groups, 0.1 ~ 0.15 " between have 19 groups, 0.15 ~ 0.2 " between have 46 groups, 0.2 ~ 0.25 " between have 82 groups, 0.25 ~ 0.3 " between have 50 groups, 0.3 ~ 0.35 " between have 18 groups, 0.35 ~ 0.4 " between have 16 groups, 0.4 ~ 0.45 " between have 15 groups, 0.45 ~ 0.5 " between have 2 groups, 0.5 ~ 0.55 " between have 5 groups, 0.55 ~ 0.6 " between have 4 groups, 0.6 ~ 0.65 " between have 1 group, it can thus be appreciated that have 214 groups in 275 groups of astronomical longitude observed reading standard deviations 0.3 " within, account for 77.8% of sum, have 61 groups 0.3 " more than, account for 22.2% of sum.

By to the observation of raw observation and observation condition and analysis, find 0.3 is exceeded for observed reading standard deviation in group " observation group; be not mainly adjusted to the best owing to seeing time keeping instrument horizontality, or observing environment is not good, wind-force is comparatively large or zenith region cloud layer is thicker.These data are screened, to choose on same survey station point, integrated GNSS and the CCD zenith telescope plumb line deviation measurement result of 10 days, and be divided into 75 groups to carry out data processing to obtained CCD star chart picture, often group comprises the CCD star chart picture that observation 15min obtains.

By the result of the difference of above-mentioned 75 groups of every group observationses and mean value, draw out group internal standard difference Butut (Fig. 4) of observed reading in group internal standard difference Butut (Fig. 3) of observed reading in astronomical latitude group and astronomical longitude group respectively.

In the measurements, if the true value of observed reading is unknown, then bessel formula (i.e. following formula (26)) is generally adopted to ask for the medial error m of observed reading, the standard as weighing accuracy of observation:

$m=\&PlusMinus;\sqrt{\frac{\left[\mathrm{vv}\right]}{n-1}}---\left(26\right);$

In above formula (27), n is observation frequency, and [vv] is the quadratic sum of each observation correction, repeatedly observes, get its arithmetic mean under identical observation condition to some amounts as or value, observation frequency is more more close to true value, each observed reading l can be calculated _icorrected value

Choose above-mentioned continuous 75 groups of observed results, draw out Fig. 2: on astronomical longitude and astronomical latitude direction or the difference v of value and every group observations _icurve map.

As shown in Figure 2, in continuous 75 group observationses, every group observations and observe for each time the difference of arithmetic mean have 74 groups ± 0.3 in astronomical latitude direction in " within, only 1 group ± 0.3 ~ 0.4 ";

73 groups are had ± 0.3 in " within, only 2 groups ± 0.3 ~ 0.4 " in astronomical longitude direction.

Can calculate according to bessel formula: in 75 group observationses in Fig. 2, the medial error of astronomical longitude observed reading is ± 0.112 ", the medial error of astronomical latitude observed reading is ± 0.124 ", to meet completely in first-class astronomical sight specification ± 0.3 " requirement.

Claims (8)

1. a high precision plumb line deviation method for fast measuring for integrated GNSS and CCD zenith telescope, is characterized in that, comprise the following steps:

The first step, survey station point sets up CCD zenith telescope, its optical axis is made to point to zenith direction, zenith region is taken at 0 ° and 180 ° of directions continuously, obtain CCD star chart picture, and the tilt data of zenith telescope is obtained by high-precision electronic level meter, the terrestrial coordinate of observation time information and survey station point is obtained by GNSS;

Second step, reads the CCD star chart picture of the FITS form of CCD zenith telescope shooting, then carries out ground unrest elimination;

3rd step, the CCD star chart picture after eliminating noise carries out the automatic search in fixed star astrology region;

4th step, uses the correction Two-Dimensional Moment method improved, calculates the image coordinate of the center of energy in the fixed star astrology region that above-mentioned automatic search goes out;

5th step, processes applicable star catalogue, the terrestrial coordinate of the temporal information obtained by GNSS and survey station point, calculates the rough astronomic coordinates of survey station point in conjunction with EGM2008;

The rough astronomic coordinates of survey station point according to calculating intercepts as corresponding region with CCD star chart in star catalogue, and the independent equatorial coordinate system in calculation exposure moment zenith region and apparent magnitude information, and independent equatorial coordinate system is projected on section be converted to planimetric coordinates;

6th step, the fixed star projected in plane the astrology in CCD image and star catalogue identifies, mates;

7th step, according to the star pair that the match is successful, utilizes the coordinate conversion coefficient in least square method computational photogrammetry projective transformation formula, and the initial equatorial coordinate of iterative computation survey station point, then be converted to the initial astronomic coordinates of survey station point;

8th step, the zenith telescope inclination correction data obtain high-precision electronic level meter and Ghandler motion correct data and are converted on meridian direction and direction at the tenth of the twelve Earthly Branches, the fourth of the twelve Earthly Branches, determine the accurate astronomic coordinates of survey station point;

9th step, in conjunction with the survey station point terrestrial coordinate that GNSS measurement data calculates, and according to Helmert formula, can calculate the plumb line deviation of survey station point.

2. the high precision plumb line deviation method for fast measuring of integrated GNSS and CCD zenith telescope according to claim 1, is characterized in that, it is 3 × 3 median filtering methods that described ground unrest eliminates the method adopted, and concrete steps are as follows:

Step one, determines 3 × 3 neighborhoods put centered by certain pixel;

Step 2, by the gray-scale value (x of each pixel in this neighborhood ₁, x ₂, x ₃..., x _n) by sorting from small to large:

x ₁≤x ₂≤x ₃...≤x _n；

Step 3, goes out the intermediate value Y of the gray-scale value of each pixel by formulae discovery:

Step 4, replaces the gray-scale value of original central point with intermediate value Y, carry out medium filtering;

Step 5, determines 3 × 3 neighborhoods put centered by one other pixel, repeats above-mentioned steps two to step 4, until all pixel filter completes.

3. the high precision plumb line deviation method for fast measuring of integrated GNSS and CCD zenith telescope according to claim 1, is characterized in that, the automatic search in described fixed star astrology region adopts algorithm of region growing to carry out.

4. the high precision plumb line deviation method for fast measuring of integrated GNSS and CCD zenith telescope according to claim 1, it is characterized in that, the correction Two-Dimensional Moment method of described improvement, is in correction Two-Dimensional Moment method in the past, use iterative calculation method to improve threshold value, its step is as follows:

1st step, extracts maximum gradation value g in CCD star chart picture _maxwith minimum gradation value g _min, order

$T^0=\frac{g_{\min }+g_{\max }}{2}---\left(2\right);$

2nd step, according to threshold value T ^kcCD star chart picture is divided into astrology region and background area, wherein, astrology regional location (x _i, y _j) gray-scale value g _k+1(x _i, y _j) > 0, and the gray-scale value of background area is 0, as shown by the equation:

Then, the average gray value in astrology region is tried to achieve by formula (4):

$g_{k+1}=\frac{\underset{g_k(x_i,y_j)>T^k}{\mathrm{\&Sigma;}}(g_k(x_i,y_i)-T^k)}{\underset{g_k(x_i,y_i)>T^k}{\mathrm{\&Sigma;}}}---\left(4\right)$

Order formula (4) is substituted into, can obtain:

$T^{k+1}=\frac{\underset{g_k(x_i,y_j)>T^k}{\mathrm{\&Sigma;}}(g_k(x_i,y_i)-T^k)}{\underset{g_k(x_i,y_i)>T^k}{\mathrm{\&Sigma;}}}---\left(5\right)$

3rd step, arranges a positive number μ, if | T ^k-T ^k+1| < μ, then termination of iterations computation process, try to achieve optimal threshold by formula (6):

T＝T ^k+B (6)

If | T ^k-T ^k+1|>=μ, then make T ^k=T ^k+1, proceed the iterative computation of second step;

4th step, by formula (7) and formula (8), obtain astrology regional center coordinate in CCD star chart picture:

$x_0=\frac{\mathrm{\&Sigma;\&Sigma;}{x}_i(g(x_i,y_j)-T)}{\mathrm{\&Sigma;\&Sigma;}(g(x_i,y_j)-T)}---\left(7\right)$

$y_0=\frac{\mathrm{\&Sigma;\&Sigma;}{y}_i(g(x_i,y_i)-T)}{\mathrm{\&Sigma;\&Sigma;}(g(x_i,y_j)-T)}---\left(8\right)$

5. the high precision plumb line deviation method for fast measuring of integrated GNSS and CCD zenith telescope according to claim 1, it is characterized in that, the described process being suitable for star catalogue is the NOVAS-F program by calling in AURIGA, and in conjunction with the temporal information that GNSS obtains, to the star place observed epoch of observation and be suitable for each influent factor that star catalogue fixed star epoch mean place there are differences and correct, to obtain the apparent place of fixed star in star catalogue;

Then, intercept region corresponding with CCD fixed star star image in star catalogue, intercepting process is divided into two steps:

The first step, before integrated CCD and GNSS numeral zenith telescope is observed, first can intercept by the fixed star information in survey station point zenith region this observation period in star catalogue, star catalogue intercepts scope so that the fixed star that may appear in visual field, zenith region during observation is all covered as principle, its survey station point equatorial coordinate (α ₂, δ ₂) tried to achieve by following formula (9):

δ ₂＝Φ ₂

α ₂＝Λ ₂+GAST (9)

In above formula (9), (Φ ₂, Λ ₂) be the initial astronomic coordinates of survey station point, this initial astronomic coordinates is the plumb line deviation adopting GNSS geodetic surveying result and EGM2008 gravity field model to calculate survey station point, then pushes away to obtain the astronomic coordinates of survey station point, gets that approximate value obtains;

Described star catalogue intercepts scope and is greater than digital zenith telescope visual field;

After intercepting scope is determined, use apparent position calculation of star model to intercept the fixed star information that may appear within the scope of this in the observation period, and these information are stored;

Second step, before carrying out the work of fixed star coupling, need intercept above star catalogue region further, and enable this region completely corresponding with CCD fixed star image-region, it is as follows that it intercepts condition:

α _min＜α ₁＜α _max

δ _min＜δ ₁＜δ _max(10)

In above formula (10):

(α ₁, δ ₁) be the equatorial coordinate of fixed star in visual field, zenith region;

α _min, α _max, δ _min, δ _maxdepend on the equatorial coordinate (α of survey station point ₂, δ ₂) and the visual field Fi of digital zenith telescope, its computing formula is as follows:

α _min＝α ₂-Fi/2

α _max＝α ₂+Fi/2

δ _min＝δ ₂-Fi/2 (11)

δ _max＝δ ₂+Fi/2

Described being projected to by independent equatorial coordinate system plane is converted to planimetric coordinates calculates by projective transformation formula (12) and (13):

$l=\frac{\tan ({\mathrm{\&alpha;}}_1-{\mathrm{\&alpha;}}_2)}{\cos (q-{\mathrm{\&delta;}}_2)}---\left(12\right)$

m＝tan(q-δ ₂) (13)

In above formula (12) and (13):

cotq＝cotδ ₁cos(α ₁-α ₂)；

(l, m) is fixed star planimetric coordinates;

(α ₁, δ ₁) be the equatorial coordinate of fixed star in visual field, zenith region.

6. the high precision plumb line deviation method for fast measuring of integrated GNSS and CCD zenith telescope according to claim 1, it is characterized in that, the identification of described fixed star, coupling adopt quadrilateral algorithm to carry out, its matching principle is that the quadrilateral of four bright stars composition in CCD star chart picture is identical with the four edges ratio of the quadrilateral that four bright stars in star catalogue form, i.e. formula (14):

$\frac{S_{\mathrm{star}1}}{S_{\mathrm{cat}1}}=\frac{S_{\mathrm{star}2}}{S_{\mathrm{cat}2}}=\frac{S_{\mathrm{star}3}}{S_{\mathrm{cat}3}}=\frac{S_{\mathrm{star}4}}{S_{\mathrm{cat}4}}---\left(14\right)$

Order (i=1,2,3,4), then side ratio SC _istandard deviation calculate by formula (15):

$\mathrm{\&sigma;}=\sqrt{\frac{\mathrm{\&Sigma;}{({\mathrm{SC}}_i-v)}^2}{3}}---\left(15\right)$

In above formula (15), ν=Σ (SC _i)/4;

Then think that four bright stars that σ is minimum in CCD star chart picture mate with the bright star of four in star catalogue;

According to four bright stars that the match is successful, calculate initial coordinate transformation model by formula (16), (17):

$X=\frac{a_{11}+a_{12}x+a_{13}y}{1+a_{31}x+a_{32}y}---\left(16\right)$

$Y=\frac{a_{21}+a_{22}x+a_{23}y}{1+a_{31}x+a_{32}y}---\left(17\right)$

In above formula (16), (17):

(X, Y) is the planimetric coordinates after bright star projection in star catalogue;

(x, y) is the image coordinate of the astrology in CCD star chart picture;

Then, utilize least square method, calculate initial coordinate transformation model by formula (18):

$\left[\begin{array}{c}{a}_{11}\\ a_{12}\\ a_{13}\\ a_{21}\\ a_{22}\\ a_{23}\\ a_{31}\\ a_{32}\end{array}\right]=\left[\begin{array}{cccccccc}n& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i& 0& 0& 0& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{X}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i{X}_i\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i^2& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i& 0& 0& 0& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i^2{X}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i{X}_i\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i^2& 0& 0& 0& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i{X}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i^2{X}_{\text{i}}\\ 0& 0& 0& n& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{Y}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i{Y}_i\\ 0& 0& 0& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i^2& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i^2{Y}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i{Y}_i\\ 0& 0& 0& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i^2& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i{Y}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i^2{Y}_i\\ -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{X}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i^2{X}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i{X}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{Y}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i^2{Y}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i{Y}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i^2{X}_i^2+x_i^2{Y}_i^2& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i{X}_i^2+x_i{y}_i{Y}_i^2\\ -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i{X}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i{X}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i^2{X}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i{Y}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i{Y}_i& -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i^2{Y}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i{X}_i^2+x_i{y}_i{Y}_i^2& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i^2{X}_i^2+y_i^2{Y}_i^2\end{array}\right]\left[\begin{array}{c}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{X}_i\\ \underset{i=1}{\overset{{}_{}n}{\mathrm{\&Sigma;}}}{X}_i{x}_i\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{X}_i{y}_i\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{Y}_i\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{Y}_i{x}_i\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{Y}_i{y}_i\\ -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{X}_i^2{x}_i+Y_i^2{x}_i\\ -\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{X}_i^2{y}_i+Y_i^2{y}_i\end{array}\right]---\left(18\right)$

According to other astrology coordinates calculated in the initial coordinate transformation model of gained, CCD star chart picture except four bright stars, and in star catalogue, the planimetric coordinates of other stars completes the coupling of institute's any stars;

Recalculate Coordinate Transformation Models according to the matching result of all stars again, after rejecting the poor fixed star of matching result, calculate final Coordinate Transformation Models.

7. the high precision plumb line deviation method for fast measuring of integrated GNSS and CCD zenith telescope according to claim 1, it is characterized in that, the calculating of described astronomic coordinates refers to, obtain fixed star in the image coordinate of the astrology in CCD star chart picture and star catalogue project after planimetric coordinates Coordinate Transformation Models after, utilize 0 ° and 180 ° of directions to photograph CCD star chart picture, carry out that survey station point astronomic coordinates calculates as follows:

(1) step, supposes that survey station point is positioned at the center of CCD star chart picture, i.e. (x, y) _z=(0,0);

(2) step, according to the Coordinate Transformation Models calculated, obtains the initial plane coordinate (X, Y) of survey station point _z, the more initial equatorial coordinate of survey station point is gone out according to the inverse operation formulae discovery of photogrammetric projective transformation formula;

(3) step, respectively according to the initial equatorial coordinate of survey station point of the CCD fixed star Image Acquisition of 0 ° and 180 ° direction shooting, calculates the initial astronomic coordinates (Φ of survey station point ₁, Λ ₁), (Φ ₂, Λ ₂), as shown in formula (19), (20):

Φ＝δ (19)

Λ＝α-GAST (20)

In above-mentioned formula (19), (20):

(Φ, Λ) is survey station point astronomic coordinates;

(δ, α) is survey station point equatorial coordinate;

GAST is the Greenwich hour angle calculated according to GNSS temporal information;

Then, (Φ is calculated by formula (21), (22) ₁, Λ ₁), (Φ ₂, Λ ₂) mean value:

$\mathrm{\&Phi;}=\frac{{\mathrm{\&Phi;}}_1+{\mathrm{\&Phi;}}_2}{2}---\left(21\right)$

$\mathrm{\&Lambda;}=\frac{{\mathrm{\&Lambda;}}_1+{\mathrm{\&Lambda;}}_2}{2}---\left(22\right);$

(4) step, by (Φ, Λ) inverse to equatorial coordinate, recalculates survey station point planimetric coordinates (X, Y) according to photogrammetric projective transformation formula _z1, (X, Y) _z2;

The CCD image coordinate (x, y) of survey station point is calculated again according to Formula of Coordinate System Transformation _z1, (x, y) _z2;

(5) step, repeats above iterative process, until the difference of the astronomic coordinates of twice adjacent calculation is less than a certain given numerical value, then computation process terminates, and gets the astronomic coordinates of average as observation station of the result of calculation of last twice astronomic coordinates;

(6) step, according to the instrument horizontal data that high-precision electronic level meter obtains, calculates the horizontal corrected value of instrument by formula (23), (24):

δΦ _l＝cos(θ+β)n ₁-sin(θ+β)n ₂(23)

δΛ _l＝sin(θ+β)n ₁+cos(θ+β)n ₂/cosΦ (24)

In above-mentioned formula (23), (24):

θ is the coordinate axis of CCD coordinate systems in image and the angle of east-west direction;

β is the coordinate axis of CCD coordinate systems in image and the angle of high-precision electronic level meter axle;

N ₁, n ₂for the instrument tilting value of high-precision electronic level meter record;

Calculate Ghandler motion corrected value by formula (25), (26) again, be subject to Ghandler motion impact with the CCD image data obtained instrument and carry out Ghandler motion correction:

δΦ _p＝-(x _pcosΛ-y _psinΛ) (25)

δΛ _p＝-(x _psinΛ+y _pcosΛ)tanΦ (26)；

In above-mentioned formula (25), (26):

(x _p, y _p) be the coordinate of instantaneous pole;

(Φ, Λ) is the astronomic coordinates of survey station point.

8. the high precision plumb line deviation method for fast measuring of integrated GNSS and CCD zenith telescope according to claim 1, is characterized in that, described Helmert formula is formula (27) and formula (28):

ξ＝Φ-φ (27)

η＝(Λ-λ)cosφ (28)；

In above-mentioned formula (27), (28):

ξ is the meridional component of plumb line deviation;

η is the fourth of the twelve Earthly Branches component at the tenth of the twelve Earthly Branches of plumb line deviation;

(Φ, Λ) is survey station point astronomic coordinates;

(φ, λ) is survey station point terrestrial coordinate.