Photoelectric image and three-dimensional virtual scene matching method based on helicopter

Abstract

The invention relates to a comprehensive visual system design technology in the field of avionics, in particular to a photoelectric image and three-dimensional virtual visual matching method based on a helicopter. The method comprises the following steps: keeping a viewport in view rendering consistent with a viewport of a photoelectric sensor; carrying out denoising operation on the photoelectric image and the three-dimensional virtual scene of the helicopter collected by the two viewports, and removing a sky area; feature point extraction is carried out on the photoelectric image after the sky area is removed and the three-dimensional virtual scene, and feature points are used as control points to be matched; finding out common characteristics of the two images from the to-be-matched control points of the photoelectric image and the to-be-matched control points of the three-dimensional virtual scene to obtain high-precision control point pairs; aligning the photoelectric image and the three-dimensional virtual scene according to the high-precision control point pair; and fusing the photoelectric and virtual visual images based on the aligned image and visual scene.

Description

Technical Field

[0001] The invention relates to a comprehensive visual system design technology in the field of avionics, and more specifically to a method for matching a helicopter-based optoelectronic image with a three-dimensional virtual visual scene. Background Art

[0002] The integrated vision system uses digital graphics and images and data processing methods to generate and display realistic external scenes based on geographic information and precise navigation information, combined with information detected by optoelectronic, radar and other sensors through multi-source fusion processing. It provides pilots with surrounding environment perception and threat warnings in complex meteorological and geographical conditions, which can effectively improve flight safety and visual navigation capabilities.

[0003] As photoelectric images have night vision capabilities and terrain data has the advantages of prior information and wide field of view, the photoelectric sensor image information captured in real time is registered and fused with the three-dimensional virtual view generated based on the digital map, and texture mapping is established to achieve enhanced display of virtual view and reduce the perceived difference between virtual view and real environment. For the matching of photoelectric images and three-dimensional virtual view, the image perspective is oblique downward. The difficulty lies in the fact that the image is distorted to varying degrees due to shooting angle, lens distortion and other reasons. The oblique deformation will change the local geometric structure of the image, making it difficult for conventional feature point detection algorithms to remain consistent under different viewing angles or shooting conditions. Such deformation causes the extracted feature points to no longer be consistent and repeatable between different images, thereby affecting the stability of the feature points. Due to the instability of feature points, traditional matching algorithms are prone to errors in oblique images. Therefore, a more robust matching algorithm is needed to cope with complex deformation situations.

[0004] Therefore, how to provide a method for matching helicopter-based optoelectronic images with three-dimensional virtual scenes is an urgent problem that those skilled in the art need to solve. Summary of the invention

[0005] Purpose of the invention: The present invention provides a method for matching optoelectronic images and three-dimensional virtual scenes based on a helicopter, which uses more robust corner point features and line features as matching feature points, uses a local transformation model to describe the transformation relationship between images, uses a feature space template matching method to obtain high-precision control point pairs, and uses a GPU+CPU joint optimization method to improve the algorithm operation efficiency, improve the matching accuracy and matching time of optoelectronic images and three-dimensional virtual scenes, and solve the current problems of low matching accuracy and poor reliability.

[0006] In order to achieve the above object, the present invention adopts the following technical solution:

[0007] The present invention provides a method for matching an optoelectronic image and a three-dimensional virtual scene based on a helicopter, comprising:

[0008] Keep the viewport when rendering the scene consistent with the viewport of the photoelectric sensor;

[0009] De-noise the helicopter's electro-optical images and 3D virtual scenes collected from two viewports and remove the sky area;

[0010] Extract feature points from the photoelectric image and the three-dimensional virtual scene after removing the sky area, and use the feature points as control points to be matched;

[0011] Find the common features of the two images from the control points to be matched in the photoelectric image and the control points to be matched in the three-dimensional virtual scene, and obtain high-precision control point pairs;

[0012] Align the photoelectric image with the 3D virtual scene based on high-precision control point pairs;

[0013] Based on the aligned images and views, optoelectronic and virtual vision images are fused.

[0014] Find the common features of the two images from the control points to be matched in the photoelectric image and the control points to be matched in the three-dimensional virtual scene, including:

[0015] The image edges of the photoelectric image and the virtual scene are extracted by edge detection operators to obtain two feature space images;

[0016] Calculate the similarity measure of the two feature space images for further adjustment and calibration to match high-precision control point pairs;

[0017] Calculate the main peak sharpness and the main-to-minor peak ratio of the two feature space images;

[0018] Correlation surface analysis is performed on the two parameters to remove mismatched points.

[0019] Alignment between images is achieved using a local weighted mean correction model;

[0020] The local correction model f(x,y) is defined as follows:

[0021]

[0022] Among them, W i (x,y) is the inverse distance weight, p i (x, y) represents a quadratic polynomial transformation through the neighborhood of each control point, and N is the number of control point pairs.

[0023] Make the viewport when rendering the scene consistent with the viewport of the photoelectric sensor, including:

[0024] According to the installation positions of the photoelectric sensor and the inertial navigation satellite antenna of the carrier aircraft, the viewport positions of the two are calculated; and the viewport during scene rendering is adjusted to be consistent with the viewport of the photoelectric sensor.

[0025] The specific steps for extracting the sky area are:

[0026] Firstly, the Sobel operator is used to extract the image gradient features, and the boundary line of the sky is extracted from the gradient feature map according to a given threshold.

[0027] An energy function is constructed to optimize the extraction of the sky boundary line. When the energy function is maximum, the sky boundary line is most accurate.

[0028] By detecting whether the boundary line has a large jump, we can determine whether the extracted sky area contains a large number of buildings. If so, we perform cluster analysis on the sky area and divide it into two categories. The real sky area is obtained by the Mahalanobis distance with the non-sky area.

[0029] Scan each column of pixels and determine whether the pixel belongs to the sky through the Mahalanobis distance between each column of pixels and the real sky area, thereby improving the sky boundary line.

[0030] Feature points are extracted from the photoelectric image and the three-dimensional virtual scene after the sky area is removed, and the feature points are used as control points to be matched, including:

[0031] Corner point detection and line segment detection are used to extract feature points from the preprocessed virtual scene and optoelectronic image respectively.

[0032] The specific steps of corner detection are:

[0033] For a two-dimensional image l(x, y), when a point at (x, y) is translated by (Δx, Δy), the autocorrelation function c(x, y; Δx, Δy) is used to represent the similarity before and after the window moves:

[0034]

[0035] Where: w(x, y) is the window centered at (x, y); w(u, v) is the weighted function of the sliding window centered at (u, v), which is a Gaussian function;

[0036] Perform a first-order Taylor expansion on I(u+Δx, v+Δy) and ignore the second-order higher-order terms:

[0037] I(u+Δx,v+Δy)=I(u,v)+Ix(u,v)Δx+Iy(u,v)Δy;

[0038] Where: Ix(x, y) and Iy(u, v) are partial derivatives of the image. Substituting into the above formula, we get:

[0039]

[0040] in:

[0041]

[0042] Where: Matrix M is the covariance matrix of the image gradient; λ1 and λ2 are the eigenvalues ​​of the matrix M;

[0043] The corner response function R(x, y) is defined by using the matrix trace Trace(M) and determinant det(M), and its expression is:

[0044] R(x, y) = det(M) - k(Trace(M)) 2 ;

[0045] Where: det(M) = λ1λ2; Trace(M) = λ1+λ2; k is an empirical constant, k∈[0.04, 0.06]; given a threshold T, when R(x, y)>T and reaches the maximum in the local window, the pixel point (x, y) is a feature point.

[0046] The specific steps of line segment detection are:

[0047] Assume that S points are close enough to the ideal straight line, then these points belong to the straight line, and the points in the remaining subset are considered to be noise points; in this model, the failure rate is:

[0048]

[0049] Among them, aS / M refers to the proportion of points in the subset that belong to the ideal straight line in the entire subset; the performance of the transformation is determined by the expected value of m*αS / M and s, as follows:

[0050] (1) If the number of points m contained in the subset is too small, the failure rate of the transformation is close to l, that is:

[0051]

[0052] (2) If m is large enough, m*αS / M ≥ 1, the transformation failure rate is low, and its failure rate is:

[0053]

[0054] At this time, when m*αS / M≥2, the larger m is, the lower the failure rate is;

[0055] After obtaining the detected line segment, the line segment control points can be obtained through line segment sampling.

[0056] Energy function J n is defined as:

[0057]

[0058] where ∑ s and∑ g are the covariance matrices of pixels described by RGB values ​​in the sky and ground regions, respectively; they are both 3×3 matrices defined as follows:

[0059]

[0060] Where N s and N g are the number of pixels in the sky and ground regions, respectively, and μ s and μ g 3×1 column vectors representing the average RGB values ​​in the sky and ground regions, respectively;

[0061] and are the eigenvalues ​​corresponding to the above two matrices, |·| represents the determinant, which measures the variance or volume of the pixel distribution in each area; where γ represents the emphasis on the uniformity of the sky area.

[0062] The beneficial effects of the present invention are as follows: the present invention solves the problem of inaccurate superposition of photoelectric images and three-dimensional virtual vision for an airborne integrated vision system, and can effectively improve the accuracy and matching time of virtual-real matching of the airborne integrated vision system, so that pilots can clearly see the runway and the surrounding environment even in a low visibility environment, thereby ensuring safe takeoff and landing of the aircraft. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings in the following description are only embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying creative work.

[0064] Figure 1 The invention is a flow chart of a method for matching helicopter-based optoelectronic images with three-dimensional virtual scenes.

[0065] Figure 2 It is a specific flow chart for implementing image preprocessing.

[0066] Figure 3 It is a specific flow chart for realizing image feature extraction.

[0067] Figure 4 It is a specific flow chart for realizing precise image registration processing. DETAILED DESCRIPTION

[0068] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0069] The present invention proposes a method for matching a helicopter-based optoelectronic image with a three-dimensional virtual scene, comprising the following steps:

[0070] S1: Calculate according to the installation position of the photoelectric sensor and the inertial navigation satellite antenna of the carrier, transform the position data in the current navigation information to the installation position of the photoelectric sensor as the coordinate origin, and ensure that the viewport during scene rendering is consistent with the viewport of the photoelectric sensor;

[0071] S2: Perform denoising on the image to reduce noise and interference in the image and retain clearer and more accurate feature information. At the same time, detect the sky area of ​​the virtual scene and remove points that are close to the sky area to improve the registration accuracy;

[0072] S3: corner point detection and line segment detection are used to extract feature points of the preprocessed virtual scene and optoelectronic image respectively, and the feature points are used as control points to be matched;

[0073] S4: Construct a feature space to find the common features of multi-source images, extract image edges through edge detection operators, and calculate the similarity measure of the two feature space images for further adjustment and calibration to obtain high-precision control point pairs; use the two parameters of main peak sharpness and main-to-secondary peak ratio in correlation surface analysis to perform correlation surface analysis to remove mismatched points as much as possible.

[0074] S5: Based on the high-precision control point pairs, a local weighted mean correction model is used to perform more accurate alignment between the real-time image and the virtual image.

[0075] The local correction model is defined as follows:

[0076]

[0077] Among them, W i (x,y) is the inverse distance weight, p i (x,y) represents a quadratic polynomial transformation through the neighborhood of each control point.

[0078] S6: A fusion algorithm for heterogeneous images based on multi-scale image block structured decomposition is used. This method performs multi-scale decomposition and reconstruction of infrared images and virtual vision images based on a multi-layer pyramid structure. Image block structured decomposition is used to extract image block features in each scale layer obtained by decomposition, and specific strategies are used at different layers to fuse optoelectronic and virtual vision image features.

[0079] Furthermore, the specific steps of extracting the sky area in S2 are:

[0080] a) First, the Sobel operator is used to extract the image gradient features, and the boundary line of the sky is extracted from the gradient feature map according to a given threshold;

[0081] b) constructing an energy function to optimize the extraction of the sky boundary line. When the energy function is maximum, the sky boundary line is most accurate;

[0082] c) By detecting whether the boundary line has a large jump, it is determined whether the extracted sky area contains a large number of buildings. If so, the sky area is clustered and divided into two categories, and the real sky area is obtained by the Mahalanobis distance with the non-sky area;

[0083] d) Scan each column of pixels and determine whether the pixel belongs to the sky by the Mahalanobis distance between each column of pixels and the real sky area, so as to improve the sky boundary line.

[0084] Furthermore, the specific steps of corner feature extraction in S3 are:

[0085] For a two-dimensional image I(x, y), when a point is translated by (Δx, Δy) at (x, y), the autocorrelation function is used to represent the similarity before and after the window moves:

[0086]

[0087] Where: w(x, y) is the window centered at (x, y); w(u, v) is the weighting function within the window, generally a Gaussian function. I(u+Δx, v+Δy) is expanded by the first order Taylor and the second order higher terms are ignored:

[0088] I(u+Δx,v+Δy)≈I(u,v)+Ix(u,v)Δx+Iy(u,v)Δy

[0089] Where: Ix(u, v) and Iy(u, v) are partial derivatives of the image. Substituting into the above formula, we get:

[0090]

[0091] in:

[0092]

[0093] Where: matrix M is the covariance matrix of image gradient; λ1 and λ2 are the eigenvalues ​​of matrix M.

[0094] The corner response function R(x, y) is defined by using the matrix trace Trace(M) and determinant det(M), and its expression is:

[0095] R(x, y) = det(M) - k(Trace(M)) 2

[0096] Where: det(M) = λ1λ2, Trace(M) = λ1+λ2; k is an empirical constant, k∈[0.04, 0.06]. Given a threshold T, when R(x, y)>T and reaches the maximum in the local window, the pixel (x, y) is a feature point.

[0097] like Figure 1 A method for matching an optoelectronic image based on a helicopter with a three-dimensional virtual scene is shown, comprising the following steps:

[0098] S1: Calculate according to the installation position of the photoelectric sensor and the inertial navigation satellite antenna of the carrier, transform the position data in the current navigation information to the installation position of the photoelectric sensor as the coordinate origin, and ensure that the viewport during scene rendering is consistent with the viewport of the photoelectric sensor;

[0099] S2: perform image preprocessing on the photoelectric image and the three-dimensional virtual scene respectively;

[0100] Specifically, the image preprocessing process is as follows: Figure 2 As shown: the image is denoised to reduce noise and interference in the image and retain clearer and more accurate feature information. At the same time, the sky area of ​​the virtual scene is detected and points that are close to the sky area are removed to improve the registration accuracy.

[0101] Among them, the specific strategy for extracting the sky area is:

[0102] e) Firstly, the Sobel operator is used to extract the image gradient features, and the boundary line of the sky is extracted from the gradient feature map according to a given threshold;

[0103] f) constructing an energy function to optimize the extraction of the sky boundary line. When the energy function is maximum, the sky boundary line is most accurate;

[0104] g) By detecting whether the boundary line has a large jump, it is determined whether the extracted sky area contains a large number of buildings. If so, the sky area is clustered and divided into two categories, and the real sky area is obtained by the Mahalanobis distance with the non-sky area;

[0105] h) Scan each column of pixels and determine whether the pixel belongs to the sky through the Mahalanobis distance between each column of pixels and the real sky area, thereby improving the sky boundary line.

[0106] The energy function is defined as:

[0107] For the original color or grayscale image, the energy function is as follows:

[0108]

[0109] where ∑ s and∑ g are the covariance matrices of pixels described by RGB values ​​in the sky and ground regions, respectively.

[0110] They are all 3×3 matrices and are defined as follows:

[0111]

[0112] Where N s and N g are the number of pixels in the sky and ground regions, respectively, and μ s and μ g 3-by-1 column vector representing the average RGB values ​​in the sky and ground regions, respectively.

[0113] and is the eigenvalue corresponding to the above two matrices, and |·| represents the determinant, which measures the variance or volume of the pixel distribution in each region. Among them, γ represents the emphasis on the uniformity of the sky region, and γ=2 is determined by experiments. Maximizing the energy function can minimize the intra-class variance of the ground and sky distributions.

[0114] S3: extracting feature points from the preprocessed virtual scene and optoelectronic image respectively, and using the feature points as matching control points;

[0115] Specifically, the image feature extraction process is as follows: Figure 3 As shown: corner point detection and line segment detection are used to detect feature points on the real-time image, and the corresponding control point pairs on the virtual image are obtained according to the rough alignment matrix.

[0116] Furthermore, the specific strategy for corner feature extraction is:

[0117] For a two-dimensional image I(x, y), when a point is translated by (Δx, Δy) at (x, y), the autocorrelation function is used to represent the similarity before and after the window moves:

[0118]

[0119] Where: w(x, y) is the window centered at (x, y); w(u, v) is the weighting function within the window, generally a Gaussian function. I(u+Δx, v+Δy) is expanded by the first order Taylor and the second order higher terms are ignored:

[0120] I(u+Δx,v+Δy)≈I(u,v)+Ix(u,v)Δx+Iy(u,v)Δy (5)

[0121] Where: Ix(u, v) and Iy(u, v) are partial derivatives of the image. Substituting into the above formula, we get:

[0122]

[0123] in:

[0124]

[0125] Where: matrix M is the covariance matrix of image gradient; λ1 and λ2 are the eigenvalues ​​of matrix M.

[0126] The corner response function R(x, y) is defined by using the matrix trace Trace(M) and determinant det(M), and its expression is:

[0127] R(x, y) = det(M) - k(Trace(M)) 2 (8)

[0128] Where: det(M) = λ1λ2; Trace(M) = λ1+λ2; k is an empirical constant, k∈[0.04, 0.06]. Given a threshold T, when R(x, y)>T and reaches the maximum in the local window, the pixel (x, y) is a feature point.

[0129] Furthermore, the specific strategy for image line segment extraction is:

[0130] If some points on the straight line in the image space deviate from the actual position, a rough model can be used: if S points are close enough to the ideal straight line, then these points belong to the straight line, and the points in the remaining subset can be considered as noise points. In this model, the failure rate is:

[0131]

[0132] Among them, αS / M refers to the proportion of points in the subset that belong to the ideal straight line in the entire subset. The performance of the transformation is determined by the expected value of m*αS / M and s, as follows:

[0133] (1) If the number of points m contained in the subset is too small, the failure rate of the transformation is close to l, that is:

[0134]

[0135] (2) If m is large enough, m*αS / M ≥ 1, the transformation failure rate is low, and its failure rate is:

[0136]

[0137] At this time, when m*αS / M≥2, the larger m is, the lower the failure rate is.

[0138] After obtaining the detected line segment, the line segment control points can be obtained through line segment sampling.

[0139] S4: performing precise registration on the image after the control point pairs are extracted;

[0140] Specifically, the image registration process is as follows: Figure 4 As shown in the figure: a feature space is constructed to find the common features of multi-source images, the image edges are extracted through edge detection operators, and the similarity measure of the two feature space images is calculated for further adjustment and calibration to obtain high-precision control point pairs; the two parameters of main peak sharpness and main-to-sub-peak ratio in correlation surface analysis are used to perform correlation surface analysis to remove mismatched points as much as possible.

[0141] The edge detection operator template is as follows:

[0142]

[0143] G in the formula X Represents the template in the horizontal direction, G Y Represents the template in the vertical direction.

[0144] Furthermore, the specific strategy of correlation analysis is:

[0145] (1) Primary and secondary peak ratio

[0146] V max Indicates the value corresponding to the highest correlation peak on the correlation surface, V sub It represents the value corresponding to the next peak on the correlation surface. The submaxratio is defined as the ratio of the highest peak to the next peak:

[0147]

[0148] When the primary and secondary peaks are relatively small, it means that there are at least two areas in the reference image that are similar to the image to be registered, and the matching result is unreliable, which may cause mismatch. The larger the value, the more unique the image to be registered is and the less likely it is to mismatch.

[0149] (2) Peak sharpness

[0150] With the maximum point of the highest peak area of ​​the correlation surface as the center, make a small circle with a radius of R1 and a large circle with a radius of R2. circle Represents the mean value of pixels in the small circle, expressed as V loop Represents the pixel mean of the area between the small circle and the large circle, and the value of sharpofmax is V circle and V loop According to the peak distribution, V loop Generally, V circle The formula for the highest peak sharpness can be expressed as follows:

[0151]

[0152] Both the primary-secondary peak ratio and the highest peak sharpness can be used as ways to remove mismatched points. The two have different focuses. The former focuses on the accuracy of the position, while the latter focuses on the precision of the matching position.

[0153] S5: Based on the high-precision control point pairs, a local weighted mean correction model is used to perform more accurate alignment between the real-time image and the virtual image.

[0154] Given N control point pairs of two images of the same scene:

[0155]

[0156] Among them, (x i ,y i ) and (X i ,Y i ) represent the control points in the i-th reference image and the real-time image respectively, and decompose the transformation model into component models f in the x- and y-directions x (x,y) and f y (x,y), then the relationship between the control points in the real-time graph and the reference graph can be expressed as follows:

[0157]

[0158] Once the transformation model is known, given any point (x, y) in the reference image, the corresponding point position (X, Y) in the real-time image can be calculated by the above formula.

[0159] Further expressed as a 3D point set:

[0160]

[0161] Among them, Xi and Yi can be regarded as two surfaces f x (x,y) and f y (x,y) in (x i ,y i), so the solution of the transformation model can be regarded as the solution of the 3D point set (x i ,y i ,X i ) and (x i ,y i ,Y i ) respectively perform surface fitting to solve the component model f x (x,y) and f y (x,y).

[0162] Given that (x i ,y i ,X i ) and (x i ,y i ,Y i ) has a similar expression, so the problem can eventually be transformed into a single-value surface fitting problem:

[0163] {(x i ,y i ,F i )|i∈[1,N]} (18)

[0164] Among them, F i =f(x i ,y i ), which can represent the transformation model in the x direction or the transformation model in the y direction. The difference lies in that the transformation coefficients in different directions are different.

[0165] The local correction model fits the transformation model of the surrounding control points through a quadratic polynomial at each control point. Therefore, the mapping position of each pixel in the real-time image is affected by the weighted average of all quadratic polynomials. The weight coefficient of the quadratic polynomial can be defined as follows:

[0166]

[0167] Among them, W i (x, y) is the inverse distance weight, n represents the neighborhood size of the control point, that is, the number of control point pairs used to calculate a quadratic polynomial model. Generally, the larger n is, the smoother the fitted surface is, but too large n will reduce the robustness of the local correction model to distortion; R n Indicates the (n-1)th nearest neighbor control point to (x i ,y i ). The f(x,y) corresponding to the local correction model is continuous and smooth. The definition of the local correction model is as follows:

[0168]

[0169] Among them, pi (x,y) represents a quadratic polynomial transformation through the n-neighborhood of each control point.

[0170] S6: A fusion algorithm for heterogeneous images based on multi-scale image block structured decomposition is used. This method performs multi-scale decomposition and reconstruction of infrared images and virtual vision images based on a multi-layer pyramid structure. Image block structured decomposition is used to extract image block features in each scale layer obtained by decomposition, and specific strategies are used to fuse infrared and virtual vision image features at different layers.

[0171] In the multi-scale decomposition structure, the average brightness information obtained by structured decomposition is used as the object of downsampling each time, and the downsampled average brightness information is used as the initial input image of the next layer.

[0172] In the multi-scale reconstruction structure, the image of the highest layer (smallest scale) is fused first. The fusion method is the same as the single-scale fusion method, and the low-frequency and high-frequency components of the highest layer are fused separately. Then the high-frequency components extracted from each layer are fused. After that, the fusion result of the highest layer is upsampled as the low-frequency component of the next highest layer, added to the high-frequency component of the next highest layer as the fusion result of the next highest layer, and then upsampled as the low-frequency component of the next lower layer... Iterate until the lowest layer is fused.

Claims

1. A method for matching an optoelectronic image and a three-dimensional virtual scene based on a helicopter, characterized in that: include: Keep the viewport when rendering the scene consistent with the viewport of the photoelectric sensor; De-noise the helicopter's electro-optical images and 3D virtual scenes collected from two viewports and remove the sky area; Extract feature points from the photoelectric image and the three-dimensional virtual scene after removing the sky area, and use the feature points as control points to be matched; Find the common features of the two images from the control points to be matched in the photoelectric image and the control points to be matched in the three-dimensional virtual scene, and obtain high-precision control point pairs; Align the photoelectric image with the 3D virtual scene based on high-precision control point pairs; Based on the aligned images and views, optoelectronic and virtual vision images are fused.

2. The method according to claim 1, characterized in that Find the common features of the two images from the control points to be matched in the photoelectric image and the control points to be matched in the three-dimensional virtual scene, including: The image edges of the photoelectric image and the virtual scene are extracted by edge detection operators to obtain two feature space images; Calculate the similarity measure of the two feature space images for further adjustment and calibration to match high-precision control point pairs; Calculate the main peak sharpness and the main-to-minor peak ratio of the two feature space images; Correlation surface analysis is performed on the two parameters to remove mismatched points.

3. The method according to claim 2, characterized in that Alignment between images is achieved using a local weighted mean correction model; The local correction model f(x,y) is defined as follows: Among them, W i (x,y) is the inverse distance weight, p i (x, y) represents a quadratic polynomial transformation through the neighborhood of each control point, and N is the number of control point pairs.

4. The method according to claim 3, characterized in that Make the viewport when rendering the scene consistent with the viewport of the photoelectric sensor, including: According to the installation positions of the photoelectric sensor and the inertial navigation satellite antenna of the carrier aircraft, the viewport positions of the two are calculated; and the viewport during scene rendering is adjusted to be consistent with the viewport of the photoelectric sensor.

5. The method according to claim 4, characterized in that The specific steps for extracting the sky area are: Firstly, the Sobel operator is used to extract the image gradient features, and the boundary line of the sky is extracted from the gradient feature map according to a given threshold. An energy function is constructed to optimize the extraction of the sky boundary line. When the energy function is maximum, the sky boundary line is most accurate. By detecting whether the boundary line has a large jump, we can determine whether the extracted sky area contains a large number of buildings. If so, we perform cluster analysis on the sky area and divide it into two categories. The real sky area is obtained by the Mahalanobis distance with the non-sky area. Scan each column of pixels and determine whether the pixel belongs to the sky through the Mahalanobis distance between each column of pixels and the real sky area, thereby improving the sky boundary line.

6. The method according to claim 5, characterized in that Feature points are extracted from the photoelectric image and the three-dimensional virtual scene after the sky area is removed, and the feature points are used as control points to be matched, including: Corner point detection and line segment detection are used to extract feature points from the preprocessed virtual scene and optoelectronic image respectively.

7. The method according to claim 6, characterized in that The specific steps of corner detection are: For a two-dimensional image I(x, y), when a point at (x, y) is translated by (Δx, Δy), the autocorrelation function c(x, y; Δx, Δy) is used to represent the similarity before and after the window moves: c(x,y;Δx,Δy)=∑ (u，v)∈W(x，y) w(u,v)(I(u,v)-I(u+Δx,v+Δy)) 2 ; Where: w(x, y) is the window centered at (x, y); w(u, v) is the weighted function of the sliding window centered at (u, v), which is a Gaussian function; Perform a first-order Taylor expansion on I(u+Δx, v+Δy) and ignore the second-order higher-order terms: I(u+Δx,v+Δy)=I(u,v)+Ix(u,v)Δx+Iy(u,v)Δy; Where: Ix(u, y) and Iy(u, v) are partial derivatives of the image; Substituting into the above formula, we get: in: Where: Matrix M is the covariance matrix of the image gradient; λ1 and λ2 are the eigenvalues ​​of the matrix M; The corner response function R(x, y) is defined by using the matrix trace Trace(M) and determinant det(M), and its expression is: R(x,y)=det(M)-k(Trace(M)) 2 ; Where: det(M) = λ1λ2; Trace(M) = λ1+λ2; k is an empirical constant, k∈[0.04, 0.06]; given a threshold T, when R(x, y)>T and reaches the maximum in the local window, the pixel point (x, y) is a feature point.

8. The method according to claim 6, characterized in that The specific steps of line segment detection are: Assume that S points are close enough to the ideal straight line, then these points belong to the straight line, and the points in the remaining subset are considered to be noise points; in this model, the failure rate is: Among them, αS / M refers to the proportion of points in the subset that belong to the ideal straight line in the entire subset; the performance of the transformation is determined by the expected value of m*αS / M and s, as follows: (1) If the number of points m contained in the subset is too small, the failure rate of the transformation is close to l, that is: (2) If m is large enough, m*αS / M ≥ 1, the transformation failure rate is low, and its failure rate is: At this time, when m*αS / M≥2, the larger m is, the lower the failure rate is; After obtaining the detected line segment, the line segment control points can be obtained through line segment sampling.

9. The method according to claim 5, characterized in that Energy function J n is defined as: where ∑ s and∑ g are the covariance matrices of pixels described by RGB values ​​in the sky and ground regions, respectively; they are both 3×3 matrices defined as follows: Where N s and N g are the number of pixels in the sky and ground regions, respectively, and μ s and μ g 3×1 column vectors representing the average RGB values ​​in the sky and ground regions, respectively; and are the eigenvalues ​​corresponding to the above two matrices, |·| represents the determinant, which measures the variance or volume of the pixel distribution in each area; where γ represents the emphasis on the uniformity of the sky area.

Images