An Image Restoration Method under Complex Optical Imaging Conditions Based on Measured Data

Abstract

The present invention discloses an image restoration method under complex optical imaging conditions based on measured data, belonging to the field of image restoration. The present invention removes the noise brought by the imaging system based on the morphological filtering method, performs edge detection on the image spectrum of the blurred image under a point light source to remove the blur brought by the imaging system, determines the motion blur parameters based on the Radon transform method, uses the method of flow field data simulation to remove the pneumatic thermal radiation noise, combines the ray tracing method and the background schlieren technology to simulate and model to remove the image blur and image distortion caused by the pneumatic optical effect, and calculates the image quality evaluation parameters to comprehensively evaluate the image restoration effect. The present invention uses two methods to model for the acquisition of flow field data, and can completely remove noise, remove blur and remove distortion and restore the degraded image under complex imaging conditions, which is convenient for subsequent target detection and positioning of the image.

Description

Technical Field

[0001] The present invention belongs to the field of image restoration, and more specifically, relates to an image restoration algorithm method under complex optical imaging conditions based on measured data. Background Art

[0002] With the increasing dependence of precision-guided weapons on the accuracy and timeliness of target information acquisition, optical imaging detection and tracking technology has become an inevitable trend in the development of precision guidance. Compared with traditional inertial guidance, command guidance, and radio terminal guidance technologies, the terminal guidance technology using a vision system has higher accuracy and stronger anti-interference ability. It uses imaging information for target detection, realizes the detection, recognition, and tracking of the target, and inputs the line-of-sight angle position and angular rate of the target to the aircraft control system, enabling the control system to control the aircraft to strike the target.

[0003] When an imaging-guided aircraft flies at high speed in the atmosphere, the flow field near the window becomes extremely complex. Its high Mach number and Reynolds number will cause a shock wave structure to be generated at the head of the aircraft, and turbulent boundary layers, shear layers, etc. will also be generated in the airflow near the window, ultimately resulting in thermal radiation interference and image transmission interference of the optical imaging detection system. This interference is called the aero-optical effect. In the engineering application of weaponry in the dense atmosphere, in order to pursue the imaging quality under the aero-optical effect, its terminal flight speed will be severely restricted, resulting in a serious decline in penetration and strike effects.

[0004] Regarding the manifestations of the aero-optical effect on image information, it is mainly divided into jitter blur, blur, phase distortion of the optical transmission path, aero-thermal radiation noise, and noise, etc. Blur is caused by the relative movement of the light beam on the imaging focal plane during a relatively long exposure time interval, ultimately resulting in attenuation of the detection object intensity, reduction of the detection distance, increase of measurement error, and even detection failure, etc.; vibration or jitter blur stems from the impact effect between the airflow and the aircraft platform, and can be characterized by amplitude, frequency, and direction. It will cause the centroid position of the detection target to shift, ultimately resulting in inaccuracy of measurement or positioning; the phase distortion of the optical transmission path is caused by the optical path difference, that is, the plane wave is advanced or lagged by different wavelength numbers at different spatial positions, ultimately resulting in the deviation of the line of sight or pointing; noise stems from the aero-thermal radiation effect, ultimately resulting in a decrease in the image signal-to-noise ratio and loss of image information. According to the reasons for its degradation, the restoration of the aero-optical effect is divided into 5 parts, which are, in order, removing the noise of the imaging system, removing the blur of the imaging system, removing the motion blur, removing the aero-thermal radiation noise, and removing the distortion and blur of the aero-optical effect.

[0005] At present, the research on restoration algorithms for degraded images under aero-optical effects at home and abroad is mainly divided into blind restoration algorithms, non-blind restoration algorithms and restoration algorithms based on deep learning. At present, the main shortcomings of aero-optical effect image restoration methods include:

[0006] (1) Since aero-optical effects can cause image degradation in many aspects, and noise, blur, and distortion can affect each other, a single algorithm cannot produce a good restored image.

[0007] (2) The aero-optical effect is time-varying and has a high frequency of change;

[0008] (3) It is difficult to obtain measured aero-optical environment data, and aero-optical effects can only be modeled and removed based on accurate data;

[0009] (4) The restoration method based on deep learning requires a large amount of data. Since the cost of obtaining measured data is high and the data is not open, it is difficult to build a complete database, making it difficult to apply deep learning methods to restore degraded images under aero-optical effects. Summary of the invention

[0010] In view of the problems that the restoration of existing aero-optical effect degraded images is easily affected by noise, it is difficult to obtain accurate flow field measured data, and it is impossible to accurately judge the image pixel displacement offset, the present invention proposes an image restoration method under complex optical imaging conditions based on measured data to solve the problems of noise, blur and distortion of degraded images under aero-optical effect.

[0011] The basic principle of the image restoration algorithm adopted in the present invention is: according to the mechanism of image degradation, the image degradation process is decomposed into five parts, namely, removing the noise caused by the imaging system, removing the blur caused by the imaging system, removing motion blur, removing aerodynamic thermal radiation noise, and removing aerodynamic optical effect distortion and blur. The degradation principle of each degraded part is studied, and partial algorithms are used in a targeted manner to reconstruct a clear image.

[0012] The method proposed in the present invention for the problem of degraded image restoration under complex optical imaging conditions can remove the noise, blur and distortion of the image under complex imaging conditions, apply certain prior information to model the degradation mechanism, improve the signal-to-noise ratio and clarity of the image, and correct the pixel displacement distortion in the image imaging process, ultimately obtaining a clearer restored target image.

[0013] To achieve the above object, the technical solution adopted by the present invention is:

[0014] The present invention first provides a method for restoring a degraded image under complex optical imaging conditions based on measured data, which comprises the following steps:

[0015] Step 1: Calculate the peak signal-to-noise ratio (PSNR), the average gradient of sharpness, and the coordinates of the maximum point of the image pixel values of the original degraded image \(g_0(x,y)\);

[0016] Step 2: Perform morphological filtering on the degraded image to remove the noise introduced by the imaging system during the imaging process and the noise generated by part of the aerodynamic heat radiation effect, obtaining a denoised image \(g_1(x,y)\); Calculate the PSNR of the image \(g_1(x,y)\), and determine whether the improvement in the PSNR value obtained in this step compared to the PSNR value obtained in Step 1 is greater than 5 dB. If so, proceed to the next step; otherwise, repeat Step 2;

[0017] Step 3: Simulate multiple point light source images passing through the imaging system, and average the multiple point light source images to reduce the influence of noise; After performing Fourier transform on the averaged image, perform edge detection, and after processing, obtain the estimated values of the blur radius \(r\) and the standard deviation \(\sigma\) of the imaging system. Through Gaussian function modeling, obtain the blur point spread function \(H_1(u,v)\) of the imaging system:

[0018]

[0019] where \(\sigma\) is the standard deviation of the normal distribution, and \((u,v)\) represents the frequency coordinates of the image; Use the blur point spread function \(H_1(u,v)\) of the imaging system to perform least-squares filtering on the image \(g_1(x,y)\) obtained in Step 2, obtaining an image \(G_2(u,v)\) that removes the blur of the imaging system

[0020]

[0021] where \(G_1(u,v)\) is the Fourier spectrum of the image \(g_1(x,y)\), \(H_1\) * (u,v) represents the complex conjugate of \(H_1(u,v)\), \(\eta_1\) represents the noise level of the image. Perform inverse Fourier transform on \(G_2(u,v)\), and then the restored image \(g_2(x,y)\) can be obtained; According to the flow field data near the detection window, calculate the radiation intensity of the aerodynamic heat radiation effect, and subtract this matrix from the image \(g_2(x,y)\) to obtain a restored image \(g_3(x,y)\) that removes the aerodynamic heat radiation;

[0022] Step 4: Perform image processing on the spectrum \(G_3(u,v)\) of the image \(g_3(x,y)\) obtained in Step 3 and then perform Radon transform; The angle corresponding to the maximum value of the Radon transform curve is the motion blur angle of the degraded image. Calculate the difference between the two sides of the main peak in the blur angle direction, which is the motion blur scale of the degraded image, and then obtain a vector line segment containing the angle and the scale as the motion blur point spread function \(H_2(u,v)\); Use the motion blur point spread function \(H_2(u,v)\) to perform least-squares filtering on the image \(g_3(x,y)\) that removes the blur of the imaging system obtained in Step 3, generating an image \(G_4(u,v)\) that removes the motion blur;

[0023]

[0024] Where \(G_3(u, v)\) is the Fourier transform of the image \(g_3(x, y)\), and \(\eta_2\) represents the noise level of the image; perform the inverse Fourier transform on \(G_4(u, v)\) to obtain the restored image \(g_4(x, y)\); calculate the average gradient of sharpness of the image \(g_4(x, y)\), and check if there is an improvement of more than 5% compared to the value in Step 1. If so, proceed to the next step; otherwise, return to Step 3 and appropriately increase the values of the noise levels \(\eta_1\) and \(\eta_2\) of the image.

[0025] Step 5: Determine whether it is possible to obtain the flow field data under complex imaging conditions in the wind tunnel test or the actual measurement scenario. For the case where the flow field data can be obtained, continue with Step 6; for the case where the flow field data cannot be obtained, directly execute Step 7.

[0026] Step 6: Simulate the incident light rays to determine the light ray incident angle. Using a fixed step size, calculate the optical path length OPL of the light rays passing through the entire flow field, and calculate the optical path difference after averaging it. Calculate the flow field point spread function \(H_3(u, v)\) based on the optical path difference; use \(H_3(u, v)\) to perform least-squares filtering on the image \(g_4(x, y)\) obtained in Step 4.

[0027]

[0028] Where \(\eta_3\) represents the noise level of the image; obtain the final clear restored image \(G_5(u, v)\) under complex optical imaging conditions, perform the inverse Fourier transform on \(G_5(u, v)\) to obtain the restored image \(g_5(x, y)\); calculate the average gradient of sharpness of the image \(g_5(x, y)\) and the coordinates of the point with the maximum pixel value of the image. Check if the average gradient of sharpness has an improvement of more than 5% compared to the value in Step 4 and if the coordinates of the point with the maximum pixel value are different from those in Step 1. If both conditions are met, proceed to the next step; otherwise, repeat Step 6 and appropriately increase the value of the noise level \(\eta_3\) of the image.

[0029] Step 7: Obtain the background random dot matrix image, and take a set of correlation images by shooting under two conditions: without a flow field and with a flow field; perform cross-correlation operations on multiple sets of correlation images, obtain the displacement change vectors of multiple sets of cross-correlation interrogation regions, and then perform time-averaging processing; perform SouthWell linear integral operation on the processed displacement vectors to obtain the wavefront optical path difference OPD, and finally calculate the flow field point spread function \(H_4(u, v)\) based on OPD; use the PSF to perform least-squares filtering on the image \(g_4(x, y)\) obtained in Step 4.

[0030]

[0031] where η4 represents the noise level of the image; the final clear restored image G6(u, v) under complex optical imaging conditions is obtained, and the inverse Fourier transform is performed on G6(u, v) to obtain the restored image g6(x, y); calculate the average gradient of the sharpness of the image g6(x, y) and the coordinates of the maximum point of the image pixel value. Whether the average gradient of the sharpness has an improvement of more than 5% compared to the value in step 4 and whether the coordinates of the maximum point of the image pixel value are different from the value in step 1. If all are satisfied, proceed to the next step; otherwise, re - execute step 7 and appropriately increase the value of the noise level η4 of the image.

[0032] Step 8: Calculate the comprehensive image quality evaluation index. When the comprehensive image quality parameter is greater than a certain threshold, it is considered that the image restoration effect is good and the image restoration is completed; otherwise, return to step 2, appropriately increase the noise level of the image in the formula, and recalculate.

[0033] Compared with the prior art, the advantages of the present invention are as follows:

[0034] (1) The present invention uses a method of obtaining the degraded image of a point light source to determine the blurring parameters of the imaging system, simplifies the process of modeling the blurring point spread function, and improves the accuracy of the point spread function model at the same time;

[0035] (2) The present invention uses two methods, the ray tracing method and the background schlieren technique, for image restoration, overcomes the problem of difficult acquisition of flow field parameters in the prior art, and improves the practicability of the technology;

[0036] (3) The present invention uses a method of decomposing the image degradation process, decomposes the image degradation into noise, imaging system blurring, motion blurring, and blurring and distortion based on the aero - optical effect, overcomes the defect of using a single model for clear image restoration in the prior art, and obtains a better image restoration effect;

[0037] (4) The present invention guides the selection of parameters in the image restoration process through the negative feedback of the image quality evaluation index, makes the image processing and image quality detection form a closed - loop, and can directly apply the restoration result to subsequent target detection, greatly improving the practicability of the technology;

[0038] (5) The present invention uses a variety of image quality evaluation indexes to evaluate the noise degree, blurring degree, and distortion degree of the image, overcomes the problem of single - dimensional evaluation of image quality in the prior art. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 is the flowchart of the image restoration method under complex optical imaging conditions provided by the embodiment of the present invention.

[0040] Figure 2(a) is the denoising result diagram of the uniformly aero - optical effect degraded image provided by the embodiment of the present invention.

[0041] Figure 2(b) is the denoising result diagram of the random pneumatic shooting effect degradation image provided by the implementation of the present invention.

[0042] Figure 3(a) is the result diagram of removing the blur of the imaging system for the uniform pneumatic effect degradation image provided by the implementation of the present invention.

[0043] Figure 3(b) is the result diagram of removing the blur of the imaging system for the random pneumatic effect degradation image provided by the implementation of the present invention.

[0044] Figure 4 It is the result diagram of the edge detection of the image spectrum provided by the implementation of the present invention.

[0045] Figure 5(a) is the result diagram of removing motion blur for the uniform pneumatic effect degradation image provided by the implementation of the present invention.

[0046] Figure 5(b) is the result diagram of removing motion blur for the random pneumatic effect degradation image provided by the implementation of the present invention.

[0047] Figure 6(a) is the result diagram of removing the blur and distortion of the pneumatic optical effect for the uniform pneumatic effect degradation image based on the ray tracing method provided by the implementation of the present invention.

[0048] Figure 6(b) is the result diagram of removing the blur and distortion of the pneumatic optical effect for the random pneumatic effect degradation image based on the ray tracing method provided by the implementation of the present invention.

[0049] Figure 7(a) is the random background dot matrix diagram provided by the implementation of the present invention under the condition of no flow field.

[0050] Figure 7(b) is the random background dot matrix diagram provided by the implementation of the present invention under the condition of having a flow field.

[0051] Figure 8(a) is the result diagram of removing the blur and distortion of the pneumatic optical effect for the uniform pneumatic effect degradation image based on the background schlieren technique provided by the implementation of the present invention.

[0052] Figure 8(b) is the result diagram of removing the blur and distortion of the pneumatic optical effect for the random pneumatic effect degradation image based on the background schlieren technique provided by the implementation of the present invention. Detailed implementation manners

[0053] In order to make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be described in detail below with reference to specific embodiments. Specific embodiments are described below to simplify the present invention. However, it should be recognized that the present invention is not limited to the described embodiments, and various modifications of the present invention are possible without departing from the basic principles, and these equivalent forms also fall within the scope defined by the appended claims of this application.

[0054] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0055] As Figure 1 shown, an image restoration algorithm under complex optical imaging conditions based on measured data provided by the present invention mainly includes the following steps:

[0056] Step 1: Calculate the peak signal-to-noise ratio (PSNR), the average gradient of sharpness, and the coordinates of the maximum point of the image pixel values of the original degraded image g0(x, y);

[0057] Step 2: Perform morphological filtering on the degraded image to remove the noise introduced by the imaging system during the imaging process and the noise generated by part of the pneumatic thermal radiation effect, and obtain the denoised image g1(x, y); calculate the peak signal-to-noise ratio (PSNR) of the image g1(x, y), and determine whether the improvement amount of this value compared with the PSNR value obtained in Step 1 is greater than 5 dB. If yes, proceed to the next step; otherwise, repeat Step 2.

[0058] Step 3: Simulate multiple point light source images passing through the imaging system, and average the multiple point light source images to reduce the influence of noise; perform Fourier transform on the averaged image and then perform edge detection. After processing, obtain the estimated values of the blur radius r and the standard deviation σ of the imaging system. Through Gaussian function modeling, obtain the blur point spread function H1(u, v) of the imaging system:

[0059]

[0060] where σ is the standard deviation of the normal distribution, and (u, v) represents the frequency coordinates of the image; use the blur point spread function of the imaging system to perform least square filtering on the image g1(x, y) obtained in Step 2 to obtain the image G2(u, v) after removing the blur of the imaging system

[0061]

[0062] where G1(u, v) is the Fourier spectrum of the image g1(x, y), H1 * (u, v) represents the complex conjugate of H1(u, v), η1 represents the noise level of the image. Perform inverse Fourier transform on G2(u, v), and then obtain the restored image g2(x, y); calculate the radiation intensity of the pneumatic thermal radiation effect according to the flow field data near the detection window, and subtract this matrix from the image g2(x, y) to obtain the restored image g3(x, y) after removing the pneumatic thermal radiation.

[0063] Step 4: Perform image processing on the spectrum G3(u, v) of the image g3(x, y) obtained in Step 3 and then perform the Radon transform; the angle corresponding to the maximum value of the Radon transform curve is the motion blur angle of the degraded image. Calculate the difference on both sides of the main peak in the blur angle direction, which is the motion blur scale of the degraded image. Then, obtain a vector line segment containing the angle and scale as the motion blur point spread function H2(u, v); use the motion blur point spread function H2(u, v) to perform least-squares filtering on the image g3(x, y) that removes the blur of the imaging system obtained in Step 3 to generate the image G4(u, v) that removes the motion blur:

[0064]

[0065] where G3(u, v) is the Fourier transform of the image g3(x, y), and η2 represents the noise level of the image; perform the inverse Fourier transform on G4(u, v) to obtain the restored image g4(x, y); calculate the average gradient of the sharpness of the image g4(x, y). Compare whether there is an increase of more than 5% compared to the value in Step 1. If so, proceed to the next step; otherwise, return to Step 3 and appropriately increase the values of the noise levels η1 and η2 of the image at the same time;

[0066] Step 5: Determine whether it is possible to obtain the flow field data under complex imaging conditions in the wind tunnel test or the actual measurement scenario. For the case where the flow field data can be obtained, continue to execute Step 6. For the case where the flow field data cannot be obtained, directly execute Step 7;

[0067] Step 6: Simulate the incident light rays to determine the light ray incident angle. Adopt a fixed-step method to calculate the optical path length OPL of the light rays passing through the entire flow field, and calculate the optical path difference after averaging it Calculate the flow field point spread function H3(u, v) based on the optical path difference; use H3(u, v) to perform least-squares filtering on the image g4(x, y) obtained in Step 4

[0068]

[0069] where η3 represents the noise level of the image; obtain the final clearly restored image G5(u, v) under complex optical imaging conditions. Perform the inverse Fourier transform on G5(u, v) to obtain the restored image g5(x, y); calculate the average gradient of the sharpness of the image g5(x, y) and the coordinates of the point with the maximum image pixel value. Compare whether the average gradient of the sharpness has an increase of more than 5% compared to the value in Step 4 and whether the coordinates of the point with the maximum image pixel value are different from the value in Step 1. If all are satisfied, proceed to the next step; otherwise, re-execute Step 6 and appropriately increase the value of the noise level η3 of the image at the same time;

[0070] Step 7: Obtain the background random dot matrix image, and take a set of correlation images under two conditions: without a flow field and with a flow field; perform cross-correlation operations on multiple sets of correlation images, obtain the displacement change vectors of multiple sets of cross-correlation interrogation regions, and then perform time-averaging processing; perform SouthWell linear integral operation on the processed displacement vectors to obtain the wavefront optical path difference OPD, and finally calculate the flow field point spread function H4(u, v) according to OPD; use the PSF to perform least-squares filtering on the image g4(x, y) obtained in Step 4

[0071]

[0072] where η4 represents the noise level of the image; obtain the final clear restored image G6(u, v) under complex optical imaging conditions, perform inverse Fourier transform on G6(u, v), and obtain the restored image g6(x, y); calculate the average gradient of sharpness of the image g6(x, y) and the coordinates of the maximum point of the image pixel value, and check whether the average gradient of sharpness has an increase of more than 5% compared to the value in Step 4 and whether the coordinates of the maximum point of the image pixel value are different from the value in Step 1. If all conditions are met, proceed to the next step; otherwise, re-execute Step 7 and appropriately increase the value of the noise level η4 of the image

[0073] Step 8: Calculate the comprehensive image quality evaluation index. When the comprehensive image quality parameter is greater than a certain threshold, it is considered that the image restoration effect is good, and the image restoration is completed; otherwise, return to Step 2, appropriately increase the noise level of the image in the formula, and recalculate

[0074] The following introduces the specific implementation methods

[0075] In a specific embodiment of the present invention, Step 1 is: calculating the image quality evaluation parameters of the original degraded image g0(x, y), and the main steps are as follows

[0076] 1-1: Calculate the peak signal-to-noise ratio PSNR of the original degraded image g0(x, y) and the reference image g0'(x, y)

[0077]

[0078]

[0079] where M and N represent the length and width of the image, g(x, y) represents the pixel value of the point (x, y) in the distorted image, g′(x, y) represents the pixel value of the reference image, and L is the maximum value of the pixel gray level, generally 255

[0080] 1-2: Calculate the average gradient of sharpness of the original degraded image g0(x, y)

[0081]

[0082] wherein k x (x, y) and k y (x, y) represents the image gradient calculated by the Sobel operator at the pixel point (x, y) in the x and y directions, and l(x, y) represents the image gradient calculated by the Laplacian operator at the pixel point (x, y).

[0083] 1 - 3: Calculate the coordinates (x max1 , y max1 ) of the point with the maximum pixel value of the original degraded image g0(x, y).

[0084] In this embodiment, step 2 is: denoise the original degraded image g0(x, y), and the main steps are as follows:

[0085] 2 - 1: Perform an erosion operation on the original degraded image g0(x, y) to smooth the noise points and make the image noise into a block of the same level;

[0086] 2 - 2: Perform morphological reconstruction on the eroded image to extract the connected regions of the features in the original image. Use the eroded image and the original image as the marker and the mask respectively to constrain the transformation process, implement the opening operation of the image, and obtain the denoised image g1(x, y). Process the uniformly aerodynamic degraded image and the randomly aerodynamic degraded image respectively, and the results are shown in Figures 2(a) and 2(b) respectively. This attached figure is the restored image after denoising the degraded image.

[0087] 2 - 3: Calculate the peak signal - to - noise ratio PSNR of the image g1(x, y), and determine whether the improvement amount of this value compared to the PSNR value obtained in step 1 - 1 is greater than 5 dB. If so, proceed to the next step; otherwise, repeat step 2.

[0088] In this embodiment, step 3 is: approximately regard the point - spread function that causes blurring in the imaging system as a Gaussian function, which is defined as:

[0089]

[0090] where σ is the standard deviation of the normal distribution, and (u, v) are the frequency coordinates of the image. The standard deviation σ determines the smoothness of the point - spread function. The larger σ is, the greater the influence of distant pixels on the central pixel, and the smoother the image. Since the point - spread function H1(u, v) of the imaging system does not vary according to different imaging scenarios, the parameters of the point - spread function H1(u, v) can be determined through experiments. Step 3 specifically includes the following steps:

[0091] 3-1: Obtain multiple degraded images of point light sources passing through the imaging system. These images contain imaging system noise and blur. First, perform weighted averaging on multiple images to suppress noise, and then perform Fourier transform on the noise-suppressed image to obtain the spectrum of the degraded image.

[0092] 3-2: Perform edge detection on the spectrum of the degraded image. Apply the Sobel operator for approximation, find the points with the maximum image gradient to locate the edges. Set the edge points to 1 and the non-edge points to 0. After determining the edge points, calculate the radius of the non-zero region to obtain the parameter r.

[0093] 3-3: The standard deviation σ generally takes values between 0 and 8. Use the parameter r and σ to perform Wiener filtering on the image g1(x,y) obtained in step 2-2, calculate the PSNR of the filtered image, and take the standard deviation σ with the optimal PSNR value as the estimated value of the Gaussian kernel standard deviation σ to obtain the point spread function H1(u,v) of the imaging system.

[0094] 3-4: Perform constrained least squares filtering on the image g1(x,y) obtained in step 2-2 to obtain the image G2(u,v) that removes the blur of the imaging system.

[0095]

[0096] where G1(u,v) is the Fourier spectrum of the image g1(x,y), H1 * (u,v) represents the complex conjugate of H1(u,v), and η1 represents the noise level of the image, which is selected according to empirical values. Perform inverse Fourier transform on G2(u,v) to obtain the restored image g2(x,y).

[0097] 3-5: According to the radiation data near the detection window, calculate the radiation intensity of the aerodynamic heat radiation effect. Calculate the radiation intensities of CO2 at 2.7μm and 4.3μm and H2O at 1.9μm, 2.7μm, and 6.3μm respectively. After interpolation, obtain a radiation matrix with the same size as the degraded image. Subtract this matrix from the image g2(x,y) obtained in step 3-4 to obtain the restored image g3(x,y) that removes the aerodynamic heat radiation. For uniformly aerodynamically degraded images, a restored image that completely filters out the aerodynamic heat radiation can be obtained; for randomly aerodynamically degraded images, since the image has been denoised in step 2, this processing also has a certain filtering effect on the random aerodynamic radiation. Therefore, no further processing is required, and the restored images are shown in Figures 3(a) and 3(b). This attached figure is the restored image after removing detector blur and aerodynamic heat radiation from the degraded image.

[0098] In this embodiment, step 4 is: Perform motion blur removal processing on the image. The main steps are as follows:

[0099] 4-1: Since the exposure time of the imaging surface CCD or the near-infrared imaging surface is extremely short, it can be approximately considered that there is a relatively uniform motion between the aircraft and the imaging detector during this exposure period. Therefore, the point spread function H2(u, v) of motion blur is a straight line. Perform Fourier transform on the image g3(x, y) to obtain the image spectrum G3(u, v), and perform logarithmic transform on the image spectrum. Since a central bright line appears around the spectrum, the image is edge-cropped. Manually move the point where u = v = 0 in the cropped spectrum image to the center position to obtain the centered spectrum image. Perform edge detection on the centered spectrum image and use the Canny operator for calculation. The detection result is as Figure 4 shown.

[0100] 4-2: The image on the imaging surface is blurred and degraded at a certain speed and angle. Its horizontal speed is a and its vertical speed is b. Then the point spread function H2(u, v) of motion blur can be expressed as

[0101]

[0102] where M and N are the image sizes, t is the image exposure time, (x, y) is the image coordinate point, u ranges from 0 to (N - 1), and v ranges from 0 to M - 1. The above formula can be expressed in discrete form

[0103]

[0104] According to the variation law of H2(u, v), it is obtained that the degraded image spectrum G3(u, v) satisfies that it shows black stripes (where black is the background, low gray level, and white is high gray level) at Finally, alternating bright and dark stripes appear. The angle α between the degraded image spectrum stripes and the horizontal is the inclination angle of the straight line and can be expressed by the formula

[0105]

[0106] After binarizing the edge detection result obtained in step 4-1 and performing Radon transform in the range of 0 - 180°, after obtaining the Radon curve, the angle corresponding to the maximum value is the motion blur angle.

[0107] 4-3: After obtaining the blur angle using the method in step 4-2, draw the Radon transform curve at this blur angle and calculate the difference between the two sides of the main peak. This difference is the blur scale of the image.

[0108] 4-4: Based on the obtained blur scale and blur angle in steps 4-2 and 4-3, draw the point spread function H2(u, v) of motion blur. Perform constrained least squares filtering on the image g3(x, y) obtained in step 3-5 to obtain the image G4(u, v).

[0109]

[0110] Where G3(u, v) is the Fourier transform of the image g3(x, y), and η2 represents the noise level of the image, which is selected according to empirical values. Perform inverse Fourier transform on G4(u, v) to obtain the restored image g4(x, y). The restoration results are shown in Figures 5(a) and 5(b). This attached figure is the restored image after removing motion blur from the degraded image.

[0111] 4-5: Calculate the average gradient of sharpness of the image g4(x, y) obtained in step 4-4. Check if there is an improvement of more than 5% compared to the value in step 1-2. If yes, proceed to the next step; otherwise, return to step 3-4, appropriately increase the values of the noise levels η1 and η2 of the image, and recalculate.

[0112] In this embodiment, step 5 is: Determine whether it is possible to obtain the flow field data under complex imaging conditions in the wind tunnel test or the actual measurement scenario. For the case where the flow field data can be obtained, continue to execute step 6; for the case where the flow field data cannot be obtained, directly execute step 7.

[0113] In this embodiment, step 6 is: For the case where the actual measurement or simulation flow field data has been obtained, use the transient refractive index of the flow field for analysis.

[0114] 6-1: Simulate the incident light. The observed target is a point light source at infinity, so the incident light can be approximated as a beam of parallel light, and the incident range is the effective window range of the optical imaging system. Determine the incident angle of the light relative to the flow field when the angle of attack of the aircraft is 0 degrees. The fixed step size is 0.001 m, and each movement of one step size is regarded as a refraction occurring. In the next step size, the air refractive index is taken as the refractive index of the midpoint of the flow field closest to this point. Let the i-th ray travel k step sizes, and each step size be set as Δl, and the refractive index in each step size be n. j , and the optical path obtained by the ray tracing method is

[0115]

[0116] 6-2: The chief ray is set as a ray directly facing the center of the optical axis of the imaging system. Let its optical path difference be The optical path difference of the i-th ray is The position at the pupil imaging is (x i , yi ) Then, the wave phase difference generated by the $i$-th ray passing through the flow field is

[0117]

[0118] where $\lambda$ represents the wavelength of light. The actual pupil function is equivalent to the product of the pupil function under ideal conditions and the complex amplitude function of the wave phase difference, that is

[0119] $A(x,y)=A_0(x,y)\times e$ jw(x,y)

[0120] where $A_0(x,y)$ is the pupil function under ideal conditions. Perform a two-dimensional Fourier transform on the pupil function to obtain the diffraction complex amplitude distribution

[0121]

[0122] Multiply the complex amplitude distribution by its conjugate function and then perform normalization processing to obtain the final point spread function $H_3(u,v)$

[0123] $H_3(u,v)=U(u,v)\times U$ * (u,v)

[0124] 6 - 3: Combine the $H_3(u,v)$ obtained in step 6 - 2 to perform constrained least - squares filtering on the image $g_4(x,y)$ obtained in step 4 - 4

[0125]

[0126] where $\eta_3$ represents the noise level of the image and is selected according to empirical values. Perform an inverse Fourier transform on $G_5(u,v)$ to obtain the restored image $g_5(x,y)$. The results are shown in Figures 6(a) and 6(b). This figure is the restored image after removing the blur and pixel distortion caused by the aero - optical effect from the degraded image.

[0127] 6 - 4: Calculate the average gradient of sharpness and the coordinates of the maximum point of the image pixel values of the image $g_5(x,y)$ obtained in step 6 - 3. Whether the average gradient of sharpness is increased by more than 5% compared to the value in step 4 - 5 and whether the coordinates of the maximum point of the image pixel values are different from the values in step 1 - 3. If all are met, proceed to the next step; otherwise, return to step 6 - 3, appropriately increase the value of the image noise level $\eta_3$, and recalculate.

[0128] In this embodiment, step 7 is as follows: For the case where measured or simulated flow field data cannot be obtained, the background schlieren technique can be used to model the point spread function. The specific steps are as follows:

[0129] 7-1: Obtain a set of random background dot matrices without a flow field and a set of random background dot matrices with a flow field. The results are shown in Figures 7(a) and 7(b). Perform enhanced contrast processing on the background dot matrix images passing through the flow field, perform cross-correlation operations on the above set of images, and use interrogation regions with windows of 32*32 and step sizes of 16*16, windows of 16*16 and step sizes of 8*8, and windows of 8*8 and step sizes of 4*4 respectively to perform cross-correlation operations based on FFT on the two images. For background dot matrix images with a size of 512*512 pixels, finally obtain 128 interrogation regions and obtain the pixel offset images in the 128 interrogation regions.

[0130] 7-2: Use the SouthWell algorithm to integrate the surrounding four points and take the average to obtain the optical path length OPL of each point, that is

[0131]

[0132] where σ represents the weight, the weight of non-existent points is taken as 0, the weights of the remaining points are taken as 1, (i,j) represents the pixel coordinates in the pixel offset image, H is the distance from the background to the center of the flow field, Δx is the offset distance of the pixels in the x direction of each interrogation region, Δy is the offset distance of the pixels in the y direction of each interrogation region, h and l are the height and width of the interrogation region respectively. Assuming that the initial value of the wavefront at the point (0,0) is 0, then the initial values of the wavefronts in each interrogation region are

[0133]

[0134]

[0135]

[0136]

[0137] After obtaining the initial values, perform iteration. The number of iterations is determined according to the relative difference between two adjacent steps. Due to the discontinuity of the difference, the wavefront data may be discontinuous at each value point. Therefore, smooth processing is performed on the obtained wavefront data.

[0138] 7-3: According to the wavefront data obtained in step 7-2, similar to step 6-2, calculate the point spread function H4(u,v) of the aerodynamic optical effect.

[0139] 7-4: Combine the H4(u,v) obtained in step 7-3 to perform constrained least squares filtering on the image G4(u,v) obtained in step 4-4

[0140]

[0141] Among them, η4 represents the noise level of the image, which is selected according to empirical values. Performing an inverse Fourier transform on G6(u, v), the restored image g6(x, y) can be obtained, and the results are shown in Figures 8(a) and 8(b). This attached figure is the restored image after removing the blur and pixel distortion caused by the aerodynamic optical effect from the degraded image.

[0142] 7-5: Calculate the average gradient of the sharpness of the image g6(x, y) obtained in step 7-4 and the coordinates of the maximum point of the image pixel values. Whether the average gradient of the sharpness has an improvement of more than 5% compared to the value in step 4-5 and whether the coordinates of the maximum point of the image pixel values are different from the values in step 1-3. If all are met, proceed to the next step; otherwise, return to step 7-4, appropriately increase the value of the noise level η4 of the image, and recalculate.

[0143] In this embodiment, step 8 is: Calculate the comprehensive image quality evaluation parameter, and the main steps are as follows:

[0144] 8-1: Calculate the SSIM index of the finally restored clear images g5(x, y) and g6(x, y)

[0145] SSIM = [L((g(x, y), g′(x, y))] α ·[C((g(x, y), g′(x, y))] β ·[S((g(x, y), g′(x, y))] γ

[0146] Among them, L((g(x, y), g′(x, y)), C((g(x, y), g′(x, y)), and S((g(x, y), g′(x, y)) are the gray similarity, contrast similarity, and structure similarity between the reference image and the distorted image respectively, and α, β, γ are weights, generally taking α = β = γ = 1.

[0147] 8-2: Calculate the MSSIM index of the finally restored clear images g5(x, y) and g6(x, y)

[0148]

[0149] Among them, M represents the number of times of low-pass filtering and 1 / 2 downsampling of the reference image and the distorted image, generally taking α M = β j = γ j And

[0150] 8-3: Calculate the GSSIM index of the finally restored clear images g5(x, y) and g6(x, y)

[0151] GSSIM = [L((g(x, y), g′(x, y))]α ·[C((g(x,y),g′(x,y))] β ·[G((g(x,y),g′(x,y))] γ

[0152] where G((g(x,y),g′(x,y)) is the structural similarity based on gradient components.

[0153] 8-4: Calculate the comprehensive image quality parameter

[0154]

[0155] where τ1, τ2, and τ3 are weights, and generally τ1 = τ2 = τ3 = 1. When the comprehensive image quality parameter is greater than a certain threshold, it is considered that the image restoration effect is good and the image restoration is completed. This threshold generally takes a value greater than 0.95; otherwise, return to step 2, appropriately increase the noise level of the image in the formula, and recalculate.

[0156] The method of the embodiment of the present invention, compared with the prior art, fully considers the mechanism of image degradation, applies multiple algorithms for different degradation methods, and can obtain a better restoration effect than the image restoration using a single algorithm. At the same time, in view of the problem of difficult acquisition of measured data, a flow field modeling method of background schlieren technology is added, which increases the practicability of the method. It can realize the whole process from image preprocessing to final image quality evaluation, and the accuracy and quality of the restoration meet the requirements.

[0157] The accompanying drawings of the embodiment of the present invention can make the purpose, technical solution and advantages of the present invention introduced more clearly and clearly. It should be noted that the specific embodiments described here are only used to explain the present invention and are not used to limit the present invention. Any equivalent replacement, improvement, etc. made within the method ideas and principles provided by the present invention should be included in the protection scope of the present invention.

Claims

1. An image restoration method under complex optical imaging conditions based on measured data, characterized in that, It includes the following steps: Step 1: Calculate the peak signal-to-noise ratio (PSNR), the average gradient of sharpness, and the coordinates of the maximum point of the image pixel values of the original degraded image \(g_0(x,y)\); Step 2: Perform morphological filtering on the degraded image to remove the noise introduced by the imaging system during the imaging process and the noise generated by part of the pneumatic thermal radiation effect, obtaining the denoised image \(g_1(x,y)\); Calculate the PSNR of the image \(g_1(x,y)\), and determine whether the improvement amount of this value compared with the PSNR value obtained in Step 1 is greater than 5 dB. If so, proceed to the next step; otherwise, repeat Step 2; Step 3: Simulate multiple point-light source images passing through the imaging system, and average the multiple point-light source images to reduce the influence of noise; After performing Fourier transform on the averaged image, perform edge detection. After processing, obtain the estimated values of the blur radius \(r\) and the standard deviation \(\sigma\) of the imaging system. Through Gaussian function modeling, obtain the blur point spread function \(H_1(u,v)\) of the imaging system: where \(\sigma\) is the standard deviation of the normal distribution, and \((u,v)\) represents the frequency coordinates of the image; Use the blur point spread function \(H_1(u,v)\) of the imaging system to perform least-squares filtering on the image \(g_1(x,y)\) obtained in Step 2, obtaining the image \(G_2(u,v)\) with the imaging system blur removed where G1(u, v) is the Fourier spectrum of the image g1(x, y), H1 * (u, v) represents the complex conjugate of H1(u, v), η1 represents the noise level of the image. Performing the inverse Fourier transform on G2(u, v) yields the restored image g2(x, y). Based on the flow field data near the detection window, calculate the radiation intensity of the aerodynamic thermal radiation effect. Subtract this matrix from the image g2(x, y) to obtain the restored image g3(x, y) with the aerodynamic thermal radiation removed; Step 4: Perform image processing on the spectrum \(G_3(u,v)\) of the image \(g_3(x,y)\) obtained in Step 3 and then perform Radon transform; The angle corresponding to the maximum value of the Radon transform curve is the motion blur angle of the degraded image. Calculate the difference between the two sides of the main peak in the blur angle direction, which is the motion blur scale of the degraded image. Then, obtain the vector line segment containing the angle and the scale as the motion blur point spread function \(H_2(u,v)\); Use the motion blur point spread function \(H_2(u,v)\) to perform least-squares filtering on the image \(g_3(x,y)\) with the imaging system blur removed obtained in Step 3, generating the image \(G_4(u,v)\) with the motion blur removed: where \(G_3(u,v)\) is the Fourier transform of the image \(g_3(x,y)\), and \(\eta_2\) represents the noise level of the image; Perform inverse Fourier transform on \(G_4(u,v)\) to obtain the restored image \(g_4(x,y)\); Calculate the average gradient of sharpness of the image \(g_4(x,y)\), and determine whether there is an improvement of more than 5% compared with the value in Step 1. If so, proceed to the next step; otherwise, return to Step 3, and at the same time appropriately increase the values of the noise levels \(\eta_1\) and \(\eta_2\) of the image; Step 5: Determine whether it is possible to obtain the flow field data under complex imaging conditions in the wind tunnel test or the actual measurement scenario. For the case where the flow field data can be obtained, continue to execute Step 6; for the case where the flow field data cannot be obtained, directly execute Step 7; Step 6: Simulate the incident light to determine the light incident angle. Adopt the method of fixed step size to calculate the optical path length OPL of the light after passing through the entire flow field, and calculate the optical path difference of the light after averaging it Calculate the point spread function H3(u, v) of the flow field according to the optical path difference; use H3(u, v) to perform least squares filtering on the image g4(x, y) obtained in Step 4 where η3 represents the noise level of the image; the final clear restored image G5(u, v) under complex optical imaging conditions is obtained, and the inverse Fourier transform is performed on G5(u, v), that is, the restored image g5(x, y) is obtained; calculate the average gradient of sharpness of the image g5(x, y) and the coordinates of the point with the maximum pixel value of the image. Whether the average gradient of sharpness is increased by more than 5% compared to the value in step 4 and whether the coordinates of the point with the maximum pixel value of the image are different from the value in step 1. If all are met, proceed to the next step; otherwise, re-execute step 6, and at the same time appropriately increase the value of the noise level η3 of the image; Step 7: Obtain the background random dot matrix image, and take a set of correlation images under two conditions of no flow field and with flow field respectively; perform cross-correlation operations on multiple sets of correlation images, and perform time-averaging processing after obtaining the displacement change vectors of multiple cross-correlation interrogation areas; Perform SouthWell linear integral operation on the processed displacement vectors to obtain the wavefront optical path difference OPD, and finally calculate the flow field point spread function H4(u, v) according to OPD; use the PSF to perform least squares filtering on the image g4(x, y) obtained in step 4 where η4 represents the noise level of the image; the final clear restored image G6(u, v) under complex optical imaging conditions is obtained, and the inverse Fourier transform is performed on G6(u, v), that is, the restored image g6(x, y) is obtained; calculate the average gradient of sharpness of the image g6(x, y) and the coordinates of the point with the maximum pixel value of the image. Whether the average gradient of sharpness is increased by more than 5% compared to the value in step 4 and whether the coordinates of the point with the maximum pixel value of the image are different from the value in step 1. If all are met, proceed to the next step; otherwise, re-execute step 7, and at the same time appropriately increase the value of the noise level η4 of the image; Step 8: Calculate the comprehensive image quality evaluation index. When the comprehensive image quality parameter is greater than a certain threshold, it is considered that the image restoration effect is good, and the image restoration is completed; otherwise, return to step 2, appropriately increase the noise level of the image in the formula, and recalculate.

2. The image restoration method under complex optical imaging conditions based on measured data according to claim 1, characterized in that, The average gradient of sharpness in step 1 is defined as Among them k x (x, y) and k y (x, y) represents the image gradient calculated by the Sobel operator at the pixel point (x, y) in the x and y directions, and l(x, y) represents the image gradient calculated by the Laplacian operator at the pixel point (x, y).

3. The image restoration method under complex optical imaging conditions based on measured data according to claim 1, characterized in that, In step 3, the point light source image passing through the imaging system is a blurred image under the point light source, which contains certain imaging system noise. Therefore, multiple images contaminated by noise and blur under the same point light source are obtained, and the multiple images are first weighted and averaged to suppress noise.

4. The image restoration method under complex optical imaging conditions based on measured data according to claim 1, characterized in that, In step 3, first use the Sobel operator for edge detection to find the point with the maximum image gradient to find the edge. The edge points are set to 1, and the non-edge points are set to 0. After determining the edge points, calculate the radius size of the non-0 region to obtain the parameter r, and then use the single variable method to determine the optimal standard deviation σ according to the maximum PSNR of the restored image.

5. The image restoration method under complex optical imaging conditions based on measured data according to claim 1, characterized in that, After obtaining the spectrum G3(u, v) of the image g3(x, y) in step 4, perform the following image processing on the image in sequence: a) Perform logarithmic transformation on G3(u, v) to increase the contrast of the spectrum image; b) Edge cropping to remove the influence of the bright lines around the image spectrum on the detection result; c) Manually move the point where u = v = 0 to the center position; d) Edge detection; e) Set the edge points to 1 and the rest of the points to 0 to obtain a binary matrix; Perform subsequent Radon transform on the processing results.

6. The image restoration algorithm under complex optical imaging conditions based on measured data according to claim 1, wherein, The specific process of Step 6 is as follows: a) Simulate the incident light rays. The target for observation is a point light source at infinity. Therefore, the incident light rays can be approximated as a beam of parallel light, and the incident range is the effective window range of the optical imaging system. Determine the angle of incidence of the light rays relative to the flow field when the angle of attack of the aircraft is 0 degrees. The fixed step size is 0.001 m. Each movement of one step size is considered as a refraction occurring, and the refractive index of air in the next step is taken as the refractive index of the midpoint of the flow field closest to this point. Assume that the i-th light ray has traveled k step sizes, each step size is set as Δl, and the refractive index in each step size is n j , and the optical path obtained by the ray tracing method is b) The chief ray is set as a ray directly facing the center of the optical axis of the imaging system, and let its optical path difference be The optical path difference of the i-th ray is The position at the pupil imaging is (x i , y i ), then the wavefront aberration generated by the i-th ray passing through the flow field is where λ represents the wavelength of light; c) The actual pupil function is equivalent to the product of the pupil function under ideal conditions and the complex amplitude function of the wave aberration, i.e., A(x,y) = A0(x,y) × e jw(x,y) where A0(x,y) is the pupil function under ideal conditions, and j is the complex number identifier; d) Perform two-dimensional Fourier transform on the pupil function to obtain the diffraction complex amplitude distribution: e) Multiply the complex amplitude distribution by its conjugate function and then perform normalization processing to obtain the final point spread function H3(u,v) H3(u, v) = U(u, v) × U * (u, v) After obtaining the point spread function H3(u,v), perform subsequent restoration on the degraded image.

7. The image restoration method under complex optical imaging conditions based on measured data according to claim 1, wherein, The specific process of Step 7 is as follows: a) Obtain the random background lattice without flow field and the random background lattice with flow field; b) Perform enhanced contrast processing on the background lattice image passing through the refractive index field, perform multiple cross-correlation calculations, and obtain pixel offset information; c) Adopt the SouthWell fitting algorithm, integrate the surrounding four points and take the average to obtain the optical path of each point, i.e., where σ represents the weight, the weight of non-existent points is taken as 0, the weights of the remaining points are taken as 1, (i,j) represents the coordinates of the image pixel points, H is the distance from the background to the center of the flow field, Δx is the offset distance of the pixel in the x direction of each interrogation area, Δy is the offset distance of the pixel in the y direction of each interrogation area, h and l are the height and width of the interrogation area respectively; assuming that the initial value of the wavefront at the point (0,0) is 0, then the initial value of the wavefront of each interrogation area is where h and l are the height and width of the interrogation area respectively; then perform iteration, and the number of iterations is determined according to the relative difference between two adjacent steps; due to the discontinuity of the difference, the wavefront data may be discontinuous at each value point, so smooth the wavefront data; d) The chief ray is set as a ray directly facing the center of the optical axis of the imaging system, and its optical path difference is The optical path difference of the i-th ray is The position at the pupil imaging is (x i , y i ), then the wavefront aberration generated by the i-th ray passing through the flow field is where λ represents the wavelength of light; e) The actual pupil function is equivalent to the product of the pupil function under ideal conditions and the complex amplitude function of the wave aberration, i.e., A(x,y) = A0(x,y) × e jw(x,y) where A0(x,y) is the pupil function under ideal conditions, and j is the complex number identifier; f) Perform two-dimensional Fourier transform on the pupil function to obtain the diffraction complex amplitude distribution g) Multiply the complex amplitude distribution by its conjugate function and then perform normalization processing to obtain the final point spread function H3(u,v) H3(u, v) = U(u, v) × U * (u, v) After obtaining the point spread function H4(u,v), perform subsequent restoration on the degraded image.

8. The image restoration method under complex optical imaging conditions based on measured data according to claim 7, characterized in that, When performing cross-correlation calculation, according to the image size, use interrogation areas with a window of 32*32 and a step size of 16*16, a window of 16*16 and a step size of 8*8, and a window of 8*8 and a step size of 4*4 to perform cross-correlation calculation. Finally, obtain 128 interrogation areas and obtain the pixel offset conditions in the 128 interrogation areas.

9. The image restoration method under complex optical imaging conditions based on measured data according to claim 1, characterized in that, The comprehensive image quality parameter in Step 8 is defined as where τ1, τ2, τ3 are weights; SSIM is defined as: SSIM = [L(g(x,y), g ′ (x,y))] α ·[C(g(x,y), g′(x,y))] β ·[S(g(x,y), g′(x,y))] γ where L(g(x,y),g′(x,y)), C(g(x,y),g′(x,y)), S(g(x,y),g′(x,y)) are the gray similarity, contrast similarity, and structure similarity between the reference image and the distorted image respectively, and α, β, γ are weights; MSSIM is defined as: where M represents the number of times of low-pass filtering and 1 / 2 downsampling on the reference image and the distorted image, α M = β j = γ j and GSSIM is defined as: GSSIM = [L(g(x, y), g′(x, y))] α ·[C(g(x, y), g′(x, y))] β ·[G(g(x, y), g′(x, y))] γ where G(g(x, y), g′(x, y)) is the structural similarity between the reference image and the distorted image based on the gradient component.

Images