CN102735216B - CCD stereoscopic camera three-line imagery data adjustment processing method - Google Patents

Abstract

The invention discloses a kind of CCD stereoscopic camera three-line imagery data adjustment processing method, the method comprising the steps of: step 1, obtains almanac data and the attitude data of lunar surface view data and correspondence; Step 2, utilizes image automatic Matching to extract the view data of same rail different visual angles and the corresponding image points of adjacent orbit view data; Step 3, utilizes above-mentioned data and corresponding image points, adopts block adjustment with independent model to carry out survey district adjustment processing; Step 4, the ball adjustment processing whole month is carried out on the basis surveying district's adjustment.Method of the present invention does not need the coordinate resolving elements of exterior orientation and pass point in adjustment processing, under the prerequisite ensureing data processing precision, significantly simplify adjustment processing process.

Description

Adjustment processing method for CCD stereo camera three-linear array image data

Technical Field

The invention relates to the field of image processing methods and application, in particular to a adjustment processing method suitable for lunar CCD stereo camera three-line array image data.

Background

The CCD stereo camera is an area array camera with only one group of optical lenses and detectors of 1024 x 1024, the 11 th line, the 512 th line and the 1013 th line in the direction perpendicular to the flight direction are read on the area array and are respectively used as image arrays with three different visual angles of front view, front view and back view, the number of pixels of each line of linear array is 512 columns, the visual angle difference between the adjacent linear arrays of the front view, the front view and the back view is 16.7 degrees, and the imaging spectral band is a visible light wave band of 0.5-0.75 mu m. According to the preset on-orbit flight parameters of the satellite, the scanning speed of the camera is set to be 11.89 frames/second, so that three two-dimensional image strips which are acquired by three linear arrays on the orbit of 200km height in a push-scanning mode in the forward direction, the downward direction and the backward direction can be used for imaging 100% of the lunar surface flown by the satellite. The three image strips are almost obtained simultaneously, the width is 60km, the satellite has 100% overlap in the flight direction, the overlap degree of adjacent orbit image strips near the equator is about 41%, the overlap degree of high-latitude areas is larger, and sufficient image information is provided for geometrical reconstruction of three-dimensional landforms of the lunar surface. Fig. 2 shows a process of image data acquisition by the CCD stereo camera.

The main purpose of adjustment processing of three-line array image data of a CCD stereo camera is to realize seamless embedding and absolute orientation of image data acquired by a satellite at different time, different positions and the like in a full-month range, and according to the principle of satellite photogrammetry, the primary task of adjustment processing is to solve external orientation elements of each line array of the three-line array image data (namely acquiring spatial position and attitude information of the satellite at the moment of the line array image). Because the CCD stereo camera only has three images of front view, front view and back view in the same scanning period, and is influenced by the change of external orientation elements, topographic relief and the like, the homonymous image points (at least 6 are needed) used for calculating the external orientation elements of the linear array cannot fall on the three images, so that the external orientation elements of the scanning period cannot be solved geometrically, and some simplified algorithms are needed for approximate processing. For satellite photography, a platform is stable, elements of exterior orientation do not change greatly, and Hofmann, king task and the like propose to adopt time with a proper large distance, namely orientation time or EFP time, discretize a course model and approximately express the course model and the elements of exterior orientation.

In the last 80 th century, Hofmann et al, Germany scholars, created a beam adjustment method for processing three-linear array CCD images, namely a directional image method, for the first time in MOMS projects. The method solves the problem that the external orientation element of each sampling period (time) of the CCD image cannot be solved only by the three-linear array image into the external orientation element of the directional image, and the external orientation element of any CCD image sampling time is obtained by interpolation from the adjacent directional image value. Therefore, the adjustment calculation only needs to solve the external orientation elements of the directional images, so that the number and the calculation amount of the external orientation elements are reduced.

The equivalent image method is a processing method proposed by Wangxiang et al, the algorithm transforms a group of linear array images with a certain time interval to generate an equivalent image of a plane center projection, and then the equivalent image is processed by the algorithm of the plane center projection. This method also reduces the number and computational effort of exterior orientation elements.

Although the number of external orientation elements in adjustment processing is obviously reduced by processing methods such as the oriented film method and the equivalent photo method, the coordinates of encrypted points (homologous image points) used for solving the external orientation elements are also unknowns to be solved in the adjustment processing, the number is huge, and the workload and the difficulty of the adjustment processing are increased.

Disclosure of Invention

The invention aims to provide a method for adjusting the balance of lunar three-linear array image data of a CCD stereo camera by using a method for approximately processing the image data of the whole flight line, such as an orientation film method, an EFP (equivalent image film) method and the like, and adopting the theory and method of area network aerial triangulation by an independent model method, thereby overcoming the problem that the orientation element outside each scanning line of the three-linear array image of the CCD stereo camera cannot be solved geometrically.

The invention discloses a CCD stereo camera three-line array image data adjustment processing method, which comprises the following steps: step 1, lunar surface image data and corresponding ephemeris data and attitude data are obtained; step 2, extracting image data of the same rail at different visual angles and homonymous image points of image data of adjacent rails by using an automatic image matching technology; step 3, using the data and the homonymous image points, and adopting an independent model method for block adjustment to carry out block adjustment processing; and 4, carrying out full moon adjustment treatment on the basis of the adjustment of the measurement area.

Preferably, the full moon surface is divided into a plurality of mapping areas, each mapping area is further divided into a plurality of measuring areas, and each measuring area is provided with a plurality of flight zone images.

Preferably, the different viewing angles include a front view, a front view and a rear view.

Preferably, the matching algorithm adopts a mode of combining Scale Invariant Feature Transform (SIFT) feature matching and least square matching.

Preferably, an SIFT feature matching algorithm is adopted to provide initial positions of features, and then least square matching is adopted to realize fine image matching.

Preferably, the independent model is constructed by: obtaining a directional slice sequence on an orthographic image in a measuring area range, forming directional slice pairs by the same-name directional slices on the foresight and the orthographic images, the orthographic images and the rear images and the foresight and the rear images, and constructing an independent model by each directional slice pair through relative directional processing, wherein the independent model is the minimum unit of the block adjustment of the independent model method.

The method of the invention does not need to solve the coordinates of the external orientation element and the encryption point in the adjustment processing, but solves another group of parameters (relative orientation parameter and absolute orientation parameter) for seamless mosaic and absolute orientation of the three-linear array image data in the whole month range, the number of unknowns to be solved is obviously less than that of the two methods in the prior art, the adjustment processing process is obviously simplified, but the data processing precision is not reduced.

Drawings

FIG. 1 is a process flow diagram of a survey area adjustment algorithm and a global adjustment algorithm;

FIG. 2 is a schematic diagram of a process for acquiring three-line array image data of a lunar surface of a CCD stereo camera;

FIG. 3 is a schematic diagram showing the division of the directional slice on the front view image of the CCD stereo camera;

FIGS. 4a and 4b are schematic diagrams showing the comparison of the results before and after the correction of the linear error of the independent model;

FIGS. 5a and 5b are schematic diagrams of the division space distribution of the adjustment processing measuring region of the CCD stereo camera;

FIGS. 6a and 6b are schematic diagrams of spatial distribution of lunar control points for survey area adjustment and global adjustment;

fig. 6c is a partial enlarged view of fig. 6 b.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to specific embodiments and the accompanying drawings.

Referring to the flowchart of fig. 1, the adjustment processing method of the lunar CCD stereo camera three-line array image data according to the present invention is described in detail. The adjustment processing method of the invention depends on the precondition that the satellite photography platform is relatively stable and the change of the elements of the exterior orientation is not large.

Referring to fig. 1, when the adjustment processing method is executed, data preparation and organization are first performed.

Taking CE-1 (goddess Chang' e I) as an example, the satellite obtains 1000 multi-orbit image data during its lifetime, and 628 (the data is merely exemplary, and any suitable number can be selected according to actual conditions) is selected from the image data which has uniform orbital tone, small difference of illumination condition and the like, and has better image quality and covers the whole moon surface. In addition, in order to perform adjustment processing, satellite ephemeris data and attitude data corresponding to the image data must be prepared, the satellite ephemeris data is satellite position data obtained by a measurement and control department performing tracking measurement on the CE-1 satellite according to ground observation, and the attitude is data such as rolling, pitching, yawing and the like of a satellite platform measured by instrument equipment carried on the satellite.

In order to ensure the consistency of the adjustment result in the global scope, the seamless mosaic and absolute orientation of the selected image data in the global scope must be realized, and the adjustment processing is carried out in the whole month scope. In order to carry out global adjustment processing, the global is divided into a plurality of measuring areas, adjustment processing of the measuring areas is carried out firstly, and then adjustment processing of the global is carried out.

The surface of the full moon is divided into a plurality of mapping areas, for example, three mapping areas of S70-N70 degrees, S70-S90 degrees and N70-N90 degrees, each mapping area is divided into a plurality of measuring areas, the whole world is divided into 202 measuring areas, and each measuring area is about 20 navigation band images. In order to ensure the connection condition between the measuring regions, the data is cut so that adjacent measuring regions overlap each other by about 4 ° in the longitudinal direction and about 2 ° in the latitudinal direction.

The results of the zoning are shown in fig. 4a and 4 b. FIG. 4a is a schematic diagram of the division of the test areas of the drawing at the middle and low latitudes of S70-N70 degrees, wherein 170 test areas are provided, the latitude span range is 14 degrees, 10 latitude bands are provided, and the longitude span range is different along with the latitude; FIG. 4b is a schematic diagram of the region division of the polar region map of S70-S90 DEG and N70-N90 DEG, and the total number of the regions is 32. The numbering information for each measurement area is given in the figure.

In this step, the acquired lunar surface image data, and the corresponding ephemeris data and attitude data, etc., are clipped according to the measurement areas and managed as input data for the subsequent image matching and adjustment processing.

And then, carrying out automatic matching processing on the three-linear array image by using the lunar surface image data acquired in the step.

The image matching mainly realizes the extraction of homonymous image points in the overlapping areas among the forward-looking, forward-looking and rear-looking images in the same orbit and between the image data of the adjacent orbits in each measuring area.

The matching algorithm adopts a mode of combining SIFT (Scale invariant feature transform) feature matching and least square matching. The feature matching based on the SIFT descriptor has been successfully applied to many fields due to invariance to image scale, brightness, contrast, rotation, translation and tiny affine, and is the best feature descriptor at present. The SIFT feature matching has invariance to the scale change and rotation of the image and strong adaptability to the illumination change and image deformation, and the DOG difference operator is mainly utilized in the calculation process, so that most of the found features are 'blobs' (round points), which is just suitable for the feature that the moon image has poor texture and round points with more tiny impact pits. Where the CCD image spatial resolution is 120m, impact pits within a diameter of 500m appear as bright circular dots on the image due to strong reflection at the edges. The image fine matching is mainly realized through least square matching, the most difficult problem of a least square matching algorithm is determination of an initial position, an SIFT feature matching algorithm is adopted to provide the initial position of features, and then least square matching is adopted to realize the image fine matching. The image matching precision is better than 0.3 pixel, and the spatial resolution of the image data of the CCD stereo camera is 120 m.

The three-linear array image automatic matching comprises automatic matching between co-rail image data (including between two forward-looking, front-looking and rear-looking images) and adjacent rail image data, and the main purpose of the automatic matching is to extract image coordinates of encrypted points (or co-name image points) for adjustment parameter (namely relative orientation parameter and absolute orientation parameter) calculation.

According to the homonymous image point matching method, homonymous image points which meet the requirements of relative orientation, adjustment of measuring areas and global adjustment are selected in quantity and distribution on the forward-looking, front-looking and rear-looking images of the homonymous zone, the adjacent zones and the adjacent measuring areas, namely the relative orientation connecting point, the absolute orientation connecting point and the global adjustment connecting point shown in figure 1. Then, the same-name image points are manually checked and corrected in a three-dimensional environment. These homonymous pixels are the input data for subsequent scout and global adjustment processing.

After the homonymous image points are acquired, the adjustment processing of the measuring region is carried out.

Referring to fig. 1, the steps within the dashed box are the survey area adjustment process. In order to realize seamless mosaic and absolute orientation of the three-linear array image data in the full-month range, the whole world is divided into a plurality of measuring areas, adjustment processing of the measuring areas is firstly carried out, and then adjustment processing of the whole world is carried out. The adjustment of the survey area and the adjustment of the global area both adopt the adjustment technology of the regional network by adopting an independent model method.

According to the principle of the oriented sheet proposed by Hofmann, Wangxiang and the like, the oriented time is selected at intervals of 50 lines on an orthophoria image strip in the measuring area range, 160 lines with the oriented time as the center are taken as the width of the oriented sheet, an oriented sheet sequence on an orthophoria image is obtained, and the boundary of the oriented sheet on the front view image and the rear view image is determined through image matching. The same name orientation sheets on the front view and the front view, the front view and the back view and the front view and the back view form orientation sheet pairs, and each orientation sheet pair is subjected to relative orientation processing to construct an independent model.

Fig. 3 is a schematic diagram of dividing the directional slice on the front view image of the CCD stereo camera. The horizontal direction in the figure represents the height of the front-view image strip of the CCD stereo camera, the height is 512 pixels, and the length is consistent with the actual image scanning line number. Black solid line in the independent modelAnd a center line (marked by a dotted line on the left side in the figure), wherein the acquisition time of the linear array data corresponding to the center line is the model orientation time, and the distance between the adjacent center lines is 50 rows. The height of the independent model is consistent with the height of the image strip, i.e. 512 pixels, centered on the independent model centerline, the model width is 160 lines. SymbolRepresenting the position of the image point of the same name on the orientation sheet, resulting from image matching. The independent model is the minimum unit of the adjustment of the area network by the independent model method, and the adjustment of the measuring area aims to connect the independent models constructed by the orientation sheet pairs in the measuring area in a seamless manner to form a ground three-dimensional model of the whole measuring area.

The adjustment processing process of the measuring area is completed by the following steps:

(1) orientation within the orientation sheet:

and converting the image coordinates of the image points with the same name into image plane coordinates with the main point of the directional image as the origin of coordinates. The transformation formula is as follows:

$\left.\begin{array}{c}{x}^{\′}=(d_x-x_0+L_0)\·\mathrm{ds}\\ y^{\′}=(y-y_0)\·\mathrm{ds}\\ d_x=x-L_i\end{array}\right\}---\left(1\right)$

wherein x 'and y' are coordinate values (in mm) of image plane, d_xIs the image coordinate value (in pixels) in the x direction in the orientation sheet; x and y are original image row and column numbers; ds is 14 μm and is the size of the detector pixel; x is the number of₀＝0.8050，y₀-0.7990 is the image plane coordinate value (unit is pixel) of the image principal point; l is₀For the position of the scanning line in the CCD area array (front view L)₀11; front view L₀512, rear view L₀＝1013)；L_iThe line number of the scan line in the image strip is the center of the slice.

And converting the coordinate values of the original image points with the same name into image plane coordinate values with the image principal point as a coordinate origin through intra-image orientation processing, wherein the coordinate values are used as input data of model relative orientation processing.

(2) Relative orientation treatment: and establishing an error equation by utilizing the coordinates of the same-name image points of the orientation sheet pair through a coplanar condition, and solving the relative orientation elements of the orientation sheet. The coplanar condition satisfies the following equation:

wherein,κ，ω ', κ' is the relative orientation element to be solved; b is the length of a base line, a forward-looking and front-looking or forward-looking and rear-looking model B is 60km, and the forward-looking and rear-looking model B is 120 km; u, v and w are left image space coordinates of the left image point, u ', v ' and w ' are right image space coordinates of the right image point, the same name image point image plane coordinates of the left and right orientation sheets are obtained through coordinate conversion of rotation matrixes of the left and right orientation sheets, and the rotation matrixes are constructed by external orientation angle elements of the orientation sheets.

The above equation (2)) is expanded by taylor series and the first order term is taken to obtain the error equation:

wherein v is_rResidual errors of the coplanar condition equation at the same-name image point are obtained;

calculating an initial value for introducing approximate values of the relative orientation elements;Δ ω ', Δ κ, Δ κ' are the corrections of the approximate values of the relative orientation elements, i.e. the corrections of the unknowns to be resolved; (x, y), (x ', y') are image coordinates of the left and right image points, and are calculated according to the formula (1); f is the focal length of the camera.

Considering that the orientation slice is an approximation, the errorThe equation takes only the first order of the taylor expansion. In order to improve the resolving precision, about 150 uniformly distributed orientation points are selected on each independent model, the orientation points are homonymous image points obtained through image matching, image coordinates of the image points are used as observation numbers, an error equation is constructed according to equation (3), a normal equation is constructed according to the least square theory, and the increment value of relative orientation parameters is resolvedΔκ，Δ ω ', Δ κ', the initial value of the relative orientation element plus the corresponding incremental value is taken as the relative orientation element value. The process requires iterative operations until the accuracy requirement is met.

And obtaining relative orientation elements of the independent model through relative orientation processing, and solving the model coordinates of the orientation points.

(3) And calculating model coordinates.

After the relative orientation, the rotation matrixes R and R' and the base length B of the left and right image pairs used for constructing the independent model are known, the homonymous light beam pairs corresponding to the homonymous image points on the left and right images are intersected, and the intersection point position is the model point of the homonymous image point on the surface of the independent model. According to the relative orientation element and the set base line length, the coordinate value of the model point in the model coordinate system can be calculated, and the calculation formula is as follows:

U＝U_s+Mu_t，V＝V_s+Mv_t，W＝W_s-Mf(4)

$\left.\begin{array}{c}{u}_t=-f\frac{a_{1}x+a_{2}y-a_{3}f}{c_{1}x+c_{2}y-c_{3}f}\\ v_t=-f\frac{b_{1}x+b_{2}y-b_{3}f}{c_{1}x+c_{2}y-c_{3}f}\\ {u^{\′}}_t=-f\frac{a_1^{\′}x+a_2^{\′}y-a_3^{\′}f}{c_1^{\′}x+c_2^{\′}y-c_3^{\′}f}\end{array}\right\}$

$M=\frac{B}{u_t-u_t^{\′}}$

wherein, U_s，V_s，W_gCoordinate values of the left photo projection center under the model coordinates; u, V and W are coordinate values of the model points under the model coordinates; a is_i，b_i，c_iAnd a'_i，b′_i，c′_i(i is 1, 2, 3) are respectively 9 matrix elements of the left and right shot rotation matrices R, R'; (x, y), (x ', y') are the image coordinates of the left and right image points; f is the focal length of the camera; u. of_t，v_t，u′_tAnd M is an intermediate variable constructed to simplify the expression of the formula.

Through model coordinate calculation, the model coordinate value of each model point is obtained and is used as input data of absolute orientation processing.

(4) And (5) correcting the coordinates of the model.

The experimental result shows that the independent model has obvious system error in the X direction, and the error is linearly changed in the X direction. To solve this problem, linear error correction in the X direction is added to the adjustment algorithm, and the results before and after error correction are shown in fig. 4a and 4 b. In the figure, a curved surface A and a curved surface B are constructed by three-dimensional coordinates of model points which are uniformly distributed on an independent model under an image space coordinate system. FIG. 4a is the adjustment result without linear error correction of the model, with a significant linear change in the X direction; FIG. 4b shows the adjustment result after linear error correction, and the systematic error is eliminated obviously. The coordinate correction equation is as follows:

U′＝U+aU(5)

wherein U is the model X coordinate. a is a coordinate correction coefficient.

The absolute orientation process will use the corrected model coordinates as input.

(5) Laser altimeter data is introduced as elevation control information.

In the adjustment processing process of the area network, the elevation information of the control point is obtained from DEM data with the spatial resolution of 3 kilometers, and the DEM data is made by using the detection data of the ChangE detector laser altimeter. The elevation control point error equation is as follows:

$v_i=Z_i-Z_i^{\mathrm{dem}}---\left(6\right)$

$Z_i^{\mathrm{dem}}=a_0+a_1{X}_i+a_2{Y}_i+a_3{X}_i^2+a_4{X}_i{Y}_i+a_5{Y}_i^2$

in the formula, v_iZ to be solved for the same-name image point_iCoordinate and DEM data acquisitionThe difference in coordinates. (X)_i，Y_i，Z_i) Is the coordinate of the lunar surface of the image point with the same name obtained by forward intersection according to (X)_i，Y_i，Z_i) Corresponding lunar surfaceObtaining longitude and latitude coordinates on DEM dataThe value is obtained. a is₀，a₁，a₂，a₃，a₄，a₅Is a coefficient of a quadratic surface, the surface is based on (X)_i，Y_i，Z_i) The number of elevation control points is i, as established by the DEM data within the 3 × 3 grid near the location on the DEM.

In the absolute orientation processing process, the error equation is used as a constraint equation to participate in adjustment, the elevation system deviation value between the absolute orientation result and the reference DEM data is eliminated, and the accuracy of absolute orientation is improved.

(6) Absolute orientation

Through the relative orientation process, the relative position and posture (in a model coordinate system) between the left photo and the right photo are calculated, but the absolute position and posture of the photo are not converted, and the relative orientation process is independent of the existence of the ground control point. In order to obtain the lunar coordinates of the model points, the independent model as a whole needs to be zoomed, translated and rotated, the connection condition between the adjacent models and the lunar control points are considered, and the absolute orientation of the independent model is realized through adjustment of the area network.

The model points are subjected to absolute directional transformation, and then the lunar surface coordinates can be obtained. The transformation adopts a spatial similarity transformation (seven-parameter transformation) formula as follows:

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\mathrm{\λR}(\mathrm{\Φ},\mathrm{\Ω},K)\left[\begin{array}{c}U+\mathrm{aU}\\ V\\ W\end{array}\right]+\left[\begin{array}{c}{X}_t\\ Y_t\\ Z_t\end{array}\right]---\left(7\right)$

wherein a is a linear correction coefficient; (X, Y, Z) are model point lunar photogrammetry coordinates; (U, V, W) are model coordinates of the model points; lambda is a model scale parameter to be solved, R (phi, omega, K) is a rotation matrix constructed by rotation angle parameters, and phi, omega and K are 3 rotation angle parameters; x_t，Y_t，Z_tIs 3 translation parameters. The lunar surface photogrammetry coordinate is a tangent plane coordinate system taking the central position of the survey area as an origin, and a conversion matrix of the tangent plane coordinate and the lunar surface photogrammetry coordinate is constructed according to the lunar fixation coordinate of the origin.

And (3) expanding the above formula by using a Taylor series and taking a first order term to obtain an error equation:

$v_X=\frac{\∂X}{\∂\mathrm{\Φ}}\mathrm{\Δ\Φ}+\frac{\∂X}{\∂\mathrm{\Ω}}\mathrm{\Δ\Ω}+\frac{\∂X}{\∂K}\mathrm{\ΔK}+\frac{\∂X}{{\∂X}_t}{\mathrm{\ΔX}}_t+\frac{\∂X}{\∂\mathrm{\λ}}\mathrm{\Δ\λ}+\frac{\∂X}{\∂a}\mathrm{\Δa}-(X-[X\left]\right)$

$v_Y=\frac{\∂X}{\∂\mathrm{\Φ}}\mathrm{\Δ\Φ}+\frac{\∂X}{\∂\mathrm{\Ω}}\mathrm{\Δ\Ω}+\frac{\∂X}{\∂K}\mathrm{\ΔK}+\frac{\∂X}{{\∂Y}_t}{\mathrm{\ΔY}}_t+\frac{\∂X}{\∂\mathrm{\λ}}\mathrm{\Δ\λ}+\frac{\∂X}{\∂a}\mathrm{\Δa}-(Y-[Y\left]\right)---\left(8\right)$

$v_Z=\frac{\∂X}{\∂\mathrm{\Φ}}\mathrm{\Δ\Φ}+\frac{\∂X}{\∂\mathrm{\Ω}}\mathrm{\Δ\Ω}+\frac{\∂X}{\∂K}\mathrm{\ΔK}+\frac{\∂X}{{\∂Z}_t}{\mathrm{\ΔZ}}_t+\frac{\∂X}{\∂\mathrm{\λ}}\mathrm{\Δ\λ}+\frac{\∂X}{\∂a}\mathrm{\Δa}-(Z-[Z\left]\right)$

in the formula, v_X，v_Y，v_ZIs the residual of the lunar coordinates, [ X ]]，[Y]，[Z]The calculation result obtained by substituting the approximate value of the unknown number for the formula (7); Δ λ, Δ Φ, Δ Ω, Δ K, Δ X_t，ΔY_t，ΔZ_tDelta a is a correction number of a pending value, namely an unknown correction number to be resolved;are error equation coefficients.

Image matching is carried out on the three-line array data, and after a certain number of uniformly distributed connection points of adjacent orientation sheets (namely homonymous image points of the overlapping area of the orientation sheets) are determined, the image matching method can be used according to the resultThe model coordinates of the connection points are calculated by the formula (4) and used as the connection conditions of adjacent independent models, then the model coordinates are substituted into an error equation (8), equations (6) and (8) are combined, and the absolute orientation parameter X can be solved by adopting a least square method_t，Y_t，Z_tΦ, Ω, K, λ and a model coordinate correction coefficient a. The connection points are chosen to ensure connection between the models in the flight band (image strips along the flight direction of the satellite) and between the models between adjacent flight bands. In the work, each model in the navigation band selects 5 rows (vertical flight direction) and 9 columns (along the flight direction) of uniformly distributed connection points as model connection points, and the overlapping area of adjacent models has 6 columns of connection points; the model overlapping area of the adjacent flight zones selects 3 columns and 4 rows of uniformly distributed connecting points as the flight zone connecting points.

Absolute orientation parameter X_t，Y_t，Z_tPhi, omega, K, lambda and the model coordinate correction coefficient a are input parameters for solving the lunar coordinates of each model point.

(7) Lunar point coordinate calculation

Using absolute orientation elements (X)_t，Y_t，Z_tΦ, Ω, K, λ) and the model coordinate correction coefficient a, are substituted into equation (7) to solve the section coordinate value (lunar photogrammetry coordinate). After the absolute orientation and the global leveling, the coordinate values under the photogrammetric coordinate system, namely the coordinate values under the tangent plane coordinate system with the central position of the measured area as the origin, are obtained. According to the exterior orientation elements of the independent model formed by the front-view and rear-view orientation sheets at the position, a conversion matrix between the tangent plane coordinates and the moon-fixed coordinates can be established, the conversion from the tangent plane coordinates to the moon-fixed coordinates is realized, and the lunar plane coordinate values of all model points are obtained.

Therefore, relative orientation parameters and absolute orientation parameters of each independent model of each independent measuring area are obtained, and the parameters can ensure seamless connection of all independent models in the measuring area, so that the construction of a ground three-dimensional model of the whole measuring area is realized. The survey area ground model is used as an independent model for global adjustment processing and is the minimum unit for global adjustment processing.

According to the above steps, the block adjustment processing is completed, and then the global adjustment processing is performed.

The survey area adjustment realizes seamless connection of each independent model in the survey area, and global adjustment is required to realize seamless mosaic of global data. The global surface is divided into a plurality of blocks (or measuring areas, see fig. 5b) according to a certain principle, and global adjustment is carried out by taking the measuring areas as units on the basis of adjustment of the measuring areas. The global adjustment also adopts the adjustment method of the area network by the independent model method, the only difference is that the independent model of the adjustment of the measuring area is a ground model constructed by the divided orientation sheet pairs, and the independent model of the global adjustment is a ground model corresponding to the whole measuring area. The input data of the survey area adjustment is connection points between adjacent independent models (constructed by directional sheets), and the connection points ensure the seamless mosaic of the data in the whole survey area; the input data of the global adjustment are the connection points between adjacent measuring areas, and the connection points ensure the seamless mosaic of all measuring areas in the global range.

Both the zonal and global adjustment must be performed with the support of control points, and the elevation control points are selected from a DEM with a spatial resolution of 3km, see equation (6) above. The plane control point coordinates are obtained by resolving a certain number and distribution of image points with the same name according to the intersection in the front. The distribution of the adjustment and global adjustment control points of the measurement area is shown in fig. 6a and 6b, fig. 6a is a distribution diagram of adjustment and global adjustment control points of the measurement area, taking the measurement area with the number of E004 as an example, a triangle represents the control points of the height (plane and elevation), and 9 control points are uniformly selected in the whole measurement area and distributed on the first model, the middle model and the last model of the first rail, the middle rail and the last rail; the round points represent elevation control points and are distributed in the center of each ground model in the measuring area; the black solid line represents the border line of the measurement area, and the number information of the adjacent measurement area is given in the figure. FIG. 6b is a distribution diagram of global adjustment lunar surface control points, wherein a triangle represents adjustment control points, and 9 control points are uniformly distributed globally; the five-pointed star represents elevation control points, the total number of the elevation control points is 202, and the elevation control points are distributed in the central position of each measuring area; the black dots represent the connection points between adjacent measurement areas, and the number of the black dots is about 23766, which are uniformly distributed in the overlapping area of the measurement areas, and fig. 6c is a partial enlarged view of the oblique line area in fig. 6b, wherein the symbol meaning in the figure is the same as that in fig. 6 b. The weight of the level control points was set to 20.0 and the elevation control points were set to 3.0, and the elevation values were obtained from global DEM data with a spatial resolution of 3km made by the "Chang' e I CE-1" altimeter.

And (3) processing by utilizing the internal orientation parameters, the block adjustment and the global adjustment of the CCD stereo camera to obtain the relative orientation parameters and the absolute orientation parameters of the independent model, and calculating the lunar surface coordinates of any image point with the same name on the three-linear array image according to the formula (1), the formula (4), the formula (5) and the formula (7) of the absolute orientation parameters of the block. In addition, enough quantity and distribution of connection points are selected in the overlapping area of adjacent independent models in the adjustment processing of the measuring area, so that seamless mosaic between the independent models in the measuring area is guaranteed, and enough quantity and distribution of connection points are also selected in the overlapping area of adjacent measuring areas in the adjustment processing of the global adjustment, so that seamless mosaic between the adjacent measuring areas in the global range is guaranteed. Therefore, the global connection problem of data is well solved through the processing of survey area adjustment and global adjustment, and processing parameters are provided for the subsequent production of a digital elevation model and an ortho-image.

The above-mentioned embodiments are intended to illustrate the objects, technical solutions and advantages of the present invention in further detail, and it should be understood that the above-mentioned embodiments are only exemplary embodiments of the present invention, and are not intended to limit the present invention, and any modifications, equivalents, improvements and the like made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (3)

1. A CCD stereo camera three-linear array image data adjustment processing method comprises the following steps:

step 1, acquiring lunar surface image data and corresponding ephemeris data and attitude data, wherein in the step, the surface of a full moon is divided into a plurality of mapping areas, each mapping area is divided into a plurality of measuring areas, and each measuring area is provided with a plurality of flight-zone images;

step 2, extracting image data of the same track at different visual angles and homonymous image points of image data of adjacent tracks by using an automatic image matching technology, wherein the different visual angles comprise a front view, a front view and a back view, in the step, a matching algorithm adopts a mode of combining Scale Invariant Feature Transform (SIFT) feature matching and least square matching, an SIFT feature matching algorithm is adopted to provide initial positions of features, and then least square matching is adopted to realize fine image matching;

step 3, using the data and the homonymous image points, and adopting an independent model method for block adjustment to carry out block adjustment processing;

step 4, carrying out full moon adjustment processing on the basis of the adjustment of the measuring area;

the independent model is constructed in the following way: obtaining a directional slice sequence on an orthographic image in a measuring area range, forming directional slice pairs by the same-name directional slices on the foresight and the orthographic image, the orthographic image and the rear image and the foresight and the rear image, and constructing an independent model by each directional slice pair through relative directional processing, wherein the independent model is a minimum unit of the adjustment of a regional network of an independent model method;

step 3 further comprises the steps of:

converting the image coordinates into image plane coordinates with the main point of the directional image as the origin of coordinates;

establishing an error equation by utilizing the coordinates of the same-name image points of the orientation sheet pair through a coplanar condition, and solving the relative orientation elements of the orientation sheet;

calculating the model coordinates of each image point with the same name on the independent model constructed by the orientation sheet pair by using relative orientation elements;

calculating the lunar coordinates of the image points with the same name,

the image coordinate is converted into an image plane coordinate with the main point of the directional slice image as the origin of coordinates, and a conversion formula is as follows:

$\left.\begin{array}{c}{x}^{\′}=(d_x-x_0+L_0)\·ds\\ y^{\′}=(y-y_0)\·ds\\ d_x=x-L_i\end{array}\right\}---\left(1\right)$

wherein, x 'and y' are coordinate values of the image plane of the directional slice, and the unit is mum; x and y are coordinate values of the original image, and the unit is a pixel; ds is the detector pixel size; x is the number of₀、y₀Is the image plane coordinate value of the image principal point, the unit is pixel, L₀The position of the scanning line in the CCD area array; l is_iScanning the line number of the line in the image strip for the center of the directional slice; d_xIn order to orient the row number of the sheet,

establishing an error equation by utilizing the coordinates of the same-name image points of the orientation sheet pair through a coplanarity condition, solving and solving the relative orientation elements of the orientation sheet, wherein the coplanarity condition satisfies the following equation:

wherein the sum of the values of ω,κ，ω ', κ' is the relative orientation element to be solved; b is the base length; u, v and w are left image space coordinates of the left image point, u ', v ' and w ' are right image space coordinates of the right image point, the coordinates of the image planes of the same name image points of the left and right orientation sheets are obtained through coordinate conversion of rotation matrixes of the coordinates, the rotation matrixes are constructed by external orientation angle elements of the orientation sheets,

and (2) expanding by using a Taylor series and taking a first order term to obtain an error equation:

in the formula,f is the focal length of the camera and,

and calculating the model coordinates of each image point with the same name on the independent model constructed by the orientation sheet pair by using the relative orientation elements, and calculating by the following formula:

U＝U_s+Mu_t，V＝V_s+Mv_t，W＝W_s-Mf(4)

$\left.\begin{array}{c}{u}_t=-f\frac{a_{1}x+a_{2}y-a_{3}f}{c_{1}x+c_{2}y-c_{3}f}\\ v_t=-f\frac{b_{1}x+b_{2}y-b_{3}f}{c_{1}x+c_{2}y-c_{3}f}\\ u_t^{\′}=-f\frac{a_1^{\′}x+a_2^{\′}y-a_3^{\′}f}{c_1^{\′}x+c_2^{\′}y-c_3^{\′}f}\end{array}\right\}$

$M=\frac{B}{u_t-u_t^{\′}}$

wherein, U_s，V_s，W_sA coordinate value of the projection center of the left directional slice under the model coordinate is obtained; u, V and W are coordinate values of the model points under the model coordinates; a is_i，b_i，c_iAnd a'_i，b′_i，c′_i9 matrix elements of left and right directional slice rotation matrices R, R', respectively, where i ═ 1, 2, 3; (x, y), (x ', y') are image plane coordinate values of the left and right orientation plates, and (x, y) is a result calculated by substituting the original image coordinate values of the same-name image points on the left orientation plate into the formula (1), and (x ', y') is a result calculated by substituting the original image coordinate values of the same-name image points on the right orientation plate into the formula (1); f is the focal length of the camera; u. of_t，v_t，u′_tAnd M is an intermediate variable constructed to simplify the expression of the formula.

2. The method of claim 1, wherein the global surface is divided into a plurality of mapping regions according to a predetermined rule, each mapping region is further divided into a plurality of measurement regions, and global adjustment is performed in units of the measurement regions on the basis of the adjustment of the measurement regions.

3. The method of claim 2, wherein the step 4 further comprises:

selecting a certain number of image points with the same name as the connecting points of the measuring area in the image data overlapping area of the adjacent measuring areas by adopting an automatic image matching technology;

laser lunar measurement data is introduced for joint adjustment, and the absolute control precision of the control network is improved;

and selecting control point data for global adjustment processing, and calculating absolute orientation parameters of the measuring area.