Automatic detection and evaluation method for time-space geometric quality of remote sensing image

Abstract

According to the automatic detection and evaluation method for the space-time geometric quality of the remote sensing image, geometric quality detection calculation is carried out by fusing the absolute geometric positioning precision, the relative positioning precision, the wave band registration precision and the inter-slice splicing precision of the remote sensing image, and the space-time geometric quality of the detection image is evaluated and expressed; automatically extracting, splicing, cutting and loading DOM (Document Object Model) and DEM (Digital Elevation Model) data of multi-source reference images with different resolutions in the ground coverage range of the to-be-evaluated image according to the geographic position and resolution information of the detection image to obtain a required reference image; on-line detection of the space-time geometric quality of the remote sensing image is realized through automatic extraction of a reference image, automatic calculation and solution of an on-line detection index of the space-time geometric quality of the remote sensing image and automatic standardized output of a geometric quality detection result; a 1A-level image product geometric quality early warning mechanism is established, the image quality problem is found in time, automatic early warning of image products is achieved, and the distribution quality of the image products is guaranteed.

Description

Technical Field

[0001] The present application relates to a big data remote sensing image detection and evaluation method, and in particular to an automated detection and evaluation method for the spatiotemporal geometric quality of remote sensing images, and belongs to the technical field of remote sensing image detection. Background Art

[0002] Remote sensing technology has been well applied in various fields such as surveying and mapping, weather forecasting, environmental monitoring, crop assessment, transportation and urban construction. The quality of remote sensing image spatiotemporal geometry is directly related to the accuracy and reliability of the information obtained, and has a great impact on various subsequent applications of remote sensing images. How to detect the geometric quality of 1A-level remote sensing images during satellite operation, analyze the changing trend of the geometric quality of 1A-level image products in time and space, and provide basic data for on-orbit testing, image quality problem analysis, product image quality improvement, and satellite imaging state parameter adjustment for the remote sensing image ground processing system during the long-term operation of the satellite in orbit, and lay the foundation for technical analysis is the focus of remote sensing image quality detection.

[0003] In the process of forming remote sensing images, the linear array CCD used by optical satellites has linear distortions such as CCD rotation, translation, and scaling, as well as nonlinear distortions such as CCD bending due to the particularity of its own structure, materials, and manufacturing process. Similarly, due to the influence of various factors such as optical distortion of camera lenses, various random and systematic errors in observations on remote sensing platforms, and shooting environment, remote sensing images will have uneven rows and columns and deformed images.

[0004] The problems that need to be solved in the prior art multi-source vector-grating remote sensing image registration and the key technical difficulties of this application include:

[0005] (1) Remote sensing images have the characteristics of high resolution, large data volume and strong real-time performance. The existing geometric quality detection methods for detecting remote sensing images have many shortcomings: 1) The detection method is not systematic, automated and intelligent enough. Most of the traditional methods for detecting the quality of remote sensing images are manual or semi-automatic, which is not only highly subjective, but also consumes a lot of manpower and material resources and is inefficient. It is difficult to conduct comprehensive detection for massive remote sensing images. 2) The detection results are not quantified. Most of the general remote sensing image quality detection methods only objectively and qualitatively describe the quality of remote sensing images, and do not conduct quantitative detection of specific indicators. 3) The detection is not comprehensive enough. It only detects certain aspects or certain indicators without comprehensive detection, and the detection results may be relatively one-sided. Therefore, it is urgent to establish an effective remote sensing image spatiotemporal geometric quality automatic detection system to detect the geometric quality of remote sensing images, ensure the reliability and geometric quality of remote sensing images, and lay a solid foundation for subsequent image applications and further processing.

[0006] (2) In order to make the spatiotemporal geometry quality of Class 1A remote sensing images intuitively displayed, it is necessary to select representative spatiotemporal geometry quality detection indicators for remote sensing images. There are many spatiotemporal geometry indicators for remote sensing images, and the existing technology lacks reliable online detection indicators for spatiotemporal geometry quality of remote sensing images. When automatically detecting the spatiotemporal geometry quality of Class 1A remote sensing images, it is necessary to refer to auxiliary data. The existing technology cannot achieve scientific management and rapid automatic extraction of reference images. The existing technology detects the spatiotemporal geometry quality of remote sensing images manually or semi-automatically, which requires human participation and is time-consuming and laborious. Although some software can automatically match control points through image matching, some matching accuracy may be insufficient and the same-name points may be inaccurate. The existing technology lacks an online automatic detection algorithm for the spatiotemporal geometry quality of remote sensing images, and cannot achieve the effect of automatic detection of the spatiotemporal geometry quality of Class 1A remote sensing images.

[0007] (3) The existing tie point matching methods only consider the image processing perspective, generally using the grayscale calculation features and associated features of the image itself, and rarely consider the imaging geometric model of the remote sensing image. As a result, the registration method cannot take into account the complex geometric deformation between images. There are theoretical defects and it is impossible to fundamentally solve the high-precision matching problem. There are radiation and geometric differences between the reference images of different phases and non-homologous sources, resulting in a lack of reliable homonymous features on the image (especially in mountainous areas and urban areas). At the same time, due to the different imaging features, time, and resolution of heterogeneous images, there is an inconsistency in the extracted features (the extracted features exist in one image but not in another), making it difficult to define an effective matching measure.

[0008] (4) The existing technology lacks a method for geometric quality detection and calculation of the absolute geometric positioning accuracy, relative positioning accuracy, band registration accuracy, and inter-slice stitching accuracy of remote sensing images. It cannot objectively evaluate the spatiotemporal geometric quality of the detected images. It cannot automatically extract, stitch, crop, and load multi-source, different-resolution reference images DOM and DEM data within the ground coverage of the evaluated image based on the geographical location and resolution information of the detected image. It cannot realize online detection of the spatiotemporal geometric quality of remote sensing images by automatically extracting reference images, automatically calculating and solving online detection indicators of the spatiotemporal geometric quality of remote sensing images, and automatically standardizing and outputting geometric quality detection results. It cannot automatically detect the quality of 1A-level image products in large quantities, quickly, and accurately. The work efficiency of remote sensing image quality detection is low, and there is a lack of an early warning mechanism for the geometric quality of 1A-level image products. Summary of the invention

[0009] This application improves the efficiency of remote sensing image quality detection by performing online detection of the spatiotemporal geometric quality of Class 1A remote sensing images, and automatically detects the quality of Class 1A image products in large quantities, quickly, and accurately, thereby realizing intelligent, automated, real-time online monitoring of image product quality. An early warning mechanism for the geometric quality of Class 1A image products has been established to issue alarms for image products with abnormal quality or suspected abnormalities, detect image quality problems in a timely manner, realize automatic early warning of image products, and ensure the distribution quality of image products. By analyzing the changing trends of the geometric quality of Class 1A image products in the temporal and spatial directions, basic data is provided for on-orbit testing, image quality problem analysis, product image quality improvement, and satellite imaging status parameter adjustment of the remote sensing image ground processing system during the long-term operation of the satellite in orbit, laying a foundation for technical analysis.

[0010] In order to achieve the above technical effects, the technical solutions adopted in this application are as follows:

[0011] The automated detection and evaluation method of the spatiotemporal geometric quality of remote sensing images integrates the absolute geometric positioning accuracy, relative positioning accuracy, band registration accuracy, and inter-slice stitching accuracy indicators of remote sensing images to perform geometric quality detection calculations, evaluate the spatiotemporal geometric quality of the detected images, and automatically extract, stitch, crop, and load multi-source, different-resolution reference images DOM and DEM data within the ground coverage of the image to be evaluated based on the geographical location and resolution information of the detected images to obtain the required reference images; the online detection of the spatiotemporal geometric quality of remote sensing images is achieved through the automatic extraction of reference images, the automatic calculation and solution of online detection indicators of the spatiotemporal geometric quality of remote sensing images, and the automatic standardized output of geometric quality detection results;

[0012] Inter-slice stitching accuracy detection: First, match the connection points of the detection image and the reference image, then calculate the object coordinates corresponding to the image points by using the DEM data and the reference image RPC parameters, and then inversely calculate the image point coordinates by using the detection image RPC parameters and the object coordinates. Finally, the mean error of the image point calculation residual is used as the inter-slice stitching accuracy of the detection image, which includes:

[0013] Step 1: High-precision matching of connection points;

[0014] Step 2: Use the image coordinates of the tie points on the reference image and the RPC of the reference image to calculate the object coordinates of the tie points. Based on Taylor's formula and the initial value (L 0 ,B 0 ,H 0) Iteratively solve (L, B, H), where B, L, and H are the geodetic latitude, geodetic longitude, and geodetic height of the geodetic coordinate system, respectively; iteratively solve (L, B), obtain H through DEM, and repeat the above steps with (L, B, H) until the two (L, B) are less than the limit difference. Finally, the geodetic coordinates (L, B, H) solved above are transformed into Gaussian plane coordinates, and then transformed into orthophoto coordinates (x, y);

[0015] Step 3: Use the object coordinates of the connection points obtained in step 2 and the RPC of the detection image to inversely calculate the image point coordinates of the connection points on the detection image; RPC is an important file in space transformation mathematics, and the image point coordinates of the control points are solved by RPC inverse calculation;

[0016] Step 4: Calculate the residuals of the image point coordinates and the matching points with the same name, and make a difference between the calculated image coordinates and the matching coordinates to get the residuals of the image point coordinates. The calculation formula is Equation 12 and Equation 13:

[0017] ΔX=X refer -X Type 12

[0018] ΔY=Y refer -Y Formula 13

[0019] The mean error of the horizontal and vertical coordinates of multiple points of the same name in the image is calculated as the inter-slice stitching accuracy. The calculation formula is:

[0020]

[0021] m and n are the number of horizontal and vertical coordinates spliced ​​respectively.

[0022] Preferably, high-precision matching of spatiotemporal geometric connection points of remote sensing images:

[0023] The first step is to detect constant feature points in key areas: establish a multi-scale space, define the scale space of the two-dimensional image, and use Gaussian difference kernels of different scales and image convolution to generate and detect constant key areas;

[0024] The second step is spatial extreme point detection: each sampling point is compared with all its neighboring points to see whether it is larger or smaller than its neighboring points in the image domain and scale domain, and the extreme points in the scale space are found. The position and scale of the key domain are accurately determined by fitting a three-dimensional quadratic function, and the low-contrast key domain and unstable edge response points are removed at the same time.

[0025] The third step is to use constant feature descriptors in key domains:

[0026] 1) Key domain direction assignment: The gradient direction distribution characteristics of the key domain neighborhood pixels are used to specify the direction parameters for each key domain. The neighborhood window centered on the key domain is sampled, and the gradient direction of the neighborhood pixels is calculated using a histogram. Each key domain has three pieces of information: location, scale, and direction, thereby determining a key domain constant feature area;

[0027] 2) Feature point descriptor generation: Rotate the coordinate axis to the direction of the key domain to ensure rotation invariance;

[0028] The fourth step is to match the constant feature vector of the key domain: first calculate the similarity, and use the Euclidean distance as the similarity of the feature; obtain the potential match between the images through the similarity, and after obtaining the constant feature vector of the key domain, use the priority kd tree to perform a priority search to search for the 2 approximate nearest neighbor feature points of each feature point. Among these two feature points, if the closest distance divided by the second closest distance is less than a certain ratio critical value, then accept this pair of matching points; eliminate mismatches, and eliminate wrong matches according to geometric restrictions and other additional constraints;

[0029] Step 5: Eliminate gross errors: Let the point set composed of all matching points be P. First, randomly select n points from P to form a subset S. Use the point set S to establish the initial model M. Let the point set excluding S from P be C. Calculate the error between the point in C and M. If the error is less than the set critical value, add the point to the point set S to form a new point set S. * , if S * If the number of midpoints is greater than N, use the point set S * Recalculate the model M based on least squares * , regenerate S and repeat the above process. If no consistent set is found, the algorithm fails. Otherwise, the maximum consistent set obtained determines the internal and external points and the algorithm ends.

[0030] Preferably, the spatiotemporal geometric quality detection of remote sensing images: first, the detection image and the reference image are connected point matched, point a and point a' are the same-name image points after matching, and the detection 1A image has no projection information. First, point a' and the reference image RPC are used. RPC is the coefficient of the rational function model. RPC associates the geodetic coordinates D (Latitude, Longitude, Height) of the ground point with its corresponding image point coordinates d (line, sample) by a ratio polynomial, and calculates the object coordinates of the ground point a' corresponding to the ground point A. Then, the object coordinates of the surface point A and the RPC of the detection image are used to inversely calculate the image coordinates of point a on the detection image, that is, point a", and the image residuals △X and △Y of points a and a" are calculated, and the mean square error of the image point residuals of all connection points is calculated.

[0031] Preferably, the reference image is automatically extracted and loaded: from the existing reference image database, according to the geographical location and resolution information of the detection image, the multi-source and different resolution reference image DOM and DEM data within the ground coverage of the image to be evaluated are automatically extracted, spliced, cropped and loaded, and the reference image database is interacted with. After the image is automatically extracted from the database, it is also necessary to mosaic, splice and crop. The specific process is:

[0032] Step 1: Read the RPC parameter file of the detection image and calculate the image center resolution according to the RPC model:

[0033] 1) Use the image point coordinates of the image center point O and the RPC parameters to calculate the corresponding object plane coordinates (X O ,Y O );

[0034] 2) Measure the image points P and Q at N pixels in the X and Y directions (the larger the value of N, the more accurate the resolution of the calculation), and then use the RPC file to solve the plane coordinates of the corresponding object point (X P ,Y P ), (X Q ,Y Q );

[0035] 3) By GSD X =|X O -X P | / N, GSD Y =|Y O -Y Q | / N calculated resolution in X and Y directions;

[0036] 4) Take the maximum value of the resolution in the X and Y directions as the resolution of the image;

[0037] Step 2: Calculate the coverage of the image based on the image point coordinates and RPC parameters of the four corner points of the detected image: measure the coordinates of the four corner points of the remote sensing image, calculate the object coordinates of the four corner points based on the RPC parameters, and the quadrilateral formed by the four points is the coverage of the image on the ground;

[0038] Step 3: Search the image resolution table of the reference image database according to the resolution of the image to be evaluated, obtain the slice size of the DOM and the corresponding slice image table, and the slice size of the DEM and the corresponding slice image table at the resolution, and extract the reference image from the DOM slice image table and the DEM slice image table respectively.

[0039] Preferably, Step 4: Find the corresponding image in the reference image library according to the image coverage range, and perform mosaicing and cutting on all retrieved sliced images to obtain the reference image data including the coverage range of the image to be evaluated; the image mosaicing process includes color matching, image edge matching and edge correction preprocessing process, image automatic mosaicing, image mosaic line mosaicing, gray-based mosaicing and color-based mosaicing;

[0040] 1) Adaptive generation of seam line network

[0041] a) Voronoi diagram considering overlap: Let a set of surfaces A = {A 1 , A 2 ,..., A n} on a plane, where no other surface is contained in any one surface, that is and overlap is allowed between different surfaces. The distance from a point to a surface in the Voronoi diagram is defined with the non-overlapping part between two surfaces as the control element. Let surface A i and A j be any two different surfaces in surface set A, and the distance d j from point p to surface A i under the constraint of surface A a (p, A i , A j ) is defined as the minimum distance from point p to A', and is calculated using Equation 1:

[0042]

[0043] Its distance d i to surface A a (p, A i ) is defined as Equation 2:

[0044]

[0045] The Voronoi polygon of any surface A i is represented using Equation 3:

[0046]

[0047] The set of Voronoi polygons of all surfaces A 1 , A 2 ,…, A n is the Voronoi diagram of surface set A:

[0048] V = {V(A 1 ), V(A 2 ), V(A 3 ),…, V(A n )} Equation 4

[0049] b) Obtain the valid range of the image

[0050] Before generating the seam line network, obtain the valid range of each scene image and exclude the invalid pixel area. First, use the boundary tracking method to obtain the set of outer contour points of the valid range of the image, then use the Hough transform method to detect the straight edges of the outer contour, and then obtain an approximate quadrilateral of the valid range of the image based on the straight edges;

[0051] The set of outer contour points of the valid range is obtained by the boundary tracking method based on the 8-neighborhood;

[0052] c) Calculate the bisector between overlapping images: Refer to the calculation method of the medial axis of a convex polygon. The specific steps are as follows:

[0053] i. Number the polygon counterclockwise as P O , …, P N , draw the angle bisectors of each vertex angle, and mark the intersection of the angle bisectors of P O and P 1 as q 1 , and so on, mark the intersection of the angle bisectors of P N and P O as q N ;

[0054] ii. Calculate the distances d i (i = 0, 1, 2, …, N) from q i to its opposite side (i.e., P i+1 P 1 ), d 2 , …, d N ;

[0055] iii. Calculate d = min(d 1 , …, d N ), let d = d 1 , reorder the vertices counterclockwise so that the distance from q 1 to the opposite side is the shortest. If there are multiple vertices with the shortest distance to the opposite side, randomly select one of them as d 1 ;

[0056] iv. Starting from vertex P O , let axis = P 0 ; m = 1, n = N; respectively find the intersection points Point_m and Point_n of the angle bisector of P i and the angle bisectors of vertex angles P m and P n ;

[0057] v. If the distance d(Point_m, axis) ≤ d(Point_n, axis), then axis = Point_m, m++; if the distance d(Point_m, axis) > d(Point_n, axis), then axis = Point_n, n--;

[0058] Record the inflection point axis of the medial axis and number it as R 1 , R 2 ,…, R N ; When recording the inflection points of the medial axis, also record the corresponding vertex numbers connected to it, including R 1 corresponding to P 1 、P 2 , R 2 to R N-1 each corresponds to a vertex, R N corresponds to two vertices;

[0059] d) Generate Voronoi polygons

[0060] After obtaining the bisectors between images, divide the effective range. All images use the bisectors to implement the division of the effective range and then form Voronoi polygons;

[0061] 2) When automatically optimizing the seam line network

[0062] Assume that the Voronoi vertex is located in the n-degree overlapping region A. There are n scenes (n ≥ 3) of images that have a common overlapping area. The pixel (x, y) is a pixel in the n-degree overlapping region. Then the difference between the n scenes of images at this pixel is defined by Equation 5:

[0063]

[0064] In the formula, D ij (x,y) is the difference between image i and image j at the pixel (x, y), and its definition is as Equation 6:

[0065]

[0066] The optimized Voronoi vertex can be calculated using Equation 7:

[0067]

[0068] The individual seam lines are the Voronoi edges. Use the shortest path algorithm to solve this problem. Let image i be the detected image and image j be the reference image. Then the cost of each path is defined by Equation 8:

[0069] f(PS) = max D ij (x, y), (x, y) ∈ PS Equation 8

[0070] Seam line optimization is to find a path that minimizes f(PS) to achieve an efficient search for the minimum-cost path. The bisection method is adopted. Let the upper and lower limit values of the search path cost be g and h respectively (for 8-bit image data, the worst-case values are 0 and 255 respectively). The current search value z is the midpoint of the search interval, that is, z = (g + h) / 2. First, according to the starting point and the ending point, check whether there is a path with a cost of z in the overlapping area of the left reference image. If it exists, the upper limit value of the search path cost becomes z; if it does not exist, the lower limit value of the search path cost becomes z + 1. The maximum number of searches does not exceed log 2 (h-g) , for 8-bit image data, the minimum-cost path can be found after no more than 8 searches.

[0071] Preferably, 3) Image mosaicking based on the seam line network:

[0072] After obtaining the optimized seam line network, perform image mosaicking processing according to the seam line network to eliminate obvious seams and obtain the final seamless mosaic image. For the reference images that have an overlapping range with the detected image, use the generated effective mosaic polygons to sequentially read the pixel values of the image points in the polygons, and then assign the read gray values to the corresponding positions in the mosaic result image. Discard the invalid pixels outside the mosaic polygons directly. Finally, perform feathering processing on the seam lines to obtain a seamless mosaic image. When writing pixels for mosaicking, perform large-area uniform illumination processing according to the parameters calculated by the uniform illumination processing at the same time, that is, integrated uniform illumination mosaicking processing;

[0073] Seam line feathering processing method: First, judge the direction of the seam line on the mosaicked image and process it separately for different directions. The method for determining and processing the line segment direction on the seam line is as follows: If the slope of the seam line is greater than 1, it is determined that the seam line is in the vertical direction. At this time, calculate the gray difference between the left and right sides of the seam line. If the slope is less than or equal to 1, it is determined to be in the horizontal direction, then calculate the gray difference between the upper and lower sides of the seam line segment. Finally, distribute the calculated gray difference within the left and right or upper and lower sides in the vertical or horizontal direction of the seam line segment;

[0074] Step 5, Image cropping: Crop the reference image data according to the coverage range of the image to be evaluated, obtain and output the DOM and DEM reference images, crop the reference image that matches the size of the detected image, and expand it by 200 pixels * 200 pixels.

[0075] Preferably, Absolute geometric positioning accuracy detection: For the detected image and the corresponding reference image, through control point matching, first obtain a large number of control points, calculate the image point residuals corresponding to the control points, and then use the image point residuals and the image spatial resolution to calculate the mean value of the positioning error as the evaluation index for the absolute geometric positioning accuracy;

[0076] First, input the detected image and the automatically extracted reference images DOM and DEM. Perform control point matching on the detected image and the reference images, and eliminate the control points with large gross errors. Then, according to the coordinates of the matched control points and the RPC parameter file of the detected image, inversely calculate the corresponding image point coordinates of the control points on the detected image, calculate the position residuals of all control points of the detected image, and then use the image point residuals and the image spatial resolution to calculate the average positioning error as the absolute geometric positioning accuracy of the image.

[0077] Preferably, the specific steps of the absolute geometric positioning accuracy detection technology are as follows:

[0078] Step 1: Select the image to be measured, calculate the geographical range of the detected image, extract the reference image according to the geographical range, and perform projection conversion, mosaicking, and cropping on the reference image;

[0079] Step 2: Load the detected image and the corresponding reference image, perform control point matching on the detected image and the corresponding reference image, and remove the gross errors. The specific process of control point matching is as follows: Control point matching supports automatic matching between the detected image and the reference DOM and DEM data, realizes automatic measurement of control points, and outputs all control point object coordinates and image point coordinate information in a standard format. It supports the pyramid matching strategy, generates pyramid images at all levels, and performs relaxation matching from coarser to finer and from top to bottom level by level. The result of the previous level is used as the constraint of the next level, reducing the matching range and uncertainty at the same time, and finally obtaining a matching result that meets the requirements;

[0080] 1) Multi-level and multi-feature relaxation matching

[0081] Adopt a multi-level matching strategy from coarser to finer. The high-reliability image matching strategy uses a method that combines features and grayscale, a strategy of increasing the search range with a pyramid, and two-dimensional adaptive relaxation method for image matching;

[0082] In the relaxation method for image feature matching, due to the local smoothness of the terrain, correctly matched points have a larger neighborhood, and incorrectly matched points have a smaller neighborhood. For each extracted feature point, according to the similarity measure, find the alternative matching points. During the relaxation iteration process, the relaxation probability value of the correct alternative matching points continuously increases, and the relaxation probability value of the incorrect alternative matching points continuously decreases. When the relaxation probability value of the correct alternative matching points converges to 1 and the incorrect alternative matching points converge to 0, an accurate matching result is finally obtained. The calculation of the neighborhood adopts a combination of the eight-neighborhood of a regular grid and adjacent nodes of a triangular mesh.

[0083] 2) Extract feature points

[0084] a) First, select a region of size n×n, and then calculate the first-order difference by differentiating all pixel points in the region to obtain their gradients g in the x and y directionsx , g y ;

[0085] b) Take σ from 0.3 to 0.9 and perform Gaussian filtering on the obtained gradient g x , g y ;

[0086] c) Calculate the intensity value M according to Equation 9, where g x is the gradient in the x direction, g y is the gradient in the y direction, and G(s) is the Gaussian template:

[0087]

[0088] d) Sort the extreme points from largest to smallest and select the required number of extreme points from largest to smallest;

[0089] 3) Coarse error rejection: Automatically and reliably detect and reject the mis-matched connection points;

[0090] Step 3: The reference image is DOM, and calculate the object coordinates according to the image coordinates of the control points;

[0091] Step 4: Use the object coordinates obtained in Step 3 and the RPC file of the detection image to inversely calculate the image coordinates of the control points in the detection image. RPC is an important file for spatial transformation mathematics, and use RPC to inversely calculate the image coordinates of the control points;

[0092] Step 5: Calculate the residuals between the calculated image coordinates and the corresponding homologous points matched on the detection image, calculate the positioning error, and calculate the average error of all control point positioning errors as the absolute geometric positioning accuracy of the image;

[0093] Use Equation 10 to calculate the difference between the calculated image coordinates and the matched coordinates to obtain the residuals of the image points, and then calculate the positioning error:

[0094] dx = ΔX * GSD

[0095] dx = △Y * GSD

[0096]

[0097] where, ΔX = X 参考 - X 图像 , ΔY = Y 参考 - Y 图像 , calculate the point position residuals of all control points in the detection image, and then use the image point residuals and the image spatial resolution to calculate the mean value of the positioning error as the absolute geometric positioning accuracy of this scene image. The calculation formula is:

[0098]

[0099] The selected area size is n×n.

[0100] Preferably, for the detection of band registration accuracy: before performing the detection of band registration accuracy, first select the band used as a reference in the multi-spectral bands, and then obtain a large number of control points through control point matching for each band and the reference band, calculate the image point residuals corresponding to the control points, and through residual analysis, calculate the mean value of the residuals as the evaluation index of band registration accuracy;

[0101] The band registration accuracy characterizes the registration accuracy between each band of the multi-spectral image and the reference band. First, select the image of one band in the multi-spectral image as the reference image, and the remaining bands as the detection band images. Perform tie point matching on the reference band image and the detection band images respectively. Then, based on the DEM data, the image coordinates of the tie points on the reference band image, and the RPC parameter file of the reference band image, calculate the corresponding object coordinates of the tie points. Next, use the RPC parameter file of the detection image and the object coordinates to inversely calculate the corresponding image coordinates of the tie points on the detection image. Finally, solve the residuals between the inversely calculated image points and the image coordinates of the matching homologous points, and calculate the root mean square error of the image point residuals as the band registration accuracy of the detection image.

[0102] The algorithm flow for the detection of band registration accuracy:

[0103] Step i: Select the multi-spectral image, perform band separation on the multi-spectral image, select the image of one band as the reference band image, and the remaining bands as the detection band images;

[0104] Step ii: Perform high-precision tie point matching on the detection band images and the reference band image;

[0105] Step iii: Use the image coordinates of the tie points on the reference band image and the RPC of the reference band image to calculate the object coordinates of the tie points through forward calculation;

[0106] Based on Taylor's formula and the initial values (L 0 , B 0 , H 0 ), iteratively solve for (L, B, H). B, L, and H are respectively the geodetic latitude, geodetic longitude, and geodetic height in the geodetic coordinate system;

[0107] The iteratively solved (L, B), obtain H through the DEM, and repeat the above steps with (L, B, H) until the difference between the previous and current (L, B) is less than the tolerance. Finally, transform the calculated geodetic coordinates (L, B, H) into Gaussian plane coordinates and then into orthoimage coordinates (x, y);

[0108] Step iv: Using the object coordinates of the connection points obtained in Step ii and the RPC of the detection image, inversely calculate the corresponding image point coordinates of the connection points on the detection image. RPC is an important document in spatial transformation mathematics. Using the RPC of the reference image, inversely calculate the corresponding image point coordinates of the control points;

[0109] Step v: Calculate the obtained image point coordinates and the residuals of the matched homologous points. Calculate the difference between the obtained image coordinates and the matched coordinates to obtain the residuals of the image point coordinates. Calculate the image point residuals, and calculate the mean square error of the horizontal and vertical coordinate errors of multiple homologous points in the image as the registration accuracy between bands.

[0110] Preferably, for relative positioning accuracy detection: Based on the multi-spectral image and the corresponding panchromatic image, first obtain a large number of control points through control point matching, calculate the residuals of the corresponding image points of the control points, and through residual analysis, calculate the mean square error of the matching as the evaluation index of relative positioning accuracy;

[0111] First, select a multi-spectral or panchromatic image as the detection image, calculate the area range of the detection image, select the corresponding panchromatic or multi-spectral image as the reference image, perform connection point matching on the detection image and the reference image, then directly calculate the object coordinates corresponding to the connection points from the DEM data and the RPC parameters of the reference image, and then inversely calculate the image point coordinates of the connection points on the detection image from the RPC parameters and the object coordinates of the detection image. Finally, calculate the residuals and evaluate the accuracy;

[0112] Relative positioning accuracy detection algorithm process:

[0113] Step a: Input the detection image and the corresponding panchromatic or multi-spectral image as the reference image;

[0114] Step b: Perform high-precision connection point matching on the detection image and the reference image;

[0115] Step c: Use the image coordinates of the connection points on the reference image and the RPC of the reference image to directly calculate the object coordinates of the connection points; Based on the Taylor formula and the initial values (L 0 , B 0 , H 0 ) Iteratively solve for (L, B, H). B, L, and H are the geodetic latitude, geodetic longitude, and geodetic height in the geodetic coordinate system respectively. The iteratively solved (L, B), obtain H through the DEM, and repeat the above steps with (L, B, H) until the difference between the previous and current (L, B) is less than the tolerance. Finally, transform the above-solved geodetic coordinates (L, B, H) into Gaussian plane coordinates and then into orthoimage coordinates (x, y);

[0116] Step d: Use the RPC of the reference image to inversely calculate the image point coordinates: RPC is an important document in spatial transformation mathematics. Using RPC inverse calculation can solve the image point coordinates of the control points;

[0117] Step e: Calculate the residuals between the calculated image point coordinates and the coordinates of the homologous points on the detected image. The difference between the calculated image coordinates and the matched coordinates is the residual of the image point coordinates.

[0118] Step f: Calculate the mean square error of the horizontal and vertical coordinate errors of multiple homologous points in the image as the relative positioning accuracy of the detected image.

[0119] Compared with the prior art, the innovation points and advantages of this application are as follows:

[0120] (1) This application integrates the absolute geometric positioning accuracy, relative positioning accuracy, band registration accuracy, and inter - tile stitching accuracy indicators of remote sensing images for geometric quality detection and calculation, evaluates and expresses the spatio - geometric quality of the detected image. According to the geographical location and resolution information of the detected image, it automatically extracts, stitches, cuts, and loads multi - source and different - resolution reference image DOM and DEM data within the ground coverage range of the image to be evaluated to obtain the required reference images; through the automatic extraction of reference images, the automatic calculation and solution of the spatio - geometric quality online detection indicators of remote sensing images, and the automatic standardized output of geometric quality detection results, it realizes the online detection of the spatio - geometric quality of remote sensing images; it establishes a geometric quality early - warning mechanism for 1A - level image products, issues alarms for image products with abnormal or suspected abnormal quality, timely discovers image quality problems, realizes the automatic early - warning of image products, and ensures the distribution quality of image products.

[0121] (2) This application first studies the influencing factors of the spatio - geometric quality of remote sensing images, clarifies the causes of spatio - geometric deformation of remote sensing images, and studies the reasons for the spatio - geometric quality of remote sensing images from the root, laying a theoretical foundation for the next step. Then it sets representative spatio - geometric quality evaluation indicators for remote sensing images, and further studies the online automatic detection algorithm for the spatio - geometric quality of remote sensing images according to the selected geometric quality evaluation indicators, objectively and quantitatively detects the spatio - geometric quality of remote sensing images, ensures the reliability and accuracy of remote sensing images, and lays a solid foundation for the further processing of subsequent remote sensing images and the analysis of the on - orbit operation status of satellites. This application can automatically detect the quality of 1A - level image products in large quantities, quickly, and accurately, improve the work efficiency of remote sensing image quality detection, and realize the intelligent and automatic real - time online monitoring of the quality of image products.

[0122] (3) This application solves two key problems: First, the online detection of the spatio-temporal geometric quality of remote sensing images. By long-term monitoring and analysis, the geometric quality of 1A-level remote sensing images can reflect the on-orbit operation status of the remote sensing platform. At the same time, the quality of the spatio-temporal geometric quality of remote sensing images is directly related to the accuracy and reliability of the acquired information, and has a great impact on various subsequent applications of remote sensing images. Good geometric quality is the basis and prerequisite for remote sensing applications. Therefore, this application establishes an effective remote sensing image quality detection system to detect the geometric quality of remote sensing images, through representative indicators of the spatio-temporal geometric quality of remote sensing images, and by calculating these detection indicators, visually and quantitatively display the geometric quality of remote sensing images. Second, the automation of the online detection of the spatio-temporal geometric quality of remote sensing images: Aiming at problems such as the need for manual participation or semi-automation in remote sensing geometric quality detection, and difficulties in batch geometric quality detection of a large number of remote sensing images, through research on online automated detection algorithms and file configuration methods for the spatio-temporal geometric quality of remote sensing images, realize the automation of accurate calculation and detection of indicators representing the geometric quality of remote sensing images from the extraction of reference images to matching, give quantitative detection results, give qualitative evaluations based on the quantitative detection results, detect image quality problems in a timely manner, and provide basic data for the ground processing system of remote sensing images during the long-term on-orbit operation of satellites, on-orbit testing, analysis of image quality problems, improvement of product image quality, and adjustment of satellite imaging state parameters, laying a technical analysis foundation. BRIEF DESCRIPTION OF THE DRAWINGS

[0123] Figure 1 It is a flow chart for connecting point matching by the key domain constant method.

[0124] Figure 2 It is a schematic diagram of automatic extraction of control points based on boundary network feature points.

[0125] Figure 3 It is a technical flow chart for automatic extraction and loading of geometric reference images.

[0126] Figure 4 It is a schematic diagram of reference image data extraction.

[0127] Figure 5 It is a schematic diagram of the 8-neighborhood of the outer contour point set outside the effective range.

[0128] Figure 6 It is a schematic diagram of solving the bisector based on the central axis.

[0129] Figure 7 It is a schematic diagram of the seam line feathering processing method.

[0130] Figure 8 It is a schematic diagram of a mixed neighborhood adapted to uneven distribution of feature points.

[0131] Figure 9 It is a reference image for the extraction of GF1 poor-quality cloudy images.

[0132] Figure 10 It is a statistical chart of the absolute geometric positioning accuracy of GF1 cloudy images with poor quality.

[0133] Figure 11 It is a statistical chart of the detection results of the inter-sheet splicing accuracy of GF1 cloudy spectral images with poor quality. Specific implementation manner

[0134] The following further describes the technical solution of the automated detection and evaluation method for the spatio-temporal geometric quality of remote sensing images provided by this application with reference to the accompanying drawings, so that those skilled in the art can better understand this application and be able to implement it.

[0135] To detect and ensure the spatio-temporal geometric quality of remote sensing images and conduct objective and quantitative detection of the spatio-temporal geometric quality of remote sensing images, this application first studies the influencing factors of the spatio-temporal geometric quality of remote sensing images, clarifies the causes of spatio-temporal geometric deformation of remote sensing images, and studies the reasons for the quality of spatio-temporal geometric quality of remote sensing images from the root, laying a theoretical foundation for the next step. Then, representative spatio-temporal geometric quality evaluation indicators of remote sensing images are set, and then an online automatic detection algorithm for the spatio-temporal geometric quality of remote sensing images is studied according to the selected geometric quality evaluation indicators to conduct objective and quantitative detection of the spatio-temporal geometric quality of remote sensing images, ensuring the reliability and accuracy of remote sensing images, and laying a solid foundation for the further processing of subsequent remote sensing images and the analysis of the on-orbit operation status of satellites.

[0136] (1) Online detection indicators for the spatio-temporal geometric quality of remote sensing images: In order to intuitively display the spatio-temporal geometric quality of 1A-level remote sensing images, it is necessary to select representative spatio-temporal geometric quality detection indicators for remote sensing images. There are many spatio-temporal geometric indicators for remote sensing images. How to select representative geometric indicators from numerous geometric indicators requires in-depth research. By analyzing and studying each geometric indicator, representative detection indicators are selected to online detect the geometric quality of 1A-level remote sensing images, making the detection results more representative and accurate. By studying and referring to the experience of spatio-temporal geometric quality detection of foreign remote sensing images, the selected geometric quality detection indicators are as follows: absolute geometric positioning accuracy, relative positioning accuracy, band registration accuracy, and inter-sheet splicing accuracy.

[0137] (2) Automatic extraction of reference images: Auxiliary data is required for the automated detection of the spatio-temporal geometric quality of 1A-level remote sensing images. How to achieve the scientific management and rapid automatic extraction of reference images is a prerequisite for the automated detection of the geometric quality of remote sensing images. This application intends to conduct in-depth research on the scientific management and automatic extraction algorithm of reference images, and realize the automated extraction, mosaicking, cutting, and loading of multi-source and different-resolution reference images (DOM and DEM data) within the ground coverage of the image to be evaluated from the existing reference image database according to information such as the geographical location and spatial resolution of the detection image, so as to achieve the rapid automatic extraction of reference images.

[0138] (3) Online automated detection algorithm for the spatio-temporal geometric quality of remote sensing images: Generally, the detection of the spatio-temporal geometric quality of remote sensing images is either manual or semi-automatic, which requires manual participation and is time-consuming and laborious. Although some software can automatically match control points through image matching, there may be problems such as insufficient matching accuracy and inaccurate homologous points. This application intends to conduct special research on the online automated detection algorithm for the spatio-temporal geometric quality of remote sensing images. Through the configuration file method and program control, it realizes the automated detection of remote sensing image spatio-temporal geometric indicators from reference image extraction to matching and finally to remote sensing image, including absolute geometric positioning accuracy, relative positioning accuracy, inter-image mosaicking accuracy, and band registration accuracy, achieving the effect of automated detection of the spatio-temporal geometric quality of 1A-level remote sensing images.

[0139] I. High-precision matching of spatio-temporal geometric connection points of remote sensing images

[0140] The existing connection point matching methods only consider the image processing perspective, generally using the gray-scale calculation features and correlation features of the image itself, and rarely considering the imaging geometric model of remote sensing images. As a result, the registration method cannot well consider the complex geometric deformations between images, and there are theoretical defects, unable to fundamentally solve the problem of high-precision matching; there are radiation and geometric differences between reference images of different temporal phases and non-homologous sources, resulting in a lack of reliable homologous features on the images (especially more serious in mountainous areas and urban areas). At the same time, due to the different imaging features, time, and resolutions of heterogeneous images, the extracted features are inconsistent (the features extracted exist in one scene image but not in another scene image), making it difficult to define an effective matching measure. Therefore, high-precision connection point matching is a prerequisite for the online detection of the spatio-temporal geometric quality of remote sensing images.

[0141] Automatically match homologous image points with a certain density and meeting certain accuracy requirements for two or more scenes of images, and output the coordinate information of all homologous image points in the standard format. The key domain constant method is used for connection point matching. The flowchart of the connection point matching algorithm is as Figure 1 shown. The connection point matching process:

[0142] Step 1: Detection of Constant Feature Points in the Key Region: Establish a multi-scale space, define the scale space of the two-dimensional image, and generate the detection of the constant key region by convolving the image with Gaussian difference kernels of different scales.

[0143] Step 2: Detection of Spatial Extreme Points: Compare each sampling point with all its adjacent points to see if it is larger or smaller than its adjacent points in the image domain and scale domain, find the extreme points in the scale space, accurately determine the position and scale of the key region by fitting a three-dimensional quadratic function, and at the same time remove the key regions with low contrast and unstable edge response points to enhance the matching stability and improve the anti-noise ability.

[0144] Step 3: Descriptor of Constant Features in the Key Region

[0145] 1) Assignment of Key Region Directions: Use the gradient direction distribution characteristics of the pixels in the key region neighborhood to assign direction parameters to each key region. Sample within the neighborhood window centered on the key region, and calculate the gradient directions of the neighborhood pixels using a histogram. Each key region has three pieces of information: position, scale, and direction, thus determining a constant feature region of the key region.

[0146] 2) Generation of Feature Point Descriptors: Rotate the coordinate axes to the direction of the key region to ensure rotational invariance.

[0147] Step 4: Matching Constant Feature Vectors of the Key Region: First, calculate the similarity, using the Euclidean distance as the similarity of the features; obtain the potential matches between images through the similarity. After obtaining the constant feature vectors of the key region, use the priority k-d tree for priority search to search for the 2 approximate nearest neighbor feature points of each feature point. Among these two feature points, if the ratio of the nearest distance to the second-nearest distance is less than a certain proportional threshold, then accept this pair of matching points; after eliminating the mismatches, eliminate the incorrect matches according to geometric constraints and other additional constraints to improve the robustness.

[0148] Step 5: Rejection of Gross Errors: Let the point set composed of all matching points be P. First, randomly select n points from P to form a subset S. Use the point set S to establish an initial model M. Let the point set in P excluding S be C. Calculate the error between the points in C and M. If the error is less than the set critical value, then add this point to the point set S to form a new point set S * , if the number of points in S * is greater than N, then use the point set S * to recalculate the model M according to the least squares method * , regenerate S, and repeat the above process. If no consistent set is found, the algorithm fails; otherwise, judge the inliers and outliers using the largest consistent set obtained, and the algorithm ends.

[0149] II. Detection of Spatiotemporal Geometric Quality of Remote Sensing Images

[0150] Such as Figure 2As shown in the figure, the left side in the figure represents the detection image, and the right side represents the reference image. First, the connection points of the detection image and the reference image are matched. Point a and point a' are the homologous image points after matching. The detection 1A image has no projection information. First, use point a' and the RPC of the reference image. RPC is the coefficient of the rational function model. RPC associates the geodetic coordinates D (Latitude, Longitude, Height) of the ground point with its corresponding image point coordinates d (line, sample) using a ratio polynomial, calculate the object space coordinates of the ground point A corresponding to the ground point a', and then use the object space coordinates of the surface point A and the RPC of the detection image to inversely calculate the image space coordinates of point a on the detection image, that is, point a'', calculate the image space residuals △X and △Y of point a and point a'', and calculate the mean square error of the image point residuals of all connection points.

[0151] III. Automatic Extraction and Loading of Reference Images

[0152] Automatically extract, splice, cut, and load the DOM and DEM data of multi-source and different-resolution reference images within the ground coverage range of the image to be evaluated according to the geographical location and resolution information of the detection image from the existing reference image database. After interacting with the reference image database and automatically extracting the image from the database, it is also necessary to mosaic, splice, and crop. The specific process is as Figure 3 shown:

[0153] Step 1: Read the RPC parameter file of the detection image and calculate the resolution of the image center point according to the RPC model:

[0154] 1) Use the image point coordinates of the image center point O and the RPC parameters to calculate its corresponding object space point plane coordinates (X O , Y O );

[0155] 2) Measure the image points P and Q at N (the larger the value of N, the more accurate the calculated resolution) pixels in the X and Y directions, and then use the RPC file to solve their corresponding object space point plane coordinates (X P , Y P ), (X Q , Y Q );

[0156] 3) Calculate the resolutions in the X and Y directions from GSD X = |X O - X P | / N, GSD Y = |Y O - Y Q | / N;

[0157] 4) Take the maximum value of the resolutions in the X and Y directions as the resolution of the image;

[0158] Step 2: Calculate the coverage range of the image based on the image point coordinates of the four corner points of the detected image and the RPC parameters: Measure the coordinates of the four corner points of the remote sensing image, calculate the object space coordinates of the four corner points according to the RPC parameters, and the quadrilateral formed by the four points is the coverage range of the image on the ground;

[0159] Step 3: Retrieve in the image resolution table of the reference image database according to the resolution of the image to be evaluated, obtain the tile size of the DOM and the corresponding tile image table, and the tile size of the DEM and the corresponding tile image table at this resolution, and extract the reference images from the DOM tile image table and the DEM tile image table respectively;

[0160] Step 4: Find the corresponding images in the reference image library according to the image coverage range, perform mosaicking and cropping on all the retrieved tile images to obtain the reference image data including the coverage range of the image to be evaluated. The schematic diagram of image extraction is as Figure 4 shown; The image mosaicking process includes color matching, image edge-matching area matching and edge-correction preprocessing process, image automatic mosaicking, image mosaic line mosaicking, gray-based mosaicking and color-based mosaicking.

[0161] 1) Adaptive generation of seam line network

[0162] a) Voronoi diagram of plane considering overlap: Let a plane face set A = {A 1 ,A 2 ,…,A n}, no other face is included by any one face, that is and overlap is allowed between different faces. The distance between a point and a face in the Voronoi diagram is defined with the non-overlapping part between two faces as the control element. Let face A i and A j be any two different faces in face set A, and the distance d j (p, A i ) from point p to face A a is defined as the minimum distance from point p to A', and is calculated using Equation 1: i ,A j ) is defined as Equation 2:

[0163]

[0164] Its distance d i to face A a (p, A i ) is defined as Equation 2:

[0165]

[0166] The Voronoi polygon of any face A i is represented using Equation 3:

[0167]

[0168] All faces A 1 , A 2 ,…, A n The set of Voronoi polygons of A

[0169] V = {V(A 1 ), V(A 2 ), V(A 3 ),…, V(A n )} Equation 4

[0170] b) Obtain the valid range of the image

[0171] Before generating the seam line network, obtain the valid range of each image and exclude the invalid pixel area. First, use the boundary tracking method to obtain the set of outer contour points of the valid range of the image, then use the Hough transform method to detect the straight edges of the outer contour, and then obtain an approximate quadrilateral of the valid range of the image based on the straight edges;

[0172] The boundary tracking method is used to obtain the set of outer contour points of the valid range, which is based on the 8-neighborhood, and the neighborhood situation is as Figure 5 shown

[0173] When using the 8-neighborhood to define the valid range, the boundary is defined as the invalid pixels in the 4-neighborhood, and the steps are as follows:

[0174] i. First, scan the image in the clockwise direction to find the untracked boundary points as the starting points for tracking. Then, take d = 5 to start tracking. If no untracked points are detected, the operation ends. d is the serial number of the 8-neighborhood pixel points, which is used to represent the direction;

[0175] ii. Start looking for 8-neighborhood pixel points in the counterclockwise direction from d. If the change from the invalid pixel point to the valid pixel point occurs at the next boundary point position d, then go to step 3) for processing. If no valid pixels are found among the 8-neighborhood pixel points, the starting point of the tracking is an isolated point, and the tracking ends;

[0176] iii. Move to the next boundary point P n . If P n-1 = P 0 , P n = P 1 , the tracking ends; otherwise, d = (d + 3) % 8 + 1, and return to step b), where % is the modulo operation;

[0177] Obtaining the outer contour point set of the effective range only involves the tracking of the outer boundary and does not involve the tracking of the inner boundary. The quadrilateral determining the effective range is performed based on the obtained outer contour point set of the effective range. For the outer contour point set, it is simplified using the Douglas-Peuker algorithm to obtain the effective range of the image;

[0178] c) Calculate the bisector between overlapping images: Referring to the calculation method of the medial axis of a convex polygon, the specific steps are as follows:

[0179] i. Number the polygon counterclockwise as P O , …, P N , and draw the angle bisectors of each vertex angle. The intersection of the angle bisectors of P O and P 1 is denoted as q 1 , and so on. The intersection of the angle bisectors of P N and P O is denoted as q N ;

[0180] ii. Calculate the distances d i (i = 0, 1, 2, …, N) from q i P i+1 ) to its opposite side (i.e., P 1 , d 2 , …, d N ;

[0181] iii. Calculate d = min(d 1 , …, d N ), let d = d 1 , and reorder the vertices counterclockwise so that the distance from q 1 to the opposite side is the shortest. If there are multiple vertices with the shortest distance to the opposite side, randomly select one of them as d 1 ;

[0182] iv. Starting from vertex P O , let axis = P 0 ; m = l, n = N; respectively find the intersection Point_m of the angle bisector of P i and the angle bisector of vertex angle P m and P n ; Point_n;

[0183] v. If the distance d(Point_m, axis) ≤ d(Point_n, axis), then axis = Point_m, m++; if the distance d(Point_m, axis) > d(Point_n, axis), then axis = Point_n, n--;

[0184] Record the inflection point axis of the central axis and number it as R 1 , R 2 ,…, R N ; When recording the inflection points of the central axis, record the corresponding vertex numbers connected to it, including R 1 corresponding to P 1 , P 2 , R 2 to R N-1 each corresponding to a vertex, and R N corresponding to two vertices;

[0185] Calculate the angular bisector teml of the included angle of the line segments P m P m+1 , P n P n-1 (if n = N, then n - 1 = 0), and calculate the intersection points of the angular bisector teml and the vertex angle p m+1 , p n and denote them as Point_m and Point_n;

[0186] i. Loop until m = n, that is, all inflection points are obtained;

[0187] ii. Denote the two intersection points between the sides of the polygon in the effective range of adjacent images as startPoint and endPoint respectively;

[0188] iii. Traverse the vertex numbers corresponding to each inflection point to find the inflection points R i , R j corresponding to startPoint and endPoint;

[0189] If i = j, then the angular bisector is the line connecting startPoint, R i , endPoint;

[0190] If i < j, then the angular bisector is the line connecting startPoint, R i , R j+1 ,…, R j , endPoint;

[0191] If i > j, then the angular bisector is the line connecting startPoint, R i , R j+1 ,…, R j , endPoint;

[0192] As Figure 6 shown, the broken line segment between startPoint and endPoint in the dotted line is the required angular bisector;

[0193] d) Generating Voronoi polygons

[0194] After obtaining the bisectors between images, divide the effective range. All images use the bisectors to implement the division of the effective range and then form Voronoi polygons. The specific processing includes:

[0195] i. Calculate the bisectors between overlapping images;

[0196] ii. For each image, use the bisectors between the images with overlapping regions to crop its effective range. The result of each cropping is used as the input data for the next cropping operation. In this way, the effective range of an image is continuously divided and finally forms a polygon, which is the Voronoi polygon to which the image belongs;

[0197] iii. Perform the operations step by step for all images to obtain the polygons of the images themselves. Finally, the effective ranges of all images are divided into non-overlapping polygons, and the Voronoi diagram is finally obtained.

[0198] 2) When automatically optimizing the seam line network

[0199] Assume that the Voronoi vertex is located in the n-degree overlapping region A. There are n images (n≥3) that have a common overlapping area. The pixel (x, y) is a pixel in the n-degree overlapping region. Then the difference between the n images at this pixel is defined by Equation 5:

[0200]

[0201] In the formula, D ij (x,y) is the difference between image i and image j at the pixel (x, y), and its definition is as Equation 6:

[0202]

[0203] The optimized Voronoi vertex can be calculated using Equation 7:

[0204]

[0205] The individual seam lines are the Voronoi edges. The shortest path algorithm is used to solve this problem. Let image i be the detected image and image j be the reference image. Then the cost of each path is defined by Equation 8:

[0206] f(PS)=maxD ij (x,y), (x,y)∈PS Equation 8

[0207] Seam line optimization is to find a path that minimizes f(PS), enabling efficient search for the minimum-cost path. The bisection method is adopted. Let the upper and lower limit values of the search path cost be g and h respectively (for 8-bit image data, the worst-case values are 0 and 255 respectively). The current search value z is the midpoint of the search interval, that is, z = (g + h) / 2. First, based on the starting point and the ending point, check whether there is a path with a cost of z in the overlapping area of the left reference image. If it exists, the upper limit value of the search path cost becomes z; if it does not exist, the lower limit value of the search path cost becomes z + 1. The maximum number of searches does not exceed log 2 (h-g) , for 8-bit image data, the minimum-cost path can be found within no more than 8 searches.

[0208] 3) Image mosaicking based on the seam line network

[0209] After obtaining the optimized seam line network, perform image mosaicking processing according to the seam line network to eliminate obvious seams and obtain the final seamless mosaic image. For the reference images extracted with overlapping ranges with the detected images, use the generated effective mosaic polygons to sequentially read the pixel values of the image points in the polygons, and then assign the read gray values to the corresponding positions in the mosaic result image. Discard the invalid pixels outside the mosaic polygons directly. Finally, perform feathering processing on the seam lines to obtain a seamless mosaic image. When writing pixels for mosaicking, perform large-area illumination equalization processing according to the parameters calculated by the illumination equalization processing at the same time, that is, integrated illumination equalization and mosaicking processing, to improve processing efficiency.

[0210] Seam line feathering processing method: First, judge the direction of the seam line on the mosaicked image and process it separately for different directions. The determination and processing method for the line segment direction on the seam line is as follows: If the slope of the seam line is greater than 1, it is determined that the seam line is in the vertical direction. At this time, calculate the gray difference between the left and right sides of the seam line. If the slope is less than or equal to 1, it is determined to be in the horizontal direction, then calculate the gray difference between the upper and lower sides of the seam line segment. Finally, distribute the calculated gray difference within the left and right or upper and lower sides in the vertical or horizontal direction of the seam line segment, as Figure 7 shown.

[0211] Step 5, Image cropping: Crop the reference image data according to the coverage range of the image to be evaluated, and obtain and output the DOM and DEM reference images. Crop the reference image that matches the size of the detected image and expand it by 200 pixels * 200 pixels.

[0212] IV. Detection of absolute geometric positioning accuracy

[0213] Detect the positioning accuracy of image pixels and control points. By analyzing the variation law of this accuracy index over a long period, it reflects the changes in the installation relationship between the camera and star sensor devices during the on-orbit operation of the satellite, as well as the time drift characteristics of the measurement errors of the star sensor, gyroscope, and GPS sensor. This is of great significance for constructing a time extrapolation model of the spatio-temporal geometric error of satellite images and improving the geometric positioning accuracy of images without control. In this application, for the detected image and the corresponding reference image, through control point matching, a large number of control points are first obtained, and the pixel residuals corresponding to the control points are calculated. Then, using the pixel residuals and the image spatial resolution, the mean positioning error is calculated as the evaluation index for absolute geometric positioning accuracy;

[0214] First, input the detected image and the automatically extracted reference image DOM, DEM. Perform control point matching on the detected image and the reference image, and eliminate the control points with large gross errors. Then, according to the coordinates of the matched control points and the RPC parameter file of the detected image, calculate the pixel coordinates corresponding to the control points on the detected image, calculate the position residuals of all control points of the detected image, and then use the pixel residuals and the image spatial resolution to calculate the mean positioning error as the absolute geometric positioning accuracy of this image.

[0215] Specific steps of the absolute geometric positioning accuracy detection technology:

[0216] Step 1: Select the image to be measured, calculate the geographic range of the detected image, extract the reference image according to the geographic range, and perform projection transformation, mosaicking, and clipping on the reference image;

[0217] Step 2: Load the detected image and the corresponding reference image, perform control point matching on the detected image and the corresponding reference image, and eliminate gross errors. The specific process of control point matching is as follows: Control point matching can support the automatic matching between the detected image and the reference DOM, DEM data, realize the automatic measurement of control points, and output all control point ground coordinates and pixel coordinate information in the standard format. It supports the pyramid matching strategy, generates pyramid images at all levels, and performs relaxation matching from coarser to finer levels from top to bottom. The result of the previous level is used as the constraint for the next level, reducing the matching range and uncertainty at the same time, and finally obtaining a matching result that meets the requirements.

[0218] 1) Multi-level and multi-feature relaxation matching

[0219] Adopt a multi-level matching strategy from coarser to finer. The high-reliability image matching strategy uses a method that combines features and grayscale, a strategy of increasing the search range with a pyramid, and two-dimensional adaptive relaxation method for image matching. Multi-view images overcome the negative impacts brought by occlusion and shadow, and at the same time enhance the reliability of the matching result.

[0220] The relaxation method for image feature matching utilizes the local smoothness of the terrain. Correctly matched points have a larger neighborhood, while incorrectly matched points have a smaller neighborhood. For each extracted feature point, according to the similarity measure, alternative matching points are found. During the relaxation iteration process, the relaxation probability value of the correct alternative matching points continuously increases, and the relaxation probability value of the incorrect alternative matching points continuously decreases. When the relaxation probability value of the correct alternative matching points converges to 1 and the incorrect alternative matching points converge to 0, an accurate matching result is finally obtained. The calculation of the neighborhood adopts a combination of the eight-neighborhood of the regular grid and the adjacent nodes of the triangular mesh to adapt to the uneven distribution of feature points. The mixed neighborhood system is as Figure 8 shown.

[0221] The process of feature matching by the relaxation method is as follows:

[0222] a) Calculate the set of matching points within the search range for each feature point, and calculate the correlation coefficient of each alternative point;

[0223] b) When the relaxation probability of the alternative point with the largest correlation with the feature point is greater than that of the second-largest correlated alternative point and greater than 0.75, the iteration ends;

[0224] The matching adopts the pyramid matching strategy, generating pyramid images at all levels. The relaxation matching is carried out step by step from coarse to fine from top to bottom. The matching result of the upper level is used as the matching constraint for the lower level, narrowing the search range and reducing the uncertainty of the matching.

[0225] 2) Extract feature points

[0226] a) First, select a region of size n×n, and then calculate the first-order difference by taking the derivative of all pixel points in the region to obtain their gradients g x , g y ;

[0227] b) Take σ from 0.3 to 0.9 to perform Gaussian filtering on the obtained gradients g x , g y ;

[0228] c) Calculate the intensity value M according to Equation 9, where g x is the gradient in the x direction, g y is the gradient in the y direction, and G(s) is the Gaussian template:

[0229]

[0230] d) Sort the extreme points from largest to smallest, and select the required number of extreme points from largest to smallest.

[0231] 3) Rejection of gross errors: Automatically and reliably detect and reject the mis-matched connection points;

[0232] Step 3: Using the DOM as the reference image, calculate the object coordinates of the control points based on their image coordinates.

[0233] Step 4: Using the object coordinates obtained in Step 3, reverse calculate the image coordinates of the control points in the detection image using the RPC file of the detection image. RPC is an important file for spatial transformation mathematics, and the image coordinates of the control points are obtained by reverse calculation using RPC.

[0234] Step 5: Calculate the residuals between the calculated image coordinates and the corresponding homologous points matched on the detection image, calculate the positioning error, and calculate the average error of all control point positioning errors as the absolute geometric positioning accuracy of the image.

[0235] The residuals of the image coordinates are obtained by performing a difference operation between the calculated image coordinates and the matched coordinates using Equation 10, and then the positioning error is calculated:

[0236] dx = ΔX * GSD

[0237] dy = ΔY * GSD

[0238]

[0239] where, ΔX = X 参考 - X 图像 , ΔY = Y 参考 - Y 图像 , calculate the point position residuals of all control points in the detection image, and then use the image point residuals and the image spatial resolution to calculate the average value of the positioning error as the absolute geometric positioning accuracy of this scene image. The calculation formula is:

[0240]

[0241] The selected area size is n × n.

[0242] V. Detection of Inter - slice Mosaic Accuracy

[0243] The inter - slice mosaic accuracy characterizes the mosaic accuracy between CCD images and is an important detection index reflecting the spatio - temporal geometric quality of the images. First, perform conjugate point matching between the detection image and the reference image, then forward calculate the object coordinates corresponding to the image points from the DEM data and the RPC parameters of the reference image, then reverse calculate the image points coordinates from the RPC parameters of the detection image and the object coordinates, and finally, the mean square error of the calculated image point residuals is used as the inter - slice mosaic accuracy of the detection image.

[0244] Step 1: High - precision conjugate point matching;

[0245] Step 2: Using the image coordinates of the conjugate points on the reference image and the RPC of the reference image, forward calculate the object coordinates of the conjugate points. Based on the Taylor formula and the initial values (L 0 , B 0 , H0 ) Iteratively solve (L, B, H), where B, L, and H are the geodetic latitude, geodetic longitude, and geodetic height in the geodetic coordinate system, respectively.

[0246] Iteratively solve for (L, B), obtain H through DEM, and repeat the above steps with (L, B, H) until the difference between the previous and current (L, B) is less than the tolerance. Finally, transform the geodetic coordinates (L, B, H) obtained from the above solution to Gauss plane coordinates, and then transform them to orthoimage coordinates (x, y);

[0247] Step 3: Use the object coordinates of the connection points obtained in Step 2 and the RPC of the detection image to inversely calculate the image coordinates of the connection points on the detection image; RPC is an important file in spatial transformation mathematics, and the image coordinates of the control points can be solved by inverse calculation using RPC;

[0248] Step 4: Statistically calculate the calculated image coordinates and the residuals of the matching homologous points. Calculate the difference between the calculated image coordinates and the matching coordinates to obtain the residuals of the image coordinates. The calculation formulas are Formulas 12 and 13:

[0249] ΔX = X refer - X Formula 12

[0250] ΔY = Y refer - Y Formula 13

[0251] Calculate the mean square error of the horizontal and vertical coordinate error values of multiple homologous points of the image as the splicing accuracy between slices. The calculation formula is:

[0252]

[0253] m and n are the numbers of horizontal and vertical coordinate splices, respectively.

[0254] VI. Detection of Band Registration Accuracy

[0255] Before performing the band registration accuracy detection, first select the band used as a reference in the multispectral band, and then obtain a large number of control points through control point matching for each band and the reference band, calculate the corresponding image residuals of the control points, and calculate the mean residual through residual analysis as the evaluation index of the band registration accuracy.

[0256] The band registration accuracy characterizes the registration accuracy between each band of the multispectral image and the reference band. First, select an image of one band in the multispectral image as the reference image, and the remaining bands as the detection band images. Perform conjugate point matching on the reference band image and the detection band images respectively. Then, calculate the corresponding object coordinates of the conjugate points by forward calculation using the DEM data, the image coordinates of the conjugate points on the reference band image, and the RPC parameter file of the reference band image. Next, obtain the corresponding image coordinates of the conjugate points on the detection image by inverse calculation using the RPC parameter file of the detection image and the object coordinates. Finally, solve the residual between the inverse-calculated image points and the image coordinates of the matching homologous points, and calculate the root mean square error of the image point residuals as the band registration accuracy of the detection image.

[0257] Band registration accuracy detection algorithm process:

[0258] Step i: Select a multispectral image, perform band separation on the multispectral image, select an image of one band as the reference band image, and the remaining bands as the detection band images;

[0259] Step ii: Perform high-precision conjugate point matching on the detection band images and the reference band image;

[0260] Step iii: Calculate the object coordinates of the conjugate points by forward calculation using the image coordinates of the conjugate points on the reference band image and the RPC of the reference band image;

[0261] Based on Taylor's formula and the initial values (L 0 , B 0 , H 0 ), iteratively solve for (L, B, H). B, L, and H are the geodetic latitude, geodetic longitude, and geodetic height in the geodetic coordinate system respectively.

[0262] The iteratively solved (L, B), obtain H through the DEM, and repeat the above steps with (L, B, H) until the difference between the previous and current (L, B) is less than the tolerance. Finally, transform the calculated geodetic coordinates (L, B, H) into Gaussian plane coordinates and then into orthoimage coordinates (x, y);

[0263] Step iv: Use the object coordinates of the conjugate points obtained in Step ii and the RPC of the detection image to inversely calculate the corresponding image point coordinates of the conjugate points on the detection image. The RPC is an important file for spatial transformation mathematics, and use the RPC of the reference image to inversely calculate the corresponding image point coordinates of the control points;

[0264] Step v: Calculate the residual between the obtained image point coordinates and the matching homologous points. Calculate the difference between the calculated image coordinates and the matching coordinates to obtain the residual of the image point coordinates, calculate the image point residuals, and calculate the root mean square error of the horizontal and vertical coordinate errors of multiple homologous points in the image as the registration accuracy between bands (sensors).

[0265] VII. Relative Positioning Accuracy Detection

[0266] The relative positioning accuracy detection characterizes the relative positioning accuracy of the multi - spectral and panchromatic images of the PMS camera, ensuring the quality of subsequent panchromatic multi - spectral image fusion.

[0267] Based on the multi - spectral image and the corresponding panchromatic image, a large number of control points are first obtained through control point matching, and the residuals of the corresponding image points of the control points are calculated. Through residual analysis, the mean square error of matching is calculated as the evaluation index of relative positioning accuracy;

[0268] First, select a multi - spectral or panchromatic image as the detection image, calculate the area range of the detection image, select the corresponding panchromatic or multi - spectral image as the reference image, perform tie - point matching on the detection image and the reference image, then calculate the object - space coordinates corresponding to the tie - points by forward calculation using the DEM data and the RPC parameters of the reference image, and then calculate the image - space coordinates of the tie - points on the detection image by inverse calculation using the RPC parameters of the detection image and the object - space coordinates. Finally, calculate the residuals and evaluate the accuracy.

[0269] The algorithm flow of relative positioning accuracy detection is as follows:

[0270] Step a: Input the detection image and the corresponding panchromatic or multi - spectral image as the reference image;

[0271] Step b: Perform high - precision tie - point matching on the detection image and the reference image;

[0272] Step c: Use the image - space coordinates of the tie - points on the reference image and the RPC of the reference image to calculate the object - space coordinates of the tie - points by forward calculation; based on Taylor's formula and the initial values (L 0 , B 0 , H 0 ) to iteratively solve (L, B, H). B, L, and H are the geodetic latitude, geodetic longitude, and geodetic height in the geodetic coordinate system respectively. The iteratively solved (L, B), together with H obtained through the DEM, are used to repeat the above steps until the difference between the (L, B) values of the previous and current iterations is less than the tolerance. Finally, the geodetic coordinates (L, B, H) obtained above are transformed into Gaussian plane coordinates and then into ortho - image coordinates (x, y).

[0273] Step d: Use the RPC of the reference image to calculate the image - space coordinates by inverse calculation: RPC is an important file in spatial transformation mathematics, and using RPC inverse calculation can solve the image - space coordinates of the control points;

[0274] Step e: Calculate the residuals between the calculated image - space coordinates and the coordinates of the homologous points on the detection image. The difference between the calculated image coordinates and the matched coordinates gives the residuals of the image - space coordinates;

[0275] Step f: Calculate the mean square error of the horizontal and vertical coordinate errors of multiple homologous points in the image as the relative positioning accuracy for detecting the image.

[0276] VIII. Standardized Output of Geometric Quality Inspection Report

[0277] According to the established standard product quality evaluation criteria, evaluate the quality of the produced standard products as a standardized product quality inspection report and output it through the report output module.

[0278] IX. Analysis of Experimental Results of GF1 Cloudy Images with Poor Quality

[0279] (1) Extract reference images

[0280] According to the geographical range of the GF1 cloudy multispectral images, automatically extract and splice and crop multi-source and different-resolution reference images (DOM, DEM data) within the ground coverage range of the images from the existing reference image database based on information such as the geographical location and resolution of the detection images to obtain the corresponding reference images as Figure 9 shown.

[0281] The left figure is the extracted DOM, and the right figure is the extracted DEM. The resolution of the extracted DOM is 8 meters, the width is 9060 pixels, and the height is 7705 pixels. The resolution of the DEM is 5 meters, the width is 14555 pixels, and the height is 12374 pixels. The extracted DOM image is clear, with good quality and reliable accuracy. There is no elevation anomaly in the DEM.

[0282] (2) Detection of absolute geometric positioning accuracy

[0283] According to the matching result map of the GF1 cloudy image with poor quality, the extracted reference image DOM, and the terrain data DEM control points, a total of 13 control points are matched. Although the number of matching points is small due to the image quality, it is sufficient to detect the absolute geometric positioning accuracy of the image, and the distribution of the matching points is uniform. After calculation, the detection results of the absolute geometric positioning accuracy of the GF1 cloudy image with poor quality are as Figure 10 shown. The experimental results show that the absolute geometric positioning of the GF1 cloudy image with poor quality has an average absolute geometric positioning error of 105.6594 meters in the x direction, an average absolute geometric positioning error of 89.4631 meters in the y direction, and an average overall absolute geometric positioning error of 115.3754 meters in the x and y directions. It can be seen from the experimental results that the absolute geometric positioning accuracy of the GF1 cloudy image with poor quality is poor.

[0284] (3) Detection of relative positioning accuracy

[0285] According to the matching result map of the GF1 multi-spectral image with poor quality under cloudy conditions and the panchromatic image, a total of 18 control points were matched. The distribution of the matching points is uniform, and the number of matching points is sufficient. After calculating the registration accuracy of the GF1 multi-spectral image with poor quality under cloudy conditions, the experimental results show that the relative positioning accuracy of the GF1 multi-spectral image with poor quality under cloudy conditions is such that the mean square error of relative positioning in the x-direction is 3.346213 pixels, and the mean square error of relative positioning in the y-direction is 4.14856 pixels. The relative positioning error of the GF2 multi-spectral image with good quality relative to the panchromatic image is relatively large, exceeding 1 pixel in both the x and y directions. The relative positioning accuracy of the GF1 multi-spectral image with poor quality under cloudy conditions relative to the panchromatic image is relatively low.

[0286] (IV) Detection of Inter-chip Mosaic Accuracy

[0287] According to the inter-chip mosaic matching result map of the GF1 multi-spectral image with poor quality under cloudy conditions, the GF1 image has a total of 3 CCDs, and 8 control points are matched between every two CCDs. After calculation, the results of the inter-chip mosaic detection of the GF1 multi-spectral image with poor quality under cloudy conditions are as Figure 11 shown. The experimental results show that the inter-chip mosaic accuracy of the GF1 multi-spectral image with poor quality under cloudy conditions is such that the mean square error in the x-direction between CCD1 and CDD2 is 0.904 pixels, and the mean square error in the y-direction is 0.891 pixels. The mean square error in the x-direction between CCD2 and CDD3 is 0.873 pixels, and the mean square error in the y-direction is 0.744 pixels. The inter-chip mosaic error between each CCD is large, and the inter-chip mosaic error in both the x and y directions between each adjacent CCD is about 1 pixel, indicating that the inter-chip mosaic accuracy of the GF1 multi-spectral image with poor quality under cloudy conditions is relatively low.

[0288] (V) Band Registration Detection

[0289] After calculation, the band registration accuracy of the GF1 multi-spectral image with poor quality under cloudy conditions is as follows: The results show that the band registration deviation between Band 1, Band 3, Band 4 of the GF1 multi-spectral image with poor quality under cloudy conditions and the reference Band 2 is relatively small, all greater than about 0.1 pixel, indicating that the band registration accuracy of the GF1 multi-spectral image with poor quality under cloudy conditions is average.

[0290] It can be seen from the above experimental results that for images with good quality, images with single texture, and images with poor quality under cloudy conditions, the program can correctly detect their absolute geometric positioning accuracy, relative positioning accuracy, inter-chip mosaic accuracy, and band registration accuracy. Moreover, through the detection results, it can be intuitively seen that the spatio-temporal geometric quality of images with good quality and images with single texture is better, while the absolute geometric positioning accuracy, relative positioning accuracy, and inter-chip mosaic accuracy of images with poor quality under cloudy conditions are relatively poor.

[0291] The registration accuracy of bands can reflect the spatio-temporal geometric quality of images, but the comparison and reflection effect is not particularly obvious. For some images, some detection indicators can significantly reflect the geometric quality of remote sensing images, while some indicators do not reflect it particularly obviously. This does not mean that these indicators cannot reflect the geometric quality of the images. It is only related to the characteristics of the images. Therefore, when detecting the geometric quality of remote sensing images, it is not possible to only detect a single indicator, but multiple indicators need to be detected simultaneously. Finally, the detection results of each indicator are integrated to determine the geometric quality of the remote sensing images. Experiments show that the detection indicators of the spatio-temporal geometric quality detection method for remote sensing images proposed in this application can effectively detect the quality of the spatio-temporal geometric quality of remote sensing images and intuitively display the detection results.

Claims

1. A method for automatically detecting and evaluating the spatiotemporal geometric quality of remote sensing images, characterized in that: The absolute geometric positioning accuracy, relative positioning accuracy, band registration accuracy, and inter-slice stitching accuracy of the integrated remote sensing image are used for geometric quality detection and calculation, and the spatiotemporal geometric quality of the detected image is evaluated. According to the geographical location and resolution information of the detected image, the multi-source and different resolution reference image DOM and DEM data within the ground coverage of the evaluated image are automatically extracted, stitched, cut, and loaded to obtain the required reference image; the online detection of the spatiotemporal geometric quality of remote sensing images is realized through the automatic extraction of reference images, the automatic calculation and solution of online detection indicators of the spatiotemporal geometric quality of remote sensing images, and the automatic standardized output of geometric quality detection results; Inter-slice stitching accuracy detection: First, match the connection points of the detection image and the reference image, then calculate the object coordinates corresponding to the image points by using the DEM data and the reference image RPC parameters, and then inversely calculate the image point coordinates by using the detection image RPC parameters and the object coordinates. Finally, the mean error of the image point calculation residual is used as the inter-slice stitching accuracy of the detection image, which includes: Step 1: High-precision matching of connection points; Step 2: Use the image coordinates of the connection points on the reference image and the RPC of the reference image to calculate the object coordinates of the connection points. Based on the Taylor formula and the initial value (L0, B0, H0), iteratively solve (L, B, H). B, L, and H are the geodetic latitude, geodetic longitude, and geodetic height of the geodetic coordinate system, respectively. Iteratively solve (L, B), get H through DEM, and repeat the above steps with (L, B, H) until the two (L, B) are less than the limit difference. Finally, the geodetic coordinates (L, B, H) solved above are transformed into Gaussian plane coordinates, and then transformed into orthophoto coordinates (x, y); Step 3: Use the object coordinates of the connection points obtained in step 2 and the RPC of the detection image to inversely calculate the image point coordinates of the connection points on the detection image; RPC is an important file in space transformation mathematics, and the image point coordinates of the control points are solved by RPC inverse calculation; Step 4: Calculate the residuals of the image point coordinates and the matching points with the same name, and make a difference between the calculated image coordinates and the matching coordinates to get the residuals of the image point coordinates. The calculation formula is Equation 12 and Equation 13: ΔX = X refer -X Equation 12 ΔY = Y refer -Y Formula 13 The mean error of the horizontal and vertical coordinates of multiple points with the same name in the image is calculated as the inter-slice stitching accuracy. The calculation formula is: m and n are the number of horizontal and vertical coordinates spliced ​​respectively.

2. The method for automatic detection and evaluation of spatiotemporal geometric quality of remote sensing images according to claim 1, characterized in that: High-precision matching of spatiotemporal geometric connection points of remote sensing images: The first step is to detect constant feature points in key areas: establish a multi-scale space, define the scale space of the two-dimensional image, and use Gaussian difference kernels of different scales and image convolution to generate and detect constant key areas; The second step is spatial extreme point detection: each sampling point is compared with all its neighboring points to see whether it is larger or smaller than its neighboring points in the image domain and scale domain, and the extreme points in the scale space are found. The position and scale of the key domain are accurately determined by fitting a three-dimensional quadratic function, and the low-contrast key domain and unstable edge response points are removed at the same time. The third step is to use constant feature descriptors in key domains: 1) Key domain direction assignment: The gradient direction distribution characteristics of the key domain neighborhood pixels are used to specify the direction parameters for each key domain. The neighborhood window centered on the key domain is sampled, and the gradient direction of the neighborhood pixels is calculated using a histogram. Each key domain has three pieces of information: location, scale, and direction, thereby determining a key domain constant feature area; 2) Feature point descriptor generation: Rotate the coordinate axis to the direction of the key domain to ensure rotation invariance; The fourth step is to match the constant feature vector of the key domain: first calculate the similarity, and use the Euclidean distance as the similarity of the feature; obtain the potential match between the images through the similarity, and after obtaining the constant feature vector of the key domain, use the priority kd tree to perform a priority search to search for the 2 approximate nearest neighbor feature points of each feature point. Among these two feature points, if the closest distance divided by the second closest distance is less than a certain ratio critical value, then accept this pair of matching points; eliminate mismatches, and eliminate wrong matches according to geometric restrictions and other additional constraints; Step 5: Eliminate gross errors: Let the point set composed of all matching points be P. First, randomly select n points from P to form a subset S. Use the point set S to establish the initial model M. Let the point set excluding S from P be C. Calculate the error between the point in C and M. If the error is less than the set critical value, add the point to the point set S to form a new point set S. * , if S * If the number of midpoints is greater than N, use the point set S * Recalculate the model M based on least squares * , regenerate S and repeat the above process. If no consistent set is found, the algorithm fails. Otherwise, the maximum consistent set obtained determines the internal and external points and the algorithm ends.

3. The method for automatic detection and evaluation of spatiotemporal geometric quality of remote sensing images according to claim 1, characterized in that: Remote sensing image spatiotemporal geometric quality detection: First, the detection image and the reference image are connected to match the points. Point a and point a' are the same-name image points after matching. The detection 1A image has no projection information. First, point a' and the reference image RPC are used. RPC is the coefficient of the rational function model. RPC associates the geodetic coordinates D (Latitude, Longitude, Height) of the ground point with its corresponding image point coordinates d (line, sample) using a ratio polynomial. The object coordinates of the ground point a' corresponding to the ground point A are calculated. Then, the object coordinates of the surface point A and the RPC of the detection image are used to inversely calculate the image coordinates of point a on the detection image, that is, point a'. The image residuals △X and △Y of points a and a' are calculated, and the mean square error of the image point residuals of all connection points is calculated.

4. The method for automatic detection and evaluation of spatiotemporal geometric quality of remote sensing images according to claim 1, characterized in that: Automatic extraction and loading of reference images: According to the geographic location and resolution information of the detection image, the multi-source and different resolution reference image DOM and DEM data within the ground coverage of the image to be evaluated are automatically extracted, spliced, cropped and loaded from the existing reference image database. After the image is automatically extracted from the database, it is also necessary to mosaic, splice and crop. The specific process is as follows: Step 1: Read the RPC parameter file of the detection image and calculate the image center resolution according to the RPC model: 1) Use the image point coordinates of the image center point O and the RPC parameters to calculate the corresponding object plane coordinates (X O ,Y O ); 2) Measure the image points P and Q at N pixels in the X and Y directions (the larger the value of N, the more accurate the resolution of the calculation), and then use the RPC file to solve the plane coordinates of the corresponding object point (X P ,Y P ), (X Q ,Y Q ); 3) By GSD X =|X O -X P | / N, GSD Y =|Y O -Y Q | / N calculated resolution in X and Y directions; 4) Take the maximum value of the resolution in the X and Y directions as the resolution of the image; Step 2: Calculate the coverage of the image based on the image point coordinates and RPC parameters of the four corner points of the detected image: measure the coordinates of the four corner points of the remote sensing image, calculate the object coordinates of the four corner points based on the RPC parameters, and the quadrilateral formed by the four points is the coverage of the image on the ground; Step 3: Search the image resolution table of the reference image database according to the resolution of the image to be evaluated, obtain the slice size of the DOM and the corresponding slice image table, and the slice size of the DEM and the corresponding slice image table at the resolution, and extract the reference image from the DOM slice image table and the DEM slice image table respectively.

5. The method for automatic detection and evaluation of spatiotemporal geometric quality of remote sensing images according to claim 4, characterized in that: Step 4: Find the corresponding image in the reference image library according to the image coverage, mosaic and crop all the retrieved slice images to obtain the reference image data containing the coverage of the image to be evaluated; the image mosaic process includes color matching, image edge area matching and edge correction preprocessing process and image automatic mosaic, image mosaic line mosaic, grayscale-based mosaic and color-based mosaic; 1) Adaptive generation of seam network a) Considering the overlapping face Voronoi diagram: Let a face set A = {A1, A2, ..., A n }, the other faces are not contained by any face, that is Different faces are allowed to overlap. The distance between a point and a face in the Voronoi diagram is defined by the non-overlapping part between the two faces as the control element. Suppose face A i and A j For any two different faces in face set A, point p is on face A j Constrain down to face A i The distance d a (p,A i , A j ) is defined as the minimum distance from point p to A', calculated using formula 1: Its to face A i The distance d a (p,A i ) is defined as Formula 2: Any face A i The Voronoi polygon is expressed using Formula 3: All faces A1, A2, ..., A n The set of Voronoi polygons is the Voronoi graph of face set A: V = {V(A1), V(A2), V(A3), …, V(A n )} Equation 4 b) Get the effective range of the image Before generating the seam line network, the effective range of each image is obtained, and the invalid pixel area is excluded. First, the outer contour point set of the effective range of the image is obtained by using the boundary tracking method, and then the straight line edge of the outer contour is detected by the Hough transform method. Then, based on the straight line edge, the approximate quadrilateral of the effective range of the image is obtained. The outer contour point set of the effective range is obtained by using the boundary tracking method based on the 8-neighborhood; c) Calculate the bisector between overlapping images: Refer to the calculation method of the medial axis of a convex polygon. The specific steps are as follows: i. Number the polygons in counterclockwise direction as P O ,…,P N , draw the angle bisectors of each vertex angle, P O The intersection point of the angle bisectors of P and P1 is denoted as q1, and so on. N , P O The intersection point of the angle bisectors is denoted by q N ; ii. Calculate q in sequence i (i=0,1,2,…,N) to its opposite side (i.e. P i P i+1 ) of the distances d1, d2, …, d N ; iii. Calculate d = min(d1,…,d N ), let d = d1, reorder the vertices in a counterclockwise direction so that the distance from q1 to the opposite side is the shortest. If there are multiple vertices with the shortest distances to the opposite sides, randomly select one of them as d1; iv. From vertex P O Initially, let axis = P0; m = l, n = N; Find P respectively i The angle bisector of the vertex angle P m With P n The intersection points of the bisectors Point_m, Point_n; v. If the distance d(Point_m,axis)≤d(Point_n,axis), then axis=Point_m,m++; if the distance d(Point_m,axis)>d(Point_n,axis), then axis=Point_n,n--; Record the turning points of the central axis and number them R1, R2, ..., R N ; When recording the inflection point of the medial axis, the numbers of the vertices connected to it are also recorded, including R1 corresponding to P1, P2, R2 to R N-1 Each corresponds to a vertex, R N Corresponding to two vertices; d) Generate Voronoi polygons After obtaining the bisector between the images, the effective range is divided. All images are divided into effective ranges using the bisector and then form Voronoi polygons. 2) When the seam network is automatically optimized Assume that the Voronoi vertex is located in the n-degree overlapping region A, there are n scenes (n ≥ 3) images, they have a common overlapping area, and the pixel (x, y) is a pixel in the n-degree overlapping region. Then the difference between the n scenes at this pixel is defined as Equation 5: Where D ij (x, y) is the difference between image i and image j at pixel (x, y), which is defined as Equation 6: The optimized Voronoi vertex can be calculated using Equation 7: A single seam line is a Voronoi edge. The shortest path algorithm is used to solve this problem. Let image i be the detection image and image j be the reference image. The cost of each path is defined as formula 8: f(PS)=max D ij (x,y),(x,y)∈PS 8 Seam line optimization is to find a path that can minimize f(PS) and achieve efficient search for the minimum cost path. The binary search method is used. The upper and lower limits of the search path cost are set to g and h respectively (for 8-bit image data, the worst case values ​​are 0 and 255 respectively). The current search value z is the midpoint of the search interval, that is, z = (g + h) / 2. First, based on the starting point and the end point, find whether a path with a cost of z exists in the overlapping area of ​​the left reference image. If it exists, the upper limit of the search path cost becomes z; if it does not exist, the lower limit of the search path cost becomes z + 1, and the maximum number of searches does not exceed log2 (h-g) , for 8-bit image data, the path with the minimum cost can be found with no more than 8 searches.

6. The method for automatic detection and evaluation of spatiotemporal geometric quality of remote sensing images according to claim 5, characterized in that: 3) Image mosaic based on seamline network: After obtaining the optimized seam line network, the image mosaicking process is performed according to the seam line network to eliminate obvious seams and obtain the final seamless mosaicking image. For the reference image that has an overlapping range with the detection image, the generated effective mosaicking polygon is used to read the pixel values ​​of the image points in the multi-deformation in turn, and then the read grayscale value is assigned to the corresponding position in the mosaicking result image. The invalid pixels outside the mosaicking polygon are directly discarded, and finally the seams are feathered to obtain a seamless mosaicking image. While mosaicking and writing pixels, the uniform light processing of a large area is performed at the same time according to the parameters calculated by the uniform light processing, that is, the uniform light mosaicking integrated processing; Seam line feathering processing method: First, determine the direction of the seam line on the mosaicked image, and process it separately for different directions. The determination and processing method of the line segment direction on the seam line is as follows: if the slope of the seam line is greater than 1, the seam line is considered to be in a vertical direction. At this time, the grayscale difference between the left and right sides of the seam line is calculated. If the slope is less than or equal to 1, it is considered to be a horizontal seam. Then, the grayscale difference between the upper and lower sides of the seam line segment is calculated. Finally, the calculated grayscale difference is distributed to the left and right or upper and lower sides according to the vertical or horizontal direction of the seam line segment. Step 5, image cropping: crop the reference image data according to the coverage of the image to be evaluated, obtain and output the DOM and DEM reference images, crop the reference image that matches the detection image size, and expand it by 200 pixels * 200 pixels.

7. The method for automatic detection and evaluation of spatiotemporal geometric quality of remote sensing images according to claim 1, characterized in that: Absolute geometric positioning accuracy detection: For the detection image and the corresponding reference image, a large number of control points are first obtained through control point matching, and the image point residuals corresponding to the control points are calculated. Then, the image point residuals and image spatial resolution are used to calculate the mean positioning error as the absolute geometric positioning accuracy evaluation index; First, the detection image and the automatically extracted reference images DOM and DEM are input, and the control points of the detection image and the reference image are matched to eliminate the control points with large gross errors. Then, the image point coordinates corresponding to the control points on the detection image are inversely calculated according to the matched control point coordinates and the RPC parameter file of the detection image. The position residuals of all control points in the detection image are calculated, and then the image point residuals and the image spatial resolution are used to calculate the mean positioning error as the absolute geometric positioning accuracy of the image.

8. The method for automatic detection and evaluation of spatiotemporal geometric quality of remote sensing images according to claim 1, characterized in that: The specific steps of absolute geometric positioning accuracy detection technology are: Step 1: Select the image to be tested, calculate the geographical range of the test image, extract the reference image according to the geographical range, and perform projection conversion, mosaicking, and cropping on the reference image; Step 2: Load the detection image and the corresponding reference image, match the control points of the detection image and the corresponding reference image, and propose gross errors. The specific process of control point matching is as follows: Control point matching can support automatic matching between the detection image and the reference DOM and DEM data, realize automatic measurement of control points, and output the object coordinates and image point coordinate information of all control points in a standard format. It supports pyramid matching strategy, generates pyramid images at all levels, and performs loose matching from coarse to fine and from top to bottom. The results of the previous level are used as constraints for the next level, narrowing the matching range and reducing uncertainty, and finally obtaining matching results that meet the requirements; 1) Multi-level multi-feature relaxed matching Adopting a coarse-to-fine, multi-level matching strategy, a high-reliability image matching strategy using a method combining features with grayscale, a pyramid strategy to increase the search range, and a two-dimensional adaptive relaxation method for image matching; The relaxation method image feature matching utilizes the local smoothness of the terrain. The correct matching point has a larger neighborhood, and the wrong matching point has a smaller neighborhood. For each extracted feature point, the candidate matching point is found according to the similarity measure. In the relaxation iteration process, the relaxation probability value of the correct candidate matching point continues to increase, and the relaxation probability value of the wrong candidate matching point continues to decrease. When the relaxation probability value of the correct candidate matching point converges to 1 and the wrong candidate matching point converges to 0, an accurate matching result is finally obtained. The neighborhood calculation adopts the combination of the eight neighborhoods of the regular grid and the adjacent nodes of the triangulated network. 2) Extract feature points a) First, select an n×n area, then take the derivative of all pixels in the area and calculate the first-order difference to obtain their gradients g in the x and y directions respectively. x , g y ; b) σ is between 0.3 and 0.9, which is equivalent to the obtained gradient g x , g y Perform Gaussian filtering; c) Calculate the intensity value M according to formula 9, where g x is the gradient in the x direction, g y is the y-direction gradient, and G(s) is the Gaussian template: d) Sort the extreme value points from large to small, and select the required number of extreme value points from large to small; 3) Gross error elimination: Automatically and reliably detect and eliminate incorrectly matched connection points; Step 3: Take the reference image as DOM and calculate the object coordinates of the control points according to their image coordinates. Step 4: Using the object coordinates obtained in step 3, use the RPC file of the detection image to inversely calculate the image point coordinates of the control point in the detection image. RPC is an important file in space transformation mathematics. Use RPC to inversely calculate the image point coordinates of the control point. Step 5: Calculate the residual of the obtained image point coordinates and the corresponding matching points on the detection image, calculate the positioning error, and calculate the average error of the positioning errors of all control points as the absolute geometric positioning accuracy of the image; The image coordinates and the matching coordinates are calculated using formula 10 to perform difference calculation to obtain the residual of the image point coordinates, and then the positioning error is calculated: dx=ΔX*GSD dx=ΔY*GSD Where ΔX = X 参考 -X 图像 , ΔY=Y 参考 -Y 图像 , calculate the position residuals of all control points in the detection image, and then use the image point residuals and image spatial resolution to calculate the mean positioning error as the absolute geometric positioning accuracy of the scene image. The calculation formula is: The size of the selected area is n×n.

9. The method for automatic detection and evaluation of spatiotemporal geometric quality of remote sensing images according to claim 1, characterized in that: Band registration accuracy test: Before the band registration accuracy test, first select the band used as a reference in the multispectral band, then obtain a large number of control points by matching the control points between each band and the reference band, and calculate the image point residuals corresponding to the control points. Through residual analysis, calculate the residual mean as the band registration accuracy evaluation index; Band registration accuracy characterizes the registration accuracy of each band of the multispectral image with the reference band. First, select an image of one band in the multispectral image as the reference image, and the remaining bands as the detection band images. Match the reference band image and the detection band image respectively. Then, the object coordinates corresponding to the connection point are calculated by forward calculation based on the DEM data, the image coordinates of the connection point on the reference band image and the RPC parameter file of the reference band image. Then, the image coordinates corresponding to the connection point on the detection image are obtained by reverse calculation based on the RPC parameter file of the detection image and the object coordinates. Finally, the residual between the inversely calculated image point and the image coordinates of the matching same-name point is solved, and the mean error of the image point residual is calculated as the band registration accuracy of the detection image. Band registration accuracy detection algorithm process: Step i: Select a multispectral image, perform band separation on the multispectral image, select one of the band images as the reference band image, and the remaining bands as the detection band images; Step ii: Perform high-precision matching of connection points between the detection band image and the reference band image; Step iii: The object coordinates of the connection point are obtained by positive calculation using the image coordinates of the connection point on the reference band image and the RPC of the reference band image; Based on Taylor's formula and initial values ​​(L0, B0, H0), iterative solution (L, B, H) is performed, where B, L, and H are the geodetic latitude, geodetic longitude, and geodetic height of the geodetic coordinate system, respectively; Iteratively solve (L, B), get H through DEM, repeat the above steps with (L, B, H), until the (L, B) twice is less than the limit difference, finally, transform the above solved geodetic coordinates (L, B, H) into Gaussian plane coordinates and then transform them into orthophoto coordinates (x, y); Step iv: Use the object coordinates of the connection points obtained in step ii and the RPC of the detection image to inversely calculate the image point coordinates corresponding to the connection points on the detection image. RPC is an important file in space transformation mathematics. Use the reference image RPC to inversely calculate the image point coordinates corresponding to the control point. Step v: Calculate the residuals of the image point coordinates and the matching same-name points. Make a difference between the calculated image coordinates and the matching coordinates to get the residuals of the image point coordinates. Calculate the image point residuals. Calculate the mean error of the horizontal and vertical coordinate error values ​​of multiple same-name points in the image as the inter-band registration accuracy.

10. The method for automatic detection and evaluation of spatiotemporal geometric quality of remote sensing images according to claim 1, characterized in that: Relative positioning accuracy detection: Based on the multispectral image and the corresponding panchromatic image, a large number of control points are first obtained through control point matching, and the residuals of the control points corresponding to the image points are calculated. Through residual analysis, the matching error is calculated as a relative positioning accuracy evaluation indicator; Firstly, a multispectral or panchromatic image is selected as the detection image, the detection image area is calculated, and the panchromatic or multispectral image corresponding to the detection image is selected as the reference image. The detection image and the reference image are matched for tie points. Then, the object coordinates corresponding to the tie points are obtained by forward calculation using the DEM data and the RPC parameters of the reference image. Then, the image point coordinates of the tie points on the detection image are obtained by reverse calculation using the RPC parameters of the detection image and the object coordinates. Finally, the residual is calculated and the accuracy is evaluated. Relative positioning accuracy detection algorithm flow: Step a: Input the detection image and the corresponding panchromatic or multispectral image as the reference image; Step b: Perform high-precision matching of connection points between the detection image and the reference image; Step c: Use the image coordinates of the connection points on the reference image and the RPC of the reference image to calculate the object coordinates of the connection points; iteratively solve (L, B, H) based on the Taylor formula and the initial value (L0, B0, H0), B, L, H are the geodetic latitude, geodetic longitude, and geodetic height of the geodetic coordinate system, iteratively solve (L, B), get H through DEM, and repeat the above steps with (L, B, H) until the two (L, B) are less than the limit difference. Finally, transform the geodetic coordinates (L, B, H) solved above into Gaussian plane coordinates and then transform them into orthophoto image coordinates (x, y); Step d: Use the reference image RPC to calculate the image point coordinates: RPC is an important file in space transformation mathematics. The image point coordinates of the control point can be solved by RPC inverse calculation; Step e: The residual of the image point coordinates and the coordinates of the same-name points on the detected image are calculated, and the residual of the image point coordinates is obtained by subtracting the calculated image coordinates and the matched coordinates; Step f: Calculate the mean error of the horizontal and vertical coordinates of multiple points of the same name in the image as the relative positioning accuracy of the detected image.

Images