CN102779333B

CN102779333B - Optical image restoration method based on Kalman filter

Info

Publication number: CN102779333B
Application number: CN201210237148.8A
Authority: CN
Inventors: 王密; 张炳先; 李德仁
Original assignee: Wuhan University WHU
Current assignee: Wuhan University WHU
Priority date: 2012-07-10
Filing date: 2012-07-10
Publication date: 2015-01-07
Anticipated expiration: 2032-07-10
Also published as: CN102779333A

Abstract

An optical image restoration method based on Kalman filter, comprising the following steps: 1. Extraction of line spread function; 2. Gaussian fitting of line spread function; 3. Generation of high-precision line spread function based on Kalman filter; 4. MTF curve and MTF matrix generation; 5. Image restoration using MTF matrix. Based on MTF theory and image quality degradation theory, the method adopts Gaussian fitting method and Kalman filter method to correct the interference components in the measured line spread function curve, so as to obtain a point spread function that is better and closer to the real imaging performance of the imaging system as an image. The basis of restoration, so as to improve the accuracy of image restoration from the source. The implementation of this method can effectively restore the image and improve the image quality, and the algorithm efficiency is also very high.

Description

An Optical Image Restoration Method Based on Kalman Filter

技术领域 technical field

本发明属于遥感图像处理领域，涉及一种基于卡尔曼滤波的光学影像复原方法。The invention belongs to the field of remote sensing image processing and relates to an optical image restoration method based on Kalman filtering.

背景技术 Background technique

卫星在轨运行期间，由于平台颤振、偏流角等因素的影响导致TDI-CCD相机在积分成像过程中会产生轨道方向和垂轨方向的像点位移，引起影像模糊现象，这严重影响了影像的后续应用，因此有必要对影像数据进行复原处理。图像复原依据的是成像模型，准确的成像模型是高精度图像复原的基础，调制传递函数(MTF)作为图像成像模型中的重要组成成分，其值直接反映了图像在不同频率的衰减情况，基于影像质量退化理论可知，如果成像系统的MTF值能够精确测量或者计算得到，那么真实影像就可以被高质量恢复。近年来，卫星在轨运行时成像系统MTF的测量方法主要是采用刀刃法或脉冲法直接从在轨获取的成像数据中提取出来，这样能够很好的解决实验室测量数据没有考虑大气和探元老化因素而导致的测量结果不准确，但是由于测量所需的样本数据是直接从成像数据中获取，因此，测量结果会受到卫星平台振动、成像数据中的噪声成分等因素的影响而导致测量结果的准确度降低，因此，采用刀刃法或脉冲法获取的MTF曲线需要剔除其中所包含的干扰频率成分，否则，复原结果会由于MTF中频率测量值的不准确而出现噪声增加的现象。目前，常用的MTF曲线干扰成分剔除的方法主要是采取整平的方法将线扩展函数(LSF)两侧由于干扰引起的波动成分直接置零或者取平均值，这种方法虽然能够去除大部分噪声成分的影响，但是也会剔除掉非干扰成分，引起MTF曲线测量的不准确，从而进一步导致复原图像中噪声的增加。During the operation of the satellite in orbit, due to the influence of factors such as platform flutter and drift angle, the TDI-CCD camera will produce image point displacement in the orbital direction and vertical orbital direction during the integral imaging process, causing image blurring, which seriously affects the image quality. Therefore, it is necessary to restore the image data. Image restoration is based on the imaging model, and an accurate imaging model is the basis of high-precision image restoration. As an important component of the image imaging model, the modulation transfer function (MTF), whose value directly reflects the attenuation of the image at different frequencies, is based on The theory of image quality degradation shows that if the MTF value of the imaging system can be accurately measured or calculated, then the real image can be restored with high quality. In recent years, the MTF measurement method of the imaging system when the satellite is in orbit is mainly to use the knife-edge method or the pulse method to directly extract the imaging data obtained on the orbit, which can solve the problem that the laboratory measurement data does not consider the atmosphere and detector elements. The measurement results are inaccurate due to aging factors, but since the sample data required for measurement are directly obtained from the imaging data, the measurement results will be affected by factors such as satellite platform vibration and noise components in the imaging data. Therefore, the MTF curve obtained by the knife-edge method or the pulse method needs to eliminate the interference frequency components contained in it, otherwise, the restoration result will increase the noise due to the inaccurate frequency measurement value in the MTF. At present, the commonly used method to remove the interference components of the MTF curve is to use the leveling method to directly set the fluctuation components caused by interference on both sides of the line spread function (LSF) to zero or take the average value. Although this method can remove most of the noise However, it will also eliminate the non-interfering components, which will cause inaccurate measurement of the MTF curve, which will further lead to the increase of noise in the restored image.

发明内容 Contents of the invention

本发明所要解决的问题是，针对国产遥感影像数据，提供一种高精度的影像复原方法。The problem to be solved by the present invention is to provide a high-precision image restoration method for domestic remote sensing image data.

本发明的技术方案为一种基于卡尔曼滤波的光学影像复原方法，其特征在于：基于含刀刃靶标的影像，执行以下步骤，The technical solution of the present invention is an optical image restoration method based on Kalman filtering, which is characterized in that: based on the image of the knife-edge target, the following steps are performed,

步骤1，进行线扩展函数提取，包括以下子步骤，Step 1, extracting the line extension function, including the following sub-steps,

步骤1.1，进行边缘检测，得到影像每一行上边缘点的子像素位置；Step 1.1, perform edge detection to obtain the sub-pixel position of the edge point on each line of the image;

步骤1.2，进行边缘拟合，根据步骤1.1所得影像每一行上边缘点的子像素位置进行线性拟合，得到拟合直线，将影像中每一行上边缘点的子像素位置调整到拟合直线上；Step 1.2, carry out edge fitting, perform linear fitting according to the sub-pixel positions of the upper edge points of each row of the image obtained in step 1.1, and obtain a fitted straight line, and adjust the sub-pixel positions of the upper edge points of each row in the image to the fitted straight line ;

步骤1.3，提取边缘扩展函数；Step 1.3, extracting the edge extension function;

步骤1.4，根据步骤1.3所得边缘扩展函数，计算获取线扩展函数；Step 1.4, according to the edge extension function obtained in step 1.3, calculate and obtain the line extension function;

步骤2，对步骤1获取的线扩展函数进行高斯拟合，得到空间状态分量；Step 2, Gaussian fitting is performed on the line spread function obtained in step 1 to obtain the space state component;

步骤3，根据步骤2所得空间状态分量和步骤1获取的线扩展函数，基于卡尔曼滤波方法拟合生成更高精度的线扩展函数；Step 3, according to the spatial state component obtained in step 2 and the line spread function obtained in step 1, a higher precision line spread function is generated based on the Kalman filtering method;

步骤4，根据步骤3所得线扩展函数生成调制传递函数矩阵；Step 4, generating a modulation transfer function matrix according to the line spread function obtained in step 3;

步骤5，根据步骤4所得调制传递函数矩阵对原始影像进行影像复原，得到复原后的影像。Step 5: Restoring the original image according to the MTF matrix obtained in Step 4 to obtain a restored image.

而且，步骤1.3中提取边缘扩展函数的实现方式如下，Moreover, the implementation of extracting the edge extension function in step 1.3 is as follows,

根据已有的采样点信息，对影像每一行，以步骤1.2调整后边缘点的子像素位置为起点分别向左和向右进行等间隔采样，采样的间隔为像元宽度，再采用三次样条插值方法以0.05像元间隔插值采样，将所得包括插值点在内的采样点灰度分布作为该行的边缘扩展函数，According to the existing sampling point information, for each line of the image, start from the sub-pixel position of the edge point adjusted in step 1.2 to sample at equal intervals to the left and right, the sampling interval is the pixel width, and then use the cubic spline The interpolation method uses 0.05 pixel interval interpolation sampling, and the obtained gray distribution of sampling points including interpolation points is used as the edge extension function of the line,

然后，对每一行的边缘扩展函数的曲线以边缘点为基准求平均值，获得影像的边缘扩展函数的曲线。Then, the curve of the edge spread function of each row is averaged based on the edge point to obtain the curve of the edge spread function of the image.

而且，步骤1.4中计算线扩展函数的实现方式为，对边缘扩展函数在每一采样点的子像素位置进行简单差分。Moreover, the implementation of calculating the line extension function in step 1.4 is to perform a simple difference on the sub-pixel positions of the edge extension function at each sampling point.

而且，步骤5中进行影像复原的实现方式为，在频率域中，利用调制传递函数矩阵采用预设的滤波器算子进行影像复原，数学模型如下，Moreover, the implementation of image restoration in step 5 is to use the modulation transfer function matrix and the preset filter operator to perform image restoration in the frequency domain. The mathematical model is as follows,

R(u,v)＝I(u,v)×P(u,v)R(u,v)=I(u,v)×P(u,v)

其中，R(u,v)为复原后的影像的频谱，I(u,v)是原始影像的频谱，P(u,v)是滤波器算子，(u,v)为影像二维频率坐标为。Among them, R(u,v) is the spectrum of the restored image, I(u,v) is the spectrum of the original image, P(u,v) is the filter operator, and (u,v) is the two-dimensional frequency of the image The coordinates are .

本发明所提供方法基于MTF理论和影像质量退化理论，并且采用高斯拟合方法和卡尔曼滤波方法校正实测线扩展函数曲线中的干扰成分，从而得到更好、更接近成像系统真实成像性能的点扩散函数作为影像复原的依据，从而从源头上提高影像复原的精度。该方法简单、易于实现，具有较强的通用性；能够有效地复原影像、改善影像质量，算法效率也很高。The method provided by the present invention is based on MTF theory and image quality degradation theory, and uses Gaussian fitting method and Kalman filtering method to correct the interference components in the measured line spread function curve, so as to obtain better and closer to the real imaging performance of the imaging system. The diffusion function is used as the basis for image restoration, thereby improving the accuracy of image restoration from the source. The method is simple, easy to implement, and has strong versatility; it can effectively restore images and improve image quality, and the algorithm efficiency is also high.

附图说明 Description of drawings

图1为本发明实施例的流程图；Fig. 1 is the flowchart of the embodiment of the present invention;

图2为本发明实施例的高精度高斯拟合线扩展函数结果示意图。Fig. 2 is a schematic diagram of the results of the high-precision Gaussian fitting line extension function according to the embodiment of the present invention.

具体实施方式 Detailed ways

本发明技术方案可采用计算机软件技术实现自动运行流程，以下结合附图和实施例详细说明本发明技术方案。参见图1，实施例的流程可以分为五个步骤：The technical solution of the present invention can use computer software technology to realize the automatic operation process, and the technical solution of the present invention will be described in detail below in conjunction with the accompanying drawings and embodiments. Referring to Figure 1, the flow of the embodiment can be divided into five steps:

步骤1，进行线扩展函数提取，包括以下子步骤：Step 1, extracting the line extension function, including the following sub-steps:

步骤1.1，进行边缘检测，得到影像每一行上边缘点的子像素位置。Step 1.1, perform edge detection to obtain the sub-pixel position of the edge point on each line of the image.

实施例先初始估计边缘的位置，对每一行灰度值作简单差分，找到最大差分点的位置即为边缘的位置，然后取该位置附近4个样本点的值作三次多项式曲线拟合，设对某一行找到的边缘的位置拟合多项式为：In the embodiment, the position of the edge is initially estimated, and a simple difference is made to each line of gray value, and the position of the maximum difference point is found to be the position of the edge, and then the values of 4 sample points near the position are used for cubic polynomial curve fitting, assuming The position fit polynomial for the edge found for a row is:

a(x)＝a₁x³+a₂x²+a₃x+a₄ (1)a(x)＝a ₁ x ³ +a ₂ x ² +a ₃ x+a ₄ (1)

其中，x为影像的列号，a₁、a₂、a₃、a₄为系数，a(x)为该位置的DN值（遥感影像像元亮度值，即记录的地物的灰度值）。由于边缘点为最大斜率点，a(x)的二阶导为0，即：Among them, x is the column number of the image, a ₁ , a ₂ , a ₃ , and a ₄ are coefficients, and a(x) is the DN value of the position (the brightness value of the remote sensing image pixel, that is, the gray value of the recorded ground object ). Since the edge point is the point of maximum slope, the second derivative of a(x) is 0, namely:

a"(x)＝6a₁x+2a₂＝0 (2)a"(x)＝6a ₁ x+2a ₂ ＝0 (2)

解得：x＝-2a₂/6a₁＝-a₂/3a₁ Solution: x＝-2a ₂ /6a ₁ ＝-a ₂ /3a ₁

x即为所求的该行灰度序列的边缘点的子像素位置，用同样的方法求出影像中每一行边缘点的子像素位置。x is the sub-pixel position of the edge point of the gray-scale sequence to be obtained, and the sub-pixel position of the edge point of each line in the image is obtained by the same method.

步骤1.2，进行边缘拟合，根据步骤1.1所得影像每一行上边缘点的子像素位置进行线性拟合，得到拟合直线，将影像中每一行上边缘点的子像素位置调整到拟合直线上。Step 1.2, carry out edge fitting, perform linear fitting according to the sub-pixel positions of the upper edge points of each row of the image obtained in step 1.1, and obtain a fitted straight line, and adjust the sub-pixel positions of the upper edge points of each row in the image to the fitted straight line .

实施例对边缘检测结果进行线性拟合，使得边缘点都位于一条直线上。The embodiment performs linear fitting on the edge detection result, so that the edge points are all located on a straight line.

线性拟合最常用的方法是最小二乘法，假设直线方程为：The most commonly used method for linear fitting is the least squares method, assuming that the equation of the line is:

y＝ax+b (3)y＝ax+b （3）

其中，x为行数，y为该行最终的边缘位置。由最小二乘得系数：Among them, x is the number of rows, and y is the final edge position of the row. Coefficients obtained by least squares:

$a a = = \frac{m m (({Σ Σ}_{i i = = 11}^{m m} {x x}_{i i} {y the y}_{i i})) - - (({Σ Σ}_{i i = = 11}^{m m} {x x}_{i i})) (({Σ Σ}_{i i = = 11}^{m m} {y the y}_{i i}))}{m m (({Σ Σ}_{i i = = 11}^{m m} {x x}_{i i}^{22})) - - {(({Σ Σ}_{i i = = 11}^{m m} {x x}_{i i}))}^{22}} - - - - - - ((44))$

$b b = = \frac{(({Σ Σ}_{i i = = 11}^{m m} {x x}_{i i}^{22})) (({Σ Σ}_{i i = = 11}^{m m} {y the y}_{i i})) - - (({Σ Σ}_{i i = = 11}^{m m} {x x}_{i i} {y the y}_{i i})) (({Σ Σ}_{i i = = 11}^{m m} {x x}_{i i}))}{m m (({Σ Σ}_{i i = = 11}^{m m} {x x}_{i i}^{22})) - - {(({Σ Σ}_{i i = = 11}^{m m} {x x}_{i i}))}^{22}} - - - - - - ((55))$

其中，m为边缘点的个数，x_i为含刀刃靶标影像的第i行，y_i为第i行边缘点的子像素位置，i的取值为1,2…m。这些数据都已在边缘检测中求出，只需代入就可得到拟合的边缘直线。然后将每一行边缘点的子像素位置调整到该拟合直线上。Among them, m is the number of edge points, x _i is the i-th line of the target image containing the blade, y _i is the sub-pixel position of the edge point in the i-th line, and the value of i is 1, 2...m. These data have been obtained in the edge detection, and only need to be substituted to get the fitted edge straight line. Then adjust the sub-pixel position of each row of edge points to the fitted straight line.

步骤1.3，提取边缘扩展函数(ESF)。Step 1.3, extract the edge extension function (ESF).

实施例提取边缘扩展函数的实现方式如下，The implementation of the embodiment to extract the edge extension function is as follows,

步骤1.4，根据步骤1.3所得边缘扩展函数，计算获取线扩展函数(LSF)。Step 1.4, according to the edge extension function obtained in step 1.3, calculate and obtain the line extension function (LSF).

线扩展函数是通过对边缘扩展函数求导得到的，这里对边缘扩展函数在每一采样点的子像素位置进行简单差分：The line extension function is obtained by deriving the edge extension function. Here, a simple difference is made to the sub-pixel position of the edge extension function at each sampling point:

LSF(n)＝ESF(n)-ESF(n-1) (6)LSF(n)＝ESF(n)-ESF(n-1)

其中，ESF(n)表示边缘扩展函数在第n个采样点处的值，ESF(n-1)表示边缘扩展函数在第n-1个采样点处的值，LSF(n)表示线扩展函数在第n个采样点处的值。Among them, ESF(n) represents the value of the edge extension function at the nth sampling point, ESF(n-1) represents the value of the edge extension function at the n-1 sampling point, LSF(n) represents the line extension function The value at the nth sample point.

通过以上公式计算获取实测的线扩展函数后，为提高精度起见，本发明采用高斯拟合和卡尔曼滤波进行去噪处理，如后续步骤2和3。After calculating and obtaining the measured line extension function through the above formula, in order to improve the accuracy, the present invention adopts Gaussian fitting and Kalman filter for denoising processing, such as subsequent steps 2 and 3.

步骤2，对步骤1获取的线扩展函数进行高斯拟合，得到空间状态分量。Step 2, Gaussian fitting is performed on the line spread function obtained in step 1 to obtain the space state component.

高斯曲线的表达式如下所示：The expression for a Gaussian curve is as follows:

$y the y = = \frac{11}{\sqrt{22 π π} σ σ} {e e}^{- - \frac{{((x x - - μ μ))}^{22}}{22 {σ σ}^{22}}} - - - - - - ((77))$

其中，μ为均值(即位置参数)，σ为标准差，x为线扩展函数的采样点的计数值(单调递增)，y为当采样点计数值等于x时线扩展函数对应的值。Among them, μ is the mean value (that is, the position parameter), σ is the standard deviation, x is the count value of the sampling point of the line extension function (monotonically increasing), and y is the value corresponding to the line extension function when the count value of the sampling point is equal to x.

由式(7)可知，高斯曲线拟合关键在于确定高斯分布的中心位置和高斯分布的离散度，但是，对于受到噪声成分影响的近似高斯分布的曲线，采用上述公式并运用迭代的方式进行拟合计算，计算量巨大且精度不高。为了简化计算复杂度，首先需要对高斯函数进行线性化处理，其处理公式如式(8)所示：It can be seen from formula (7) that the key to Gaussian curve fitting is to determine the center position of the Gaussian distribution and the dispersion of the Gaussian distribution. Combined calculation, the calculation is huge and the accuracy is not high. In order to simplify the computational complexity, the Gaussian function needs to be linearized first, and the processing formula is shown in formula (8):

$In In ((y the y)) = = - - \frac{{((x x - - μ μ))}^{22}}{22 {σ σ}^{22}} + + In In ((\frac{11}{\sqrt{22 π π} σ σ})) - - - - - - ((88))$

令In(y)＝Y，x＝X则公式(8)可以简化为如下形式：Let In(y)=Y, x=X then formula (8) can be simplified to the following form:

Y＝aX²+bX+c (9)Y＝ ^aX2 +bX+c (9)

公式(9)是典型的二次曲线方程，那么对于高斯分布的拟合就变成了对二次曲线的拟合，依据最小二乘曲线拟合原理求取参数a,b,c的公式如下所示：Formula (9) is a typical quadratic curve equation, so the fitting of Gaussian distribution becomes the fitting of quadratic curve, and the formulas for obtaining parameters a, b, and c according to the principle of least squares curve fitting are as follows Shown:

$\{\begin{matrix} a a = = \frac{((n no {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22} - - {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i} {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i})) ((n no {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22} {y the y}_{i i} - - {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22} {Σ Σ}_{i i = = 11}^{n no} {y the y}_{i i})) - - ((n no {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i} {y the y}_{i i} - - {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i} {Σ Σ}_{i i = = 11}^{n no} {y the y}_{i i})) ((n no {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{33} - - {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22} {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}))}{((n no {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{44} - - {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22} {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22})) ((n no {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22} - - {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i} {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i})) - - ((n no {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{33} - - {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22} {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i})) ((n no {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{33} - - {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22} {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}))} \\ b b = = \frac{((n no {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i} {y the y}_{i i} - - {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i} {Σ Σ}_{i i = = 11}^{n no} {y the y}_{i i})) - - ((n no {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{33} - - {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22} {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i})) * * a a}{((n no {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22} - - {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i} {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}))} \\ c c = = \frac{11}{n no} (({Σ Σ}_{i i = = 11}^{n no} {y the y}_{i i} - - a a * * {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i}^{22} - - b b * * {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i})) \end{matrix} - - - - - - ((1010))$

其中，x_i指线扩展函数的采样点的计数值等于i，y_i为当采样点计数值等于x_i时线扩展函数对应的值，n表示采样点的个数。由于解是二次曲线的拟合结果，因此需要对其进行逆变换来获取所需的高斯曲线的拟合结果，其求解公式如下所示：Among them, _xi refers to the count value of the sampling point of the line extension function equal to i, y _i is the value corresponding to the line extension function when the count value of the sampling point is equal to x _i , and n represents the number of sampling points. Since the solution is the fitting result of the quadratic curve, it needs to be inversely transformed to obtain the required fitting result of the Gaussian curve. The solution formula is as follows:

$\{\begin{matrix} μ μ = = - - \frac{b b}{22 a a} \\ σ σ = = \sqrt{- - \frac{11}{22 a a}} \end{matrix} - - - - - - ((1111))$

为了使得高斯拟合结果更加精确，本发明建议采用步距拟合的方式对二次曲线拟合结果进行微量调整，即小尺度的改变高斯曲线的中心值和离散度，寻找与实测曲线差异最小的高斯拟合结果，实测曲线及其拟合结果如附图的图2所示，高斯拟合结果与实测曲线的差异很小。In order to make the Gaussian fitting result more accurate, the present invention proposes to use the method of step fitting to slightly adjust the fitting result of the quadratic curve, that is, to change the central value and dispersion of the Gaussian curve on a small scale, and to find the smallest difference from the measured curve The Gaussian fitting result, the measured curve and its fitting result are shown in Figure 2 of the accompanying drawings, and the difference between the Gaussian fitting result and the measured curve is very small.

步骤3，根据步骤2所得空间状态分量和步骤1获取的线扩展函数，基于卡尔曼滤波方法拟合生成更高精度的线扩展函数。Step 3: According to the spatial state component obtained in step 2 and the line spread function obtained in step 1, a higher-precision line spread function is generated based on the Kalman filter method.

卡尔曼滤波（KALMAN滤波）实质上是一个最优化自回归数据处理算法，其基本原理是根据某一时刻得到的系统状态值和实际测量值来拟合出该时刻系统的真实估计，本发明中，即步骤2所得高精度高斯拟合结果——空间状态分量作为系统状态值，而步骤1所得线扩展函数为实际测量值。Kalman filter (KALMAN filter) is essentially an optimized autoregressive data processing algorithm, and its basic principle is to fit the real estimate of the system at that time based on the system state value and actual measurement value obtained at a certain time. In the present invention , that is, the high-precision Gaussian fitting result obtained in step 2—the space state component is used as the system state value, and the line extension function obtained in step 1 is the actual measured value.

KALMAN滤波具体实现为现有技术，为便于实施参考起见，提供说明如下：The specific implementation of KALMAN filtering is an existing technology. For the convenience of implementation and reference, the description is provided as follows:

要了解KALMAN滤波的基本原理，需要引入一个离散控制过程的系统，该系统可以用一个线性随机微分方程来描述：To understand the basic principles of KALMAN filtering, it is necessary to introduce a system of discrete control processes, which can be described by a linear stochastic differential equation:

X(k)＝AX(k-1)+BU(k)+W(k) (12)X(k)＝AX(k-1)+BU(k)+W(k) (12)

再引入系统的实际测量值：Then introduce the actual measured value of the system:

Z(k)＝HX(k)+V(k) (13)Z(k)＝HX(k)+V(k)

式(12)、(13)中，X(k)是k时刻的系统状态，U(k)是k时刻对系统的控制量。A、B为系统参数，对于多模型系统，它们表现为矩阵形式。Z(k)是k时刻的测量值，H为测量系统的参数，对于多模型系统，它表现为矩阵形式。W(k)、V(k)分别表示过程和测量中的噪声，令它们的方差为Q、R。In formulas (12) and (13), X(k) is the state of the system at time k, and U(k) is the control amount of the system at time k. A and B are system parameters, and for multi-model systems, they are expressed in matrix form. Z(k) is the measured value at time k, H is the parameter of the measurement system, and for a multi-model system, it is in the form of a matrix. W(k), V(k) represent the noise in the process and measurement respectively, let their variance be Q, R.

首先，可以通过系统的过程模型预测下一状态的系统，令当前的系统状态为k，则：First, the system in the next state can be predicted through the process model of the system, let the current system state be k, then:

X(k|k-1)＝AX(k-1|k-1)+BU(k) (14)X(k|k-1)＝AX(k-1|k-1)+BU(k) (14)

式(14)中，X(k|k-1)是利用上一状态预测的结果，X(k-1|k-1)是上一状态的最优结果。U(k)是现在状态的控制量，如果没有，可以认定U(k)为0。In formula (14), X(k|k-1) is the result predicted by the previous state, and X(k-1|k-1) is the optimal result of the previous state. U(k) is the control quantity of the current state, if not, U(k) can be considered as 0.

在更新系统过程模型时，其对应的系统噪声方差也需要做同步更新，其更新公式为：When updating the system process model, the corresponding system noise variance also needs to be updated synchronously, and the update formula is:

P(k|k-1)＝AP(k-1|k-1)A'+Q (15)P(k|k-1)＝AP(k-1|k-1)A'+Q

在获得现有状态的预测结果后，结合现在状态下的实测值便可以得到当前状态的最优化估计值X(k|k)：After obtaining the prediction results of the current state, combined with the measured value in the current state, the optimal estimated value X(k|k) of the current state can be obtained:

X(k|k)＝X(k|k-1)+Gain(k)(Z(k)-HX(k|k-1)) (16)X(k|k)＝X(k|k-1)+Gain(k)(Z(k)-HX(k|k-1)) (16)

其中Gain(k)为卡尔曼增益：Where Gain(k) is the Kalman gain:

Gain(k)＝P(k|k-1)H'/(HP(k-1|k-1)H'+R) (17)Gain(k)＝P(k|k-1)H'/(HP(k-1|k-1)H'+R) (17)

其中H'为H的转置矩阵。在数学计算中，为了得到数据在函数收敛时的最优估计，需要令卡尔曼滤波不断的运行下去直至系统收敛，并不断更新k时刻的系统状态X(k)的噪声方差，使得噪声方差趋于最小，其方差更新公式如下所示:where H' is the transpose matrix of H. In mathematical calculations, in order to obtain the optimal estimate of the data when the function converges, it is necessary to make the Kalman filter run continuously until the system converges, and continuously update the noise variance of the system state X(k) at time k, so that the noise variance tends to At the minimum, its variance update formula is as follows:

P(k|k)＝(I-Gain(k)H)P(k|k-1) (18)P(k|k)＝(I-Gain(k)H)P(k|k-1)

其中，P(k|k)是系统状态X(k)对应的噪声方差，I为单位矩阵。Among them, P(k|k) is the noise variance corresponding to the system state X(k), and I is the identity matrix.

步骤4，根据步骤3所得线扩展函数生成调制传递函数矩阵。Step 4, generating a modulation transfer function matrix according to the line spread function obtained in step 3.

实施例在获得高精度的线扩展函数后，对线扩展函数曲线做傅立叶变换，取傅立叶变换后结果的模为各频率的MTF值，并以第一个MTF值为基准作归一化处理，即得到MTF曲线。然而利用刀刃边缘影像（即含刀刃靶标的影像）求得的MTF曲线只是水平或垂直方向的MTF曲线，而建立MTF复原模型必须要构建二维的MTF矩阵，常规的处理方式是将水平MTF列向量乘以垂直MTF列向量，即Embodiment After obtaining the high-precision line spread function, do Fourier transform to the line spread function curve, get the modulus of the result after Fourier transform to be the MTF value of each frequency, and do normalization process with the first MTF value as a benchmark, That is, the MTF curve is obtained. However, the MTF curve obtained by using the knife-edge image (that is, the image containing the knife-edge target) is only the MTF curve in the horizontal or vertical direction, and the establishment of the MTF restoration model must build a two-dimensional MTF matrix. The conventional processing method is the horizontal MTF. vector multiplied by the vertical MTF column vector, i.e.

MTF_(u,v)＝MTF_u×MTF_v (19)MTF _(u,v) = MTF _u × MTF _v (19)

其中，MTF_u为在频率u处水平的MTF值，MTF_v为在频率v处垂直的MTF值，MTF(u,v)为二维频率坐标为(u,v)处的MTF值。Among them, MTF _u is the horizontal MTF value at frequency u, MTF _v is the vertical MTF value at frequency v, and MTF(u,v) is the MTF value at the two-dimensional frequency coordinates (u,v).

这种方法求得的45°方向的MTF值与水平或垂直方向差别很大，因此，实施例中采用顾行发等人提出的矩阵构造方法，在求得水平、垂直方向的MTF曲线后，取水平与垂直方向0.5频率处MTF值的平均值再衰减90%作为45°方向0.5频率处的MTF值；再根据水平与垂直方向的MTF向量之间的比例关系进行插值，得到二维插值MTF矩阵。由于模的对称性，只需求出0~0.5频率处的MTF值，根据对称性即可得到-0.5~0频率处的MTF值(0.5频率即截止频率的一半)。The MTF value in the 45° direction obtained by this method is very different from that in the horizontal or vertical direction. Therefore, in the embodiment, the matrix construction method proposed by Gu Xingfa and others is adopted. After obtaining the MTF curves in the horizontal and vertical directions, take the horizontal The average value of the MTF value at 0.5 frequency in the vertical direction is attenuated by 90% as the MTF value at 0.5 frequency in the 45° direction; then interpolation is performed according to the proportional relationship between the MTF vectors in the horizontal and vertical directions to obtain a two-dimensional interpolation MTF matrix. Due to the symmetry of the mode, only the MTF value at the frequency of 0~0.5 is required, and the MTF value at the frequency of -0.5~0 can be obtained according to the symmetry (the frequency of 0.5 is half of the cut-off frequency).

在频率域中，可直接利用MTF矩阵采用各种滤波器算子进行影像复原。其数学模型表示为：In the frequency domain, the MTF matrix can be directly used for image restoration with various filter operators. Its mathematical model is expressed as:

R(u,v)＝I(u,v)×P(u,v) (20)R(u,v)＝I(u,v)×P(u,v)

其中，R(u,v)为复原后的影像的频谱，I(u,v)是原始影像的频谱，P(u,v)是选取的滤波器算子。Among them, R(u,v) is the spectrum of the restored image, I(u,v) is the spectrum of the original image, and P(u,v) is the selected filter operator.

实施例采用维纳滤波器算子。维纳滤波法考虑了噪声的影响，对噪声起到了抑制和减少的作用，维纳滤波器算子为：An embodiment employs a Wiener filter operator. The Wiener filtering method takes into account the influence of noise, and plays a role in suppressing and reducing noise. The Wiener filter operator is:

$P P ((u u,, v v)) = = \frac{11}{MTF MTF ((u u,, v v))} \times \times [[\frac{MTF MTF {((u u,, v v))}^{22}}{MTF MTF {((u u,, v v))}^{22} + + {k k}_{w w}}]],, - - - - - - ((21 twenty one))$

其中，k_w是与影像信噪比有关的先验常数，实施例取值0.02带入公式。MTF(u,v)为二维的MTF矩阵。采用传统方法和本实施例所提供方法进行影像复原的效果对比，其中基于传统方法的复原图像中的噪声信息明显增多。Wherein, k _w is a priori constant related to the image signal-to-noise ratio, and the value of 0.02 in the embodiment is brought into the formula. MTF(u,v) is a two-dimensional MTF matrix. The effect of image restoration using the traditional method and the method provided in this embodiment is compared, and the noise information in the restored image based on the traditional method is obviously increased.

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which the present invention belongs can make various modifications or supplements to the described specific embodiments or adopt similar methods to replace them, but they will not deviate from the spirit of the present invention or go beyond the definition of the appended claims range.

Claims

1. An optical image restoration method based on Kalman filtering, characterized in that: based on the image of the knife-edge target, the following steps are performed,

Step 1, extracting the line extension function, including the following sub-steps,

Step 1.1, perform edge detection to obtain the sub-pixel position of the edge point on each line of the image;

Step 1.2, carry out edge fitting, perform linear fitting according to the sub-pixel positions of the upper edge points of each row of the image obtained in step 1.1, and obtain a fitted straight line, and adjust the sub-pixel positions of the upper edge points of each row in the image to the fitted straight line ;

Step 1.3, extracting the edge extension function;

Step 1.4, according to the edge extension function obtained in step 1.3, calculate and obtain the line extension function;

Step 2, Gaussian fitting is performed on the line spread function obtained in step 1 to obtain the space state component;

Step 3, according to the spatial state component obtained in step 2 and the line spread function obtained in step 1, a higher precision line spread function is generated based on the Kalman filtering method;

Step 4, generating a modulation transfer function matrix according to the line spread function obtained in step 3;

Step 5: Restoring the original image according to the MTF matrix obtained in Step 4 to obtain a restored image.

2. the optical image restoration method based on Kalman filter as claimed in claim 1, is characterized in that: the implementation mode of extracting edge spread function in step 1.3 is as follows,

According to the existing sampling point information, for each line of the image, start from the sub-pixel position of the edge point adjusted in step 1.2 to sample at equal intervals to the left and right, the sampling interval is the pixel width, and then use the cubic spline The interpolation method uses 0.05 pixel interval interpolation sampling, and the obtained gray distribution of sampling points including interpolation points is used as the edge extension function of the line,

Then, the curve of the edge spread function of each row is averaged based on the edge point to obtain the curve of the edge spread function of the image.

3. The optical image restoration method based on Kalman filter as claimed in claim 2, characterized in that: the implementation of calculating the line spread function in step 1.4 is to carry out simple difference to the edge spread function at the sub-pixel position of each sampling point .

4. The optical image restoration method based on Kalman filtering as claimed in claim 1 or 2 or 3, characterized in that: when the line spread function obtained in step 1 is carried out to Gaussian fitting in step 2, first the Gaussian function is linearized Processing, the fitting of the Gaussian distribution is changed to the fitting of the quadratic curve, and the solution is solved according to the least squares curve fitting principle and the inverse transformation is performed to obtain the required Gaussian curve fitting result.

5. The optical image restoration method based on Kalman filtering as claimed in claim 4, characterized in that: the step fitting method is used to slightly adjust the fitting result of the quadratic curve, including changing the central value of the Gaussian curve on a small scale and dispersion, looking for the Gaussian fitting result with the smallest difference from the measured curve.

6. The optical image restoration method based on Kalman filter as claimed in claim 1 or 2 or 3, characterized in that: in step 4, the realization of generating the modulation transfer function matrix according to the line spread function obtained in step 3 is, in obtaining the level , After the MTF curve in the vertical direction, take the average value of the MTF value at the frequency of 0.5 in the horizontal and vertical directions and then attenuate 90% as the MTF value at the frequency of 0.5 in the 45° direction; then according to the proportional relationship between the MTF vectors in the horizontal and vertical directions Perform interpolation to obtain a two-dimensional interpolation MTF matrix.

7. The optical image restoration method based on Kalman filter as claimed in claim 1 or 2 or 3, characterized in that: the implementation of image restoration in step 5 is, in the frequency domain, using the modulation transfer function matrix to adopt the preset The filter operator for image restoration, the mathematical model is as follows,

R(u,v)=I(u,v)×P(u,v)

Among them, R(u,v) is the spectrum of the restored image, I(u,v) is the spectrum of the original image, P(u,v) is the filter operator, and (u,v) is the two-dimensional frequency of the image coordinate.