Method for compensating motion information by using gyroscope

Abstract

The invention relates to a method for compensating motion information by using a gyroscope, which comprises the following steps of: modeling a gyroscope error by using an autoregressive moving average model, estimating model parameters, and establishing a gyroscope error model; processing gyroscope data by using a gyroscope error model, removing errors, calculating a camera rotation matrix, and carrying out reverse compensation on an image to obtain an image without rotation influence, namely a racemization image; extracting a feature point vector of the de-rotation image by adopting an acceleration feature algorithm, and carrying out rough matching and fine matching in sequence so as to estimate a motion vector sequence of the inter-frame image; processing the image motion vector sequence through Kalman filtering to obtain random jitter of the image, and performing motion compensation on the racemization image according to the random jitter to obtain a stable image and a subjective motion vector of the image; and detecting a target in the stable image by using a moving target detection algorithm in a difference image based on a large image coordinate system. According to the invention, the error of the gyroscope can be better inhibited, and the cost is reduced on the premise of realizing high precision.

Description

Technical Field

[0001] This invention relates to a method for compensating motion information using a gyroscope, belonging to the field of moving target detection in motion scenes. Background Technology

[0002] Random video jitter can severely impact image quality, causing irregular changes in the target's position within the image. This significantly interferes with subsequent target identification and tracking. In weapon system target identification and tracking systems, platform jitter can prevent the system from accurately tracking targets, potentially leading to target loss or misjudgment, resulting in false alarms and missed detections, severely impacting the entire system.

[0003] Commonly used image stabilization technologies are divided into three categories: mechanical stabilization, optical stabilization, and electronic stabilization. Mechanical stabilization (electromechanical stabilization) uses a servo motor system and gyroscope sensors to compensate for the relative motion of the imaging system, thereby achieving image stabilization. This technology is typically used on large platforms, such as artillery and tanks, for correcting observation and aiming systems and fire control systems. Mechanical stabilization offers good real-time performance and stability, but its large platform size, complex structure, and high cost make it 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. This technology 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 image stabilization technology, electronic image stabilization technology integrates computer, microelectronics, and digital signal processing technologies, and has advantages such as directness, accuracy, and low cost, with broad development prospects.

[0004] Mainstream algorithms for moving target detection include statistical learning methods such as inter-frame difference, optical flow, and background modeling. Inter-frame difference uses the subtraction of adjacent frames to obtain a difference image, and then uses a threshold to select moving target information. Although this method is computationally simple and easy to implement, it often struggles to obtain the complete outline of the moving target, easily resulting in "double images" and "holes," leading to inaccurate target detection. Background subtraction first constructs a background model, then subtracts the current frame image from the background to detect moving targets. 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. One research focus is on finding efficient and practical background models and updating them. Optical flow segmentation uses the vector features of moving targets that change over time to detect moving regions in image sequences, performing particularly well when the camera is moving. However, this algorithm is complex, requiring significant computation and specialized hardware support, making it 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] The technical problem solved by the present invention is to overcome the above-mentioned shortcomings of the prior art and provide a method for compensating motion information using a gyroscope.

[0007] The technical solution of this invention is:

[0008] A method for compensating motion information using a gyroscope includes:

[0009] The gyroscope error is modeled using an autoregressive moving average model. The model parameters are estimated by analyzing the data variation patterns, thereby establishing the gyroscope error model.

[0010] The gyroscope data is processed using a gyroscope error model. After removing the error, the camera rotation matrix is ​​calculated. The camera rotation matrix is ​​then used to perform reverse compensation on the image to obtain an image that eliminates the effects of rotation, i.e., a derotated image.

[0011] The feature point vectors of the derotated image are extracted using a robust accelerated feature algorithm. The feature point vectors are coarsely matched using the nearest neighbor and second nearest neighbor method, and fine matching is performed using an improved matching algorithm, thereby estimating the motion vector sequence of the inter-frame image.

[0012] By processing the image motion vector sequence through Kalman filtering, the random jitter of the image is obtained, and motion compensation is performed on the derotated image based on the random jitter to obtain a stable image and the subjective motion vector of the image.

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

[0014] Furthermore, the autoregressive moving average model is used to model the gyroscope error. The model parameters are estimated by analyzing the data variation patterns, thereby establishing the gyroscope error model. Specifically:

[0015] (1.1) Establish an autoregressive moving average model of order A(2,2).

[0016] h(k)=m1h(k-1)+m2h(k-2)+n1ε(k-1)+n2ε(k-2)+ε(k)

[0017] Where h(k) is the error data of the k-th frame after first-order differencing, h(k-1) is the error data of the (k-1)-th frame after first-order differencing, and h(k-2) is the error data of the (k-2)-th frame after first-order differencing; ε(k) is white noise with a mean of 0 and a constant variance; m1, m2, n1, and n2 are model parameters; ε(k-1) and ε(k-2) are the white noise of the (k-1)-th and (k-2)-th frames, respectively.

[0018] (1.2) Parameter estimation of the autoregressive moving average model using Kalman filtering

[0019] Parameter estimation is performed using Kalman filtering, and the corresponding state variables are X = [m1(k), m2(k), n1(k), n2(k)] T The observation matrix is

[0020] Furthermore, the step of processing gyroscope data using a gyroscope error model, removing errors, and then calculating the camera rotation matrix specifically involves:

[0021] Let the rotation angle of the camera along the three axes between the k-th frame and the reference frame be vector θ. c =[xyz], find vector θ c The modulus δ:

[0022]

[0023] 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, denoted by vector θ. c Dividing by δ, we obtain the normalized rotation vector β = [β1 β2 β3], as shown in the following formula:

[0024]

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

[0026]

[0027] According to the Rodriguez rotation matrix principle, the rotation matrix D is obtained as follows:

[0028]

[0029] Furthermore, the method of using the camera rotation matrix to perform inverse compensation on the image to obtain an image that eliminates the effects of rotation, i.e., a derotated image, specifically involves:

[0030] Let the coordinates of a point in the world coordinate system be T. w The coordinates in the i-th frame are T i The coordinates in the j-th frame are T j According to the camera imaging model, we have:

[0031] T i =KRTw

[0032] Where K is the intrinsic parameter matrix and R is the rigid body transformation matrix;

[0033] Assuming the camera only undergoes rotational motion from frame i to frame j, and the rotation matrix is ​​Z(α,β,γ), then we have

[0034] T j =KZ(α,β,γ)RT w

[0035] =KZ(α,β,γ)R(KR)-1T i

[0036] =KZ(α,β,γ)K -1 T i

[0037] 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:

[0038] T j_new =KZ -1 (α,β,γ)K -1 T j

[0039] T j_new That is, the derotation image.

[0040] Furthermore, the method employs a robust accelerated feature extraction algorithm to extract feature point vectors from the derotation image, uses the nearest neighbor and second nearest neighbor method for coarse matching of feature point vectors, and utilizes an improved matching algorithm for fine matching to estimate the motion vector sequence of the inter-frame images. Specifically:

[0041] (3.1) Calculate the integral image based on the derotation image to obtain the image function;

[0042] (3.2) Feature points are detected by the Hessian matrix, which is composed of the second derivative of the Laplacian function of Gauss. The Laplacian function of Gauss is the result of the convolution of the image function and the Gaussian function. The filtered image is obtained by replacing the second derivative of the Gaussian function with a box filter.

[0043] (3.3) Establish scale space; the scale space is divided into four groups, each group includes four layers. Each group represents a series of response maps of the same input image filtered by progressively enlarged filter templates.

[0044] (3.4) Perform feature point localization. Perform non-maximum suppression on the three adjacent matrix images in each group of images. Perform non-maximum suppression on the 3×3×3 stereo neighborhood in the scale space. If the middle point is an extreme value, it is considered a feature point. At the same time, record the scale space and position of the point.

[0045] (3.5) Perform feature point vector description to obtain feature point vectors;

[0046] (3.6) Perform coarse matching of feature points: use feature point vectors to pair feature points, use Euclidean distance to measure the similarity of feature points, and obtain a set of matching point pairs;

[0047] (3.7) Perform fine matching of feature points to obtain a set of matching point pairs after removing mismatched points;

[0048] (3.8) Estimate the motion vector sequence based on the set of matching point pairs after removing mismatched points.

[0049] Furthermore, step (3.5) describes the feature point vectors, specifically as follows:

[0050] 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 ​​of the image in the x and y directions 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 greater weights for closer points. The cumulative Haar wavelet response values ​​are calculated within a sliding sector region with a rotation center of 60° to determine the principal direction of the feature point.

[0051] After determining the principal direction of the feature point, calculate the feature description vector within a 20α×20α rectangular region centered on the feature point along this principal direction. Divide this matrix into 4×4 sub-regions, and calculate the response values ​​in the x and y directions within each sub-region using a 2×2 Haar wavelet template, calculating ∑dx, ∑dy, ∑|dx|, and Σ|dy| respectively. Each sub-region is represented by a vector as follows:

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

[0053] A complete feature point description consists of vectors from 4×4 sub-regions, i.e., a 4×4×4 = 64-dimensional vector representation. This vector is then normalized to obtain the feature point vector.

[0054] Furthermore, step (3.7) involves fine-tuning feature point matching, specifically as follows:

[0055] Determine M in the reference image M i With all the feature points M in itself j The slope relationship, determining N in the image N to be matched i With all the feature points N in itself j Based on the slope relationship, and evaluating the correspondence between two points according to the similarity of the two slope relationships, the following evaluation function is proposed:

[0056]

[0057] K(i,j)=[k(M i N i )+k(M j N j )] / 2

[0058] Where I(i,j) represents M i N i The relative difference in slope between each interest point and the image; K(i,j) represents M i N i The average slope of each interest point in the respective image; removing the evaluation value T(i) in Points outside the range, among which Th is the average value of the evaluation value T(i), and th is the elimination threshold.

[0059] k(Mi,Ni) represents the slope between feature point Mi of the reference image and feature point Ni of the image to be matched;

[0060] k(Mj,Nj) represents the slope between feature point Mj of the reference image and feature point Nj of the image to be matched;

[0061] k(Mi,Mj) represents the slope between feature points Mi of the reference image M and feature points Mj of the image to be matched.

[0062] k(Ni,Nj) represents the slope between feature points Ni of the reference image and feature points Nj of the image to be matched in the reference image N.

[0063] Furthermore, step (3.8) estimates the motion vector sequence, specifically as follows:

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

[0065]

[0066] Where q0 is a point in the reference image M, q1 is the corresponding point in the current image N, and d x It is the horizontal motion vector, d y It is a vertical motion vector;

[0067] By removing incorrectly matched points from the set of correctly matched points, the motion model parameters between frames, i.e., the motion vector d, are estimated. x d y :

[0068]

[0069] In the formula, (x i ,y i (X) represents the coordinates of the matching point in the reference image M in the i-th correctly matched feature point pair, where (X) i ,Y i ) represents the coordinates of the N matching points in the current image of the i-th correctly matched feature point pair.

[0070] Furthermore, the process of processing the image motion vector sequence through Kalman filtering to obtain random jitter in the image, and then performing motion compensation on the derotated image based on the random jitter to obtain a stable image and the image's subjective motion vector, specifically involves:

[0071] (4.1) The motion vector (d) of the derotated image in the current frame x (k), d 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: That is, the subjective motion vector of the image;

[0072] (4.2) Calculate the random jitter amount and perform motion compensation, specifically as follows:

[0073]

[0074] 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;

[0075]

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

[0077] Furthermore, the method of using a moving target detection algorithm based on a large-image coordinate system in a difference map to detect targets in a stable image specifically involves:

[0078] (5.1) Perform static background estimation: Based on the subjective motion vector of the image, perform a translation operation on the static background image obtained in the previous frame, and then accumulate it with the current frame image to obtain an updated static background image.

[0079] (5.2) A gradient-based method is used for background suppression of the difference map;

[0080]

[0081] Where a is the background suppression coefficient; T(i,j) is the static background image; C(i,j) is the difference map; th is the adjustment threshold; G(i,j) represents the gradient value at coordinate (i,j) in the image; and a(th) is the background suppression coefficient related to the threshold th.

[0082] (5.3) Extract curvature points from the difference map after background suppression;

[0083] (5.4) Extract the dot information;

[0084] Let t denote a curvature point i as t i This is also the collection of information about this point, namely:

[0085] t i ={x i ,y i gray i ,cur i ,sot i}

[0086] x i x-axis pixel coordinates, y-axis pixel coordinates i Here, y represents the pixel coordinate in the y-direction, gray represents the grayscale value, cur represents the curvature value, and sot represents the target class index.

[0087] Suppose an image contains d curvature points, then the set t of valid point information carried by the image is represented as:

[0088] t = {t1, t2, t3, ..., t d};

[0089] (5.5) Perform track association;

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

[0091]

[0092] 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;

[0093] use Let be the description vector of the i-th point detected in the image at time t on scale s, and it is expressed as:

[0094]

[0095] 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 that point;

[0096] The similarity between the two descriptive vectors is assessed by calculating the cosine of the angle between them, i.e., calculating their cosine similarity. The formula is as follows:

[0097]

[0098] 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 track detected in the image at scale s at time t.

[0099] In the process of matching tracks and points within a gate, two conditions must be met simultaneously: first, the similarity must be greater than a similarity threshold; second, the point with the highest similarity must be selected for association. This can be expressed mathematically as follows:

[0100]

[0101] In the formula, ch represents the similarity threshold. This means that similarity calculation is performed only on the dots within the gate area; 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 track has achieved a match with the current point within the gate;

[0102] (5.6) Track clustering to extract moving targets;

[0103] The clustering expression for the tracks is:

[0104]

[0105] Where T represents the feature point track. In each frame, all tracks are re-clustered and a decision is made. Only when a track is identified as a suspected target for multiple consecutive frames is it finally identified as a moving target.

[0106] Among them, the upper limit of the average speed range of the background track is V. g _Tu, lower limit of range is V g _Tl, average speed of the trajectory is

[0107] The advantages of this invention compared to the prior art are:

[0108] (1) This invention has low requirements for gyroscope accuracy. The model parameters can be estimated by using the data variation patterns through software, and a corresponding mathematical model can be established, thereby better suppressing gyroscope errors. Under the premise of achieving high accuracy, the cost is reduced because an excessively high-precision gyroscope is not required, thus saving costs.

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

[0110] (3) The present invention has strong target detection capability, which can effectively suppress the interference of background on the target, improve the image signal-to-noise ratio, and make it possible to screen out moving targets in complex moving backgrounds.

[0111] (4) The present invention has a high sensitivity detection capability, which can quickly detect the target and its trajectory through the correlation of a small number of frames, thereby improving the efficiency and accuracy of target detection. Attached Figure Description

[0112] Figure 1 This is a flowchart of a hardware method for compensating motion information using a gyroscope.

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

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

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

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

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

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

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

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

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

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

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

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

[0125] Figure 14 Curvature background suppression.

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

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

[0128] The specific embodiments of the present invention will now be described in further detail with reference to the accompanying drawings.

[0129] like Figure 1 As shown, this invention proposes a hardware method for compensating motion information using a gyroscope, comprising the following steps:

[0130] Step 1: Gyroscope Error Analysis and Modeling: An autoregressive moving average model is used to model the obtained gyroscope error. The model parameters are estimated using the data variation patterns to establish the gyroscope error model. The specific steps are as follows:

[0131] (1) Establish an autoregressive moving average model of order A(2,2).

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

[0133]

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

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

[0136] h(k)=m1h(k-1)+m2h(k-2)+n1ε(k-1)+n2ε(k-2)+ε(k)

[0137] h(k) is the error data of the kth frame after first-order differencing, z(k-1) is the error data of the (k-1)th frame after first-order differencing, h(k-2) is the error data of the (k-2)th frame after first-order differencing, ε(k) is white noise with a mean of 0 and a constant variance, and m1, m2, n1, and n2 are model parameters.

[0138] (2) Parameter estimation of the autoregressive moving average model using Kalman filtering

[0139] Parameter estimation is performed using Kalman filtering, and the corresponding state variables are X = [m1(k), m2(k), n1(k), n2(k)] T The observation matrix is

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

[0141] Use MATLAB's built-in functions to perform A(2,2) modeling on the collected data.

[0142] h(k) = 0.205h(k-1) - 0.106h(k-2) - 0.635ε(k-1) - 0.715ε(k-2) + ε(k) The autocorrelation and partial correlation of the residuals are as follows: Figure 3 As shown:

[0143] 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.

[0144]

[0145] The A(2,2) model obtained using Kalman filtering is:

[0146] h(k)=0.197h(k-1)-0.082h(k-2)-0.642ε(k-1)-0.731ε(k-2)+ε(k)

[0147] Step 2: 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. The rotation matrix is ​​then used to perform inverse compensation on the rotated image to obtain a derotated image. The specific steps are as follows:

[0148] (1) Calculate the rotation matrix

[0149] By integrating the angular velocity output by the gyroscope, the rotation angle around the three axes can be obtained. Let ω i To measure the angular velocity, the gyroscope sampling is discrete, and Δt is the gyroscope's sampling period. Therefore, the rotation angle during the time interval from t1 to t2 is:

[0150] 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 δ:

[0151]

[0152] 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 β = [β1 β2 β3], as shown in the following formula:

[0153]

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

[0155]

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

[0157]

[0158] (2) Reverse compensation

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

[0160] T i =KRT w

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

[0162] Assuming the camera only undergoes rotational motion from frame i to frame j, and the rotation matrix is ​​Z(α,β,γ), then we have

[0163] T j =KZ(α,β,γ)RT w

[0164] Substituting, we get

[0165] T j =KZ(α,β,γ)RT w

[0166] =KZ(α,β,γ)R(KR)-1T i

[0167] =KZ(α,β,γ)K-1T i

[0168] 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:

[0169] T j_new =KZ -1 (α,β,γ)K -1 T j

[0170] T j_new The original image and the derotated image are shown below. Figure 4 As shown.

[0171] Step 3: Translation estimation and compensation based on image matching: A robust accelerated feature extraction algorithm is used to extract feature point vectors from the derotated image. Coarse matching of the feature point vectors is performed using the nearest neighbor / second nearest neighbor method, followed by fine matching using an improved matching algorithm, thereby estimating the motion vector sequence between frames. The specific steps are as follows:

[0172] (3.1) Find the integral graph:

[0173] To calculate the sum of pixel gray levels in a region of an image, the gray levels of all pixels within that region need to be accumulated. If the sum of gray levels of sub-blocks in an image is frequently calculated, an integral image can be used for quick computation. An integral image is defined as follows:

[0174]

[0175] In the above formula, σ(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:

[0176]

[0177] In the formula, Iσ(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.

[0178] (3.2) Feature point detection based on the Hessian matrix:

[0179] 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:

[0180] L(x,y,α)=G(x,y,α)*σ(x,y)

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

[0182]

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

[0184]

[0185] In the formula, M xx (x,y,α),M xy (x,y,α),M yy (x,y,α) is the second derivative of the function M(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.

[0186] 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:

[0187] The derotated image is filtered using three different box filters (b1, b2, b3) of varying sizes to obtain the filtered image. The specific method is as follows:

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

[0189]

[0190] Where G is the derotation image, b1, b2, and b3 are three box filters, m is the number of rows in the derotation image, n is the number of columns in the derotation image, Iσ(x,y) is the value at the (x,y) coordinate position of the integral image, and Z... xx Z yy Z xy These are images filtered by three different box filters;

[0191] (3.3) Establishment of scale space:

[0192] 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.

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

[0194]

[0195] 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.

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

[0197] (3.4) Feature point localization:

[0198] Each scale-space box filter has three filter templates, b1, b2, and b3. Each template is convolved with the image to obtain the corresponding Z-value. xx Z yy Z xy The values ​​are then used to generate a determinant image, and the calculation process is as follows: Figure 8 As shown.

[0199] 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.

[0200] 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:

[0201]

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

[0203] det(H approx ) = Z xx Z yy -(λZ xy ) 2

[0204] 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.

[0205] (3.5) Feature point vector description:

[0206] 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 responses of the image in the x and y directions 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 is 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.

[0207] 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.

[0208] 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:

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

[0210] 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.

[0211] (3.6) Coarse matching of feature points:

[0212] 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.

[0213] 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:

[0214]

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

[0216] Suppose there are images M and N, P M Let P be the set of feature points in image M. N Let M be the set of feature points in image N. i For P M A feature vector, N m N n For P N In and M 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 N < T, then N 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 M to image N.

[0217] Perform the same operation on all feature points in image N to obtain the matching point pair set D from image N to image M. The intersection of point pair set G and point pair set D is the final matching point pair set.

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

[0219] (3.7) Fine-grained feature point matching:

[0220] 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:

[0221]

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

[0223]

[0224] As stated above, at least eight equations, or four non-collinear point pairs, are needed to determine the eight parameters in the projection transformation matrix P. Since there are mismatched point pairs, a specific method is required to eliminate them.

[0225] The random sample consensus algorithm is used to refine the correspondences of feature points, thereby reducing the error of the transformation matrix. This algorithm assumes a set of data points is a correct set of inliers. It fits model parameters to this set of data and then determines how many data points within a certain error threshold satisfy the model. Those satisfying the threshold are called inliers; otherwise, they are called outliers. The more inliers, the more data points satisfy the model, and the higher the model's reliability. The algorithm terminates under two conditions: the maximum number of iterations and the maximum number of inliers. When the number of inliers exceeds the maximum number or the number of iterations reaches its maximum, the algorithm stops. At this point, the set of inliers with the maximum number of inliers is considered the correct set of data points. Taking projection transformation as an example, the steps of the random sample consensus algorithm are as follows:

[0226] 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 P.

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

[0228] 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.

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

[0230] 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.

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

[0232]

[0233] K(i,j)=[k(M i N i )+k(M j N j )] / 2

[0234] Where I(i,j) represents M i N i The relative difference in slope between each interest point and the image; K(i,j) represents M i N 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.

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

[0236] (3.8) Estimating the motion vector sequence:

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

[0238]

[0239] Where q0 is a point in the reference image M, q1 is the corresponding point in the current image N, and d x It is the horizontal motion vector, d y It is a vertical motion vector;

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

[0241]

[0242] In the formula, (x i ,y i (X) represents the coordinates of the matching point in the reference image M in the i-th correctly matched feature point pair, where (X) i ,Y i ) represents the coordinates of the N matching points in the current image of the i-th correctly matched feature point pair.

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

[0244]

[0245] In the formula, x is the motion vector [x,y], e is the change in x, and P(0,μ) is Gaussian white noise.

[0246] Expanding the above equation, we get:

[0247]

[0248] In the above formula, P.0, μ x / 、P.0,μ y / represents independent Gaussian white noise. μ 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.

[0249] The system measurement equation is:

[0250]

[0251] 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.

[0252] Step 4: Process the image motion vector sequence through Kalman filtering to obtain the random jitter of the image, and perform motion compensation on the derotated image based on the random jitter to obtain the stable image and the subjective motion vector of the image.

[0253] (4.1) The motion vector (d) of the derotated image in the current frame x (k), d 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:

[0254] (4.2) Perform motion compensation and calculate the random jitter amount, specifically:

[0255]

[0256] 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;

[0257]

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

[0259] Step 5: 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-map coordinate system in the difference map. Specific steps include:

[0260] (5.1) Static background estimation:

[0261] Due to 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:

[0262] T n (x,y)=γ·T n-1 (x,y)+(1-γ)·I n (x,y)

[0263] Among them, T n (x,y) represents the updated background at the current moment; T 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.

[0264] For static background estimation, we need to perform a translation operation on the static background image obtained in the previous frame based on the image's subjective motion vectors, and then accumulate it with the current frame image to obtain an updated static background image. Therefore, to obtain an accurate static background image, it is necessary to accurately calculate and extract the image's subjective motion vectors, and then generate the final static background image based on these extracted image subjective motion vectors.

[0265] (5.2) Differential plot background suppression:

[0266] Background subtraction is a method for detecting moving targets. It works by subtracting each frame of a video sequence from a currently constructed background image. This method is commonly used for object segmentation due to its high algorithm efficiency and ease of implementation. However, in practical applications, sudden changes in lighting and other factors such as the movement of interfering objects can complicate target detection. The formula for calculating the difference map is:

[0267] C k (x,y)=|I k (x,y)-T k (x,y)|

[0268] Among them, C k (x,y) is the difference graph, I k (x,y) represents the current frame image, T k (x,y) is the background image.

[0269] For a noise-free image containing both a stationary background and a moving target, the stationary background can be estimated by accumulating the image using tracks (trajectories). Theoretically, by subtracting the stationary background image from the original image, the resulting difference image should contain only information about the moving target.

[0270] In practice, background information is difficult to completely remove, especially in the edge regions, leaving a significant amount of information in the difference image. Therefore, background suppression is necessary. In small target detection and tracking, the diversity of both target and background clutter, particularly in complex detection environments, directly impacts the performance of subsequent target detection. To address the presence of abundant edge information in the difference image, edge extraction followed by edge suppression can yield a difference image that highlights the target signal.

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

[0272]

[0273] Where a is the background suppression coefficient; T(i,j) is the background image; C(i,j) is the difference image; and th is the adjustment threshold.

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

[0275] 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 becomes weaker. This effectively eliminates most of the background edge information in the image, thereby reducing the false alarm rate of target detection. This algorithm can effectively suppress the background, highlight the target, improve the signal-to-noise ratio of the image, and has real-time and efficient processing capabilities due to its low computational complexity.

[0276] (5.3) Curvature point extraction:

[0277] In different background environments, small targets typically exhibit a Gaussian-like distribution in the 3D response map of an image. Reflected in the image's grayscale, the curvature values ​​of the grayscale feature curves at the target point are relatively large in any direction. Therefore, potential small targets can be detected by searching for regions with large curvature values ​​in different directions. In the difference map after background suppression, curvature points—points with large curvature in the image—can be extracted. These points may correspond to target portions in the image, thus aiding in further localization and identification of small targets.

[0278] According to the point target imaging model, an ideal point target forms a light spot on the imaging surface after passing through the optical system, and the intensity distribution of this light spot can usually be 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:

[0279] k = |y”|

[0280] 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.

[0281] The fitting coefficients can be obtained by template operations on image pixel values, and can then be expressed as:

[0282] T i =f(x+u,y+v)*W i

[0283] In the formula, "*" represents the convolution operator. It indicates the coefficient T. 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 {Q...} i (u,v)} into W i From the expression, we can obtain

[0284]

[0285]

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

[0287]

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

[0289]

[0290] 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:

[0291] I c =m(I d 0 ·I d 90 )+n(I d 45 ·I d 135 )

[0292] In the formula, I c Let m represent the curvature plot after weighted summation, where m and n are the weight values. In this paper, we set m = n = 0.5.

[0293] 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:

[0294] U=λ+aσ 2

[0295] In the formula, λ is the mean of the curvature diagram, a is a coefficient, and σ is... 2 Let represent the variance of the curvature map. After thresholding and binarization, pixels with values ​​greater than the threshold U in the resulting curvature map are initialized to 255, and pixels with values ​​less than the threshold U are initialized to 0. This yields the curvature background suppression map of the image.

[0296] Curvature background suppression, such as Figure 14 As shown:

[0297] (5.4) Extraction of dot information:

[0298] 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 t denote a curvature point i. i This is also the collection of information about this point, namely:

[0299] t i ={x i ,y i gray i ,cur i ,sot i}

[0300] Assuming an image contains d curvature points, the set t of valid point information carried by the image can be represented as:

[0301] t = {t1, t2, t3, ..., t d}

[0302] (5.5) Track correlation:

[0303] 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.

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

[0305]

[0306] 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.

[0307] 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:

[0308]

[0309] 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.

[0310] 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:

[0311]

[0312] 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. In the field of target tracking, a track can be associated with at most one point at any given time; therefore, associating with the point with the highest similarity is considered the optimal choice. Furthermore, for all points to be matched, a similarity threshold must be met; only points with similarity greater than the threshold can be considered true target points. In summary, the matching process between tracks and points within the gate requires two conditions to be met simultaneously: first, the similarity must be greater than the similarity threshold; second, the point with the highest similarity must be selected for association. This can be expressed mathematically as:

[0313]

[0314] In the formula, ch 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.

[0315] (5.6) Track clustering to extract moving targets:

[0316] 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.

[0317] 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.

[0318] 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:

[0319]

[0320] 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.

[0321] 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:

[0322]

[0323] 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.

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

[0325] 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. g _Tu、Lower limit of range V g _Tl and average speed of the track It can be represented as follows:

[0326]

[0327] The clustering expression for the tracks is:

[0328]

[0329] Here, T represents the feature point track. After global motion clustering segments background tracks and suspected target tracks, due to potential errors and incorrect associations, a few background track points may be misclassified as suspected target tracks. However, such misclassifications are generally short-lived and accidental, rather than long-term and continuous. Therefore, a re-clustering decision is made for all tracks in each frame, and a track is only ultimately classified as a moving target if it has been identified as a suspected target for multiple consecutive frames.

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

[0331] 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.

[0332] The parts of this invention not described in detail are common knowledge to those skilled in the art.

Claims

1. A method for compensating motion information using a gyroscope, characterized in that... include: The gyroscope error is modeled using an autoregressive moving average model. The model parameters are estimated by analyzing the data variation patterns, thereby establishing the gyroscope error model. The gyroscope data is processed using a gyroscope error model. After removing the error, the camera rotation matrix is ​​calculated. The camera rotation matrix is ​​then used to perform reverse compensation on the image to obtain an image that eliminates the effects of rotation, i.e., a derotated image. The feature point vectors of the derotated image are extracted using a robust accelerated feature algorithm. The feature point vectors are coarsely matched using the nearest neighbor and second nearest neighbor method, and fine matching is performed using an improved matching algorithm, thereby estimating the motion vector sequence of the inter-frame image. By processing the image motion vector sequence through Kalman filtering, the random jitter of the image is obtained, and motion compensation is 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 compensating motion information using a gyroscope according to claim 1, characterized in that: The method uses an autoregressive moving average model to model the gyroscope error, estimating model parameters by analyzing data variation patterns, thereby establishing the gyroscope error model. Specifically: (1.1) Establish an autoregressive moving average model of order A(2,2). h(k)=m1h(k-1)+m2h(k-2)+n1ε(k-1)+n2ε(k-2)+ε(k) Where h(k) is the error data of the k-th frame after first-order differencing, h(k-1) is the error data of the (k-1)-th frame after first-order differencing, and h(k-2) is the error data of the (k-2)-th frame after first-order differencing; ε(k) is white noise with a mean of 0 and a constant variance; m1, m2, n1, and n2 are model parameters; ε(k-1) and ε(k-2) are the white noise of the (k-1)-th and (k-2)-th frames, respectively. (1.2) Parameter estimation of the autoregressive moving average model using Kalman filtering Parameter estimation is performed using Kalman filtering, with the corresponding state variables being X = [m1(k), m2(k), n1(k), n2(k)]T and the observation matrix being...

3. The method for compensating motion information using a gyroscope according to claim 1, characterized in that: The process of processing gyroscope data using a gyroscope error model, removing errors, and then calculating the camera rotation matrix is ​​as follows: Let the rotation angle of the camera along the three axes between the k-th frame and the reference frame be vector θ. c =[xyz], find vector θ c The modulus δ: 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, denoted by vector θ. c Dividing by δ, we obtain the normalized rotation vector β = [β1 β2 β3], as shown in the following formula: 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. A method for compensating motion information using a gyroscope according to claim 3, characterized in that: The process of using the camera rotation matrix to perform inverse compensation on the image to obtain an image free from the effects of rotation, i.e., a derotated image, is as follows: Let the coordinates of a point in the world coordinate system be T. w The coordinates in the i-th frame are T i The coordinates in the j-th frame are T j According to the camera imaging model, we have: T i =KRT w Where K is the intrinsic parameter matrix and R is the rigid body transformation matrix; Assuming the camera only undergoes rotational motion from frame i to frame j, and the rotation matrix is ​​Z(α,β,γ), then we have T j =KZ(α,β,γ)RT w =KZ(α,β,γ)R(KR) -1 T i =KZ(α,β,γ)K -1 T i 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: T j_new =KZ -1 (a,b,c)K -1 T j T j_new That is, the derotation image.

5. A method for compensating motion information using a gyroscope according to claim 3, characterized in that: The method employs a robust accelerated feature extraction algorithm to extract feature point vectors from the derotated image, uses the nearest neighbor and second nearest neighbor method for coarse matching of feature point vectors, and then uses an improved matching algorithm for fine matching to estimate the motion vector sequence of the inter-frame image. Specifically: (3.1) Calculate the integral image based on the derotation image to obtain the image function; (3.2) Feature points are detected by the Hessian matrix, which is composed of the second derivative of the Laplacian function of Gauss. The Laplacian function of Gauss is the result of the convolution of the image function and the Gaussian function. The filtered image is obtained by replacing the second derivative of the Gaussian function with a box filter. (3.3) Establish scale space; the scale space is divided into four groups, each group includes four layers. Each group represents a series of response maps of the same input image filtered by progressively enlarged filter templates. (3.4) Perform feature point localization. Perform non-maximum suppression on the three adjacent matrix images in each group of images. Perform non-maximum suppression on the 3×3×3 three-dimensional neighborhood in the scale space. If the middle point is an extreme value, it is considered a feature point. At the same time, record the scale space and position of the point. (3.5) Perform feature point vector description to obtain feature point vectors; (3.6) Perform coarse matching of feature points: use feature point vectors to pair feature points, use Euclidean distance to measure the similarity of feature points, and obtain a set of matching point pairs; (3.7) Perform fine matching of feature points to obtain a set of matching point pairs after removing mismatched points; (3.8) Estimate the motion vector sequence based on the set of matching point pairs after removing mismatched points.

6. A method for compensating motion information using a gyroscope according to claim 5, characterized in that: Step (3.5) describes the feature point vectors, specifically as follows: 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 ​​of the image in the x and y directions 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 greater weights for closer points. The cumulative Haar wavelet response values ​​are calculated within a sliding sector region with a rotation center of 60° to determine the principal direction of the feature point. After determining the principal direction of the feature point, calculate the feature description vector within a 20α×20α rectangular region centered on the feature point along this principal direction. Divide this matrix into 4×4 sub-regions, and calculate the response values ​​in the x and y directions within each sub-region using a 2×2 Haar wavelet template, calculating ∑dx, ∑dy, ∑|dx|, and ∑|dy| respectively. Each sub-region is represented by a vector as follows: v=[∑dx,∑|dx|,∑dy,∑|dy|] A complete feature point description consists of vectors from 4×4 sub-regions, i.e., a 4×4×4 = 64-dimensional vector representation. This vector is then normalized to obtain the feature point vector.

7. A method for compensating motion information using a gyroscope according to claim 5, characterized in that: Step (3.7) involves fine-tuning feature points, specifically as follows: Determine M in the reference image M i With all the feature points M in itself j The slope relationship, determining N in the image N to be matched i With all the feature points N in itself j Based on the slope relationship, and evaluating the correspondence between two points according to the similarity of the two slope relationships, the following evaluation function is proposed: K(i,j)=[k(M i ,N i )+k(M j ,N j )] / 2 Where I(i,j) represents M i N i The relative difference in slope between each interest point and the image; K(i,j) represents M i N i The average slope of each interest point in the respective image; removing the evaluation value T(i) in Points outside the range, among which Th is the average value of the evaluation value T(i), and th is the elimination threshold. k(Mi,Ni) represents the slope between feature point Mi of the reference image and feature point Ni of the image to be matched; k(Mj,Nj) represents the slope between feature point Mj of the reference image and feature point Nj of the image to be matched; k(Mi,Mj) represents the slope between feature points Mi of the reference image M and feature points Mj of the image to be matched. k(Ni,Nj) represents the slope between feature points Ni of the reference image and feature points Nj of the image to be matched in the reference image N.

8. A method for compensating motion information using a gyroscope according to claim 5, characterized in that: Step (3.8) estimates the motion vector sequence, specifically as follows: Constructing the derotated image motion model is as follows: Where q0 is a point in the reference image M, q1 is the corresponding point in the current image N, and d x It is the horizontal motion vector, d y It is a vertical motion vector; By removing incorrectly matched points from the set of correctly matched points, the motion model parameters between frames, i.e., the motion vector d, are estimated. x d y : In the formula, (x i ,y i (X) represents the coordinates of the matching point in the reference image M in the i-th correctly matched feature point pair, where (X) i ,Y i ) represents the coordinates of the N matching points in the current image of the i-th correctly matched feature point pair.

9. A method for compensating motion information using a gyroscope according to claim 5, characterized in that: The process involves processing the image motion vector sequence using Kalman filtering to obtain random jitter in the image, and then performing motion compensation on the derotated image based on this random jitter to obtain a stable image and the image's subjective motion vector. Specifically: (4.1) The motion vector (d) of the derotated image in the current frame x (k), d 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: That is, the subjective motion vector of the image; (4.2) Calculate the random jitter amount and perform motion compensation, specifically as follows: 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; Where Q1 is the stable image, Q0 is the derotated image of the current frame, and Δx(k) and Δy(k) are the random jitter amounts.

10. A method for compensating motion information using a gyroscope according to claim 5, characterized in that: The method of using a moving target detection algorithm based on a large-image coordinate system in a difference map to detect targets in a stable image specifically involves: (5.1) Perform static background estimation: Based on the subjective motion vector of the image, perform a translation operation on the static background image obtained in the previous frame, and then accumulate it with the current frame image to obtain an updated static background image. (5.2) A gradient-based method is used for background suppression of the difference map; Where a is the background suppression coefficient; T(i,j) is the static background image; C(i,j) is the difference map; th is the adjustment threshold; G(i,j) represents the gradient value at coordinate (i,j) in the image; and a(th) is the background suppression coefficient related to the threshold th. (5.3) Extract curvature points from the difference map after background suppression; (5.4) Extract the dot information; Let t denote a curvature point i as t i This is also the collection of information about this point, namely: t i ={x i ,y i ,gray i ,cur i ,sot i } x i x-axis pixel coordinates, y-axis pixel coordinates i Here, y represents the pixel coordinate in the y-direction, gray represents the grayscale value, cur represents the curvature value, and sot represents the target class index. Suppose an image contains d curvature points, then the set t of valid point information carried by the image is represented as: t={t1,t2,t3,...,t d }; (5.5) Perform track association; use Let represent the description vector of the j-th track at time t-1 on scale s, which is expressed as: 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; use Let be the description vector of the i-th point detected in the image at time t on scale s, and it is expressed as: In the formula, This represents the angular position of the i-th point detected in the image at time t on scale s. This represents the curvature value of that point; The similarity between the two descriptive vectors is assessed by calculating the cosine of the angle between them, i.e., calculating their cosine similarity. The formula is as follows: 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 track detected in the image at scale s at time t. In the process of matching tracks and points within a gate, two conditions must be met simultaneously: first, the similarity must be greater than a similarity threshold; second, the point with the highest similarity must be selected for association. This can be expressed mathematically as follows: In the formula, ch represents the similarity threshold. This means that similarity calculation is performed only on the dots within the gate area; 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 track has achieved a match with the current point within the gate; (5.6) Track clustering to extract moving targets; The clustering expression for the tracks is: Where T represents the feature point track. In each frame, all tracks are re-clustered and a decision is made. Only when a track is identified as a suspected target for multiple consecutive frames is it finally identified as a moving target. Among them, the upper limit of the average speed range of the background track is V. g _Tu, lower limit of range is V g _Tl, average speed of the track is

Images