Method for temporal filtering of images based on motion compensation

Abstract

The application discloses a kind of based on motion compensation time domain filtering method for image, comprising: using autoregressive moving average model to the obtained gyroscope error is modeled, the change rule of data is used to estimate model parameter, and gyroscope error model is established;Using the gyroscope data of removing error obtained by gyroscope error model, the rotation matrix of camera is calculated, and the rotation image is compensated in reverse using rotation matrix, and the de-rotation image is obtained.Through kalman filtering to image motion vector sequence, image random jitter is obtained, and the de-rotation image is compensated according to random jitter, and the stable image and the subjective motion vector of image are obtained;Using the motion target detection algorithm in difference map based on large map coordinate system to detect target in stable image.The application has high detection sensitivity, and can detect target and its trajectory quickly through a small amount of frame track correlation.

Description

Technical Field

[0001] This invention belongs to the field of moving target detection in motion scenes, specifically a method for temporal filtering of images based on motion compensation. Background Technology

[0002] Random video jitter can severely affect image quality, and irregular changes in the position of targets within the image can greatly interfere with subsequent target recognition and tracking. In weapon system target recognition and tracking systems, platform jitter can lead to the loss or misjudgment of tracked targets, ultimately resulting in false alarms and missed detections, causing significant impact.

[0003] Commonly used image stabilization technologies are divided into three categories: mechanical stabilization, optical stabilization, and electronic stabilization. Mechanical stabilization, also known as electromechanical stabilization, uses a servo motor system and gyroscope sensors to compensate for the relative motion of the imaging system, thereby achieving image stabilization. Mechanical stabilization is mainly used on large platforms such as artillery and tanks, for example, to correct the tank's observation and aiming systems and fire control systems. It offers good real-time performance and stability, but mechanical stabilization platforms are large, complex in structure, expensive, and difficult to miniaturize. Optical stabilization uses optical elements such as variable wedges to adjust the position and angle according to jitter, keeping the optical path stable. Optical stabilization is relatively complex to control and maintain, therefore it is mainly used in small systems with small jitter amplitudes. Electronic stabilization technology processes the acquired image sequence, using image processing algorithms or inertial navigation devices to calculate the amount of jitter in the image sequence. Compared to traditional stabilization technologies, electronic stabilization technology integrates computer, microelectronics, and digital signal processing technologies, possessing advantages such as directness, accuracy, and low cost, and has broad development prospects.

[0004] The mainstream algorithms for moving target detection mainly include statistical learning methods such as inter-frame difference, optical flow, and background modeling. Inter-frame difference uses the subtraction of two or three adjacent frames in a time-series image to obtain a difference image, which is then used to select the moving target information through thresholding. This method is computationally simple and easy to implement, but two-frame difference methods generally struggle to obtain the complete outline of the moving target, easily producing "double images" and "holes" within the target, leading to inaccurate target detection. Background subtraction first constructs a background image, then subtracts the current frame image from the background, and detects the moving target based on the difference image. This method is simple and easy to implement, and can extract target feature data well, but it is very sensitive to changes in the external environment and is suitable for situations where the background is known. Finding an efficient and practical background model and updating it is the focus of research for this method. Optical flow segmentation uses the vector features of moving targets that change over time to detect moving regions in image sequences. It performs well when the camera is moving, but the algorithm is complex, computationally intensive, requires specialized hardware support, and is difficult to meet the requirements of real-time video stream processing.

[0005] In summary, the problems with existing technologies can be summarized as follows: high cost; large computational load; poor detection performance against complex background targets; and difficulty in hardware miniaturization. Summary of the Invention

[0006] This invention proposes a method for temporal filtering of images based on motion compensation.

[0007] The technical solution to achieve the objective of this invention is as follows: a method for temporal filtering of images based on motion compensation, comprising the following steps:

[0008] An autoregressive moving average model is used to model the obtained gyroscope error. The model parameters are estimated by utilizing the variation patterns of the data, thus establishing a gyroscope error model.

[0009] Using the gyroscope error model to obtain error-removed gyroscope data, the camera rotation matrix is ​​calculated. The rotation matrix is ​​then used to perform reverse compensation on the rotated image to obtain a derotated image.

[0010] A robust feature extraction algorithm based on an accelerated version is used to extract feature point vectors from the despin image after gyroscope rotation compensation. Coarse matching of the feature point vectors is performed using the nearest neighbor-second nearest neighbor method, followed by fine matching of the coarsely matched point pairs using an improved nearest neighbor-second nearest neighbor matching algorithm. The motion vectors of the inter-frame images are estimated using the correctly matched feature points. Kalman filtering of the image motion vector sequence yields random image jitter, which is then used to perform motion compensation on the despin image, resulting in a stable image and the image's subjective motion vector.

[0011] A moving target detection algorithm based on a large graph coordinate system is used to detect targets in stable images.

[0012] Compared with the prior art, the significant advantages of this invention are:

[0013] 1. This invention has low requirements for gyroscope accuracy. It can use software to estimate model parameters by utilizing the changing patterns of data, and finally establish a corresponding mathematical model, which better suppresses gyroscope errors and greatly reduces costs while achieving high accuracy.

[0014] 2. This invention improves the traditional image registration method, greatly reduces the amount of iterative calculation, obtains a higher feature point matching accuracy, has low computational complexity, and has real-time and efficient processing capabilities.

[0015] 3. This invention has strong target detection capabilities, effectively suppressing background prominence of targets, improving image signal-to-noise ratio, and enabling the screening of moving targets in complex moving backgrounds.

[0016] 4. This invention has high detection sensitivity and can quickly detect targets and their trajectories by associating a small number of frames of track data.

[0017] The present invention will now be described in further detail with reference to the accompanying drawings. Attached Figure Description

[0018] Figure 1 This is a flowchart of a method for temporal filtering of images based on motion compensation.

[0019] Figure 2 The result is the estimation result of random noise ARMA(2,2).

[0020] Figure 3 The images show the residual autocorrelation and partial autocorrelation plots obtained from running the MATLAB code.

[0021] Figure 4 The images are the original diagram and the racemic diagram.

[0022] Figure 5 This is a schematic diagram of an integral graph.

[0023] Figure 6 This is a box filter template.

[0024] Figure 7 This is a diagram showing the size distribution of the box filter in scale space.

[0025] Figure 8 This is a diagram illustrating the process of generating a Hesse determinant image.

[0026] Figure 9 A schematic diagram illustrating the determination of the main vector.

[0027] Figure 10 This is a Haar wavelet template.

[0028] Figure 11 This is a schematic diagram of nearest neighbor matching.

[0029] Figure 12 The effect diagram of the improved random sampling consensus algorithm.

[0030] Figure 13 This is a schematic diagram of difference plot suppression.

[0031] Figure 14 Curvature background suppression.

[0032] Figure 15 The average speed of the background trajectory is shown in the diagram.

[0033] Figure 16 This includes a flight path effect diagram and a flight path speed diagram. Detailed Implementation

[0034] like Figure 1 As shown, a method for temporal filtering of images based on motion compensation includes the following steps:

[0035] (10) Gyroscope Error Analysis and Modeling: The obtained gyroscope error is modeled using an autoregressive moving average model. The model parameters are estimated using the data variation patterns to establish the gyroscope error model. The specific steps are as follows:

[0036] (11) Establish an autoregressive moving average model of order ARMA(2,2).

[0037] Traditional autoregressive moving average (ARMA) model parameter estimation, when the data volume is large enough, uses the MATLAB function "ARMAX" to model the collected gyroscope random noise data using ARMA(p,q), provided that p and q are known. The general mathematical expression of the ARMA model is:

[0038]

[0039] In the formula, ε(k) is white noise with a mean of 0 and a constant variance, and p and q are the model order.

[0040] Change the order of the autoregressive moving average model to ARMA(2,2):

[0041] z(k)=a1z(k-1)+a2z(k-2)+b1ε(k-1)+b2ε(k-2)+ε(k)

[0042] z(k) is the error data of the kth frame after first-order difference processing, z(k-1) is the error data of the (k-1)th frame after first-order difference processing, z(k-2) is the error data of the (k-2)th frame after first-order difference processing, ε(k) is white noise with a mean of 0 and a constant variance, and a1, a2, b1, b2 are model parameters.

[0043] (12) Parameter estimation of autoregressive moving average model using Kalman filtering

[0044] Parameter estimation is performed using Kalman filtering, and the corresponding state variables are X = [a1(k), a2(k), b1(k), b2(k)]. T The observation matrix is

[0045] The estimation results of random noise ARMA(2,2) are shown in the figure below. Figure 2 As shown.

[0046] The collected data were modeled using ARMA(2,2) using MATLAB's built-in functions.

[0047] z(k)=0.205z(k-1)-0.106z(k-2)-0.635ε(k-1)-0.715ε(k-2)+ε(k)

[0048] The results show the autocorrelation and partial correlation of the residuals as follows: Figure 3 As shown.

[0049] In the actual iteration process, the iteration termination condition is set to the difference between the continuous parameter estimates being less than the threshold th. With th = 0.001, the ideal estimate can be reached after 2350 iterations on the previous data.

[0050]

[0051] The ARMA(2,2) model obtained using Kalman filtering is:

[0052] z(k)=0.197z(k-1)-0.082z(k-2)-0.642ε(k-1)-0.731ε(k-2)+ε(k)

[0053] (20) Rotation matrix estimation and compensation based on gyroscope information: Using the error-removed gyroscope data obtained from the gyroscope error model, the camera rotation matrix is ​​calculated, and the rotation matrix is ​​used to perform inverse compensation on the rotated image to obtain the derotated image. The specific steps are as follows:

[0054] (21) Calculate the rotation matrix

[0055] By integrating the angular velocity output by the gyroscope, the rotation angles around the three axes can be obtained. Let ωi be the measured angular velocity, the gyroscope sampling be discrete, and Δt be the gyroscope sampling period. Then, the rotation angles from t1 to t2 are:

[0056] The camera rotates at an angle θ along its three axes between the k-th frame and the reference frame. c Let vector θ c =[xyz], find vector θ c The modulus ξ:

[0057]

[0058] Where x is the camera's horizontal axis rotation angle, y is the camera's vertical axis rotation angle, z is the camera's optical axis rotation angle, and ξ is the normalized length, vector θ c Dividing by ξ, we obtain the normalized rotation vector e = [e1 e2 e3], as shown in the following formula:

[0059]

[0060] Rewrite the normalized rotation vector as a skew-symmetric matrix:

[0061]

[0062] According to the Rodriguez rotation matrix principle, the rotation matrix D can be obtained, as shown in the following formula:

[0063]

[0064] (22) Reverse compensation

[0065] Assume the coordinates of a point P in the world coordinate system w The coordinates of the i-th frame image are P i The coordinates of the j-th frame are P j According to the camera imaging model, we have:

[0066] P i =KRP w

[0067] Where K is the intrinsic parameter matrix and R is the rigid body transformation matrix.

[0068] Assume that the camera only undergoes rotational motion from frame i to frame j, and the rotation matrix is: Then there is

[0069]

[0070] Substituting the expression into the equation, we get...

[0071]

[0072] The above formula represents the pixel coordinate change when the camera rotates. To achieve reverse compensation, we only need to use the jitter component of the current frame to perform reverse compensation on the current frame, that is:

[0073]

[0074] Original image and racemogram as follows Figure 4 As shown.

[0075] (30) Translation estimation and compensation based on image matching: Feature point vectors of the derotated image after gyroscope rotation compensation are extracted based on an accelerated, robust feature algorithm. Coarse matching of the feature point vectors is performed using the nearest neighbor-second nearest neighbor method. A refined nearest neighbor-second nearest neighbor matching algorithm is then used to fine-match the coarsely matched point pairs. The motion vectors of the inter-frame image are estimated using the correctly matched feature points. Kalman filtering is applied to the image motion vector sequence to obtain random image jitter. Motion compensation is then performed on the derotated image based on the random jitter to obtain a stable image and the image's subjective motion vector. The specific steps are as follows:

[0076] (31) Find the integral graph:

[0077] To calculate the sum of pixel gray levels in a region of an image, you need to accumulate the gray levels of all pixels within that region. If you frequently need to calculate the sum of gray levels for sub-blocks of an image, you can use an integral image for quick calculation. An integral image is defined as follows:

[0078]

[0079] In the above formula, I(x, y) represents the grayscale value of pixel (x, y). The integral image is introduced to simplify the convolution operation of the Hessian determinant. The convolution calculation requires calculating the sum of pixel grayscale values ​​in different regions. Using an integral image, the complex accumulation operation can be transformed into simple arithmetic operations, thus reducing the computational load. Figure 5 The rectangle contains four points A, B, C, and D. The sum of the pixel grayscale values ​​within this matrix region is:

[0080]

[0081] In the formula, II(X) is the integral value of a point X in the integral image. From the above formula, it can be seen that as long as the integral image is known, calculating the sum of pixel grayscale values ​​within any rectangular region of the image is easy.

[0082] (32) Feature point detection based on the Hessian matrix:

[0083] An accelerated, robust feature detection algorithm detects feature points using the Hessian matrix. The Hessian matrix is ​​composed of the second derivative of the Laplacian function of Gauss, which is the result of convolving an image function with a Gaussian function, as shown in the equation:

[0084] L(x,y,σ)=G(x,y,σ)* / (x,y)

[0085] The two-dimensional Gaussian function G(x, y, σ) is defined as follows:

[0086]

[0087] Where σ is the scaling factor of the Gaussian function. The formula for the Hessian matrix is ​​as follows:

[0088]

[0089] In the formula, L xx (x, y, σ), L xy (x, y, σ), L yy (x, y, σ) is the second derivative of the function L(x, y, σ). In practical applications, the Gaussian function G(x, y, σ) is discretized, and a robust feature point based on an accelerated version is used to replace the Gaussian second-order differential operator with a box filter.

[0090] A Gaussian filter template with a scaling factor σ of 1.2 is approximately a 9x9 box filter template in both the yy and xx directions, as shown below. Figure 6 As shown.

[0091] The derotated image is filtered using three box filters of different sizes (box1, box2, and box3) to obtain the filtered image. The specific method is as follows:

[0092] The filtering formulas for three types of 9x9 box filters are as follows:

[0093]

[0094]

[0095]

[0096]

[0097] Where IMG is the derotation image, box1, box2, and box3 are three types of box filters, m and n are the number of rows and columns of the derotation image, II(x,y) is the value at the (x,y) coordinate position of the integral image, and D... xx D yy D xy These are images filtered by three different box filters;

[0098] (33) Establishment of scale space:

[0099] To achieve scale invariance in feature point detection, detecting features in only one scale space is insufficient. An accelerated, robust feature detection algorithm establishes a multi-scale space by convolving the same image with box filters of different sizes. If the output image of the filter template with a size factor σ of 1.2 and a size of 9×9 is used as the initial layer, the remaining scale layers can be obtained by expanding the initial box filter.

[0100] The relationship between σ for each layer and the filter template size L is as follows:

[0101]

[0102] The scale space is divided into four groups, each containing four layers. Each group represents a series of response maps of the same input image filtered by progressively enlarged filter templates. Each group consists of several fixed layers. Due to the discretization of the integral image, the minimum scale change between two layers is determined by the response length L of the second-order Gaussian differential filter to positive and negative speckles in the differential direction, which is 1 / 3 of the box filter template size. For a 9×9 template, L = 3. The response length of the next layer should be increased by at least 2 pixels based on L to ensure one pixel on each side, i.e., L = 5. Thus, the template size is 15×15. Following this pattern, we can obtain a sequence of templates with increasing sizes: 9×9, 15×15, 21×21, and 27×27. The length of the black and white regions is increased by an even number of pixels to ensure the existence of a center pixel.

[0103] The size distribution diagram of the box filter in scale space is shown below. Figure 7 As shown.

[0104] (34) Feature point localization:

[0105] Each scale-space box filter has three filter templates: box1, box2, and box3. Each template is convolved with the image to obtain the corresponding D. xx D yy D xy The value is then used to generate a determinant image using a formula. The calculation process is as follows: Figure 8 As shown.

[0106] When locating feature points, non-maximum suppression (NMS) needs to be applied to the three adjacent matrix images in each image set. NMS is also applied to the 3×3×3 stereo neighborhood in the multi-scale space. If the middle point is an extremum, it is considered a candidate feature point, and its scale space and location are recorded. Clearly, the top and bottom layers in a matrix image set cannot be used as intermediate layers for NMS. Typically, a threshold needs to be set for NMS to filter out feature points with small response values. A larger threshold results in fewer feature points; a suitable threshold can improve the robustness of the feature points.

[0107] By using the result of convolving the box filter with the image to replace the second derivative of the Gaussian Laplacian function, an approximate matrix of the Hessian matrix is ​​obtained, as shown in the equation:

[0108]

[0109] The determinant of the Hessian matrix at pixel (x, y) can be approximated as:

[0110] det(H approx ) = D xx D yy -(ωD xy ) 2

[0111] Because the box filter uses an approximation, a weighting factor ω is needed to balance the effects of the approximation when calculating the determinant; ω is typically set to 0.9. Using an integral plot can significantly improve the computational speed when calculating the components of the Hessian matrix.

[0112] (35) Feature point vector description:

[0113] To ensure rotation invariance of the feature vectors, the principal orientation of the feature points needs to be determined. Let σ be the scale of the feature point. Within a circular region centered on the feature point and with a radius of 6σ, the Haar wavelet response values ​​in the x and y directions of the image are statistically analyzed. The Haar wavelet template size is 4σ × 4σ. During the statistical analysis, weights are assigned based on the distance from the feature point to the center of the circle, with closer points receiving larger weights. Finally, the cumulative Haar wavelet response values ​​are calculated within a sliding sector region with a rotation center of 60° to determine the principal orientation of the feature point. Figure 9 As shown.

[0114] The feature point vector utilizes Haar wavelet features. Haar wavelets, such as... Figure 10 As shown, the left figure shows the calculation of the Haar wavelet response value in the horizontal direction, and the right figure shows the Haar wavelet response value in the vertical direction.

[0115] After determining the principal direction of the feature point, the feature description vector is calculated within a 20σ×20σ rectangular region centered on the feature point along this principal direction. This matrix is ​​divided into 4×4 sub-regions. Within each sub-region, the response values ​​in the x and y directions are calculated using a 2×2 Haar wavelet template, and ∑dx, ∑dy, ∑|dx|, and ∑|dy| are calculated respectively. Each sub-region can be represented by a vector:

[0116] v=[∑dx,∑|dx|,∑dy,∑|dy|]

[0117] Therefore, a complete feature point descriptor can be composed of vectors from 4×4 sub-regions, i.e., a 4×4×4 = 64-dimensional vector representation. Finally, this vector is normalized to obtain the feature vector.

[0118] (36) Coarse matching of feature points:

[0119] Feature point pairing is performed using feature point vectors. Two images taken in the same scene have similar feature descriptors for corresponding feature point pairs. A similarity metric can be used to determine the similarity between the two feature vectors. Euclidean distance is typically used to measure the similarity of feature points for feature point matching.

[0120] Suppose that the reference image and the image to be matched have feature points X and Y, respectively, and their feature vectors are X = (x1, x2, x3, ..., x...). n ), Y=(y1, y2, y3,...,y n If X and Y are denoted as X, then the Euclidean distance between them is:

[0121]

[0122] The bidirectional matching algorithm based on the nearest neighbor and second nearest neighbor method has the following specific process:

[0123] (1) Assume images A and B, P A Let P be the set of feature points in image A. B Let A be the set of feature points in image B. i For P A One of the feature vectors, B m B n For P B In and A i The eigenvectors of the nearest and second nearest neighbors in Euclidean distance, d m d n For its distance, if (d m / d n If ) < T, then B m These are the matching feature point pairs, typically T is set to 0.6. This yields the set of matching point pairs G from image A to image B.

[0124] (2) Perform the same operation as in step 1 on all feature points in image B to obtain the matching point pair set H from image B to image A. The intersection of point pair set G and point pair set H is the final matching point pair set.

[0125] The effect diagram based on nearest neighbor matching is shown below. Figure 11 As shown.

[0126] (37) Fine-grained feature point matching:

[0127] As can be seen from the camera imaging model, in the absence of camera distortion, the image space transformation can be represented by a projection transformation matrix, with the general expression as follows:

[0128]

[0129] The formula for transforming the coordinates of corresponding points in two images is:

[0130]

[0131]

[0132] As shown in the above equation, to find the 8 parameters in the projection transformation matrix M, at least 8 equations are needed, which are equivalent to 4 non-collinear point pairs. Due to the existence of mismatched point pairs, the solutions to the above equations will be inaccurate. Therefore, it is necessary to remove mismatched point pairs through certain methods. The random sampling consensus algorithm can be used to purify the point pairs of feature points to reduce the error of the transformation matrix. The random sampling consensus algorithm assumes that a certain set of data is the correct set. It fits the model parameters through this set of data and then judges how many data points within a certain error threshold can satisfy the model. Those that satisfy the threshold condition are called inliers, and those that do not are called outliers. The more inliers there are, the more data points satisfy the model, and the higher the reliability of the model. The filtering algorithm has two termination conditions: the maximum number of iterations and the maximum number of inliers. When the number of inliers is greater than the maximum number of inliers or the number of iterations is greater than the maximum number of iterations, the algorithm ends, and the set of inliers with the largest number of inliers is the correct data point. Taking projection transformation as an example, the steps of the random sampling consensus algorithm are as follows:

[0133] 1. Assuming there are k pairs of candidate feature points, randomly select 4 pairs of points, establish a system of equations, and solve for the 8 parameters of the transformation matrix M.

[0134] 2. Calculate the point obtained by transforming the remaining k-4 feature points through the transformation matrix M, and calculate the distance between this point and the candidate matching point.

[0135] 3. Determine if the distance is less than the threshold. If it is, the point is an interior point; otherwise, it is an exterior point.

[0136] 4. Count the number of interior points that satisfy the transformation matrix.

[0137] 5. Randomly select four more point pairs and repeat steps 1-4 several times, selecting the set with the largest number of interior points. Use this set to re-estimate the transformation model and obtain the final result.

[0138] An improved random sampling consensus algorithm is proposed, utilizing the relationship of equal relative slopes among perfectly matched points to eliminate mismatched points, and then calculating the homography matrix. Let (A i A j ) and (B i B j Let A and B be a pair of correct matches between the reference image A and the image to be matched B. Then, A... i and B i The absolute slope k(A) i B i It should be similar to A j and B j The absolute slope k(A) j B j Based on the above principles, this paper utilizes A in reference image A. i With all the feature points A in itself j The slope relationship and B in the image to be matched i With all the feature points B in itself j To evaluate the similarity of the slope relationship between two points, the following evaluation function is proposed:

[0139]

[0140]

[0141] K(i,j)=[k(A i B i )+k(A j B j )] / 2

[0142] Where R(i, j) represents A i B i The relative difference in slope between each interest point and the image; K(i,j) represents A i B i The average slope of each interest point in the respective image. The evaluation value T(i) is removed in... Points outside the range, among which is the average value of the evaluation value T(i), and th is the elimination threshold.

[0143] The improved random sampling consensus algorithm is shown in the following figure. Figure 12 As shown:

[0144] (38) Estimate the motion vector:

[0145] Constructing the derotated image motion model is as follows:

[0146]

[0147] Where p0 is a point in reference image A, p1 is the corresponding point in the current image B, and d x It is the horizontal motion vector, d y It is a vertical motion vector;

[0148] The motion model parameters between frames, i.e., the motion vector d, are estimated by correctly matching point pairs. x d y :

[0149]

[0150] In the formula, (x i y i Let (X) be the coordinates of the matching point in reference image A in the i-th correctly matched feature point pair, and (X) be the coordinates of the matching point in reference image A. i Y i ) represents the coordinates of the matching point B in the current image of the i-th correctly matched feature point pair.

[0151] Using Kalman filtering, the horizontal motion vector d is processed respectively. x Vertical motion vector d y Filtering is performed separately. The system's state variables are:

[0152]

[0153] In the formula, x is the motion vector [x, y], v is the change in x, and N(0, σ) is Gaussian white noise.

[0154] Expanding the above equation, we get:

[0155]

[0156] In the above formula, N(0, σ) x ), N(0, σ y ), which are independent Gaussian white noises. σ x σ y The variances are the process noise variances. Their magnitudes determine the smoothness of the image. The larger the variance, the more variable the state variables are, which increases the randomness of the estimated values ​​and leads to instability in the output image sequence. Therefore, a smaller variance value should be chosen.

[0157] The system measurement equation is:

[0158]

[0159] In the above formula, These parameters describe the variability of jitter, representing the observation noise variance. A smaller observation noise variance leads to faster Kalman filter convergence; however, if the value is too small, it can result in unstable images.

[0160] (39) Exercise compensation:

[0161] The motion vector (dx(k), d) of the derotated image in the current frame y (k)) Perform Kalman filtering, d x (k), d y (k) is the motion vector d x d y The time series and the filtered results are as follows: This refers to the subjective motion vector of the image:

[0162] The random jitter is calculated as follows:

[0163]

[0164]

[0165] The derotated image of the current frame is transformed based on the jitter amounts Δx(k) and Δy(k) to obtain a stable image sequence;

[0166]

[0167] Where P1 is the stable image, P0 is the derotated image of the current frame, and Δx(k) and Δy(k) are the random jitter amounts.

[0168] (40) Target detection based on difference map and track association: Target detection in stable images is performed using a moving target detection algorithm based on a large image coordinate system in the difference map. The specific steps include:

[0169] (41) Static background estimation:

[0170] Due to various factors such as lighting and scene changes, the background cannot remain unchanged after it is created. To enable the background model to adapt to environmental changes, the initial model must be continuously updated. The essence of background updating is to use the current frame to correct the background model established in the past. A common adaptive background update algorithm is as follows:

[0171] B n (x, y) = α·B n-1 (x, y) + (1-α)·I n (x, y)

[0172] Among them, B n(x, y) represents the updated background at the current moment; B n-1 (x, y) represents the background established in the previous moment; I n (x, y) represents the current frame; α is the weighting coefficient, which determines how fast the background is updated. When α is small, the background weight in the previous moment is small, and the weight of the current frame is large, indicating that the current frame image is trusted more. The background image is mainly updated based on the current frame image, and the background image is updated faster. Conversely, when α is large, the background image is updated slower.

[0173] For static background estimation, the static background image obtained from the previous frame is translated based on the image's subjective motion vector, and then accumulated with the current frame image. Therefore, to obtain an accurate large static background image, it is necessary to accurately calculate and extract the image's subjective motion vector, and then obtain the final large static background image based on the extracted image's subjective motion vector.

[0174] (42) Difference map background suppression:

[0175] Background subtraction is a method for detecting moving targets by subtracting the background image from each frame in a video sequence. Background subtraction is used for target segmentation, and the algorithm is efficient and easy to implement. However, in practical applications, sudden changes in lighting and interference from object movement can make target detection difficult. The formula for calculating the difference image is:

[0176] D k (x, y) = |I k (x, y)-B k (x, y)|

[0177] Among them, D k (x, y) is the difference graph, I k (x, y) represents the current frame image, B k (x, y) is the background image.

[0178] For a noise-free image, its signal includes a stationary background and moving targets. By accumulating images using flight tracks, the stationary background can be estimated. Theoretically, the difference image obtained by subtracting the stationary background image from the original image only contains the moving targets.

[0179] However, in practice, since the background cannot be completely removed, especially in the edge regions, the difference image still retains a significant amount of information. Therefore, it is necessary to suppress the background information in the difference image. In the detection and tracking of small targets, due to the dual diversity of target and background clutter, especially in complex detection backgrounds, background suppression techniques directly affect the performance of subsequent target detection. Since the difference image contains a large amount of edge information, edges can be extracted and then suppressed to obtain a difference image where the target signal is prominent.

[0180] Background suppression of the difference map is performed using a gradient-based method, and the calculation formula is as follows:

[0181]

[0182] Where a is the background suppression coefficient; B(i,j) is the background image; D(i,j) is the difference image; and th is the adjustment threshold.

[0183] A schematic diagram of difference plot suppression is shown below. Figure 13 As shown:

[0184] The left image is the background image, and the right image is the difference image after background suppression. It can be seen that after background edge suppression, the moving target signal is enhanced, while the background edge signal is weakened, effectively eliminating most of the background edge information in the image and reducing the false alarm rate of target detection. This algorithm can effectively suppress background-prominent targets, improve the image signal-to-noise ratio, and has low computational complexity, providing real-time and efficient processing capabilities.

[0185] (43) Curvature point extraction:

[0186] In different background environments, small targets typically exhibit a Gaussian-like distribution in the 3D response map of an image. This is reflected in the image's grayscale map, where the curvature value of the grayscale feature curve in any direction at the target point is relatively large. Therefore, regions with large curvature values ​​in different directions can be searched to find potential small targets. Curvature points are extracted in the post-suppression difference map.

[0187] According to the point target imaging model, an ideal point target will generally form a light spot on the imaging surface after passing through the optical system, and the intensity distribution of this light spot is usually approximated as a Gaussian distribution. Therefore, we can assume that there is a maximum or minimum value at the target point on the imaging surface, that is, the first derivative at the target point is zero. Thus, the formula for the curvature of the target point curve can be simplified to:

[0188] k=|y″|

[0189] The spatial curvature of a point on a discrete three-dimensional surface is characterized by calculating the curvature of the curves at a given pixel in four directions, with angles of 0°, 45°, 90°, and 135° to the x-axis.

[0190] The fitting coefficients can be obtained by template operations on image pixel values, and the formula can be expressed as follows:

[0191] K i = f(x+r, y+c)*W i

[0192] In the formula, "*" represents the convolution operator. The formula represents the coefficient K. i The acquisition can be further simplified to obtaining the image and the kernel function W. i Perform convolution operations. This is done by using the discrete orthogonal polynomial function {P}. i Substituting (r, c) into W i From the expression, we can obtain

[0193]

[0194]

[0195] With the kernel function described above, the fitting coefficient K can be obtained. i The combined formula yields the pixel (x) o y o The formula for the first derivative in the direction is:

[0196]

[0197] The formula for the second derivative in this direction is:

[0198]

[0199] To comprehensively characterize image information, this paper calculates the second derivatives of image pixels in four directions (0°, 45°, 90°, and 135°) and then sums these second derivatives in four directions using weighted summation to obtain the final curvature result. The three-dimensional curvature plot of the image can be represented as:

[0200] I c =a(I d 0 ·I d 90 )+b(I d 45 ·I d 135 )

[0201] In the formula, I c Let represent the curvature plot after weighted summation, where a and b are the weight values, and in this paper, we set a = b = 0.5.

[0202] Next, the obtained curvature map is used to replace the original grayscale image, and adaptive segmentation is performed on the curvature map. The segmentation threshold is as follows:

[0203] T=μ+kσ 2

[0204] In the formula, μ is the mean of the curvature diagram, k is a coefficient, and σ 2 Let T be the variance of the curvature map. After thresholding and binarization, pixels with values ​​greater than the threshold T in the resulting curvature map are initialized to 255, and pixels with values ​​less than the threshold T are initialized to 0. This will give us the curvature background suppression map of the image.

[0205] Curvature background suppression, such as Figure 14 As shown.

[0206] (44) Extraction of dot information:

[0207] A curvature point contains the following information: pixel coordinates in the x-direction, pixel coordinates in the y-direction, azimuth angle (az), pitch angle (el), grayscale value (gray), curvature value (cur), and sot value. Therefore, within the scope of small object detection, an image will be represented by the information contained in a series of such curvature points. Furthermore, in subsequent trajectory association, only the curvature point information is used, not the original image information. Let a curvature point i be denoted as cdot. i This is also the collection of information about this point, namely:

[0208] cdot i ={x i y i gray i ,cur i , sot i}

[0209] Assuming an image contains dotnum curvature points, then the set of valid point information cdot carried by the image can be represented as:

[0210] cdot={cdot1,cdot2,cdot3,...,cdot dotnum}

[0211] (45) Track correlation:

[0212] A flight path is composed of individual points. Matching flight paths with points essentially means matching points representing the same target at different times. Since angle values ​​already include pixel coordinates and provide a more accurate description of the target's position, and since grayscale values ​​are more susceptible to changes in camera settings or the environment, while curvature values ​​represent a relative value between the target's grayscale and the background grayscale, they are less affected by external factors. Therefore, this paper selects the azimuth angle (az), pitch angle (el), and curvature value (cur) to construct a descriptive vector. The more similar these three values ​​are for two points at different times, the higher the similarity between the two points, and the greater the probability that they represent the same target.

[0213] use Let represent the description vector of the j-th track at time t-1 on scale s, which can be expressed as:

[0214]

[0215] In the formula, This indicates the predicted gate center position at time t based on the flight path information. This represents the curvature value of the track associated with the point at time t-1.

[0216] use Let be the description vector representing the i-th point trace detected in the image at time t on scale s, which can be expressed as:

[0217]

[0218] In the formula, Cur represents the angular position of the i-th point detected in the image at time t on scale s. i (s，t) This represents the curvature value of the point.

[0219] The similarity between two descriptive vectors is assessed by calculating the cosine of the angle between them, which is called cosine similarity. The formula is as follows:

[0220]

[0221] In the formula, · represents the dot product operation, and ||·|| represents the modulus operation of the vector. This represents the similarity between the j-th track at scale s at time t and the i-th point detected in the image at time t at scale s. Since a track can be associated with at most one point at any given time, associating it with the point with the highest similarity is the optimal choice. Furthermore, for all points to be matched, the similarity must be greater than a certain threshold; only points with similarity greater than the threshold can be considered true target points. In summary, matching tracks and points within a gate requires two conditions: first, the similarity must be greater than a similarity threshold; second, the point with the highest similarity must be selected. This can be expressed mathematically as:

[0222]

[0223] In the formula, csim_thre represents the similarity threshold, and the value selected in this paper is 0.95. This means that similarity calculation is performed only on the dots within the gate range. This represents the maximum similarity between the j-th track at scale s at time t and the in-gate point traces detected in the image at scale s at time t. If the value is not 0, it means that the trajectory has achieved a match with the current point within the gate.

[0224] (46) Track clustering to extract moving targets:

[0225] By statistically analyzing the average speed of a large number of tracks at the same moment, the range of the average speed of the background tracks at this time is estimated, and motion clustering of the background tracks is achieved, initially screening out tracks with large differences in average speed.

[0226] Let t1 represent the current time. Taking the azimuth angle as an example, the average velocity of the j-th track on the s-scale during the time interval t0 to t1 can be expressed as:

[0227]

[0228] In the formula, Let represent the azimuth angle of the j-th track at time t on scale s. To ensure more accurate calculation of the track's average velocity, this paper stipulates that the track must have more than 10 associated points between t0 and t1. The total number of tracks satisfying this condition at time t1 is denoted as .

[0229] After obtaining the average velocity of all eligible tracks within the same time period, the range of the average velocity of background tracks is estimated based on the average velocity of a large number of tracks. This paper uses a fixed multiple of the mean of the average velocities of all eligible tracks plus or minus their standard deviation as the range of the velocity of background track points. Tracks whose average velocity falls within this range can be identified as background tracks; tracks outside the range may be moving target tracks or background tracks.

[0230] Therefore, the mean and standard deviation of the average speed of all tracks were further calculated. Similarly, to ensure the representativeness and reasonableness of this mean, the total number of tracks meeting the conditions was determined. It must be a large value. Let s be the mean of the average velocity of the conditional track within the time interval t0 to t1 on the s-scale. Let s be the standard deviation of the average velocity of the conditional track at scale s during the time interval t0 to t1:

[0231]

[0232]

[0233] In the formula, This represents the number of all tracks at time t1 on scale s. This represents the total number of tracks that meet the condition that the number of associated points is greater than 10 during the time period t0 to t1.

[0234] After obtaining the mean and standard deviation of the average velocity of the track at time t1, the estimated range of the average velocity of the background track can be expressed as:

[0235]

[0236] In the formula, k is a coefficient that determines the size of the range. When the coefficient is large, the calculated mean and standard deviation are more accurate. In this case, the coefficient k can be appropriately increased. When the coefficient k is small, the calculated mean and standard deviation will contain a large error. In this case, the coefficient k can be appropriately reduced.

[0237] Background trajectory average velocity diagram as follows Figure 15 As shown:

[0238] Tracks whose average speed falls within the estimated range of the average speed of background tracks are classified as background point tracks; those whose average speed falls outside the estimated range are classified as suspected moving target tracks. The upper limit of the average speed range for background tracks is V. global -Tu, lower limit of range V global -Tl and average speed of the track It can be represented as follows:

[0239]

[0240] The clustering expression for the tracks is:

[0241]

[0242] Here, Trace represents the feature point track. After segmenting background tracks and suspected target tracks through global motion clustering, due to potential errors and incorrect associations, a few background track points may be incorrectly classified as suspected target tracks. However, such incorrect track classifications are generally short-lived and accidental, rather than long-term and continuous. Therefore, the solution to this problem is to re-cluster all tracks in each frame. Only when a track is classified as a suspected target for multiple consecutive frames will it be ultimately classified as a moving target.

[0243] Flight path renderings and flight path speed diagrams are as follows: Figure 16 As shown:

[0244] In the trajectory rendering, the yellow asterisks represent the background curvature points at the current moment, and the green dots represent the background curvature points at historical moments. A yellow asterisk and a series of green dots form a background point trajectory. The red boxes and red dots represent detected moving targets and their trajectories. In the trajectory velocity graph, the vertical axis represents the average velocity of the trajectory, the blue circles represent the average velocity of the background trajectory, the green dots represent the threshold, and the red circles represent the average velocity of the moving target's trajectory.

Claims

1. A method for temporal filtering of images based on motion compensation, characterized in that, Includes the following steps: An autoregressive moving average model was used to model the obtained gyroscope error, and the model parameters were estimated by utilizing the data variation patterns to establish a gyroscope error model. Using the gyroscope error model to obtain error-removed gyroscope data, the camera rotation matrix is ​​calculated. The rotation matrix is ​​then used to perform reverse compensation on the rotated image to obtain a derotated image. Feature point vectors of the despin image after gyroscope rotation compensation are extracted based on an accelerated, robust feature extraction algorithm; the specific method is as follows: Constructing an integral image from the racemic image, specifically: The grayscale value of the pixel (x, y); The derotated image is filtered using three box filters (box1, box2, and box3) of different sizes. The specific formula is as follows: in, For the derotation image, There are three types of box filters, where m and n are the number of rows and columns of the derotation image, and II(x,y) is the value at the (x,y) coordinate position of the integral image. , , These are images filtered by three different box filters; The Hessian determinant image is calculated for filtered images of different sizes, and a multi-scale image pyramid is constructed using the Hessian determinant image. Specifically, the Hessian determinant image is: In the formula, Let be the determinant of the Hessian matrix; Strong feature points are located by performing extreme value detection on the middle two layers of each group of the image pyramid, and the main direction of the strong feature points is obtained by statistically analyzing the cumulative value of the Haar wavelet response of the strong feature points. The feature point vector is obtained by performing Haar wavelet response on strong feature points along the main direction. The feature point vectors are coarsely matched using the nearest neighbor and second nearest neighbor method. The coarsely matched point pairs are then finely matched using an improved nearest neighbor and second nearest neighbor matching algorithm. The motion vectors of the inter-frame images are estimated using the correctly matched feature points. The random jitter of the image is obtained by performing Kalman filtering on the motion vector sequence. Motion compensation is then performed on the derotated image based on the random jitter to obtain a stable image and the subjective motion vector of the image. A moving target detection algorithm based on a large graph coordinate system is used to detect targets in stable images.

2. The method for temporal filtering of images based on motion compensation according to claim 1, characterized in that, An autoregressive moving average model is used to model the obtained gyroscope error. The model parameters are estimated using the variation patterns of the data to establish the gyroscope error model. The specific method is as follows: Establish an autoregressive moving average model of order ARMA(2,2). It is the error data of the kth frame after first-order difference processing. It is the error data of the (k-1)th frame after first-order difference processing. It is the error data of the (k-2)th frame after first-order difference processing. It is white noise with a mean of 0 and a constant variance. , , , These are model parameters; The parameters of the autoregressive moving average model were estimated by Kalman filtering to obtain an ARMA(2,2) order autoregressive moving average model.

3. The method for temporal filtering of images based on motion compensation according to claim 1, characterized in that, The specific method for calculating the camera rotation matrix using error-removed gyroscope data obtained from the gyroscope error model is as follows: Let the rotation angle of the camera along the three axes between the k-th frame and the reference frame be a vector. Find the vector model : Where x is the camera's horizontal axis rotation angle, y is the camera's vertical axis rotation angle, and z is the camera's optical axis rotation angle. To normalize the length, the vector Divide by To obtain the normalized rotation vector The formula is as follows: Rewrite the normalized rotation vector as a skew-symmetric matrix: According to the Rodriguez rotation matrix principle, the rotation matrix D is obtained as follows: 。 4. The method for temporal filtering of images based on motion compensation according to claim 1, characterized in that, The rotated image is inversely compensated using a rotation matrix to obtain a derotated image. The inverse compensation method is as follows: Assume the coordinates of a point in the world coordinate system are The coordinates in the i-th frame image are The coordinates in the j-th frame are And from frame i to frame j, the camera only undergoes rotational motion, with the rotation matrix being... The pixel coordinate change formula during machine rotation is as follows: K is the intrinsic parameter matrix. These represent the camera's roll, pitch, and yaw angles, respectively. The formula for inverse compensation of a rotated image is: in, The coordinates of the derotated image in the j-th frame are .

5. The method for temporal filtering of images based on motion compensation according to claim 1, characterized in that, The specific method for coarse matching of feature point vectors using the nearest neighbor and second nearest neighbor method, followed by fine matching of the coarsely matched points using an improved nearest neighbor and second nearest neighbor matching algorithm, is as follows: The current frame derotated image is denoted as current frame image B, and the previous frame derotated image is denoted as reference image A; A bidirectional matching algorithm based on the nearest neighbor distance ratio is adopted. The feature point vector is used to perform coarse matching of feature points between image A and image B to obtain the coarse matching point set Ai and Bi of image AB. An improved random sampling consensus algorithm is used to perform fine matching on the coarse matching point sets Ai and Bi. The specific process is as follows: Calculate the evaluation function value between matching points and discard the evaluation values. exist[ , Points outside the range, among which Evaluation value The average value of , where th is the elimination threshold, and the formula for calculating the evaluation function value is: in This represents the relative difference between Ai, Bi, and the slope of each interest point in their respective images; This represents the average slope of Ai, Bi, and each interest point in their respective images.

6. The method for temporal filtering of images based on motion compensation according to claim 1, characterized in that, Motion vectors between frames are estimated using correctly matched feature points. Constructing the derotated image motion model is as follows: in It is a point in reference image A. It is the corresponding point in the current image B. It is a horizontal motion vector. It is a vertical motion vector; The motion model parameters, i.e., motion vectors, between frames are estimated by correctly matching point pairs. , : In the formula, Let A be the coordinates of the matching point in reference image A in the i-th correctly matched feature point pair. Let be the coordinates of the matching point B in the current image of the i-th correctly matched feature point pair.

7. The method for temporal filtering of images based on motion compensation according to claim 1, characterized in that, Kalman filtering is used to obtain random jitter in the image. Motion compensation is then performed on the derotated image based on the random jitter to obtain a stable image and the image's subjective motion vector. Motion vector of the derotated image in the current frame Perform Kalman filtering. Motion vector , The time series and the filtered results are as follows: This is the subjective motion vector of the image: The random jitter is calculated as follows: Based on the amount of shaking , Transform the derotated image of the current frame to obtain a stable image sequence; in, To stabilize the image, For the current frame derotation image, , This represents the random jitter.

8. The method for temporal filtering of images based on motion compensation according to claim 1, characterized in that, The specific method for detecting targets in stable images using a moving target detection algorithm based on a large image coordinate system in a difference map is as follows: The static background is estimated using subjective motion vectors, and the specific formula is as follows: in, This is the background image obtained after the update at the current moment; The background image created in the previous moment; For the current frame; These are the weighting coefficients; The difference image is obtained by subtracting the stable image of the current frame from the background image, and background suppression is performed on the difference image; Curvature point detection is performed on the difference map after background suppression. Based on the difference in motion characteristics between target and background points, the curvature difference map is mapped to the large image coordinate system for track association to filter out targets. The specific method is as follows: Calculate the difference image pixels after background suppression. The second derivatives in the four directions are calculated, and the weighted sum of the second derivatives in the four directions is used to obtain the final curvature result, specifically: In the formula, This represents the curvature plot after weighted summation, where a and b are the weight values, respectively. The obtained curvature map is used to replace the original grayscale image. The curvature map is then adaptively segmented and binarized to obtain the curvature background suppression map of the image. Point trace information is extracted from all points in the curvature background suppression map. Directional pixel coordinates Direction pixel coordinates, azimuth angle Pitch angle Grayscale value curvature value , Values ​​are used to represent; The points are associated with flight paths. The cosine similarity between the suspected points and the latest points on the established flight paths is calculated. The points with the highest similarity are then associated. The calculation formula is as follows: The first to represent The flight path is The time-time description vector is as follows: This indicates a prediction based on flight track information. At the center position of the time gate, This indicates that the flight path is in The curvature value is constantly correlated with the point trace; Indicates in The first time detected in the time image The description vector of each point trace is as follows: Clustering flight paths to extract moving targets; 。

Images