CN101604018B

CN101604018B - Method and system for processing high-definition remote sensing image data

Info

Publication number: CN101604018B
Application number: CN2009101616382A
Authority: CN
Inventors: 张力; 张继贤; 燕琴; 谢文寒; 程春泉; 高武俊; 许骥; 袁枫; 孙钰珊
Original assignee: Chinese Academy of Surveying and Mapping
Current assignee: Chinese Academy of Surveying and Mapping
Priority date: 2009-07-24
Filing date: 2009-07-24
Publication date: 2011-09-21
Anticipated expiration: 2029-07-24
Also published as: CN101604018A

Abstract

The invention provides a method and a system for processing high-definition remote sensing image data by applying a rational function model, a satellite image block adjustment and auto-match and distributed parallel processing technology. The method comprises the following steps: analyzing and processing input primary image data and computing parameters of the corresponding rational function model so as to obtain a stereo image pair by applying the computed rational function model; performing measurement and singlechip orientation of a control point and a connecting point on the primary image data, and performing block adjustment processing on an error between stereo image pixel data and the primary image data so as to improve the fit precision between the stereo image pixel data and the primary image data; generating an image similar to an epipolar ray, defining a matching template, measuring matched seed points and lines, performing automatic image matching and generating a digital elevation model; and respectively generating a digital orthophoto map and a digital line graph, and finishing manufacturing a product by applying the obtained data.

Description

High-resolution remote sensing image data processing method and system

Technical Field

The invention relates to a method and a system for processing high-resolution remote sensing image data and a corresponding computer product program, in particular to a method and a system for processing the high-resolution remote sensing image data by utilizing a rational function model, satellite image block adjustment with rare large-range area or without ground control, automatic matching based on multi-baseline multi-matching characteristics and a distributed parallelization processing technology based on a high-speed network and a corresponding computer product program.

Background

Since the twenty-first century, spatial technology and information technology have been further developed, demand for geographic information resources has rapidly increased, and geographic information resources are used as important strategic information resources in all countries of the world. With the accumulation of a large amount of geographic information data and the change of the demand of national economy and social informatization development, providing geographic information service for the whole society becomes a main task of surveying and mapping development in a new period, and meanwhile, a surveying and mapping system is also developing informatization characterized by the spatialization and real-time acquisition of basic geographic information, the automation and intellectualization of processing and the socialization of service networking.

In the aspect of aerial remote sensing systems, in recent years, with the rapid development of digital imaging technology, large-breadth frame-type digital aerial imaging systems represented by DMC (digital media controller), UltraCam-D, UltraCamX of Vexcel and SWDC (Single-dimensional space digital imaging) system of China surveying and mapping sciences research institute and three-linear-array imaging systems represented by ADS40 of Leicageosystem are rapidly popularized and applied. Various small and medium-breadth digital imaging systems based on low-altitude airplanes and unmanned aerial vehicle platforms are developed more and more rapidly, so that a full digital data chain consisting of data acquisition, data processing, intelligent information extraction, basic geographic information database building, updating and the like is formed by a photogrammetry technology and data processing thereof which are important acquisition means of spatial geographic information.

The aviation digital imaging system is a system integration of a multispectral panchromatic digital photographic system and a GPS/inertial navigation system, all the new characteristics are full-automatic processing of digital images, and the possibility is provided for modification and optimization of a traditional digital photogrammetric data processing algorithm. It is predicted that the aerial digital imaging system will gradually replace the traditional aerial camera and become the mainstream of aerial photographic instruments, which will change the traditional aerial photogrammetric data acquisition mode and data processing mode.

In the aspect of aerospace remote sensing systems, high-resolution remote sensing satellites have been vigorously developed since the first civil remote sensing satellite emission in 1972. With the rapid development of decades, and especially the breakthrough of two key technologies of high-precision navigation, positioning, attitude and time measurement system on high-resolution earth observation satellite and sensor platform, satellite remote sensing data has formed a global, various-resolution, multi-temporal, high-signal-to-noise ratio, high-geometric and radiation resolution (greater than 8 bits/pixel) image sequence covering from Landsat series (80 m) -TM (30 m) -SPOT1/2/4(10 m) -SPOT5HRS (10 × 5 m) -SPOT5HRG (2.5 m) -ikonoos (1.0 m) -QuickBird (0.6 m), which become necessary geospatial information for resource promotion, ecological environment protection and social economy planning. Earth observation satellites with topographic features, such as SPOT-5, IKONOS, etc., can acquire extremely large-range co-orbital or hetero-orbital stereoscopic images, and in addition, integration with high-precision navigation, positioning, attitude and time measurement systems (GPS/star/inertial integrated navigation system) greatly affects topographic measurement and topographic techniques.

High resolution optical satellite imaging systems are increasingly being used in the fields of remote sensing and photogrammetry. Such imaging systems are capable of providing not only high resolution panchromatic, multispectral imagery, but also stereoscopic imaging capabilities. The satellite image data is important data guarantee of surveying and mapping, and large-scale potential geographic information can be rapidly acquired. The method for extracting the large-scale basic geographic information by using the high-resolution remote sensing data plays an important role in digital region/digital earth construction, and has a huge prospect.

At present, the processing technology of the remote sensing information with medium and low resolution is mature, and the processing technology and the method system for extracting the basic geographic information with large scale by using the high-resolution remote sensing image are lagged, so that the technical problem which is needed to be broken through in the development of the surveying and mapping industry is formed. Therefore, it is necessary to conduct research on mapping techniques based on integration of high-resolution remote sensing data processing techniques. In particular, the development of high-precision block adjustment orientation of no control points or a small number of control points of satellite images, the generation of a full-automatic digital ground model, the manufacture of an orthophoto image with high resolution ratio, the stereo measurement of a digital line drawing and the like are key technologies of high-resolution satellite mapping.

In addition, in terms of remote sensing data processing algorithms and systems, in remote sensing image mapping and monitoring application systems, development of photogrammetry has entered the digital photogrammetry era in the last 90 years through the times of analog photogrammetry and analytic photogrammetry. The digital photogrammetry system has entered a stage of commercialization; the digital photogrammetry system is combined with a geographic information system, so that the digitization and automation of the surveying and mapping production process are promoted; and the GPS is utilized to determine the exterior orientation elements of aerial photography, so that aerial photogrammetry without ground control points or with few ground control points is realized, and heavy field control measurement work is avoided.

With the rapid development of computer hardware and software, Digital Photogrammetry Workstations (DPW) have gradually developed from original dedicated hardware-based devices and UNIX systems to more practical, more modular, and more automated systems using PC/Window system platforms, and have gained widespread use in the fields of daily topographic map production, three-dimensional geographic information extraction, remote sensing, and the like. At the present stage, there are already DPWs based on PC/Windows systems that differ in price, performance, complexity and functionality, appearing in the market. For example, SOCET by BAE System

Of Autometric

Of Z/I Imaging

Of PCI Geometics

Provided is a system. The basic functions of DPW can be summarized as: storage, processing and management of maximum or mass image data, support of stereo image observation, stereo vector data superposition and support of stereo vector dataExtracting and editing three-dimensional information; automatic or semi-automatic image orientation, including internal orientation, relative orientation and absolute orientation; automatic or semi-automatic image connection point extraction and measurement, automatic and high-precision area network adjustment; automatic Digital Elevation Model (DEM), Digital Orthophoto Map (DOM) generation, and Digital Line Graph (DLG) data acquisition, among others.

Currently, with the rapid development of the aerospace Digital Imaging (Digital Imaging) technology, the Active Sensing technology (Active Sensing), the sensor autonomous positioning (direct georgeffering) technology based on the DGPS/IMU combined system, and the automation/intelligent data processing technology, the photogrammetry and remote Sensing software system is also developing towards a new direction, and the software system is also developing in the form of the ISTARPixel-factor of infotera corporation of france^TM("pixel factory") software systems are examples, which are embodied in: (1) the processing of the aerospace remote sensing data is integrated, namely various aerospace images, optical and radar image data and LIDAR data can be processed by adopting a general algorithm, and the tradition that different types of remote sensing data are processed by adopting different professional modules in the prior art is broken through; (2) starting to adopt an algorithm based on multi-image processing, and simultaneously processing more than two images to acquire reliable and high-precision three-dimensional information of an imaging area; (3) the data processing is more intelligent and automatic; (4) and a distributed data processing mode is adopted, so that the efficiency of processing the mass remote sensing data is greatly improved. JX-4 digital photogrammetry workstation of four-dimensional information technology company and VirtuoZo full digital photogrammetry system of software company Limited are mainly used for the production of Digital Elevation Models (DEM), digital ortho images (DOM) and Digital Line Graphs (DLG) of various scales. At present, JX-4 and VirtuoZo digital photogrammetry workstations expand a module for mapping high-resolution satellite stereoscopic images on the basis of the original mature aerial image mapping function, and provide a technical means for mapping high-resolution satellite images.

The cluster system implementation technology is gradually mature, so that the large-scale parallel system becomes a computing resource which is easy to obtain. If a parallel computing technology can be applied, a plurality of computers and remote sensing image parallel processing software break through the traditional photogrammetry and remote sensing image processing procedures, do not necessarily take photos and image pairs as units, do not take operators as units, integrate production, quality detection and management into a whole, reasonably arrange the work of people and machines, and further improve the efficiency of digital photogrammetry, remote sensing image data processing and even spatial information extraction. Aiming at the problems, Wuhan university combines computer network technology, parallel processing technology, high-performance computing technology and digital photogrammetry processing technology based on research results in digital photogrammetry for many years, and researches and provides a set of high-performance new-generation aerospace digital photogrammetry processing platform: digital Photogrammetry Grid (DPgrid).

Disclosure of Invention

The invention provides a remote sensing data processing system which aims at the requirement of large-scale topographic map mapping by using high-resolution remote sensing images. Specifically, the invention provides and develops new technologies such as large-range remote sensing image sparse/uncontrolled area network adjustment based on an RFM (radio frequency mass spectrometry) general imaging model, high-precision DEM/DSM (digital elevation model/digital surface model) automatic extraction based on multi-baseline/multiple matching characteristics, vector contour line data semi-automatic acquisition, high-precision image map making and splicing and the like by fusing a computer technology and a network communication technology on the basis of modern photogrammetry and a remote sensing scientific technology theory, and realizes cluster distributed parallel computing which is based on a loose coupling parallel service middleware and can run in a high-speed local area network.

According to one aspect of the invention, a method for processing high-resolution remote sensing image data by utilizing a rational function model, satellite image block adjustment with a large range of area rare or without ground control, automatic matching based on multi-baseline and multi-matching characteristics and a distributed parallelization processing technology based on a high-speed network is provided, and the method comprises the following steps: analyzing and processing the input original image data, and calculating parameters of a corresponding rational function model so as to obtain a stereoscopic image pair by using the calculated rational function model; measuring control points and connection points of original image data, orienting a single chip, and performing block adjustment processing on errors between stereo image pixel data and the original image data so as to improve the fitting precision between the stereo image pixel data and the original image data; generating an approximate epipolar line image, defining a matching template, measuring a matching seed point line, automatically matching the images and generating a digital elevation model; and respectively generating a digital orthographic projection image and a digital line drawing, and finishing the manufacture of the product by using the obtained data.

According to a preferred embodiment of the present invention, calculating the parameters of the rational function model comprises: dividing input image data by a preset grid to obtain coordinates of grid points, and layering fluctuation ranges of image coverage areas to obtain preset height planes; constructing a strict geometric imaging model, and calculating the corresponding ground coordinates of each grid point on each equal altitude plane according to the strict geometric imaging model to obtain three-dimensional virtual object grid points; and iteratively solving coefficients of the rational function model by using the obtained coordinates of the three-dimensional virtual object grid points and adopting a ridge estimation method, thereby obtaining the rational function model.

According to a preferred embodiment of the present invention, calculating the parameters of the rational function model further comprises selecting a grid point in the input image as a check point; calculating the three-dimensional ground coordinates of each inspection point to each contour plane by using a strict geometric imaging model; calculating the corresponding image coordinate of the three-dimensional ground coordinate by using the obtained rational function model; and counting the fitting errors of the original image coordinates on the check points and the calculated corresponding image pixel coordinates to obtain the fitting precision of the rational function model to the strict geometric imaging model.

According to a preferred embodiment of the present invention, the area network adjustment process includes the steps of: converting the ground coordinates of the image grid points to image pixel coordinates using the following formula and rational function model

Wherein q is_i，0，a_i，1，a_i，2And b_i，0，b_i，1，b_i，2Is directed to 6 orientation parameters of image i, and (x)_k，y_k) And

is the image and ground coordinates of the equal altitude plane with the label k; and using the image pixel coordinate obtained by calculation and its actual measurement coordinate to make 2 translation parameters a_i，0And b_i， ₀Or all 6 parameters a_i，0，a_i，1，a_i，2And b_i，0，b_i，1，b_i，2The adjustment estimation is performed according to the following formula

v＝AΔ+l；P

v＝[v_x，v_y]_k ^T；

Wherein P is a weight matrix defined according to the image coordinate measurement accuracy.

According to a preferred embodiment of the present invention, the process of generating an approximate epipolar line image, defining a matching template, measuring a matching seed point line, automatically image matching and generating a digital elevation model comprises the steps of: preprocessing an image to be matched to reduce image noise and other image defects, adaptively enhancing image textures and linear features, and adaptively enhancing image shadows and image contrast of a region with less textures so as to perform feature extraction and image matching afterwards; based on a multi-scene image correlation matching algorithm, combining matching results of multiple matching primitives, utilizing local and global information of an image, and adopting a rough-to-fine image matching strategy to carry out influence matching; and performing fine matching on the matching result obtained by the multi-scene image multi-matching primitive matching module by adopting a least square matching algorithm and a matching algorithm capable of performing local fine matching on the landform details so as to improve the precision of the matching result and position and eliminate small matching gross errors.

According to a preferred embodiment of the present invention, the matching primitive includes feature points, feature lines, and lattice points.

According to a preferred embodiment of the present invention, the multi-scene image correlation matching algorithm comprises: selecting one scene image of the multi-scene images as a reference image I₀Using the rest of the images as search images I_iForming a corresponding image pair with the reference image; for a given feature point p on the reference image₀Obtaining the point p by the collinearity equation₀Three-dimensional coordinate P of₀(X₀，Y₀，Z₀) Wherein Z is₀Represents a point p₀A rough elevation value of (a); corresponding to the elevation Z₀- Δ Z and Z₀Two object points P of + Δ Z_minAnd P_maxIs taken of_minP_maxUsing the precise orientation element of the image to image the line segment P_minP_maxProjected on the search image to obtain a given point p₀A quasi-epipolar line on each image; at point p₀Defining a matching window W for the center, obtaining an automatic digital elevation model of the imaged area, and projecting the window onto the approximate automatic digital elevation model using the orientation elements of the reference image, such that the deformation caused by topography, image imaging geometry and imaging scale is atAutomatically obtaining compensation in the image correlation process; and obtaining the point p by the following formula₀In a reference image I₀And searching for images I_iNormalized correlation coefficient NCC between upper corresponding matching windows W_i

<math><mrow><msub><mi>NCC</mi><mi>i</mi></msub><mrow><mo>(</mo><msub><mi>p</mi><mn>0</mn></msub><mo>,</mo><mi>Z</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><munder><mi>Σ</mi><mrow><mi>i</mi><mo>&Element;</mo><mi>W</mi></mrow></munder><mrow><mo>(</mo><msub><mi>I</mi><mn>0</mn></msub><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>I</mi><mn>0</mn></msub><mo>)</mo></mrow><mo>×</mo><mrow><mo>(</mo><msub><mi>I</mi><mi>i</mi></msub><mrow><mo>(</mo><msub><mi>s</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>Z</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>I</mi><mi>i</mi></msub><mo>)</mo></mrow><mo>)</mo></mrow></mrow><mrow><msqrt><munder><mi>Σ</mi><mrow><mi>i</mi><mo>&Element;</mo><mi>W</mi></mrow></munder><msup><mrow><mo>(</mo><msub><mi>I</mi><mn>0</mn></msub><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>I</mi><mn>0</mn></msub><mo>)</mo></mrow><mn>2</mn></msup></msqrt><msqrt><munder><mi>Σ</mi><mrow><mi>i</mi><mo>&Element;</mo><mi>W</mi></mrow></munder><msup><mrow><mo>(</mo><msub><mi>I</mi><mi>i</mi></msub><mrow><mo>(</mo><msub><mi>s</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>Z</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>-</mo><msub><mi>I</mi><mi>i</mi></msub><mo>)</mo></mrow><mn>2</mn></msup></msqrt></mrow></mfrac></mrow></math>

<math><mrow><mi>Here</mi><mo>,</mo><msub><mi>I</mi><mn>0</mn></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mi>m</mi><mo>×</mo><mi>n</mi></mrow></mfrac><munder><mi>Σ</mi><mrow><mi>i</mi><mo>&Element;</mo><mi>W</mi></mrow></munder><msub><mi>I</mi><mn>0</mn></msub><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow><mo>;</mo></mrow></math>

<math><mrow><msub><mi>I</mi><mi>i</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mi>m</mi><mo>×</mo><mi>n</mi></mrow></mfrac><munder><mi>Σ</mi><mrow><mi>i</mi><mo>&Element;</mo><mi>W</mi></mrow></munder><msub><mi>I</mi><mi>i</mi></msub><mrow><mo>(</mo><msub><mi>s</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>Z</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>,</mo><mi>i</mi><mo>=</mo><mn>1,2</mn></mrow></math>

And are provided with

Performing image matching, wherein W and s are a matching window on the reference image and an image point in the window respectively; m and n are the size of the matching window W; s_iAnd (Z) is the image point corresponding to s on the search image.

According to the preferred embodiment of the present invention, the processing of matching images by using a rough-to-fine image matching strategy based on the multi-scene image correlation matching algorithm, in combination with the matching results of multiple matching primitives, using local and global information of the images, comprises the steps of: initializing a digital elevation model for matching in an iterative process; after the image is subjected to self-adaptive preprocessing and a multi-level pyramid image structure is generated, based on the matching result of different matching primitives on the original resolution image, a multi-scene image related matching algorithm is utilized to perform rough-to-fine step-by-step extraction and matching from a low-resolution image of the multi-level pyramid image structure; for the non-texture or less-texture area, taking the grid points as main matching features, and matching through a probability relaxation matching algorithm based on an object space for dense grid points; after the pyramid image registration of each level is finished, all the matched features form an automatic digital elevation model of an irregular triangular grid; and taking the automatic digital elevation model as a matching initial value for the high-level pyramid image matching, performing adaptive adjustment on matching parameters, and entering the conditions into a high-resolution image for registration to generate an automatic digital elevation model for fine matching.

According to another aspect of the present invention, there is provided a system for processing high resolution remote sensing image data using a rational function model, a block adjustment of satellite images with a large range of sparse regions or without ground control, automatic matching based on multi-baseline multi-matching features, and a distributed parallelization processing technique based on a high-speed network, comprising: the image preprocessing part is used for analyzing and processing the input original image data and calculating the parameters of the corresponding rational function model so as to obtain a stereoscopic image pair by using the calculated rational function model; the image orientation and area network adjustment part is used for measuring control points and connection points of original image data, orienting a single chip and performing area network adjustment processing on errors between stereoscopic image pixel data and the original image data so as to improve the fitting precision between the stereoscopic image pixel data and the original image data; an automatic digital elevation model extraction part used for generating an approximate epipolar line image, defining a matching template, measuring a matching seed point line, automatically matching the image and generating a digital elevation model; and an ortho-image and line map measuring part for generating an ortho-image, generating a contour line, performing single-sheet mapping and stereo measurement, and completing map product manufacturing by using the obtained data.

According to still another aspect of the present invention, there is provided a computer product for causing a computer to implement a method for processing high resolution remote sensing image data using a rational function model, satellite image block adjustment with a wide range of areas rare or without ground control, automatic matching based on multi-baseline multi-matching features, and distributed parallelization processing based on a high speed network, wherein the method comprises the steps of: analyzing and processing the input original image data, and calculating parameters of a corresponding rational function model so as to obtain a stereoscopic image pair by using the calculated rational function model; measuring control points and connection points of original image data, orienting a single chip, and performing block adjustment processing on errors between stereo image pixel data and the original image data so as to improve the fitting precision between the stereo image pixel data and the original image data; generating an approximate epipolar line image, defining a matching template, measuring a matching seed point line, automatically matching the images and generating a digital elevation model; and respectively generating a digital orthographic projection image and a digital line drawing, and finishing the manufacture of the product by using the obtained data.

According to the method and the system and the corresponding computer product, the photogrammetry processing of satellite images and aerial images such as IKONOS, WorldView-I/QuickBird, IRS-P5, SPOT-5, ALOS/PRISM and the like can be realized, and the production tasks of surveying and mapping products such as DLG, DEM/DSM, DOM and the like of the remote sensing images from aerial triangulation to corresponding scales can be completed.

The system and the method are based on the multi-baseline and multi-feature automatic matching technology, adopt automatic matching and fusion of multi-scene images (more than two scene images are matched simultaneously), for example, the automatic extraction of a digital elevation model in a terrain relief area or a terrain broken area can be efficiently solved by utilizing the technology of simultaneously matching three stereopair consisting of a SPOT-510 multiplied by 5 m resolution HRS and a 5 m resolution HRG image, and the algorithm and the strategy of various image features (such as feature points, grid points and linear features) are taken into consideration for local image information and global information, and adopt a global image matching algorithm based on a probability relaxation method to organically combine point features and linear features. The workload of manual editing of the automatically extracted three-dimensional information of the earth surface is reduced, and the operation efficiency is improved.

According to the system and the method, the key technology and the algorithm for mapping the novel large-breadth frame type aerial digital image (such as DMC, UltraCam-D, SWDC) based on multi-baseline processing can simultaneously ensure the precision and the reliability of image matching, greatly reduce the matching gross error, effectively solve the full-automatic extraction of DEM/DSM and the semi-automatic extraction of contour lines under complex terrain conditions (large-area arid and semi-arid broken landforms, deep cutting mountains, large-area deserts/deserts, city central areas and the like), and finally form a new data processing method and a software system which take the areas or the image frames as basic operation units, thereby fully utilizing the characteristics of high geometric/radiation resolution, high ground coverage overlapping rate and the like of the aerial digital image, solving the difficulties and the problems of the traditional aerial digital image mapping based on the stereopair image, such as the reduction of elevation measurement precision caused by a shorter baseline, The number of models is increased due to large ground coverage overlapping rate, the automation degree of aerial survey internal data processing is further improved, full-automatic DEM/DSM and DOM and semi-automatic extraction of contour line data and ground feature elements are realized, and the collection and the update of medium and large scale basic geographic information (including novel products such as true orthographic images and the like) are efficiently carried out.

According to the system and the method, key technologies such as contour line data adaptive filtering and smoothing based on terrain gradient and altitude difference analysis and important landform characteristics are adopted, semi-automatic extraction of contour lines of part of hills, most of mountains and mountain areas in the mapping area is carried out, and accordingly workload of interior data acquisition is reduced.

The high-speed network parallel distributed processing part comprises a high-speed storage network and related services, a parallel cluster computing system and a user operating environment.

The system supports distributed parallel processing based on a high-speed network, and is based on the research of the key algorithm for automated and intelligent processing of data such as the adjustment of a cluster computer system IBMHS21-8853(8 computing nodes, each computing node comprises 8-core CPU and 2 72GB hard disks) and the near-two-year no/rare control area network, multiple baselines, multiple matching characteristics matching and the like in order to enable a software system to have large-scale parallel processing capacity and larger data processing throughputBasically, distributed parallel computing based on loosely coupled parallel service middleware is realized, i.e. all the workstations (including PC and high-performance cluster computer-blade computer) interconnected in the local area network communicate and cooperate in a software mode, and the distributed processing work of the image data is completed together by a certain task scheduling strategy. The distributed parallel processing can not only reduce the workload of personnel, but also realize the high automation of the operation steps of image preprocessing, epipolar image generation, image matching, image orthorectification and the like, thereby improving the efficiency. The software system of the current version can reach 400,000-600,000Km for SPOT-5HRS/HRG high-resolution satellite remote sensing image data²Production capacity of digital elevation/surface model (DEM/DSM) and Orthophoto Map (DOM) of/month; for aerial imagery, taking the IBM HS21-8853 cluster computer as an example, which only uses 8 computing nodes, the processing speed of DEM/DSM and DOM can reach 150-.

Drawings

FIG. 1 is a flow chart of high resolution remote sensing image data processing;

FIG. 2 is a flow chart of solving RFM parameters using a SPOT-5 satellite imagery rigorous geometric imaging model;

FIG. 3 is a flow chart of solving RFM parameters using a satellite imagery rigorous geometric imaging model;

FIG. 4 is a flow diagram of DEM automatic extraction based on multi-baseline, multi-match feature image matching;

FIG. 5 is a diagram of a multi-view image correlation matching algorithm based on object-side geometric constraints;

fig. 6 shows a user interface of a cluster distributed parallel processing module in a software system of the high resolution remote sensing image data processing method according to the present invention; and

fig. 7 shows a high resolution remote sensing image data processing system according to the present invention.

Detailed Description

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

First, a high-resolution remote sensing image data processing method of the present invention is described below with reference to fig. 1. Fig. 1 shows a basic flow of high-resolution remote sensing image data processing. Referring to fig. 1, the high resolution remote sensing image data processing method of the present invention calculates Rational Function Model (RFM) parameters for an input image by data analysis of an original remote sensing image (pair) directly acquired by a satellite in step S101.

Using RFM, parameters of RFM, which is a substitute model of the strict geometric imaging model, must be calculated according to the imaging mechanism of the satellite images and the ephemeris parameters and attitude angle data of the satellite images. The calculated RFM must replace the original rigid imaging model with extremely high accuracy and should have a fit error of no more than 0.05 pixels (the polynomial described herein, and the method of calculating the RFM parameters, will be described in detail below with reference to fig. 2).

For satellite imagery that directly provides RFM model parameters (e.g., IKONOS, QuickBird, IRS-P5, etc.), the necessary data format conversion (consistent with the RFM model data format of IKONOS satellite imagery) should be performed. The system adopts an automatic matching technology based on multiple base lines and multiple matching features (feature points, grid points and feature lines), so that in order to conveniently extract the feature lines, facilitate image matching, point location measurement and block adjustment calculation, the original satellite image is generally subjected to image self-adaptive enhancement. Meanwhile, the color tone of the acquired original satellite image is dark in general, the contrast of the image is small, the interpretation and the subsequent data processing of the image are greatly influenced, and the original image must be enhanced by adopting a proper preprocessing algorithm before the digital orthophoto map is manufactured. The processed stereo image pair, ground control data, and auxiliary data (calculated RFM parameters) are used as system inputs.

The rfm (rational Function model) is a commonly used model for restoring the imaging geometry of remote sensing images, and maps an object in a three-dimensional space into pixels in a two-dimensional image space by a rational polynomial ratio Function. Here, the RFM is used for processing wide-range remote sensing image data. In RFM, the relationship between the image pixel coordinates and their ground three-dimensional coordinates is expressed using a rational function model, where all coordinates are normalized. For convenience of description, assume now that x_nAnd y_nIn order to normalize the coordinates of the image pixels,to normalize geographic coordinates (longitude, latitude, and ellipsoidal height), the mathematical expression for RFM is:

wherein,

(t is 1, 2, 3, 4) is a general polynomial, and each coordinate component in the polynomial isλ_n，h_nThe power of (a) does not exceed 3 at maximum, nor does the sum of the powers of each coordinate component exceed 3.

From equation (1), it can be seen that the co-linear equation familiar from the preferred embodiment of the present invention is, in fact, only a specific example of an RFM. Generally, in RFM, distortion caused by optical perspective projection is expressed as a first order polynomial, and imaging distortion caused by intrinsic errors such as atmospheric refraction, earth curvature, lens distortion, satellite-borne GPS/IMU can be approximated by a second order polynomial, and unknown distortion of other higher order parts can be partially simulated by a third order polynomial.

Some high resolution commercial remote sensing satellites such as IKONOS, IRS-P5, etc. only provide RFM model coefficients to the user, while SPOT-5HRS/HRG provides the position-to-earth data of the image via metadata files (DIMAP format files), which include CCD scanning frequency, CCD instantaneous imaging time, orbit and attitude parameters, linear array internal orientation parameters, etc. of the image, and these data can be used to construct a rigorous geometric model of the SPOT-5 image.

Typically, SPOT-5 satellite imagery products have a nominal accuracy (without control) of 50-150 meters. For the calculation of the SPOT-5 satellite imagery RFM parameters, a terrain-independent solution proposed by Tao., C.V. (see Tao, C.V. (2000), Semi-Automated Object Measurement Using Multi-Image Matching from Mobile mapping Image Sequences, PE & RS, Vol.67, No.12, pp.1347-1357) is used in accordance with a preferred embodiment of the present invention. The method is based on the idea that a group of virtual three-dimensional object grid points are established by utilizing an orbit attitude model of the SPOT-5 image to serve as control points, and then the parameters of the RFM, namely RPCs, are solved by adopting a ridge estimation method. By adopting the scheme, the RFM can realize high-precision fitting of the SPOT-5 image rigorous geometric imaging model, and further the subsequent photogrammetry processing of the image is completed.

FIG. 2 is a schematic diagram of the solution of RFM parameters using a SPOT-5 satellite imagery rigorous geometric imaging model. Referring now to FIG. 2, a method of resolving the RFM parameters is described, the method flow being shown in FIG. 3:

in step S301, the entire scene image is divided into m lines and n columns of grid images, and coordinates of (m +1) × (n +1) grid points are obtained. Here, the values of m and n may be determined by factors such as the size of the image. For example, the row or column pitch is typically around 200 pixels. Further, for example, m rows and n columns may all be equidistant.

In step S302, a rigorous geometric imaging model, i.e., an orbit attitude model, is constructed using the SPOT-5 orbit and attitude angle parameters.

In step S303, to ensure the fitting accuracy of RFM, the fluctuation range (Z) of the image coverage area is determined_min～Z_max) Dividing into k layers (generally k is more than 5) to obtain (k +1) equal height planes. And calculating the corresponding ground coordinate point of each grid point on each layer of equal altitude plane by using the SPOT-5 satellite rigorous geometric imaging model. Thus, the (m +1) × (n +1) × (k +1) three-dimensional virtual object grid points which are uniformly distributed in space are obtained. Similarly, the fluctuation range (Z) can be considered here_min～Z_max) Are uniformly divided into k layers.

In step S304, the RFM coefficients are solved, and RPC is iteratively solved by using the virtual control point coordinate pairs with sufficient number of generated spatial distributions and using a ridge estimation method.

In step S305, checking points (px, py) are selected from the image for accuracy check, and these checking points are generally located at the center of the grid for calculating the RFM. And calculating the three-dimensional ground coordinates of each check point to each elevation surface by using the SPOT-5 rigorous geometric imaging model, then calculating the corresponding image coordinates (px ', py') of the RFM according to the three-dimensional ground coordinates by using the RFM obtained by calculation, and counting the fitting errors (px-px ', py-py') on each check point to obtain the fitting accuracy of the RFM to the rigorous geometric imaging model.

Referring back to fig. 1, step S102 will be described. In step S102, block adjustment of the video is performed. According to the preferred embodiment of the invention, the high-resolution satellite image block adjustment technology with rare or no ground control in a large area, namely RFM-based satellite image orientation and block adjustment, is adopted.

The rational function model RFM is high-precision fitting of a strict geometric imaging model, and the block adjustment model can be realized by image space translation and affine transformation on the basis of RFM and model transformation. In practice, the RFM can be considered as being a result of recombining the parameters of its corresponding rigid imaging geometry model. Various errors present in the inner and outer orientation elements of the sensor platform can also cause errors in the RFM. If the resulting RFM is calculated based on the initial ephemeris parameters, attitude angle data and internal orientation parameters, the RFM must be refined using a number of control points, i.e., based on the RFM's block adjustment, in accordance with a preferred embodiment of the present invention.

Theoretically, the RFM can be applied to object coordinates in different coordinate systems, such as geocentric coordinates, geographic coordinates or any map projection coordinate system, and considering that the coverage area of satellite images is generally large, and especially the coverage area of images is larger when the integral regional net adjustment of multi-view SPOT-5 stereo images is performed, the geographic coordinates (i.e. longitude, latitude and geodetic height) under WGS84 are selected as the object coordinates (as shown in formula (1)) in the RFM according to the preferred embodiment of the present invention, which is beneficial to solving the problems such as the correction of the earth curvature and the crossing of the projection zone in the regional net adjustment process. In view of the above, a basic method for adjusting the block size of a large-scale satellite image is provided according to the preferred embodiment of the present invention.

The specific block adjustment method comprises the following steps: according to the study of Grodecki and Dial (Grodecki, J., Dial, G. (2003), Block Adjustment of High-Resolution software imagedescription by random Synthesis, PE & RS, Vol.69, No.1, pp.59-68), image orientation for RFM-based imaging models can be performed in image space as well as in object space. A preferred embodiment according to the invention employs an orientation method in image space, which can be described as follows (see equation (2)):

wherein, a_i，0，a_i，1，a_i，2And b_i，0，b_i，1，b_i，2Is for 6 orientation parameters of image i; (x)_k，y_k) And

is the image of the contour plane labeled k and the ground coordinates. In addition, the high-resolution linear array CCD sensor has the characteristics of high flying height, narrow imaging light beam and approximate parallel projection, and has strong correlation between orientation parameters, for example, the instantaneous imaging position of the sensor in the flying direction has strong correlation with the pitching angle of the sensor, and the mode and the size of an image coordinate error caused by the error of the parameters are basically consistent. Therefore, using the orientation model, the adjustment parameter b_i，0Will absorb the errors in the image line direction caused by the position and attitude errors of all the satellite-borne sensors in the flight direction, the adjustment parameter a_i， ₀Errors in the image column direction caused by position and attitude errors in the scanning direction of all the satellite-borne sensors are absorbed; since the rows of the image generally correspond to the flight direction of the satellite-borne sensor, the rows of the image are related to the instantaneous imaging time of each CCD linear array, and the adjustment parameter b_i，1And a_i，2Will absorb the image error caused by the drift error of the satellite-borne GPS and inertial navigation system, and the parameter a_i，1And b_i，2The image error caused by the error of the internal orientation parameter is absorbed.

According to the preferred embodiment of the present invention, the ground coordinates of the points are converted to image pixel coordinates using the RFM model using equation (2), and the calculated image pixel coordinates and its actual measured coordinates are used to pair 2 translation parameters a_i，0And b_i，0(orientation method IMG-2) or all 6 parameters a_i，0，a_i， ₁，a_i，2And b_i，0，b_i，1，b_i，2(orientation method IMG-6) adjustment estimation is performed. Mathematical models of these orientation methodsSee the following equation:

v＝AΔ+1；P

v＝[v_x，v_y]_k ^T；

wherein, P is a weight matrix defined according to the image coordinate measurement precision.

The adjustment model formed by the formulas (2) and (3) can be applied to single-image orientation of high-resolution satellite images, and can also be used for forming a regional network by using a plurality of satellite images so as to perform combined regional network adjustment on a plurality of images in an imaging area or image data of a plurality of different sensor platforms. It should be noted that the image orientation model based on image space constructed by the equations (2) and (3) is only suitable for high-resolution satellite images with narrow imaging beams. For satellite images with 1-2.5 m ground resolution such as IKONOS, Quickbird, SPOT-5HRG and the like, as the instantaneous imaging light beam is extremely narrow, for shorter flight paths (such as about 50-100 km), higher orientation precision can be achieved by using orientation methods IMG-2 and 1-4 well-defined ground control points; for SPOT-5HRS or longer IKONOS, Quickbird image strips, the orientation method IMG-6 is used for orientation, and more than 4-6 ground control points are needed on each image.

In step S103, an approximate epipolar line image is generated, and in step S104, image matching is performed and a digital elevation model DEM is generated. Aiming at the characteristic that almost all high-resolution satellite imaging systems are CCD linear array imaging, a set of high-precision matching algorithm specially aiming at linear array images is provided for automatic DEM extraction according to the preferred embodiment of the invention. The algorithm can generate dense, accurate and reliable image matching results. The algorithm adopts a Coarse-to-fine (Coarse-to-fine) multistage image matching strategy (step-by-step matching from a low-resolution image to a high-resolution image), integrates multiple mature image matching algorithms with complementary performance comprehensively, performs quality control among all sub-modules of the matching algorithm, automatically performs matching gross error positioning and elimination, and obtains a high-precision DEM (digital elevation model) of an imaging area by utilizing new characteristics (high signal-to-noise ratio, high-contrast image, high ground coverage overlapping rate and the like) provided by the high-resolution satellite image. The method can also be applied to traditional aerial photos and novel digital aerial images after being expanded.

The DEM automatic extraction method based on multi-baseline matching mainly comprises three mutually connected steps: self-adaptive preprocessing of images; matching multi-baseline and multi-matching primitive images; and high-precision landform detail refinement matching.

The operations of the above-described step S103 and step S104 in fig. 1 will now be described with reference to fig. 4. Fig. 4 illustrates a specific flow of step 103 and step 104 shown in fig. 1. In FIG. 4, the inputs to the method are the high resolution satellite imagery and its precise orientation parameters (obtained directly from the DGPS/IMU autonomous positioning system or from the block adjustment calculation).

In step S410, the step as the image adaptive preprocessing section is applied to both 8-bit and higher than 8-bit digital images, and it mainly preprocesses the images participating in matching, reduces image noise and other image defects, adaptively enhances image texture and linear features crucial to matching, adaptively enhances image contrast of image shadow and less-texture regions, and provides good image data for the subsequent feature extraction and image matching stages.

In step S420, the multi-scene image, multi-feature matching part is a core part of the algorithm, wherein the image matching adopts a matching algorithm based on multi-scene images (equal to or greater than two images) in combination with multiple matching algorithmsMatching results of the matching elements (feature points, feature lines and grid points) utilize local and global information of the images and adopt an image matching strategy from coarse to fine. Here, the conventional image Correlation matching algorithm (Cross-Correlation) is extended according to the preferred embodiment of the present invention, and a unique multi-scene image Correlation matching algorithm (GC: geometric Constrained Cross-Correlation) based on object-space geometric constraint, which can process multi-scene images simultaneously, is proposed³). Meanwhile, aiming at the analysis of the matching algorithm based on the surface element and the image characteristic point and line, the invention provides a combined matching algorithm which can combine the characteristic point, the characteristic line and the grid point for matching, can simultaneously utilize the local information and the global information of the image, so that the automatically matched characteristic point line can fully express the important characteristics of the landform and the landform, and can also carry out effective transition on the matching of the dense grid point in the non-texture or less-texture area by the probability relaxation method based on the image global geometric constraint.

The step S420 may be divided into sub-modules, such as feature point-based matching, feature line-based matching, mesh point matching based on image global geometric constraint, and matching feature fusion, where each sub-module provides an initial matching value and support, and provides a check to locate and reject a large matching gross error. Wherein the matching of the characteristic points and the characteristic lines is realized by GC³The algorithm is completed, and the matching of grid points adopts GC³The algorithm is carried out with an improved probability relaxation method, wherein the improved probability relaxation method is a probability relaxation matching algorithm based on an object space, the algorithm is an extension of the traditional probability relaxation method, and the algorithm can simultaneously process more than two images.

Here, the related matching algorithm-GC of multi-scene images based on the geometric constraint of object space³The basic principle of the method is as follows: GC³The algorithm is an extension of the traditional image correlation matching algorithm, abandons the traditional matching strategy based on image space, adopts the matching strategy based on image-object space relationship, applies the concept of matching multi-scene images guided by the constraint of object space geometric conditions, abandons the traditional concept of simultaneous matching based on stereo imageThe double-image matching algorithm of the image pair directly obtains the three-dimensional information of the features by simultaneously matching the multi-scene images, so that the reliability and the precision of the algorithm are simultaneously improved.

Fig. 5 is a schematic diagram of a multi-scene image correlation matching algorithm based on object-space geometric constraint. As shown in fig. 5, it is assumed that there are three IKONOS images according to the preferred embodiment of the present invention, which constitute an IKONOS three-dimensional image pair (Triplet Stereo). According to the preferred embodiment of the present invention, an intermediate scene image is selected as the reference image I₀The other two images are search images I_iAnd i is 1 and 2. They form two stereopair, i.e. image pair I₀-I₁And I₀-I₂. For a given feature point p on the reference image₀The corresponding point of the object space is positioned at the passing point p₀Of the photographing light Cp₀Upper (for CCD line image, C is the corresponding point p₀Instantaneous photography center).

First, suppose Z₀Is a point p₀According to a preferred embodiment of the present invention, the approximate elevation value of (A) may be calculated by a collinearity equation (corresponding to the camera light Cp)₀With the elevation plane Z₀Intersection) of the three-dimensional coordinates P of the acquisition points₀(X₀，Y₀，Z₀) Then, assume that a given rough elevation value Z is₀Is Δ Z, a tolerance corresponding to the elevation Z may be obtained according to a preferred embodiment of the invention₀- Δ Z and Z₀Two object points P of + Δ Z_minAnd P_maxAnd point p₀The corresponding point in the object space is positioned on the photographic light line segment P_minP_maxAbove. Then, using the precisely oriented elements of the image, the photographing light line segment P is divided according to the preferred embodiment of the present invention_minP_maxRespectively projected to the search image I_iTo obtain a given point p₀In image I_iQuasi-epipolar line (segment) above, and point p₀The matching point of (A) is also positioned in the I of the search image_iOn the upper quasi-epipolar line. Thus, GC³The algorithm implicitly utilizes the epipolar constraint conditions between the images。

Definition I₀(p) and I_i(p) is a gray-scale matrix of the image. At a given point p on the reference image₀A matching window W is defined for the center. Assuming that an approximate DEM of the imaged area (which may be an average elevation plane or a sketch DEM interpolated from the matching results of the higher order pyramid images) has been obtained according to the preferred embodiment of the present invention, the window is projected onto the sketch DEM using the orientation element of the reference image. In object space a small bin corresponding to the matching window W can be obtained according to a preferred embodiment of the invention.

Also, the small facet may be projected onto the search image to obtain a given point p according to a preferred embodiment of the present invention₀A corresponding matching window on the search image. This process is called Image-restoration. With Image-restoration, a rectangular matching window on the reference Image can be correlated with a non-rectangular window on the search Image. In this way, the deformations caused by topography, image imaging geometry and imaging scale can be automatically compensated for in the image correlation process.

Defining a point p₀(for a certain elevation value Z) in the reference image I₀And searching for images I_iNormalized correlation coefficient between the corresponding matching windows is NCC_i(formula 4). Wherein, W and s are respectively a matching window on the reference image and an image point in the window; m and n are the size of the matching window W; s_i(Z) is the Image point corresponding to s on the search Image, which can be obtained by Image-restoration process, and the gray value can be obtained by matching the search Image I_iAnd carrying out bilinear interpolation to obtain the target.

<math><mrow><msub><mi>I</mi><mi>i</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mi>m</mi><mo>×</mo><mi>n</mi></mrow></mfrac><munder><mi>Σ</mi><mrow><mi>i</mi><mo>&Element;</mo><mi>W</mi></mrow></munder><msub><mi>I</mi><mi>i</mi></msub><mrow><mo>(</mo><msub><mi>s</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>Z</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>,</mo><mi>i</mi><mo>=</mo><mn>1,2</mn><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow></math>

As can be seen from equation (4), unlike the conventional definition of normalized correlation coefficient, here NCC_iIs a picture point p₀And its elevation Z: z belongs to [ Z ]₀-ΔZ，Z₀+ΔZ]As a function of (c). Thus, a feature point p on the reference image is given₀And its corresponding approximate elevation Z₀And a tolerance Δ Z, for all stereo pairs (referred to herein as I)₀-I₁And I₀-I₂)，NCC_iAre defined in the object space, so that the preferred embodiment according to the invention adopts the literature^[6]The method as set forth in (1), normalized correlation coefficient NCC of all stereo pairs_iThe addition and averaging are performed, thereby defining SNCC, i.e.:

as can be seen from equation (5), SNCC is also the pixel p₀And Z: z belongs to [ Z ]₀-ΔZ，Z₀+ΔZ]Is the elevation Z at which the SNCC reaches a maximum, is exactly the given point p₀Of the correct elevation value, wherein point p₀The search range on the search image is determined by the tolerance Δ Z of the initial elevation.

Briefly, the algorithm exploits the fact that: a correct match exhibits a peak in the NCC curves of all image pairs, and corresponds to the true elevation value at a given point; whereas a false match (e.g., a false match caused by similar texture) may also exhibit a peak in the NCC curve of a single image pairHowever, due to differences in imaging geometry of each image pair, the peak positions in the NCC curves of different image pairs are different. Thus, by adding and averaging the NCC curves for each image pair, the correct match is further enhanced in the SNCC curve, while the peaks corresponding to the false matches cancel each other and are suppressed in the SNCC curve. Thus, by defining SNCC, it is possible to match multiple views simultaneously and directly obtain a given feature point p according to a preferred embodiment of the present invention₀Three-dimensional coordinates of (a). Equation (4) can also be easily extended to the case of n +1(n ≧ 1) scene images (equation (6)):

referring back to fig. 4, the specific flow of step S420 is shown in fig. 4. Step S421 is to match the initial digital elevation model in the iterative process.

In step S422, after adaptive preprocessing of the image and generation of a multi-level pyramid image structure (the image pyramid structure is stored and displayed at different resolutions according to the needs of the user under a uniform spatial reference to form a pyramid structure with resolution from coarse to fine and data size from small to large), three different types of matching elements, i.e., feature points, grid points and feature lines, are used to match the result on the image with the original resolution by using GC³The algorithm starts to carry out rough-to-fine step-by-step extraction and matching from a low-resolution image of a multi-step pyramid image structure.

In step S423, for the non-textured or less-textured region, the grid points are used as the main matching features, and matching is performed by using an object space-based probability relaxation matching algorithm for dense grid points.

In step S424, after the pyramid image registration of each level is completed, all the matching features form DEM data of a TIN (irregular triangular grid) structure (in the TIN, all the point features form nodes, and the line features form structural edges of the TIN).

In step S425, the DEM is further entered into the higher resolution image for registration as a matching initial value for the higher-level pyramid image matching and a condition for adaptive adjustment of matching parameters. After iteration of the matching process, the generated DEM is used for the fine matching part.

Next, step S430, which is a step of high-precision refinement matching, is described. The step adopts a simplified high-precision least square Matching algorithm and a special Matching algorithm (GRM) capable of carrying out local fine Matching on the geomorphic details, carries out fine Matching on a Matching result obtained by a multi-scene image multi-Matching element Matching module, improves the precision of the Matching result, and positions and eliminates smaller Matching gross errors. The part is characterized by having the function of local fine matching of the landform details. In general, general matching algorithms have the function of locally smoothing matching results, and GRM has an important meaning for maintaining fine features such as small gullies and linear features.

And the final digital elevation model is formed by interpolation on the basis of fusing all matched feature primitives.

Referring back to fig. 1, in step S105, a digital ortho-image (DOM) and a Digital Line Graph (DLG) are generated, respectively. The invention utilizes the DEM generated in the steps to carry out orthorectification on the panchromatic image and the multispectral image, selects multispectral wave band combination through tests according to the characteristics of the operation area, and adopts a proper algorithm to carry out fusion to form the image with the optimal color effect which truly reflects the earth surface characteristics. Because all images in the measuring area are unified to the same coordinate system after adjustment by the area network, the orthorectified multispectral image and the panchromatic image can be directly fused (or only a few registration points are needed). Finally, proper enhancement preprocessing (including deblurring, cloud and fog removal and the like) is respectively carried out on the panchromatic image and the multispectral image on the basis of not damaging the original color gradation, so that the fused image is bright in color and clear in detail texture characteristics.

The method adopts key technologies such as contour line data adaptive filtering, smoothing and the like based on terrain gradient and altitude difference analysis and important landform characteristics maintenance to perform semi-automatic extraction of contour lines in partial hills, most mountainous regions and mountain regions of the mapping region.

And (4) finishing the production of the topographic map, the image map, the thematic map and other thematic products in the step (S106) by using the obtained data.

In addition, in order to enable the invention to have large-scale parallel processing capacity and larger data processing throughput, the inventor realizes distributed parallel computing based on loosely coupled parallel service middleware on the basis of the research and the proposal of an automatic and intelligent processing key algorithm of a purchased cluster computer system IBM HS21-8853(8 computing nodes, each computing node comprises 8-core CPUs and 2 72GB hard disks) and data of no/rare control area network adjustment, multiple baselines, multiple matching characteristic matching and the like, namely all workstations (including PCs and high-performance cluster computers-blade computers) interconnected in a local area network communicate and cooperate in a software mode, and the image data processing work is completed together by a certain task scheduling strategy. The distributed parallel processing can not only reduce the workload of personnel, but also realize the high automation of the image orthorectification, thereby improving the efficiency.

The software system of the current version can reach 400,000-600,000Km for SPOT-5HRS/HRG high-resolution satellite remote sensing image data²Monthly digital elevation/surface model (DEM/DSM) and Orthophoto Map (DOM).

Fig. 6 shows a user interface of cluster distributed parallel processing in a software system for high resolution remote sensing image data processing according to the present invention. Referring to fig. 6, through the interface, a user may check resource configuration conditions of computing nodes (including a high-performance cluster computer (blade) and a general PC) operating in a high-speed network environment, including information of CPU usage, memory usage, and a limited space of a hard disk of each node; according to the resource allocation condition of each computing node, a user can send batch processing commands (including satellite image RFM model calculation, image format conversion and adaptive enhancement, epipolar line image generation, matching, orthoimage correction and the like) capable of being processed automatically, and the progress of each processing task can be monitored.

The remote sensing data processing is a type of data processing which is intensive in calculation and storage, therefore, according to the size of the data volume to be transmitted, according to the preferred embodiment of the invention, the remote sensing data processing task is divided into the calculation intensive processing (such as satellite image RFM model calculation and the like) and the storage intensive processing (such as image format conversion and adaptive enhancement, epipolar image generation, matching, orthoimage correction and the like), and the software system automatically distributes the data processing task according to the size of the data volume and the transmission speed of the network so as to achieve the highest data processing speed. Taking the example of an IBM HS21-8853 clustered computer using only 8 compute nodes, its processing speed can be increased 56 times for compute-intensive processing tasks (7 core CPUs for each compute node), and 16 times for storage-intensive processing tasks at gigabit speeds (2 core CPUs for each compute node).

A high-resolution remote sensing image data processing system according to the present invention is described below with reference to fig. 7. Fig. 7 shows a schematic configuration of a high-resolution remote sensing image data processing system according to the present invention. As shown in fig. 7, the system includes an image preprocessing section 701, an image orientation and block adjustment section 702, an automatic digital elevation model extraction section 703, a digital orthophoto and wiredrawing measurement section 704, and a high-speed network parallel distributed processing section 705.

The main image sources processed by the system are high-resolution (stereo) satellite images such as SPOT, IKONOS, Quickbird, IRS-P5, and the like. The relatively common tasks comprise SPOT-5HRS stereoscopic images (or strips), SPOT-5HRG images (panchromatic and multispectral waveband images), IKONOS co-rail/off-rail stereoscopic images or three stereoscopic image pairs, IKONOS and Quickbird full-color and multispectral images, IRS-P52.5 m resolution panchromatic waveband stereoscopic image pairs, ALOS/PRISM three-linear array stereoscopic images and the like, and the method can be respectively suitable for the tasks of manufacturing a 1: 50000 Digital Elevation Model (DEM) and a Digital Linear Graph (DLG), adjusting and drawing images and manufacturing a digital orthographic image (DOM).

Referring to fig. 7, the image preprocessing section 701 is used to perform data analysis and sorting on the original image, and perform parameter calculation of a rational function model to perform format conversion and image adaptive enhancement on the image data.

The satellite imaging system generally adopts a CCD linear array imaging technology, and a unified geometric imaging model, namely a Rational Function Model (RFM), can be calculated and adopted according to the characteristics of high flying height, narrow imaging light beam and approximate parallel projection of a high-resolution linear array CCD sensor. RFMs are applicable to a variety of sensors, including the latest aerospace and aeronautical sensors. The RFM is beneficial to the adjustment and orientation of the united area network of satellite images of various different sensor platforms and different ground resolutions. The system innovatively uses the RFM for processing the data of the remote sensing images in a large range, and supports the processing of mainstream optical sensor data including high-resolution satellite images (including Quickbird/WorldView-I, IKONOS, SPOT1-4, SPOT-5HRS/HRG, IRS-P5, Orbview, ALOS/PRISM and the like), traditional aerial image data and novel digital aerial images such as DMC and UCD; and can process satellite remote sensing images larger than 2G, such as Quickbird and WorldView-I.

The image orientation and area network adjustment portion 702 is used to perform control point and connection point measurements on the image, single slice orientation, and area network adjustment of the image.

The satellite imaging system generally flies at high altitude, the image coverage is wide, a large area can be controlled by a small number of control points on the ground, and the imaging requirement can be met. The system adopts a high-resolution satellite image area network adjustment technology with rare or no ground control in a large-range area. Tests show that for a long-strip SPOT-5HRS three-dimensional image with the coverage area of 7.2 ten thousand square kilometers, the requirement of 1: 50000 topographic map surveying and mapping can be met only by using 8-12 ground control points; for the IKONOS image with large area coverage, the requirement of 1: 10000 topographic map mapping only needs to meet the condition that a control point is arranged in the coverage range of each image, which undoubtedly greatly reduces the workload of field control.

The system also adopts a scheme of controlling free net adjustment or one side under the condition of no control in the border area, effectively solves the problem of difficulty in border area control and mapping, and also makes the mapping of the overseas topographic map possible. Experiments show that for SPOT-5HRS images, a small number of control points are utilized to carry out unilateral control, 600 kilometers can be overridden in the flight direction of a satellite, and the mapping requirements of hilly land, mountain land and high mountain land areas are still met by 1: 50000.

The automatic digital elevation model extraction part 703 is used to generate an approximate epipolar line image, define a matching template, measure a matching seed point line, perform automatic image matching, and generate a digital elevation model.

The system provides and develops a unique multi-baseline and multi-feature-based automatic matching technology, adopts the technology of automatically matching and fusing multi-scene images (simultaneously matching more than two scenes; for example, the technology of simultaneously matching three stereopair consisting of SPOT-510 x 5 m resolution HRS and 5 m resolution HRG images can effectively solve the automatic extraction of a digital elevation model in a terrain relief area or a terrain breakage area), and algorithms and strategies of various image features (such as feature points, grid points and linear features), gives consideration to image local information and global information, and adopts a probability relaxation-based global image matching algorithm to organically combine point features and linear features. The workload of manual editing of the automatically extracted three-dimensional information of the earth surface is reduced, and the operation efficiency is improved.

The key technology and algorithm for mapping the novel large-breadth frame type aerial digital image (such as DMC and UltraCam-D, SWDC) based on multi-baseline processing can simultaneously ensure the precision and reliability of image matching, greatly reduce the matching gross error, effectively solve the full-automatic extraction of DEM/DSM and the semi-automatic extraction of contour lines under complex terrain conditions (large-area arid and semi-arid broken landforms, deep cutting mountains, large-area deserts/deserts, urban central areas and the like), and finally form a new data processing method and a software system which take the areas or the image frames as basic operation units, thereby fully utilizing the characteristics of high geometric/radiation resolution ratio, high ground coverage overlapping rate and the like of the aerial digital image, solving the difficulties and problems of the traditional aerial digital image mapping based on stereopair image, such as the reduction of elevation measurement precision caused by a shorter baseline, The number of models is increased due to large ground coverage overlapping rate, the automation degree of aerial survey internal data processing is further improved, full-automatic DEM/DSM and DOM and semi-automatic extraction of contour line data and ground feature elements are realized, and the collection and the update of medium and large scale basic geographic information (including novel products such as true orthographic images and the like) are efficiently carried out.

The orthographic image and linear map measuring part 704 is used for generating an orthographic image, generating a contour line, performing single-chip mapping and stereo measurement, and using the obtained data to complete the production of a topographic map, an image map, a thematic map and other thematic products.

The system adopts key technologies such as contour line data adaptive filtering and smoothing based on terrain gradient and altitude difference analysis and important landform characteristics maintenance to semi-automatically extract contour lines in partial hills, most mountains and mountain areas of the mapping area, so that the workload of interior data acquisition is reduced.

The high-speed network parallel distributed processing portion 705 includes a high-speed storage network and associated services, a parallel cluster computing system, and a user operating environment.

The system supports distributed parallelization processing based on a high-speed network, and combines automation and intelligent of cluster computer systems IBMHS21-8853(8 computing nodes, each computing node comprises 8-core CPU and 2 72GB hard disks) and data of near-two-year no/rare control area network adjustment, multi-baseline and multiple matching characteristic matching and the like in order to enable a software system to have large-scale parallelization processing capacity and larger data processing throughputOn the basis of the research of a chemical processing key algorithm, distributed parallel computing based on loosely coupled parallel service middleware is realized, namely all workstations (including a PC (personal computer) and a high-performance cluster computer-blade computer) interconnected in a local area network are communicated and cooperated in a software mode, and the distributed processing work of image data is completed together by a certain task scheduling strategy. The distributed parallel processing can not only reduce the workload of personnel, but also realize the high automation of the operation steps of image preprocessing, epipolar image generation, image matching, image orthorectification and the like, thereby improving the efficiency. The software system of the current version can reach 400,000-600,000Km for SPOT-5HRS/HRG high-resolution satellite remote sensing image data²Production capacity of digital elevation/surface model (DEM/DSM) and Orthophoto Map (DOM) of/month; for aerial imagery, taking the IBM HS21-8853 cluster computer as an example, which only uses 8 computing nodes, the processing speed of DEM/DSM and DOM can reach 150-.

Finally, the user can use the obtained data to complete production and production of topographic maps, image maps, thematic maps and other thematic products according to requirements.

The method and the system thereof are mainly used for completing the work of multi-source satellite image preprocessing, field control/quick drawing preparation, adjustment of a large-range satellite image sparse control area network, automatic DEM generation, semi-automatic contour line data extraction, image fusion, high-precision orthoimage generation and the like.

Those skilled in the art will appreciate that the method can be embodied as computer readable codes recorded on a computer readable recording medium. The computer readable recording medium is any data storage device that can store data which can be read by a computer system. Examples of the computer readable recording medium include read-only memory (ROM), random-access memory (RAM), CD-ROMs, magnetic tapes, floppy disks, optical data storage devices, and carrier waves (such as data transmission through the internet). The computer readable recording medium can also be distributed over network coupled computer systems so that the computer readable code is stored and executed in a distributed fashion.

While the present invention has been described with reference to exemplary embodiments, it will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the spirit and scope of the invention as defined by the appended claims. The preferred embodiments should be considered in all respects only as illustrative and not restrictive. Therefore, the detailed description of the present invention does not limit the scope of the present invention, which should be defined by the appended claims, and all the distinguishing technical features within the scope of the present invention should be construed as being included in the present invention.

Claims

1. A method for processing high-resolution remote sensing image data by utilizing a rational function model, satellite image block adjustment with rare or no ground control in a large-range area, automatic matching based on multi-baseline and multi-matching characteristics and a distributed parallelization processing technology based on a high-speed network comprises the following steps:

1) analyzing and processing the input original image data, and calculating parameters of a corresponding rational function model so as to obtain a stereoscopic image pair by using the calculated rational function model;

2) measuring control points and connection points of original image data, orienting a single chip, and performing block adjustment processing on errors between stereo image pixel data and the original image data so as to improve the fitting precision between the stereo image pixel data and the original image data;

3) generating an approximate epipolar line image, defining a matching template, measuring a matching seed point line, automatically matching the images and generating a digital elevation model; and

4) and respectively generating a digital orthophoto map and a digital line map, and finishing the manufacturing of a map product by using the obtained data.

2. The method of claim 1, wherein computing the parameters of the rational function model comprises:

dividing input image data by a preset grid to obtain coordinates of grid points, and layering fluctuation ranges of image coverage areas to obtain preset height planes;

constructing a strict geometric imaging model, and calculating the corresponding ground coordinates of each grid point on each equal altitude plane according to the strict geometric imaging model to obtain three-dimensional virtual object grid points; and

and iteratively calculating coefficients of the rational function model by using the obtained coordinates of the three-dimensional virtual object grid points and adopting a ridge estimation method, thereby obtaining the rational function model.

3. The method of claim 2, further comprising the step of

Selecting a lattice point in an input image as a check point;

calculating the three-dimensional ground coordinates of each inspection point to each contour plane by using a strict geometric imaging model;

calculating the corresponding image coordinate of the three-dimensional ground coordinate by using the obtained rational function model; and

and (4) counting the fitting error between the original image coordinate at each check point and the calculated corresponding image coordinate to obtain the fitting precision of the rational function model to the strict geometric imaging model.

4. The method of claim 1, wherein the step 2) comprises:

converting the ground coordinates of the image grid points to image coordinates using the following formula and rational function model

Wherein a is_i，0，a_i，1，a_i，2And b_i，0，b_i，1，b_i，2Is directed to 6 orientation parameters of image i, and (x)_k，y_k) And (a)

λ_k，h_k) Is the image and ground coordinates of the equal altitude plane with the label k; and

using the image coordinate obtained by calculation and its actual measurement coordinate to make 2 translation parameters a_i，0And b_i， ₀Or all 6 parameters a_i，0，a_i，1，a_i，2And b_i，0，b_i，1，b_i，2The adjustment estimation is performed according to the following formula

v＝AΔ+l；P

v＝[v_x，v_y]_k ^T；

where

5. The method of claim 1, wherein the step 3) comprises:

3-1) preprocessing the image to be matched to reduce image noise and other image defects, adaptively enhancing image texture and linear features, and adaptively enhancing image shadow and image contrast of a less-texture area so as to perform feature extraction and image matching afterwards;

3-2) based on a multi-scene image correlation matching algorithm, combining matching results of multiple matching primitives, utilizing local and global information of the image, and adopting a rough-to-fine image matching strategy to perform image matching; and

and 3-3) performing fine matching on the matching result obtained by the multi-scene image multi-matching element matching module by adopting a least square matching algorithm and a matching algorithm capable of performing local fine matching on the landform details so as to improve the precision of the matching result and position and eliminate small matching gross errors.

6. The method of claim 5, wherein the matching primitives comprise feature points, feature lines, and grid points.

7. The method of claim 5, wherein the multi-view image correlation matching algorithm comprises:

selecting one scene image of the multi-scene images as a reference image I₀Using the rest of the images as search images I_iForming a corresponding image pair with the reference image;

for a given feature point p on the reference image₀Obtaining the point p by the collinearity equation₀Three-dimensional coordinate P of₀(X₀，Y₀，Z₀) Wherein Z is₀Represents a point p₀A rough elevation value of (a);

corresponding to the elevation Z₀- Δ Z and Z₀Two object points P of + Δ Z_minAnd P_maxIs taken of_minP_maxUsing the precise orientation element of the image to image the line segment P_minP_maxProjected on the search image to obtain a given point p₀A quasi-epipolar line on each image;

at point p₀Defining a matching window W for the center, obtaining an automatic digital elevation model of an imaging area, and projecting the window onto the approximate automatic digital elevation model by using the orientation elements of a reference image, so that the deformation caused by topographic relief, image imaging geometry and an imaging scale is automatically compensated in the image correlation process; and

the point p is obtained by the following formula₀In a reference image I₀And searching for images I_iNormalized correlation coefficient NCC between upper corresponding matching windows W_i

<math><mrow><msub><mi>NCC</mi><mi>i</mi></msub><mrow><mo>(</mo><msub><mi>p</mi><mn>0</mn></msub><mo>,</mo><mi>Z</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><munder><mi>Σ</mi><mrow><mi>s</mi><mo>&Element;</mo><mi>W</mi></mrow></munder><mrow><mo>(</mo><msub><mi>I</mi><mn>0</mn></msub><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>I</mi><mn>0</mn></msub><mo>)</mo></mrow><mo>×</mo><mrow><mo>(</mo><msub><mi>I</mi><mi>i</mi></msub><mrow><mo>(</mo><msub><mi>s</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>Z</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>-</mo><msub><mi>I</mi><mi>i</mi></msub><mo>)</mo></mrow></mrow><mrow><msqrt><munder><mi>Σ</mi><mrow><mi>s</mi><mo>&Element;</mo><mi>W</mi></mrow></munder><msup><mrow><mo>(</mo><msub><mi>I</mi><mn>0</mn></msub><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>I</mi><mn>0</mn></msub><mo>)</mo></mrow><mn>2</mn></msup></msqrt><msqrt><munder><mi>Σ</mi><mrow><mi>s</mi><mo>&Element;</mo><mi>W</mi></mrow></munder><msup><mrow><mo>(</mo><msub><mi>I</mi><mi>i</mi></msub><mrow><mo>(</mo><msub><mi>s</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>Z</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>-</mo><msub><mi>I</mi><mi>i</mi></msub><mo>)</mo></mrow><mn>2</mn></msup></msqrt></mrow></mfrac></mrow></math>

<math><mrow><mi>Here</mi><mo>,</mo><mi></mi><msub><mi>I</mi><mn>0</mn></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mi>m</mi><mo>×</mo><mi>n</mi></mrow></mfrac><munder><mi>Σ</mi><mrow><mi>s</mi><mo>&Element;</mo><mi>W</mi></mrow></munder><msub><mi>I</mi><mn>0</mn></msub><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow><mo>;</mo></mrow></math>

<math><mrow><msub><mi>I</mi><mi>i</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mi>m</mi><mo>×</mo><mi>n</mi></mrow></mfrac><munder><mi>Σ</mi><mrow><mi>s</mi><mo>&Element;</mo><mi>W</mi></mrow></munder><msub><mi>I</mi><mi>i</mi></msub><mrow><mo>(</mo><msub><mi>s</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>Z</mi><mo>)</mo></mrow><mo>)</mo></mrow></mrow></math>

i＝1，2

And are provided with

8. The method as claimed in claim 7, wherein the step 3-2) comprises:

initializing a digital elevation model in a matching manner in an iterative process;

after the image is subjected to self-adaptive preprocessing and a multi-level pyramid image structure is generated, based on the matching result of different matching primitives on the original resolution image, a multi-scene image related matching algorithm is utilized to perform rough-to-fine step-by-step extraction and matching from a low-resolution image of the multi-level pyramid image structure;

for the non-texture or less-texture area, taking the grid points as main matching features, and matching through a probability relaxation matching algorithm based on an object space for dense grid points;

after the pyramid image registration of each level is finished, all the matched features form an automatic digital elevation model of an irregular triangular grid; and

and taking the automatic digital elevation model as a matching initial value for the high-level pyramid image matching and performing adaptive adjustment on matching parameters, and entering the automatic digital elevation model into a higher-resolution image for registration to generate an automatic digital elevation model for fine matching.

9. A system for processing high-resolution remote sensing image data by utilizing a rational function model, satellite image block adjustment with rare or no ground control in a large range, automatic matching based on multi-baseline and multi-matching characteristics and a distributed parallelization processing technology based on a high-speed network comprises the following steps:

the image preprocessing part is used for analyzing and processing the input original image data and calculating the parameters of the corresponding rational function model so as to obtain a stereoscopic image pair serving as an input item by using the calculated rational function model;

the image orientation and area network adjustment part is used for measuring control points and connection points of original image data, orienting a single chip and performing area network adjustment processing on errors between stereoscopic image pixel data and the original image data so as to improve the fitting precision between the stereoscopic image pixel data and the original image data;

an automatic digital elevation model extraction part used for generating an approximate epipolar line image, defining a matching template, measuring a matching seed point line, automatically matching the image and generating a digital elevation model; and

and the orthographic image and linearized map measuring part is used for generating an orthographic image, generating a contour line, carrying out single-chip mapping and three-dimensional measurement and completing map product manufacturing by utilizing the obtained data.

10. The system of claim 9, wherein the image preprocessing section comprises:

a device for dividing the input image data by a predetermined grid to obtain the coordinates of each grid point and layering the fluctuation range of the image coverage area to obtain a predetermined number of altitude planes;

the device is used for constructing a strict geometric imaging model and calculating the corresponding ground coordinates of each grid point on each equal altitude plane according to the strict geometric imaging model so as to obtain the grid points of the three-dimensional virtual object; and

and the device is used for iteratively solving the coefficient of the rational function model by using the obtained coordinates of the three-dimensional virtual object grid points and adopting a ridge estimation method.

11. The system of claim 10, wherein the image preprocessing section further comprises:

means for calculating three-dimensional ground coordinates of each inspection point for each contour plane using a rigorous geometric imaging model based on the inspection points selected from the input images;

a device for calculating the corresponding image coordinate from the three-dimensional ground coordinate by using the obtained rational function model; and

and the device is used for counting the fitting errors of the original image coordinates on each check point and the calculated corresponding image coordinates so as to obtain the fitting precision of the rational function model to the strict geometric imaging model.

12. The system of claim 9, wherein the image orientation and area network adjustment portion comprises:

device for converting ground coordinates of image grid points into image coordinates by using following formula and rational function model

using the image coordinate obtained by calculation and its actual measurement coordinate to make 2 translation parameters a_i，0And b_i， ₀Or all 6 parameters a_i，0，a_i，1，a_i，2And b_i，0，b_i，1，b_i，2Means for performing adjustment estimation according to the following formula

v＝AΔ+l；P

v＝[v_x，v_y]_k ^T；

where

13. The system of claim 9, wherein the automatic digital elevation model extraction section comprises:

means for preprocessing the image to be matched to reduce image noise and other image artifacts, adaptively enhance image texture and linear features, and adaptively enhance image contrast in image shadows and less-textured areas for subsequent feature extraction and image matching;

the device is used for matching images by combining the matching results of various matching primitives and utilizing local and global information of the images and adopting a rough-to-fine image matching strategy based on a multi-scene image correlation matching algorithm; and

and the device is used for precisely matching the matching result obtained by the multi-scene image multi-matching element matching module by adopting a least square matching algorithm and a matching algorithm capable of performing local precise matching on the landform details so as to improve the precision of the matching result and positioning and eliminating small matching gross errors.

14. The system of claim 13, wherein the matching primitives comprise feature points, feature lines, and grid points.