CN112669459A - Satellite image optimal mosaic line generation method based on feature library intelligent decision - Google Patents

Abstract

The invention discloses a satellite image optimal mosaic line generation method based on feature library intelligent decision, which specifically comprises the steps of obtaining an effective overlapping area of a satellite image; firstly, acquiring an image effective area of a satellite image, and then acquiring an overlapping area between adjacent images; generating an initial mosaic reference line; generating an initial mosaic reference line by adopting a mosaic line automatic extraction method based on pixel difference feathering; an overlapping area segmentation method based on a primary feature library; automatically generating segment embedded lines based on a secondary feature library; connecting all the sections by embedding wires; after generating a complete mosaic line, storing the mosaic line into a feature library, and taking the mosaic line as a mosaic line for next mosaic of the area under the condition that the variation difference of the ground features is small; and under the condition that the ground feature variation difference is large, the obtained complete mosaic line is used as a mosaic reference line, and the steps are repeated. The method can better meet the requirements of automatic, rapid and high-quality global DOM production in large areas and relevant research requirements.

Description

Satellite image optimal mosaic line generation method based on feature library intelligent decision

Technical Field

The invention belongs to the technical field of remote sensing data processing and image mosaic, and particularly relates to a satellite image optimal mosaic line generation method based on intelligent decision of a feature library.

Background

The current remote sensing satellite image large area mosaic field lacks the mosaic system of independent intellectual property right, and is lacked to the quick intelligent mosaic processing of multisource satellite data and data update ability, and along with means increase, resolution ratio is higher, the platform that the function is stronger constantly develops, and the remote sensing satellite constantly develops, and the data bulk of acquireing is just increasing progressively with TB magnitude. For example, in addition to the conventional photogrammetry method, LiDAR and InSAR technologies are new means for acquiring a digital elevation model, and the development of various remote sensing data acquisition technologies accelerates the speed of acquiring high-resolution earth observation data and acquiring digital ortho-image (DOM) data, but the problems of low efficiency, uneven quality, strong manual intervention, low automatic acquisition degree, low embedding speed and the like still exist in the current image embedding line selection, so that how to quickly, reasonably and automatically acquire digital geographic information products with industrial significance, especially the automatic high-quality production of digital ortho-images in a large area range, becomes a problem to be solved by urgent needs in the field of satellite remote sensing image updating embedding production and research.

Over the past decade, researchers have developed over 30 methods of damascene line automatic generation and optimization. These mosaic methods can be roughly classified into 4 types, including methods based on image difference in the overlapping region, methods based on auxiliary data, methods based on homologous points, and methods based on morphology. In the method for generating mosaic lines by image mosaic, the following limitations are mainly existed at present: firstly, the automatic mosaic line generation method based on the image difference of the overlapping area only considers the pixel difference and cannot consider the geometric difference problem of the same-name pixels of actual ground objects, so that the generation result is not real; secondly, image radiation information is not considered when the mosaic lines are generated based on the same-name point method, the mosaic lines obtained by the method easily pass through building areas, and the obtained mosaic results are not ideal; in addition, when point cloud data is adopted as auxiliary data, the method relates to a Digital Surface Model (DSM) and optical image precise matching technology, and has greater technical difficulty; and fourthly, the geometric registration precision of the remote sensing image is high by the morphology-based method, but the accurate judgment of the classification of the ground objects and the boundary of the ground objects is difficult, and the obtained mosaic line still needs to be continuously optimized by the conventional point cloud-based auxiliary data method. Therefore, the existing methods cannot completely avoid the characteristics of high-altitude objects such as buildings, cannot ensure the integrity of the target elements of the objects, and have more parts needing manual editing in the DOM generation process, so that the improvement of the generation quality of the mosaic lines and the automation degree of image updating mosaic are problems to be solved urgently at present.

In the existing mosaic line automatic generation method, for example, the mosaic line automatic extraction method based on the image difference of the overlapping region mainly comprises the following steps: 1) calculating pixel difference, local texture difference and projection difference; 2) and selecting the inlaid line according to a search strategy, wherein the search strategy comprises a line-by-line search method, a dynamic programming method, a shortest path method, an optimal spanning tree method, a dynamic contour model method and the like, and the inlaid line is obtained according to the search strategy. However, such search strategies have several major problems: 1) searching line by line to be incapable of specifying the starting point and the end point of the seam line; 2) the dynamically planned path trend is not globally optimal, and only mosaic lines can be selected according to the horizontal and vertical trends; 3) the exhaustive method of the shortest path method has extremely large calculation amount; 4) the optimal spanning tree method cannot guarantee the starting point and the end point of the initial mosaic line and can avoid the problem of falling on a building; 5) the dynamic profile selection strategy may lose the global minimum.

In summary, although the current mosaic line generation methods are various, most of the mosaic line generation methods have a slow search strategy speed, and multiple repeated feature classifications generate a great amount of computation for large-area mosaic, most of the methods have no intelligent decision effect, are weak in current situation, and currently, no mosaic line generation method based on feature library intelligent decision is studied. At present, no satellite image region mosaic method based on feature library intelligent decision is applied to large-range multi-feature region fast intelligent mosaic, and particularly, the optimization of the most important mosaic line fast generation method in the mosaic process needs to be researched urgently. The method disclosed by the invention is used for researching and generating the optimal mosaic line by adopting a search strategy based on an intelligent decision-making method of a feature library, has a quick automatic selection function, and is favorable for improving the mosaic efficiency of the next region when the generated mosaic line is put in storage, and under the condition that the variation difference of the ground features in the region is small.

Disclosure of Invention

Aiming at the problems that the search strategy speed of the mosaic line generating method of the existing remote sensing satellite image is low, repeated ground feature classification for many times can generate great calculated amount for large-area mosaic, the method has no intelligent decision effect and is poor in situational performance, the invention provides a method for generating an optimal mosaic line based on segmentation selection of a primary feature library and mosaic line selection strategy of a secondary feature library, and aims to improve the calculation efficiency of mosaic line generation, improve the search speed, enable the mosaic line to be generated reasonably to the maximum extent, enable the mosaic line to follow road vector data in the secondary feature library as far as possible and avoid crossing buildings and areas with large geometric differences to the maximum extent. Compared with the traditional mosaic line generation method, the method is a quick intelligent mosaic method which fully utilizes feature data of ground feature classification, the calculated amount in the mosaic process is small, meanwhile, the mosaic quality is guaranteed to be high, the strategy of the feature library can be suitable for various ground feature type image segments, and the method has good universality and expansibility.

The invention discloses a satellite image optimal mosaic line generating method based on feature library intelligent decision, which comprises the following steps,

s1, acquiring an effective overlapping area of the satellite images; firstly, acquiring an image effective area of a satellite image, and then acquiring an overlapping area between adjacent images;

for the image effective area acquisition of the satellite images, the method comprises the following steps:

if the obtained satellite image is a multiband image, adding multiband pixel values of each pixel point in the satellite image and storing the added value, and if the added value is 0, judging that the added value is an invalid value;

if the obtained satellite image is a single-band image, sequentially circulating from the first line to the subsequent lines according to the number of rows and columns of the image, judging each pixel point value from left to right in each line, continuously judging the next pixel point value if the pixel point value is 0, recording the corresponding row number and column number of the pixel point when the first non-0 pixel point is judged, setting the point as an upper left corner C point, acquiring the upper left corner point, sequentially circulating from the last line to the upper line, judging whether each pixel point value is 0 from right to left in each line, and setting the point as a lower right corner D point when the first non-0 pixel point value is met; and continuously circulating each row of the image from left to right, judging each pixel point value from bottom to top in each row, continuously judging the next pixel point value if the value is 0, recording the corresponding row number column number of the first non-0 pixel point, setting the point as a lower left corner B point, circularly reading the image rows from the rightmost row from right to left after obtaining the lower left corner point, judging each pixel point value from top to bottom in each row, and setting the point as an upper right corner A point when encountering the first non-0 pixel point value.

And creating a vector file according to the obtained row and column numbers of the four corner points to finish the acquisition of the effective area of the satellite image.

And acquiring the overlapping area between the adjacent images, acquiring the corresponding amplitude data vectors of the two adjacent images, respectively superposing the respective vector data of the two images, calculating the overlapping edge of the two images, and acquiring the effective overlapping area range of ACBD.

And S2, generating an initial mosaic reference line.

According to the overlapping area between the adjacent images acquired in the step S1, generating an initial mosaic reference line by adopting a mosaic line automatic extraction method based on pixel difference feathering, wherein the initial mosaic reference line is represented by adopting a coordinate sequence in the overlapping area;

the mosaic line automatic extraction method based on pixel difference eclosion is realized by adopting a Dijkstra method, the Dijkstra method is used for solving the shortest path of minimum cost (cost) between two points, wherein the cost is replaced by a pixel path with a weight, a normalized pixel value s (I, j) is distinguished by using the brightness difference delta I (I, j) and the gradient of an overlapping area as the pixel difference value of the position of the overlapping area, the gray gradient difference delta g is the absolute value of the gradient value difference, and the gradient value is the image gray level in the horizontal and vertical directionsWith a value g of magnitude of change in the direction of the gradient

The specific calculation process of the luminance difference Δ I (I, j) is:

ΔIΔI(i，j)＝|imgl(i，j)-img2(i，j)|，

wherein img1(i, j) and img2(i, j) respectively represent the brightness values of the points (i, j) of the overlapping area on two adjacent images to be mosaiced; magnitude g (i, j) of gradient value and gradient direction of point (i, j) of overlapping region

The calculation formula of (2) is as follows:

g(i，j)＝[(img(i，j)-img(i+1，j)²+(img(i，j)-img(i，j+1))²)]^1/2

img (i, j) represents the luminance value of the point (i, j) of the overlap region, and the gray gradient difference Δ g is represented as:

wherein Δ g (i, j) represents the gray gradient difference value at the pixel (i, j), g1(i, j) represents the gradient value of the pixel in one of the two adjacent images, g2(i, j) represents the gradient value of the pixel in the other of the two adjacent images,

respectively, the gradient directions corresponding to the two pixel points. The maximum value of the normalized brightness difference or gray gradient difference represents the maximum difference value of the pixels in the overlapping area, the values of the maximum difference value are max (delta I (I, j), delta g (I, j)), the corresponding pixel points are called pixel difference points, the set of the pixel difference points is CC, the upper right corner point of the boundary of the image in the overlapping area is used as the initial point A of the overlapping area, and the initial point A is used as the initial point A of the overlapping areaThe coordinate of the point A is (x0, y0), the starting point A reaches the difference point a of all pixels in the overlapping area_iBelongs to CC, and the shortest path in the overlap region is recorded as Dist [ A, a ]_i]And i is a subscript, then:

Dist[A，a_i]≥short[A，a_i]，

wherein short [ A, a_i]Representing the difference point a from the start point A of the overlap region to all pixels_iBelongs to the shortest path of CC; by continuously calculating Dist [ A, a ]_i]Up to a_iTo the lower left corner B of the overlap region in step S1, at which time Dist [ A, a ]_i]＝short[A，a_i]And obtaining an initial mosaic reference line according to the mosaic line generation method based on the pixel feather difference.

S3, an overlapping area segmentation method based on a primary feature library, comprising the following steps:

s31, selecting and grading a feature library;

and (3) carrying out regional object feature classification on the overlapping region in the acquired satellite image, and establishing a primary feature library and a secondary feature library, wherein the primary feature library is used for determining the segmentation limit, and the secondary feature library is used for determining the selection of the mosaic line. Selecting a primary characteristic library and a secondary characteristic library from the characteristic library data, wherein:

the primary feature library comprises dense boundary vector data and sparse boundary vector data;

the secondary feature library comprises road vector data, water system vector data, scarp data and the like.

According to the ground feature type and feature library data established by various remote sensing geographic information resources, dividing an overlap region into a dense section, a sparse section and a mixed section according to dense boundary vector data and sparse boundary vector data in a primary feature library, obtaining the dense section in the overlap region according to step S32 according to the dense boundary vector data, obtaining the sparse section in the overlap region according to step S33 according to the sparse boundary vector data, wherein the regions except the dense section and the sparse section are collectively called as the mixed section, and the specific segmentation method is shown in steps S32 and S33; the dense boundary vector data establishes a feature library according to data (first level and second level) of buildings, high-level roads and the like, noctilucent remote sensing data, water system vectors and the like, and the sparse boundary vector data establishes a feature library according to desert vectors, DEM data, expressways, low-level road data (third level and fourth level), water system vectors and the like.

S32, dense segment segmentation, which specifically includes:

s321, firstly, judging whether dense boundary vector data exists in an overlapping area, if so, segmenting according to the step S322, if not, segmenting, and directly selecting an initial mosaic reference line in the step S2;

s322, constructing a minimum circumscribed rectangle according to the dense boundary vector data in the primary feature library, wherein the specific construction process is as follows:

drawing the minimum circumscribed rectangle by utilizing a minAreaRect () function meter in opencv, wherein the minimum circumscribed rectangle comprises four end points and four edges of the minimum circumscribed rectangle: firstly, defining a circumscribed rectangle set to be a circumscribed rectangle set formed by the outlines of the dense boundary vectors and the sparse boundary vectors in the overlapping area ABCD in the step S1, and then defining a minimum circumscribed rectangle set to be a set formed by the minimum circumscribed rectangle of the outlines of each dense boundary vector and each sparse boundary vector; the outline here refers to a set of points within the overlap area ABCD obtained in step S1 corresponding to the boundary vector; then drawing the central point of the minimum circumscribed rectangle to obtain four end points of the minimum circumscribed rectangle; finally, drawing four sides of the minimum circumscribed rectangle to finish drawing the minimum circumscribed rectangle;

rotating a polygon formed by original dense boundary vector data end points, wherein the rotation degree is-90 degrees, obtaining a simple external rectangle of the polygon after rotating each degree, recording the area, vertex coordinates and rotation degree of each simple external rectangle, obtaining the external rectangle with the minimum area, and finally rotating the external rectangle to the rotation angle with the same direction and the same size as the recorded rotation degree, thus obtaining the minimum external rectangle.

S323, in the case where the minimum bounding rectangle obtained in step S322 intersects with or includes the bounding rectangles of the remaining boundary vectors, the overlap region in step S1 is segmented according to the dense boundary vectors.

S324, according to the upper and lower boundaries of the minimum bounding rectangle of the dense boundary vector in step S322, the minimum bounding tangent equation is expressed as: and y is kx + b, and the equations of the upper minimum external tangent line and the lower minimum external tangent line are solved according to the coordinates of the four end points of the minimum external rectangle:

y₁=k₁x+b₁

y₂=k₂x+b₂，

wherein the slope k₁=k₂，b₁And b₂Is a constant term of the two smallest circumscribed tangents.

S33, sparse segment segmentation, which specifically includes:

firstly, judging whether the overlap region has sparse boundary vector data or not, if so, constructing a minimum circumscribed rectangle according to the step S322, obtaining a lower boundary intersection point S by intersecting the overlap region boundary and the known sparse segment boundary, and regarding the sparse segment segmentation boundary line y30, the expression is y₃=k₃x+b₃Slope k of₃=k₁=k₂Substituting the coordinate of the lower boundary intersection point S into the expression of y30 to calculate b₃Values, resulting in an expression for sparse segment segmentation line y 30; if no sparse boundary vector data exists, the initial mosaic reference line in the step S2 is directly selected without segmentation, and the other segments except the dense segment and the sparse segment are collectively called as a mixed segment.

S4, matching the road vector, the water-based vector data, and the image pixel data in the secondary feature library, that is, the pixel matrix of the overlap area image, specifically includes:

the three-parameter method is that translation of X, Y and Z axes is carried out between two reference surfaces, and according to three parameters of translation and reference points, the road vector and water system vector coordinates of the input overlapping area image are converted by the three-parameter method to obtain output vector geographic coordinate data, and the calculation is as follows:

wherein [ X, Y, Z]_newTo obtain vector data, [ X, Y, Z]_orginIs originalVector coordinates, [ dX, dY, dZ]Three parameters for translation.

And matching road vector data and water system vector data with geographic entity coordinate points to the overlapping area image according to the prior characteristic data of the secondary characteristic library, and matching entities with the same name. The position of the pixel point in the image is represented by the pixel coordinate, when the road is matched with the water system vector data and the image pixel of the overlapping area, the corresponding geographic coordinate needs to be converted into the pixel coordinate, and the conversion method comprises the following steps:

first, the geographic coordinate [ GEO ] of the upper left corner of the overlap region image_x，GEO_y]Known, and the resolution R of the image_x，R_yIt is known to compute the geographic coordinates of an arbitrary pixel location from the geographic coordinates of the upper left corner of the image and the resolution of the image. The formula for converting the geographic coordinates to pixel coordinates is:

the formula for converting the pixel coordinates into the geographic coordinates is as follows:

X＝GEO_x+col*R_x

Y＝GEO_y+row*R_y，

wherein X and Y in the above two coordinate transformation formula represent geographical coordinates, PIX_x，PIX_yRepresenting pixel coordinates, GEO_x、GEO_yRepresenting the geographic coordinates of the upper left corner of the image, col, row representing the number of columns and rows of pixels of the image, R_xAnd R_yRepresenting the row and column resolution of the image. And matching the road vector with the A point coordinate as a starting point to the image pixel to prepare as data generated by a subsequent mosaic line.

S5, automatically generating the segment embedded line based on the secondary feature library, which specifically comprises:

s51, aiming at the dense segment, adopting an automatic mosaic line generating method based on a water system vector and a road vector feature vector, which comprises the following specific steps:

when there is no secondary feature library data, the mosaic line is the initial mosaic reference line in step S2;

when secondary feature library data exist, an automatic mosaic line selection method based on a water system vector and a road vector is adopted to obtain mosaic lines:

the starting point of the optimized mosaic line in the dense section is A₂And fitting the water system or road vector coordinates of the secondary feature library by selecting N lower intersection points of the water system or road vector curve in the overlapping area and the lower boundary straight line y20 of the dense section, namely the dense section road end point to obtain the water system or road vector curve, wherein the water system vector curve is as follows:

where ρ is₀，ρ₁，…，

Polynomial coefficient of water-based vector curve, N₀The order of the water system or the road vector curve is as follows:

wherein, ω is₀，ω₁，…，

Polynomial coefficients of the road vector curve;

s511, if the secondary feature library has water system vector data, collecting a water system vector curve f₁(x)，f₂(x)，f₃(x)，...，f_N(x) Forming an equation set in a simultaneous manner with the straight line y20, and solving the intersection point set to obtain a water system vector terminal set in the dense section;

s512, if the secondary feature library has no water system vector data, the road vector curve h is processed₁(x)，h₂(x)，...，h_N(x) And (4) combining the curve set and the straight line y20 to form an equation set, and solving to obtain a road vector terminal set in the dense section.

S513, taking the intersection point of the reference mosaic line of the overlapping area obtained in the step S2 and the upper border line y10 of the dense section as a starting point, taking the water system or road vector point which is continuous on the lower border line y20 of the dense section and penetrates through the upper and lower borders of the dense section as an end point, searching for M water systems or roads according to constraint conditions, wherein N water systems or roads (N is less than or equal to M) are water systems or road vectors containing end points until the water system or road vector line which is continuous on the lower border line y20 and penetrates through the upper and lower borders of the dense section is searched, selecting the searched continuous water systems or road vector lines which penetrate through the upper and lower borders of the dense section as optimal mosaic lines, if the continuous water systems or roads are not searched, discarding the road or water system vector data, and selecting the initial mosaic line in the step S2 as the optimal mosaic line.

S514, selecting the optimal mosaic line under the condition of running through of the water system or road vector characteristic data;

the process of judging the penetration of the water system or road vector characteristic data comprises the steps of fitting the coordinates of discrete points of each water system or road into a water system or road vector curve equation, substituting the coordinates of the starting point and the end point of the boundary point of each water system or road into y10 and y20 equations of a boundary straight line, judging that the water system or road vector penetrates through the upper boundary and the lower boundary of a dense section or a sparse section when the coordinates of the end point of the water system or road vector are on the boundary straight line, selecting the penetrated road or water system vector data, discarding the road or water system vector data if the coordinates of the end point of the water system or road vector are not on the boundary straight line, and selecting the initial inlaid line in the step S2 as the inlaid line in the section.

The selection process of the inlaid wire comprises the following steps: selecting a vector of a shortest penetrating path according to the sequence of a water system, a high-level road and a low-level road; then according to the distance from the lower boundary end point of the water system or road vector line in the dense section to the intersection point of the reference mosaic line in the overlapping area in the lower boundary of the dense section, according to the principle of the shortest distance priority, selecting the water system or road vector corresponding to the shortest distance end point as the mosaic line after optimization, and for the intersection of the lower boundary end point of the water system or road vector line in the dense section to the lower boundary of the reference mosaic line in the dense sectionDistance between points, selecting

Constraint of proximity of water system or road vector to reference line, d_iRepresents the distance between the terminal point of the ith water system or road vector line at the lower boundary of the dense section and the intersection point of the reference mosaic line at the lower boundary of the dense section, N₁The number of the searched water penetrating systems or road vectors is shown. The method for calculating the length of the water system or the road path comprises the following steps: equally dividing the water system or the road curve penetrated in the step S51 into n equal parts according to the actual road vector length, and solving the total curve length of the equally divided curve by adopting a short straight line method, wherein the method specifically comprises the following steps:

wherein L represents the total curve length, Δ x_iThe arc length in the ith subinterval corresponds to the transverse axis length, Δ y_iThe arc length in the ith sub-interval corresponds to the longitudinal axis length, x_iIs the abscissa of any point in the ith subinterval, i is the curve dividing interval number, wherein i is 1,2, …, n, f'²(x_i) Is the equation f (x) at x for the water system or road curve_iThe derivative of the point is squared, and finally n equal curve lengths are accumulated to obtain the total curve length L, i.e. the length of the corresponding water system or road path penetrating the upper and lower boundaries of the dense section, and the curve lengths of all water systems or road vectors in the dense section are calculated to make L_min＝mm{L₁，L₂，...，L_NSelecting L_minAs a water system or road vector shortest path constraint. Obtaining the optimal mosaic line segment A according to the method in step S514₂A₅。

S515, the coordinate of the end point of the shortest water system or road vector line is coincident with the coordinate of the end point of the lower boundary of the reference mosaic line in the dense section or d_minUnder the condition that the distance is less than or equal to 10m, directly connecting the end point of the shortest water system or road vector line with the end point of the embedded reference line in the step S2 in the lower boundary y20 of the dense section, and connecting the end point of the shortest water system or road vector line with the end point of the embedded reference line according to the shortest water system or road vector line and the embedded reference lineThe distance between two intersection points of the line on the lower boundary of the dense section is used as the diameter of a circle, the middle point of the two points is used as the center of the circle to establish a buffer circle, and the optimal mosaic line section A in the step S514 is used in the circle₂A₅Intersection A with lower boundary y20 of dense segment₅(xx5, yy5) as the starting point, searching for a water system or road vector, and selecting a mosaic line A in the section according to the steps from S511 to S514₅A₃。

And S516, obtaining the dense section inlaid wire as an optimal inlaid wire according to the steps from S511 to S515.

S52, aiming at the sparse segment, when there is no secondary feature library vector, keeping the original mosaic reference line in the step S2 as the mosaic line segment of the region;

when a secondary feature library vector exists, an automatic mosaic line generating method based on a road and a water system vector is adopted, and the method specifically comprises the following steps:

s521, obtaining the optimal mosaic line AA with priority on the road by taking the point A as the starting point according to the selection mode of the through road of the first road and the second water system in the sparse terrain section and according to the steps S511 to S515₄If there is no road or water system that satisfies the condition in step S514, the initial mosaic reference line in step S2 is selected.

S522, the shortest road or water system vector line terminal point and the reference mosaic line are coincided in the boundary coordinates of the sparse section or d_minUnder the condition that the distance is less than or equal to 10m, directly connecting the end point of the shortest road or the water system vector line in the lower boundary y30 of the sparse section with the end point of the embedded reference line in the step S2; distance d between terminal and reference line_minWhen the distance between two intersection points of the shortest water system or the road vector line and the reference mosaic line on the lower boundary of the sparse section is larger than 10m, the diameter of the circle is taken as the distance between the two intersection points, the midpoint of the two intersection points is taken as the center of the circle to establish a buffer circle, and the optimal mosaic line segment AA in the sparse section is taken as the center of the circle₄Intersection A with sparse segment boundary y30₄Searching water system or road vector as starting point, and selecting mosaic line A according to steps S511-S514₄A₁。

S523, obtaining AA by taking the mosaic line penetrating through the sparse segment obtained in the steps from S521 to S522 as an optimal mosaic line₄、A₄A₁As the optimal mosaic line segment.

S53, aiming at the generation of mosaic lines in the mixed sections except the dense sections and the sparse sections, respectively selecting the intersection points of the initial mosaic reference lines and the sparse sections and the intersection points of the initial mosaic reference lines or the boundary lines of the dense sections in the step S2 as the starting points and the ending points of the mosaic line sections, and adopting the initial mosaic reference lines in the step S2 as the mosaic lines.

And S6, connecting the embedded wires among the sections.

The hybrid section uses the initial mosaic reference line; the dense segment uses the feature vector optimal mosaic line a acquired in step S51₂A₅、A₅A₃(ii) a The sparse segment uses the optimal mosaic line AA obtained in step S52₄、A₄A₁. And finally, combining the endpoints of the segment-to-segment mosaic lines on the segment boundary lines of the mixed segment, the sparse segment and the dense segment into a complete optimal mosaic line vector in an overlapping region.

S7, after generating a complete mosaic line, storing the mosaic line into a feature library, and taking the mosaic line as the mosaic line of the next mosaic area under the condition that the variation difference of the ground features is small; and under the condition of large terrain variation difference, taking the obtained complete mosaic line as a mosaic reference line, and continuing to perform the mosaic line optimization algorithm of the steps from S1 to S6, thereby improving the mosaic line generation efficiency in the region mosaic.

The invention has the beneficial effects that:

1. compared with other mosaic line generation methods, the method disclosed by the invention has the advantages that the calculation amount of mosaic line generation is reduced, the mosaic quality is improved, the method has strong current performance, the integrated geographic information data resources and the remote sensing data resources are fully and reasonably utilized, the advantage of higher information utilization degree is presented, and the method can better meet the requirements of large-area DOM production and related research in an automatic, rapid and high-quality global scope.

2. Compared with the traditional method, the method has higher accuracy and stronger universality when the mosaic line generated by the method is used for carrying out regional mosaic on the satellite image, for example, the method can provide an intelligent strategy for the existing satellite image to be mosaic, and carries out intelligent decision selection on the mosaic line according to the characteristic library data of the object region, thereby being beneficial to solving the problem that the mosaic line passes through high-rise buildings and other high-rise buildings, ensuring the quality precision of the selected mosaic line and further having certain practical significance for generating a full sphere, particularly a high-quality and automatic digital orthographic product (DOM) in a large regional range.

3. The method of the invention aims at the secondary region mosaic of the same region of different source data, and can directly utilize the existing mosaic line in the library to carry out region mosaic or be used as the mosaic line reference line of the region mosaic. After the mosaic line is searched, the mosaic line is put into a warehouse, the region mosaic is carried out on the same region for different source data, the mosaic line generated by the method can be used as the reference mosaic line of the next region, or the existing mosaic line in the warehouse can be directly utilized to carry out rapid intelligent region mosaic under the conditions that the change of the terrain in the region is not large and the image overlapping regions are basically consistent.

Drawings

FIG. 1 is a flow chart of an optimal mosaic line selection method based on a feature library;

FIG. 2 is a schematic diagram of an image to be embedded and an overlapping area;

FIG. 3 is a schematic diagram of an overlap reference damascene line;

FIG. 4 is a schematic view of a minimum circumscribed rectangle;

FIG. 5 is a schematic view of an overlapping section;

FIG. 6 is a schematic diagram of the relationship between two reference planes between three parameters;

FIG. 7 is a schematic diagram of a geographic coordinate to pixel coordinate transformation;

FIG. 8 is a schematic diagram of a dense segment and sparse segment damascene line generation process;

FIG. 9 is a diagram of the effect of optimal damascene line merging.

Detailed Description

For a better understanding of the present disclosure, an example is given here.

The embodiment discloses a method for generating an optimal mosaic line of a satellite image based on intelligent decision of a feature library, the flow chart of which is shown in fig. 1, the method comprises the following specific steps,

s1, acquiring an effective overlapping area of the satellite images; firstly, an image effective area of a satellite image is obtained, and then an overlapping area between adjacent images is obtained, wherein a schematic diagram of an image to be embedded and the overlapping area is shown in fig. 2.

For the image effective area acquisition of the satellite images, the method comprises the following steps:

if the obtained satellite image is a multiband image, adding multiband pixel values of each pixel point in the satellite image and storing the added value, and if the added value is 0, judging that the added value is an invalid value;

if the obtained satellite image is a single-band image, sequentially circulating from the first line to the subsequent lines according to the number of rows and columns of the image, judging each pixel point value from left to right in each line, continuously judging the next pixel point value if the pixel point value is 0, recording the corresponding row number and column number of the pixel point when the first non-0 pixel point is judged, setting the point as an upper left corner C point, acquiring the upper left corner point, sequentially circulating from the last line to the upper line, judging whether each pixel point value is 0 from right to left in each line, and setting the point as a lower right corner D point when the first non-0 pixel point value is met; and continuously circulating each row of the image from left to right, judging each pixel point value from bottom to top in each row, continuously judging the next pixel point value if the value is 0, recording the corresponding row number column number of the first non-0 pixel point, setting the point as a lower left corner B point, circularly reading the image rows from the rightmost row from right to left after obtaining the lower left corner point, judging each pixel point value from top to bottom in each row, and setting the point as an upper right corner A point when encountering the first non-0 pixel point value.

And creating a vector file according to the obtained row and column numbers of the four corner points to finish the acquisition of the effective area of the satellite image.

For the acquisition of the overlapping area between the adjacent images, the corresponding amplitude data vectors of the two adjacent images are acquired, the respective vector data of the two images are respectively superposed, the overlapping edge of the two images is calculated, and the range of the effective overlapping area is acquired as ACBD, as shown in fig. 2.

And S2, generating an initial mosaic reference line.

According to the overlapping area between the adjacent images acquired in the step S1, an initial mosaic reference line is generated by using an automatic mosaic line extraction method based on pixel difference feathering, and the initial mosaic reference line is represented in the overlapping area by using a coordinate sequence, wherein the coordinate sequence is as follows:

(x1,y1),(x2,y2),(x3,y3),...,(xi,yi)，

wherein, (x1, y1) is the value of the first coordinate point in the coordinate sequence, the values of the subsequent coordinate points are analogized in turn, and (xi, yi) is the value of the ith coordinate point in the coordinate sequence;

the mosaic line automatic extraction method based on the pixel difference eclosion is realized by adopting a Dijkstra method, the Dijkstra method is used for solving the shortest path of the minimum cost (cost) between two points, wherein the cost is replaced by a pixel path with a weight, a normalized pixel value s (I, j) is distinguished by using the brightness difference delta I (I, j) and the gradient of an overlapping area as the pixel difference value of the position of the overlapping area, the gray gradient difference delta g is the absolute value of the gradient value difference, the gradient value is the change value g of the image gray in the horizontal and vertical directions, and the gradient direction is the gradient direction

The specific calculation process of the luminance difference Δ I (I, j) is:

ΔI(i，j)＝|img1(i，j)-img2(i，j)|，

wherein img1(i, j) and img2(i, j) respectively represent the brightness values of the points (i, j) of the overlapping area on two adjacent images to be mosaiced; magnitude g (i, j) of gradient value and gradient direction of point (i, j) of overlapping region

The calculation formula of (2) is as follows:

g(i，j)＝[(img(i，j)-img(i+1，j)²+(img(i，j)-img(i，j+1))²)]^1/2

img (i, j) represents the luminance value of the point (i, j) of the overlap region, and the gray gradient difference Δ g is represented as:

wherein Δ g (i, j) represents the gray gradient difference value at the pixel (i, j), g1(i, j) represents the gradient value of the pixel in one of the two adjacent images, g2(i, j) represents the gradient value of the pixel in the other of the two adjacent images,

respectively, the gradient directions corresponding to the two pixel points. The maximum value of the normalized brightness difference or gray gradient difference represents the maximum difference value of the pixels in the overlapping area, the values of the maximum difference value are max (delta I (I, j), delta g (I, j)), the corresponding pixel points are called pixel difference points, the set of the pixel difference points is CC, the upper right corner point of the boundary of the image in the overlapping area is used as the initial point A of the overlapping area, the coordinates of the initial point A are (x0, y0), and the initial point A reaches all the pixel difference points a in the overlapping area_iBelongs to CC, and the shortest path in the overlap region is recorded as Dist [ A, a ]_iI is a subscript, then:

Dist[A，a_i]≥short[A，a_i]，

wherein short [ A, a_i]Representing the difference point a from the start point A of the overlap region to all pixels_iBelongs to the shortest path of CC; by continuously calculating Dist [ A, a ]_i]Up to a_iTo the lower left corner B of the overlap region in step S1, at which time Dist [ A, a ]_i]＝short[A，a_i]An initial mosaic reference line is obtained according to the above mosaic line generation method based on the pixel feathering difference, and a schematic diagram of the initial mosaic reference line obtained by using the overlapping region obtained in step S1 as an input is shown in fig. 3.

S3, an overlapping area segmentation method based on a primary feature library, comprising the following steps:

s31, selecting and grading a feature library;

and (3) carrying out regional object feature classification on the overlapping region in the acquired satellite image, and establishing a primary feature library and a secondary feature library, wherein the primary feature library is used for determining the segmentation limit, and the secondary feature library is used for determining the selection of the mosaic line. Selecting a primary characteristic library and a secondary characteristic library from the characteristic library data, wherein:

the primary feature library comprises dense boundary vector data and sparse boundary vector data;

the secondary feature library comprises road vector data, water system vector data, scarp data and the like.

According to the ground feature type and feature library data established by various remote sensing geographic information resources, dividing an overlap region into a dense section, a sparse section and a mixed section according to dense boundary vector data and sparse boundary vector data in a primary feature library, obtaining the dense section in the overlap region according to step S32 according to the dense boundary vector data, obtaining the sparse section in the overlap region according to step S33 according to the sparse boundary vector data, wherein the regions except the dense section and the sparse section are collectively called as the mixed section, and the specific segmentation method is shown in steps S32 and S33; the dense boundary vector data establishes a feature library according to data (first level and second level) of buildings, high-level roads and the like, noctilucent remote sensing data, water system vectors and the like, and the sparse boundary vector data establishes a feature library according to desert vectors, DEM data, expressways, low-level road data (third level and fourth level), water system vectors and the like.

S32, dense segment segmentation, which specifically includes:

s321, firstly, judging whether dense boundary vector data exists in an overlapping area, if so, segmenting according to the step S322, if not, segmenting, and directly selecting an initial mosaic reference line in the step S2;

s322, constructing a minimum circumscribed rectangle according to the dense boundary vector data in the primary feature library, wherein the specific construction process is as follows:

the minimum bounding rectangle, as shown in fig. 4, is drawn by using the minAreaRect () function in opencv, and includes four end points and four sides of the minimum bounding rectangle: firstly, defining a circumscribed rectangle set to be a circumscribed rectangle set formed by the outlines of the dense boundary vectors and the sparse boundary vectors in the overlapping area ABCD in the step S1, and then defining a minimum circumscribed rectangle set to be a set formed by the minimum circumscribed rectangle of the outlines of each dense boundary vector and each sparse boundary vector; the outline here refers to a set of points within the overlap area ABCD obtained in step S1 corresponding to the boundary vector; then drawing the central point of the minimum circumscribed rectangle to obtain four end points of the minimum circumscribed rectangle; finally, drawing four sides of the minimum circumscribed rectangle to finish drawing the minimum circumscribed rectangle; the concrete formula is as follows: when a point (x1, y1) on the plane of the overlap region is rotated counterclockwise by an angle θ around another point (x0, y0) is (x2, y2), the following formula is given:

x2＝(x1-x0)cosθ-(y1-y0)sinθ+x0

y2＝(x1-x0)sinθ+(yl-y0)cosθ+y0，

rotating a polygon formed by original dense boundary vector data end points, wherein the rotation degree is-90 degrees, obtaining a simple external rectangle of the polygon after rotating each degree, recording the area, vertex coordinates and rotation degree of each simple external rectangle, obtaining the external rectangle with the minimum area, and finally rotating the external rectangle to the rotation angle with the same direction and the same size as the recorded rotation degree, thus obtaining the minimum external rectangle.

S323, in the case where the minimum bounding rectangle obtained in step S322 intersects with or includes the bounding rectangles of the remaining boundary vectors, the overlap region in step S1 is segmented according to the dense boundary vectors. FIG. 5 is a schematic view of an overlapping section.

S324, according to the upper and lower boundaries of the minimum bounding rectangle of the dense boundary vector in step S322, the minimum bounding tangent equation is expressed as: and y is kx + b, and the equations of the upper minimum external tangent line and the lower minimum external tangent line are solved according to the coordinates of the four end points of the minimum external rectangle:

y₁＝k₁x+b₁

y₂＝k₂x+b₂，

wherein the slope k₁＝k₂，b₁And b₂Is a constant term of the two smallest circumscribed tangents.

S33, sparse segment segmentation, which specifically includes:

firstly, judging whether the overlap region has sparse boundary vector data or not, if so, constructing a minimum circumscribed rectangle according to the step S322, obtaining a lower boundary intersection point S by intersecting the overlap region boundary and the known sparse segment boundary, and regarding the sparse segment segmentation boundary line y30, the expression is y₃＝k₃x+b₃Slope k of₃＝k₁＝k₂Substituting the coordinate of the lower boundary intersection point S into the expression of y30 to calculate b₃Values, resulting in an expression for sparse segment segmentation line y 30; if no sparse boundary vector data exists, the initial mosaic reference line in the step S2 is directly selected without segmentation, and the other segments except the dense segment and the sparse segment are collectively called as a mixed segment.

S4, matching the road vector, the water-based vector data, and the image pixel data in the secondary feature library, that is, the pixel matrix of the overlap area image, specifically includes:

the coordinates in different spatial reference coordinate systems are different on different reference planes, the reference plane is a part of the coordinate system, because the reference plane involves translation or rotation relative to the geocentric during positioning, so we cannot directly perform such conversion, and a conversion parameter is needed, and these parameters are based on different models, the three-parameter method in this embodiment is to perform translation of X, Y, Z axes between two reference planes, fig. 6 is a relation schematic diagram of two reference planes between three parameters, and according to the three parameters of translation and the reference point, the road vector and water system vector coordinates of the input overlapping area image are converted by using the three-parameter method to obtain the output vector geographic coordinate data, which is calculated as:

wherein [ X, Y, Z]_newTo obtain vector data, [ X, Y, Z]_orginAs original vector coordinates, [ dX, dY, dZ]Three parameters for translation.

And matching road vector data and water system vector data with geographic entity coordinate points to the overlapping area image according to the prior characteristic data of the secondary characteristic library, and matching entities with the same name. The position of the pixel point in the image is represented by the pixel coordinate, when the road is matched with the water system vector data and the image pixel of the overlapping area, the corresponding geographic coordinate needs to be converted into the pixel coordinate, and the conversion method comprises the following steps:

first, the geographic coordinate [ GEO ] of the upper left corner of the overlap region image_x，GEO_y]Known, and the resolution R of the image_x，R_yIt is known to compute the geographic coordinates of an arbitrary pixel location from the geographic coordinates of the upper left corner of the image and the resolution of the image. If a geographical coordinate is given, the position of the geographical coordinate on the image can also be calculated, and fig. 7 is a schematic diagram of the conversion of the geographical coordinate and the pixel coordinate.

The formula for converting the geographic coordinates to pixel coordinates is:

the formula for converting the pixel coordinates into the geographic coordinates is as follows:

X＝GEO_x+col*R_x

Y＝GEO_y+row*R_y，

wherein X and Y in the above two coordinate transformation formula represent geographic coordinates, PIPIX_x，PIX_yRepresenting pixel coordinates, GEO_x、GEO_yRepresenting the geographic coordinates of the upper left corner of the image, col, row representing the number of columns and rows of pixels of the image, R_xAnd R_yRepresenting the row and column resolution of the image. And matching the road vector with the A point coordinate as a starting point to the image pixel to prepare as data generated by a subsequent mosaic line.

S5, automatically generating the segment embedded line based on the secondary feature library, which specifically comprises:

s51, aiming at the dense segment, adopting an automatic mosaic line generating method based on a water system vector and a road vector feature vector, which comprises the following specific steps:

when there is no secondary feature library data, the mosaic line is the initial mosaic reference line in step S2;

when secondary feature library data exist, an automatic mosaic line selection method based on a water system vector and a road vector is adopted to obtain mosaic lines:

the starting point of the optimized mosaic line in the dense section is A₂And fitting the water system or road vector coordinates of the secondary feature library by selecting N lower intersection points of the water system or road vector curve in the overlapping area and the lower boundary straight line y20 of the dense section, namely the dense section road end point to obtain the water system or road vector curve, wherein the water system vector curve is as follows:

where ρ is₀，ρ₁，…，

Polynomial coefficient of water-based vector curve, N₀The order of the water system or the road vector curve is as follows:

wherein, ω is₀，ω₁，…，

Polynomial coefficients of the road vector curve;

s511, if the secondary feature library has water system vector data, collecting a water system vector curve f₁(x)，f₂(x)，f₃(x)，...，f_N(x) Forming an equation set in a simultaneous manner with the straight line y20, and solving the intersection point set to obtain a water system vector terminal set in the dense section;

s512, if the secondary feature library has no water system vector dataThe road vector curve h₁(x)，h₂(x)，...，h_N(x) And (4) combining the curve set and the straight line y20 to form an equation set, and solving to obtain a road vector terminal set in the dense section.

S513, taking the intersection point of the reference mosaic line of the overlapping area obtained in the step S2 and the upper border line y10 of the dense section as a starting point, taking the water system or road vector point which is continuous on the lower border line y20 of the dense section and penetrates through the upper and lower borders of the dense section as an end point, searching for M water systems or roads according to constraint conditions, wherein N water systems or roads (N is less than or equal to M) are water systems or road vectors containing end points until the water system or road vector line which is continuous on the lower border line y20 and penetrates through the upper and lower borders of the dense section is searched, selecting the searched continuous water systems or road vector lines which penetrate through the upper and lower borders of the dense section as optimal mosaic lines, if the continuous water systems or roads are not searched, discarding the road or water system vector data, and selecting the initial mosaic line in the step S2 as the optimal mosaic line.

S514, selecting the optimal mosaic line under the condition of running through of the water system or road vector characteristic data;

the process of judging the penetration of the water system or road vector characteristic data comprises the steps of fitting the coordinates of discrete points of each water system or road into a water system or road vector curve equation, substituting the coordinates of the starting point and the end point of the boundary point of each water system or road into y10 and y20 equations of a boundary straight line, judging that the water system or road vector penetrates through the upper boundary and the lower boundary of a dense section or a sparse section when the coordinates of the end point of the water system or road vector are on the boundary straight line, selecting the penetrated road or water system vector data, discarding the road or water system vector data if the coordinates of the end point of the water system or road vector are not on the boundary straight line, and selecting the initial inlaid line in the step S2 as the inlaid line in the section.

The selection process of the inlaid wire comprises the following steps: selecting a vector of a shortest penetrating path according to the sequence of a water system, a high-level road and a low-level road; then according to the distance from the terminal point of the lower boundary of the dense section of the water system or the road vector line to the intersection point of the lower boundary of the reference mosaic line in the overlapping area and according to the principle of priority of the shortest distance,selecting a water system or road vector corresponding to the shortest distance end point as an optimized mosaic line, and selecting the distance from the end point of the lower boundary of the water system or road vector line in the dense section to the intersection point of the reference mosaic line in the lower boundary of the dense section

Constraint of proximity of water system or road vector to reference line, d_iRepresents the distance between the terminal point of the ith water system or road vector line at the lower boundary of the dense section and the intersection point of the reference mosaic line at the lower boundary of the dense section, N₁The number of the searched water penetrating systems or road vectors is shown. The method for calculating the length of the water system or the road path comprises the following steps: equally dividing the water system or the road curve penetrated in the step S51 into n equal parts according to the actual road vector length, and solving the total curve length of the equally divided curve by adopting a short straight line method, wherein the method specifically comprises the following steps:

wherein L represents the total curve length, Δ x_iThe arc length in the ith subinterval corresponds to the transverse axis length, Δ y_iThe arc length in the ith sub-interval corresponds to the longitudinal axis length, x_iIs the abscissa of any point in the ith subinterval, i is the curve dividing interval number, wherein i is 1,2, …, n, f'²(x_i) Is the equation f (x) at x for the water system or road curve_iThe derivative of the point is squared, and finally n equal curve lengths are accumulated to obtain the total curve length L, i.e. the length of the corresponding water system or road path penetrating the upper and lower boundaries of the dense section, and the curve lengths of all water systems or road vectors in the dense section are calculated to make L_min＝mm{L₁，L₂，...，L_NSelecting L_minAs a water system or road vector shortest path constraint. Obtaining the optimal mosaic line segment A according to the method in step S514₂A₅。

S515, the coordinate of the end point of the shortest water system or road vector line is coincident with the coordinate of the end point of the lower boundary of the reference mosaic line in the dense section or d_minLess than or equal to 10mUnder the condition, the shortest water system or road vector line end point and the mosaic reference line end point in the step S2 are directly connected with the end point in the lower boundary y20 of the dense section, the distance between two intersection points of the shortest water system or road vector line and the mosaic reference line on the lower boundary of the dense section is taken as the diameter of a circle, the middle points of the two points are taken as the circle centers to establish a buffer circle, and the optimal mosaic line segment A in the step S514 is taken as the circle center in the circle₂A₅Intersection A with lower boundary y20 of dense segment₅(xx5, yy5) as the starting point, searching for a water system or road vector, and selecting a mosaic line A in the section according to the steps from S511 to S514₅A₃。

S516, the dense section inlaid line is obtained as the optimal inlaid line according to the steps from S511 to S515, in this embodiment, A is obtained₂A₅、A₅A₃As an optimal mosaic line segment, fig. 8 is a schematic diagram of a mosaic line generation process for dense segments and sparse segments.

S52, for the sparse segment, when there is no secondary feature library vector, the original mosaic reference line in the step S2 is reserved as the mosaic line segment of the region, that is, the intersection point A of the top right corner A of the start of the overlap region and the sparse segment boundary in the step S1₁Constructed curve AA₁；

When a secondary feature library vector exists, an automatic mosaic line generating method based on a road and a water system vector is adopted, and the method specifically comprises the following steps:

s521, obtaining the optimal mosaic line AA with priority on the road by taking the point A as the starting point according to the selection mode of the through road of the first road and the second water system in the sparse terrain section and according to the steps S511 to S515₄If there is no road or water system that satisfies the condition in step S514, the initial mosaic reference line in step S2 is selected.

S522, the shortest road or water system vector line terminal point and the reference mosaic line are coincided in the boundary coordinates of the sparse section or d_minUnder the condition that the distance is less than or equal to 10m, directly connecting the end point of the shortest road or the water system vector line in the lower boundary y30 of the sparse section with the end point of the embedded reference line in the step S2; distance d between terminal and reference line_minWhen the distance is more than 10m, the distance between two intersection points of the shortest water system or the road vector line and the reference inlaid line on the lower boundary of the sparse section is used as the diameter of a circle, the midpoint of the two intersection points is used as the center of the circle to establish a buffer circle, and the buffer circle is arranged on the circleOptimal mosaic line segment AA in medium sparse segment₄Intersection A with sparse segment boundary y30₄Searching water system or road vector as starting point, and selecting mosaic line A according to steps S511-S514₄A₁。

S523, obtaining AA by taking the mosaic line penetrating through the sparse segment obtained in the steps from S521 to S522 as an optimal mosaic line₄、A₄A₁As an optimal mosaic line segment, as shown in fig. 8.

S53, aiming at the generation of mosaic lines in the mixed sections except the dense sections and the sparse sections, respectively selecting the intersection points of the initial mosaic reference lines and the sparse sections and the intersection points of the initial mosaic reference lines or the boundary lines of the dense sections in the step S2 as the starting points and the ending points of the mosaic line sections, and adopting the initial mosaic reference lines in the step S2 as the mosaic lines.

And S6, connecting the embedded wires among the sections.

The hybrid section uses the initial mosaic reference line; the dense segment uses the feature vector optimal mosaic line a acquired in step S51₂A₅、A₅A₃(ii) a The sparse segment uses the optimal mosaic line AA obtained in step S52₄、A₄A₁. Finally, the endpoints of the segment-to-segment mosaic lines on the segment boundary lines of the mixed segment, the sparse segment and the dense segment are combined into a complete optimal mosaic line vector in the overlapping region, as shown in fig. 9, i.e. the complete mosaic line AA in the embodiment₄A₁A₂A₅A₃B。

S7, after generating a complete mosaic line, storing the mosaic line into a feature library, and taking the mosaic line as the mosaic line of the next mosaic area under the condition that the variation difference of the ground features is small; under the condition of large terrain variation difference, the obtained complete mosaic line is taken as a mosaic reference line, and the mosaic line optimization algorithm of the steps S1 to S6 in the embodiment is continued, so that the mosaic line generation efficiency in the region mosaic is improved.

The above description is only an example of the present application and is not intended to limit the present application. Various modifications and changes may occur to those skilled in the art. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should be included in the scope of the claims of the present application.

Claims (7)

1. A method for generating optimal mosaic lines of satellite images based on intelligent decision of a feature library is characterized by comprising the following specific steps,

s1, acquiring an effective overlapping area of the satellite images; firstly, acquiring an image effective area of a satellite image, and then acquiring an overlapping area between adjacent images;

s2, generating an initial mosaic reference line;

according to the overlapping area between the adjacent images acquired in the step S1, generating an initial mosaic reference line by adopting a mosaic line automatic extraction method based on pixel difference feathering, wherein the initial mosaic reference line is represented by adopting a coordinate sequence in the overlapping area;

s3, an overlapping area segmentation method based on a primary feature library, comprising the following steps:

s31, selecting and grading a feature library;

s32, dense segment segmentation:

s33, sparse segment segmentation;

s4, matching the road vector, the water system vector data and the image pixel data in the secondary feature library,

s5, automatically generating the segment embedded line based on the secondary feature library, which specifically comprises:

s51, aiming at the dense segment, adopting an automatic mosaic line generating method based on a water system vector and a road vector feature vector, which comprises the following specific steps:

when there is no secondary feature library data, the mosaic line is the initial mosaic reference line in step S2;

when secondary feature library data exist, an automatic mosaic line selection method based on a water system vector and a road vector is adopted to obtain mosaic lines:

the starting point of the optimized mosaic line in the dense section is A₂And fitting the water system or road vector coordinates of the secondary feature library by selecting N lower intersection points of the water system or road vector curve in the overlapping area and the lower boundary straight line y20 of the dense section, namely the dense section road end point to obtain a water system or road vectorA quantity curve, wherein the water system vector curve is:

wherein the content of the first and second substances,

polynomial coefficient of water-based vector curve, N₀The order of the water system or the road vector curve is as follows:

wherein the content of the first and second substances,

polynomial coefficients of the road vector curve;

s511, if the secondary feature library has water system vector data, collecting a water system vector curve f₁(x)，f₂(x)，f₃(x),...,f_N(x) Forming an equation set in a simultaneous manner with the straight line y20, and solving the intersection point set to obtain a water system vector terminal set in the dense section;

s512, if the secondary feature library has no water system vector data, the road vector curve h is processed₁(x),h₂(x),...,h_N(x) The curve set and the straight line y20 are combined to form an equation set, and a road vector terminal set in the dense section is obtained through solution;

s513, taking the intersection point of the reference mosaic line of the overlapping area obtained in the step S2 and the upper border line y10 of the dense section as a starting point, taking the continuous water system or road vector point which penetrates through the upper border and the lower border of the dense section on the lower border line y20 of the dense section as an end point, searching the penetrated water system or road according to constraint conditions, searching M water systems or roads, wherein N water systems or roads (N is less than or equal to M) are water systems or road vectors containing end points, selecting the searched continuous water system or road vector line which penetrates through the upper border and the lower border of the dense section as an optimal mosaic line until the continuous water system or road vector line which penetrates through the upper border and the lower border of the dense section is searched at the end point on the lower border line y20 of the dense section, if the penetrated water system or road is not searched, discarding the road or water system vector data, and selecting the initial mosaic reference line in the step S2 as the optimal mosaic;

s514, selecting the optimal mosaic line under the condition of running through of the water system or road vector characteristic data;

the process of judging the penetration of the water system or road vector characteristic data comprises the steps of fitting the coordinates of discrete points of each water system or road into a water system or road vector curve equation, substituting the coordinates of the starting point and the end point of the boundary point of each water system or road into y10 and y20 equations of a boundary straight line, judging that the water system or road vector penetrates through the upper and lower boundaries of a dense or sparse section when the coordinates of the end point of the water system or road vector are on the boundary straight line, selecting the penetrated road or water system vector data, discarding the road or water system vector data if the coordinates of the end point of the water system or road vector are not on the boundary straight line, and selecting an initial embedded line in the step S2 as an embedded line in the section;

the selection process of the inlaid wire comprises the following steps: selecting a vector of a shortest penetrating path according to the sequence of a water system, a high-level road and a low-level road; then according to the distance from the lower boundary end point of the water system or road vector line in the dense section to the intersection point of the reference mosaic line in the overlapping area in the lower boundary of the dense section, according to the principle of the shortest distance priority, selecting the water system or road vector corresponding to the shortest distance end point as the mosaic line after optimization, and selecting the distance from the lower boundary end point of the water system or road vector line in the dense section to the intersection point of the reference mosaic line in the lower boundary of the dense section

Constraint of proximity of water system or road vector to reference line, d_iRepresents the distance between the terminal point of the ith water system or road vector line at the lower boundary of the dense section and the intersection point of the reference mosaic line at the lower boundary of the dense section, N₁The number of the searched water penetrating systems or road vectors is set; the method for calculating the length of the water system or the road path comprises the following steps: equally dividing the water system or the road curve penetrated in the step S51 into n equal parts according to the actual road vector length, and solving the total curve length of the equally divided curve by adopting a short straight line method, wherein the method specifically comprises the following steps:

wherein L represents the total curve length, Δ x_iThe arc length in the ith subinterval corresponds to the transverse axis length, Δ y_iThe arc length in the ith sub-interval corresponds to the longitudinal axis length, x_iIs the abscissa of any point in the ith subinterval, i is the curve dividing interval number, wherein i is 1,2, …, n, f'²(x_i) Is the equation f (x) at x for the water system or road curve_iThe derivative of the point is squared, and finally n equal curve lengths are accumulated to obtain the total curve length L, i.e. the length of the corresponding water system or road path penetrating the upper and lower boundaries of the dense section, and the curve lengths of all water systems or road vectors in the dense section are calculated to make L_min＝min{L₁,L₂,...,L_NSelecting L_minAs a water system or road vector shortest path constraint; obtaining the optimal mosaic line segment A according to the method in step S514₂A₅；

S515, the coordinate of the end point of the shortest water system or road vector line is coincident with the coordinate of the end point of the lower boundary of the reference mosaic line in the dense section or d_minUnder the condition that the distance between the shortest water system or road vector line end point and the mosaic reference line end point in the lower boundary of the dense section is less than or equal to 10m, directly connecting the shortest water system or road vector line end point with the mosaic reference line end point in the step S2, taking the distance between two intersection points of the shortest water system or road vector line and the mosaic reference line in the lower boundary of the dense section as the diameter of a circle, taking the middle points of the two points as the circle centers to establish a buffer circle, and taking the optimal mosaic line segment A in the step₂A₅Intersection A with lower boundary y20 of dense segment₅(xx5, yy5) as the starting point, searching for a water system or road vector, and selecting a mosaic line A in the section according to the steps from S511 to S514₅A₃；

S516, obtaining the embedded line in the dense section as an optimal embedded line according to the steps from S511 to S515;

s52, aiming at the sparse segment, when there is no secondary feature library vector, keeping the original mosaic reference line in the step S2 as the mosaic line segment of the region;

when a secondary feature library vector exists, an automatic mosaic line generating method based on a road and a water system vector is adopted, and the method specifically comprises the following steps:

s521, obtaining the optimal mosaic line AA with priority on the road by taking the point A as the starting point according to the selection mode of the through road of the first road and the second water system in the sparse terrain section and according to the steps S511 to S515₄If there is no road or water system satisfying the condition in step S514, selecting the initial mosaic reference line in step S2;

s522, the shortest road or water system vector line terminal point and the reference mosaic line are coincided in the boundary coordinates of the sparse section or d_minUnder the condition that the distance is less than or equal to 10m, directly connecting the end point of the shortest road or the water system vector line in the lower boundary y30 of the sparse section with the end point of the embedded reference line in the step S2; distance d between terminal and reference line_minWhen the distance between two intersection points of the shortest water system or the road vector line and the reference mosaic line on the lower boundary of the sparse section is larger than 10m, the diameter of the circle is taken as the distance between the two intersection points, the midpoint of the two intersection points is taken as the center of the circle to establish a buffer circle, and the optimal mosaic line segment AA in the sparse section is taken as the center of the circle₄Intersection A with sparse segment boundary y30₄Searching water system or road vector as starting point, and selecting mosaic line A according to steps S511-S514₄A₁；

S523, obtaining AA by taking the mosaic line penetrating through the sparse segment obtained in the steps from S521 to S522 as an optimal mosaic line₄、A₄A₁As an optimal mosaic line segment;

s53, aiming at the generation of mosaic lines in the mixed sections except the dense sections and the sparse sections, respectively selecting the intersection points of the initial mosaic reference lines and the sparse sections and the intersection points of the initial mosaic reference lines or the boundary of the dense sections in the step S2 as the starting points and the ending points of the mosaic line sections, and adopting the initial mosaic reference lines in the step S2 as mosaic lines;

s6, connecting the embedded wires among the sections;

the hybrid section uses the initial mosaic reference line; the dense segment uses the feature vector optimal mosaic line a acquired in step S51₂A₅、A₅A₃(ii) a The sparse segment uses the optimal mosaic line AA obtained in step S52₄、A₄A₁(ii) a Finally, embedding lines between the sections on the section edges of the mixed section, the sparse section and the dense sectionThe upper end points of the boundary lines are combined into a complete optimal mosaic line vector in the overlapping area;

s7, after generating a complete mosaic line, storing the mosaic line into a feature library, and taking the mosaic line as the mosaic line of the next mosaic area under the condition that the variation difference of the ground features is small; and under the condition of large terrain variation difference, taking the obtained complete mosaic line as a mosaic reference line, and continuing the mosaic line optimization method of the steps from S1 to S6, thereby improving the mosaic line generation efficiency in the region mosaic.

2. The method for generating optimal mosaic lines for satellite imagery based on feature library intelligent decision making according to claim 1, wherein the step of obtaining the image effective area of the satellite imagery in step S1 comprises:

if the obtained satellite image is a multiband image, adding multiband pixel values of each pixel point in the satellite image and storing the added value, and if the added value is 0, judging that the added value is an invalid value;

if the obtained satellite image is a single-band image, sequentially circulating from the first line to the subsequent lines according to the number of rows and columns of the image, judging each pixel point value from left to right in each line, continuously judging the next pixel point value if the pixel point value is 0, recording the corresponding row number and column number of the pixel point when the first non-0 pixel point is judged, setting the point as an upper left corner C point, acquiring the upper left corner point, sequentially circulating from the last line to the upper line, judging whether each pixel point value is 0 from right to left in each line, and setting the point as a lower right corner D point when the first non-0 pixel point value is met; continuously circulating each row of the image from left to right, judging each pixel point value from bottom to top in each row, continuously judging the next pixel point value if the value is 0, recording the corresponding row number column number of the first non-0 pixel point when the first non-0 pixel point is judged, setting the point as a lower left corner B point, circularly reading the image rows from right to left from the rightmost row after the lower left corner point is obtained, judging each pixel point value from top to bottom in each row, and setting the point as an upper right corner A point when the first non-0 pixel point value is met;

creating a vector file according to the obtained row and column numbers of the four corner points to finish the acquisition of the effective area of the satellite image;

and acquiring the overlapping area between the adjacent images, acquiring the corresponding amplitude data vectors of the two adjacent images, respectively superposing the respective vector data of the two images, calculating the overlapping edge of the two images, and acquiring the effective overlapping area range of ACBD.

3. The method as claimed in claim 1 or 2, wherein the automatic extraction method of mosaic line based on pixel difference feather in step S2 is implemented by Dijkstra' S method for solving the shortest path of minimum cost (cost) between two points, wherein the cost is replaced by weighted pixel path, the normalized pixel value S (I, j) is distinguished by using the brightness difference Δ I (I, j) and the gradient of the overlapping region as the pixel difference value of the position of the overlapping region, the gray gradient difference Δ g is the absolute value of the gradient value difference, the gradient value is g, the change of the image gray in the horizontal and vertical directions, and the gradient direction is g

The specific calculation process of the luminance difference Δ I (I, j) is:

ΔI(i,j)＝|img1(i,j)-img2(i,j)|，

wherein img1(i, j) and img2(i, j) respectively represent the brightness values of the points (i, j) of the overlapping area on two adjacent images to be mosaiced; magnitude g (i, j) of gradient value and gradient direction of point (i, j) of overlapping region

The calculation formula of (2) is as follows:

g(i,j)＝[(img(i,j)-img(i+1,j)²+(img(i,j)-img(i,j+1))²)]^1/2

img (i, j) represents the luminance value of the point (i, j) of the overlap region, and the gray gradient difference Δ g is represented as:

wherein Δ g (i, j) represents the gray gradient difference value at the pixel (i, j), g1(i, j) represents the gradient value of the pixel in one of the two adjacent images, g2(i, j) represents the gradient value of the pixel in the other of the two adjacent images,

respectively representing the gradient directions corresponding to the two pixel points; the maximum value of the normalized brightness difference or gray gradient difference represents the maximum difference value of the pixels in the overlapping area, the values of the maximum difference value are max (delta I (I, j), delta g (I, j)), the corresponding pixel points are called pixel difference points, the set of the pixel difference points is CC, the upper right corner point of the boundary of the image in the overlapping area is used as the initial point A of the overlapping area, the coordinates of the initial point A are (x0, y0), and the initial point A reaches all the pixel difference points a in the overlapping area_iBelongs to CC, and the shortest path in the overlap region is recorded as Dist [ A, a ]_i]And i is a subscript, then:

Dist[A,a_i]≥short[A,a_i]，

wherein short [ A, a_i]Representing the difference point a from the start point A of the overlap region to all pixels_iBelongs to the shortest path of CC; by continuously calculating Dist [ A, a ]_i]Up to a_iTo the lower left corner B of the overlap region in step S1, at which time Dist [ A, a ]_i]＝short[A,a_i]And obtaining an initial mosaic reference line according to the mosaic line generation method based on the pixel feather difference.

4. The method for generating optimal mosaic lines for satellite images based on intelligent decision of feature library as claimed in claim 1, wherein said step S31, feature library selection and classification, specifically comprises:

carrying out regional object feature classification on an overlapping region in the obtained satellite image, and establishing a primary feature library and a secondary feature library, wherein the primary feature library is used for determining a segmentation boundary, and the secondary feature library is used for determining the selection of a mosaic line; selecting a primary characteristic library and a secondary characteristic library from the characteristic library data, wherein:

the primary feature library comprises dense boundary vector data and sparse boundary vector data;

the secondary feature library comprises road vector data, water system vector data and scarp data;

according to the ground feature type and feature library data established by various remote sensing geographic information resources, dividing an overlap region into a dense section, a sparse section and a mixed section according to dense boundary vector data and sparse boundary vector data in a primary feature library, obtaining the dense section in the overlap region according to step S32 according to the dense boundary vector data, obtaining the sparse section in the overlap region according to step S33 according to the sparse boundary vector data, wherein the regions except the dense section and the sparse section are collectively called as the mixed section, and the specific segmentation method is shown in steps S32 and S33; the dense boundary vector data establishes a feature library according to buildings, high-grade road data, noctilucent remote sensing data, water system vectors and the like, and the sparse boundary vector data establishes a feature library according to desert vectors, DEM data, expressways, low-grade road data, water system vectors and the like.

5. The method for generating optimal mosaic lines for satellite images based on intelligent decision of feature library as claimed in claim 1, wherein said step S32, dense segment segmentation, specifically comprises:

s321, firstly, judging whether dense boundary vector data exists in an overlapping area, if so, segmenting according to the step S322, if not, segmenting, and directly selecting an initial mosaic reference line in the step S2;

s322, constructing a minimum circumscribed rectangle according to the dense boundary vector data in the primary feature library, wherein the specific construction process is as follows:

drawing the minimum circumscribed rectangle by utilizing a minAreaRect () function meter in opencv, wherein the minimum circumscribed rectangle comprises four end points and four edges of the minimum circumscribed rectangle: firstly, defining a circumscribed rectangle set to be a circumscribed rectangle set formed by the outlines of the dense boundary vectors and the sparse boundary vectors in the overlapping area ABCD in the step S1, and then defining a minimum circumscribed rectangle set to be a set formed by the minimum circumscribed rectangle of the outlines of each dense boundary vector and each sparse boundary vector; the outline here refers to a set of points within the overlap area ABCD obtained in step S1 corresponding to the boundary vector; then drawing the central point of the minimum circumscribed rectangle to obtain four end points of the minimum circumscribed rectangle; finally, drawing four sides of the minimum circumscribed rectangle to finish drawing the minimum circumscribed rectangle;

rotating a polygon formed by original dense boundary vector data end points, wherein the rotation degree is-90 degrees, obtaining a simple external rectangle of the polygon after rotating each degree, recording the area, vertex coordinates and rotation degree of each simple external rectangle, obtaining the external rectangle with the minimum area, and finally rotating the external rectangle to a rotation angle with the same size and the direction opposite to the recorded rotation degree, thus obtaining the minimum external rectangle;

s323, in the case where the minimum bounding rectangle obtained in step S322 intersects with or includes the bounding rectangles of the remaining boundary vectors, segmenting the overlap region in step S1 according to the dense boundary vectors;

s324, according to the upper and lower boundaries of the minimum bounding rectangle of the dense boundary vector in step S322, the minimum bounding tangent equation is expressed as: and y is kx + b, and the equations of the upper minimum external tangent line and the lower minimum external tangent line are solved according to the coordinates of the four end points of the minimum external rectangle:

y₁＝k₁x+b₁

y₂＝k₂x+b₂，

wherein the slope k₁＝k₂，b₁And b₂Is a constant term of the two smallest circumscribed tangents.

6. The method for generating optimal mosaic lines for satellite images based on intelligent decision of feature library as claimed in claim 1, wherein said step S33, sparse segment segmentation, specifically comprises:

first, it is determined whether the overlap region has sparse boundary vector data, e.g.If sparse boundary vector data exists, constructing a minimum circumscribed rectangle according to step S322, obtaining a lower boundary intersection point S by intersecting the boundary of the overlapping region with the boundary of the known sparse segment, and regarding the sparse segment segmentation boundary line y30, the expression is y₃＝k₃x+b₃Slope k of₃＝k₁＝k₂Substituting the coordinate of the lower boundary intersection point S into the expression of y30 to calculate b₃Values, resulting in an expression for sparse segment segmentation line y 30; if no sparse boundary vector data exists, the initial mosaic reference line in the step S2 is directly selected without segmentation, and the other segments except the dense segment and the sparse segment are collectively called as a mixed segment.

7. The method for generating optimal mosaic lines for satellite images based on intelligent decision of feature library as claimed in claim 1, wherein said step S4 specifically comprises:

the three-parameter method is that translation of X, Y and Z axes is carried out between two reference surfaces, and according to three parameters of translation and reference points, the road vector and water system vector coordinates of the input overlapping area image are converted by the three-parameter method to obtain output vector geographic coordinate data, and the calculation is as follows:

wherein [ X, Y, Z]_newTo obtain vector data, [ X, Y, Z]_orginAs original vector coordinates, [ dX, dY, dZ]Three parameters for translation;

matching road vector data and water system vector data with geographic entity coordinate points to an overlapping area image according to prior feature data of a secondary feature library, and matching entities with the same name; the position of the pixel point in the image is represented by the pixel coordinate, when the road is matched with the water system vector data and the image pixel of the overlapping area, the corresponding geographic coordinate needs to be converted into the pixel coordinate, and the conversion method comprises the following steps:

first, the geographic coordinate [ GEO ] of the upper left corner of the overlap region image_x,GEO_y]Known, and the resolution R of the image_x,R_yCalculating the geographic coordinate of any pixel position according to the geographic coordinate of the upper left corner of the image and the resolution of the image; the formula for converting the geographic coordinates to pixel coordinates is:

the formula for converting the pixel coordinates into the geographic coordinates is as follows:

X＝GEO_x+col*R_x

Y＝GEO_y+row*R_y，

wherein X and Y in the above two coordinate transformation formula represent geographical coordinates, PIX_x，PIX_yRepresenting pixel coordinates, GEO_x、GEO_yRepresenting the geographic coordinates of the upper left corner of the image, col, row representing the number of columns and rows of pixels of the image, R_xAnd R_yA row and column resolution representing an image; and matching the road vector with the A point coordinate as a starting point to the image pixel to prepare as data generated by a subsequent mosaic line.

Images