Color adaptive equalization method for arbitrary orientation seam based on spatiotemporal kalman filter

CN123033127APending Publication Date: 2026-11-03自然资源部第六地形测量队
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202611548140.1
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2026-10-08
Publication Date
2026-11-03

AI Technical Summary

Technical Problem

Existing technologies struggle to achieve smooth parameter estimation and uniform adaptation of arbitrary-direction polygonal seams along the seam direction in remote sensing image stitching, resulting in color patches and segment boundaries. Furthermore, the computational load is high, making it difficult to meet the engineering production needs of large-format remote sensing images.

Method used

A spatiotemporal Kalman filter-based approach is adopted. Through a two-stage architecture of global linear pre-alignment and local smoothing correction, Kalman filtering is used to perform forward and reverse recursive smoothing of local correction parameters along the arc length of the joint. Combined with unified arc length parameterization and anisotropic spatial weight interpolation, noise-resistant smoothing estimation and unified parameter adaptation along the joint direction are achieved.

Benefits of technology

It effectively suppresses random noise and parameter fluctuations caused by abnormal features in the local statistical process, improves the visual consistency of color balance and parameter stability at the seams, reduces computational overhead, adapts to arbitrary orientations and complex seam shapes, and meets the requirements of large-format remote sensing image engineering production.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN123033127A_ABST
    Figure CN123033127A_ABST
Patent Text Reader

Abstract

The application discloses a color adaptive equalization method for arbitrary trend joint based on space-time Kalman filtering, and belongs to the field of remote sensing image processing. In order to solve the problem that the existing joint color uniform method is difficult to adapt to arbitrary trend fold line joint and the parameter fluctuation caused by local statistical noise, the application comprises the following steps: globally pre-aligning buffer pixels; adopting a unified arc length parameterization setting control point along the joint and mapping into a one-dimensional sequence; calculating local correction parameters of each point as an observation sequence, and establishing a state space model, and recursively smoothing through bidirectional Kalman filtering; and interpolating the smoothed parameters to each pixel in the buffer area by using anisotropic spatial weight, and combining a hybrid coefficient controlled by the distance attenuation along the joint normal to control the correction strength. The application is used for color radiation consistency correction of remote sensing image mosaic splicing, and has the advantages of noise resistance smoothing, adaptation to arbitrary fold line joint and high calculation efficiency.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing, specifically to an adaptive color equalization method for arbitrary orientation seams based on spatiotemporal Kalman filtering. Background Technology

[0002] With the rapid development of remote sensing technology, the mosaicking of multi-source, multi-temporal remote sensing images has become a core component of the production of Digital Orthophoto Maps (DOM) products, and is widely used in many fields such as natural resource surveys, geographic information updates, and ecological environment monitoring. However, due to multiple factors such as differences in imaging time phases, variations in sunlight angle, sensor response characteristic drift, and seasonal changes in the Earth's surface, there are generally significant systematic differences in radiometric brightness and color between remote sensing images of adjacent map sheets. If mosaicking is performed directly, visible color abrupt changes will appear at the seams, which not only reduces the overall visual quality of remote sensing products but also seriously interferes with subsequent applications such as ground feature interpretation and quantitative remote sensing inversion.

[0003] Currently, mainstream color radiometric correction methods can be broadly categorized into two types: general color transfer methods and seam color fusion schemes specifically designed for remote sensing mosaicking. In the general color transfer approach, global linear color transfer methods, such as Reinhard color transfer and Adaptive Instance Normalization (AdaIN), and histogram matching methods based on the cumulative distribution function alignment of the entire image, typically achieve color correction by statistically analyzing the mean, standard deviation, or histogram features of the entire image and then performing a global linear transformation. While these methods can quickly eliminate the overall radiometric deviation between two images, their inherent limitation—based on a unified mapping parameter constructed using a full-image statistical algorithm—means the correction process lacks spatial awareness and cannot adapt to local scenes where the radiometric difference on both sides of the seam continuously changes along the seam direction. This results in residual tonal variations still remaining in the local seam area of ​​the globally corrected image. Local window adaptive methods, such as Wallis filtering, attempt to overcome this deficiency. These methods establish local statistical windows centered on pixels and estimate local correction parameters window by window, which can respond to local radiometric differences to some extent. However, in scenarios with long seams, curved seams, and large-format images, these methods need to perform window statistics pixel by pixel. The computational load increases geometrically with the image size and seam length. Furthermore, the correction parameters of adjacent windows are estimated independently, lacking global consistency constraints. Random noise and abnormal ground object interference in the local statistical process are directly transmitted to the parameter sequence, causing the correction parameters to jitter along the seam direction. This can easily generate new local color discontinuities and color spot artifacts along the seam in the corrected image.

[0004] In seam color fusion schemes specifically designed for remote sensing mosaicking, the weighted transition method with seam buffers masks visual differences at the seam by setting transition bands on both sides of the seam and weighting the two images according to the distance from the pixel to the seam line. However, this method is essentially a pixel-level numerical mixing rather than a true correction at the radiometric parameter level. When the radiometric differences between the two sides are large, the weighted mixing is prone to ghosting, blurring, and even information loss. Furthermore, this method cannot output reusable explicit correction parameters, which is not conducive to traceability and batch reuse in subsequent production processes. Although the multi-resolution pyramid fusion method alleviates the abrupt transition problem under large radiometric differences to some extent by using different transition widths at different frequency levels, it is still essentially a pixel mixing method. When there are large differences in low-frequency components, it can lead to large-scale color bleeding. At the same time, constructing a full image pyramid requires a huge amount of memory and computing resources, making it difficult to adapt to the engineering production needs of large-format remote sensing images. The gradient domain Poisson fusion method reconstructs images by solving the Poisson equation using seam boundary conditions, achieving gradient continuity at ideal boundaries. However, solving large-scale sparse linear equations incurs extremely high memory and time complexity. Furthermore, this method diffuses radiometric differences at the seam throughout the entire area to be corrected, causing an overall tonal shift in the image, which does not meet the stringent requirements of faithful reproduction of original radiometric features in remote sensing image production. The piecewise linear stretching method along the seam divides the transition zone into several segments and performs linear stretching segment by segment based on statistics from both sides. However, the correction parameters for each segment are also independently estimated and unconstrained, leading to visible segmental boundaries due to parameter jumps between segments. Moreover, the segmented statistical results are highly sensitive to local noise and anomalous features, exhibiting poor parameter stability.

[0005] Overall, both the local window adaptive method and the seam segmented linear stretching method have insurmountable technical defects in practical applications. On the one hand, the correction parameters of each window or segment are estimated independently, lacking smooth consistency constraints along the seam direction. This leads to random noise and anomalous features from local statistics directly entering the correction parameter sequence, causing the parameters to fluctuate irregularly along the seam direction, ultimately resulting in bright and dark spots and segment boundaries distributed along the seam. On the other hand, both methods implicitly assume that the seam is approximately horizontal or vertical, or rely excessively on pixel-by-pixel dense window calculations. For irregular polygonal seams with arbitrary orientations and multiple turns, existing methods either cannot achieve unified parameterized modeling along the seam direction, or require splitting into multiple segments for separate processing, with parameters prone to abrupt changes at the segment points. At the same time, the computational load increases dramatically with the seam length and image size.

[0006] Therefore, how to simultaneously achieve noise-resistant smooth estimation of seam parameters and uniformly adapt to arbitrary orientation polyline seams has become a technical problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0007] To alleviate or partially alleviate the above-mentioned technical problems, the solution of the present invention is as follows:

[0008] This invention discloses an adaptive color equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering, comprising:

[0009] Step S1: Obtain pixel data of the reference image and the image to be corrected within the seam buffer, and perform global linear pre-alignment on the image to be corrected based on the global statistics within the seam buffer;

[0010] Step S2: Set control points along the seam line at a fixed step size, map the control points into a one-dimensional sequence ordered along the seam arc length using a unified arc length parameterization method, and determine the tangent and normal at each control point;

[0011] Step S3: Perform local statistics within the preset neighborhood of each control point, calculate the local correction parameters at each control point, and obtain the observation parameter sequence corresponding to the arc length order of the one-dimensional sequence;

[0012] Step S4: Establish a Kalman filter state space model, use the observation parameter sequence as the observation input of the state space model, and sequentially perform forward Kalman filter recursion and reverse Kalman filter recursion along the arc length order of the one-dimensional sequence.

[0013] The arithmetic mean of the forward Kalman filter recursion result and the reverse Kalman filter recursion result is used to obtain the smoothed correction parameter sequence at each control point.

[0014] Based on the tangent and normal directions at each control point, the smoothing correction parameters of each control point are interpolated to each pixel in the seam buffer using anisotropic spatial weights, thereby obtaining the correction parameters of each pixel in the seam buffer.

[0015] Step S5: Use the correction parameters of each pixel to perform a linear transformation on the pixels of the image to be corrected, and combine the mixing coefficient of attenuation along the normal distance of the seam to control the correction intensity, and output the color-corrected image.

[0016] In one implementation example, step S1, the global linear pre-alignment specifically includes:

[0017] The pixels of the reference image and the image to be corrected located in the seam buffer are converted to the Lab color space, and the global mean and global standard deviation of the L, a, and b channels are calculated respectively.

[0018] A linear transformation is performed on the pixel values ​​of each channel of the image to be corrected, so that the mean and standard deviation of the transformed image in the L, a, and b channels are aligned with the mean and standard deviation of the L, a, and b channels of the reference image, respectively.

[0019] In one implementation example, in step S2, the unified arc length parameterization method is specifically as follows:

[0020] For joints with a broken or curved shape, the cumulative arc length is calculated along the broken line of the joint, and the control points are sorted in ascending order of arc length to form the one-dimensional sequence.

[0021] The local tangent direction of each control point in the middle of the polyline is taken as the bisector of the angle between the directions of the polyline segments before and after the control point, and the direction of the first and last control points is taken as the direction of the polyline segment in which they are located; and the tangent direction of adjacent control points is smoothed by moving average, and the local normal direction of each control point is obtained by rotating the local tangent direction of the control point counterclockwise by 90°.

[0022] In one implementation example, in step S2, for seams with near-straight-line shapes, the tangential and normal directions are determined by principal component analysis:

[0023] Principal component analysis is performed on the spatial coordinates of pixels within the seam buffer to solve for the principal eigenvectors of the coordinate covariance matrix, thereby obtaining the overall tangent direction vector and normal direction vector of the seam line.

[0024] The criteria for determining a near-straight line shape are: the ratio of the secondary eigenvalue to the principal eigenvalue of the covariance matrix is ​​less than 0.05, and the maximum turning angle of adjacent line segments of the seam polygonal line is less than 15°.

[0025] In one implementation example, step S3 specifically includes:

[0026] Within a circular neighborhood centered on the control point and with a neighborhood radius proportional to the sampling step size, the mean and standard deviation of the Lab three channels are statistically analyzed for the corresponding pixels of the reference image and the image to be corrected, and the local correction parameters are calculated.

[0027] During the local statistical process, the following processing is performed: invalid pixel removal, saturated pixel removal, and outlier removal are performed sequentially, and observation confidence scores are recorded simultaneously.

[0028] The outlier removal method employs a multi-iterative removal approach based on the initial mean and standard deviation. After removing outliers whose pixel values ​​deviate from the initial mean by more than a preset multiple, the statistical analysis is performed again.

[0029] The observation confidence record is to record the number of valid pixels that ultimately participate in the statistics within the neighborhood, and control points with a number of valid pixels lower than a preset threshold are marked as low-confidence observations.

[0030] In one implementation example, step S4 specifically includes:

[0031] The forward Kalman filter recursion is performed sequentially from the first control point to the last control point: calculate the predicted covariance based on the state covariance and process noise of the previous control point; calculate the Kalman gain based on the predicted covariance and observation noise; update the current state estimate based on the Kalman gain, the previous state estimate, and the current observation; and update the current state covariance based on the Kalman gain and the predicted covariance.

[0032] The reverse Kalman filter recursion, from the last control point to the first control point, uses the state estimate and covariance at the end of the forward Kalman filter recursion as initial values ​​and performs the same prediction and update steps.

[0033] In one implementation example, step S4 also includes an engineering parameter tuning step:

[0034] The observation noise covariance is estimated based on the robust residuals of the observation sequence. First, median filtering is performed on the observation sequence to obtain the trend term, and then the observation noise covariance is calculated based on the deviation between the observations and the trend term.

[0035] The process noise covariance is tuned by setting the ratio of process noise to observation noise, which is controlled within the range of 0.01 to 0.5.

[0036] For control points marked as low-confidence observations, the corresponding observation noise covariance is adaptively amplified based on the ratio of the actual effective pixel count to the reference effective pixel count.

[0037] In one implementation example, the specific method for performing interpolation using anisotropic spatial weights in step S4 is as follows:

[0038] Calculate the tangential and normal distances from any pixel within the seam buffer to each control point;

[0039] The interpolation weights for each control point corresponding to the pixel are calculated using an anisotropic Gaussian kernel function. The interpolation weights are obtained by performing an exponential function operation based on the ratio of the tangential distance to the tangential diffusion radius and the ratio of the normal distance to the normal diffusion radius.

[0040] Wherein, the tangential diffusion radius is greater than the normal diffusion radius;

[0041] The final correction parameter at the pixel is equal to the normalized weighted sum of the smooth correction parameters at all control points.

[0042] In one implementation example, step S5, where the mixing coefficient that attenuates along the seam normal distance controls the correction intensity, specifically includes: obtaining the absolute value of the normal distance from the pixel to the seam line; when the absolute value of the normal distance is not greater than half the width of the seam buffer, the mixing coefficient attenuates with the normal distance according to the square cosine law; when the absolute value of the normal distance is greater than half the width of the seam buffer, the mixing coefficient is 0; and,

[0043] The original pixel value and the linearly corrected pixel value are weighted and summed using the mixing coefficients to obtain the final output pixel value.

[0044] In one type of implementation example, the joint segment connection process is also included:

[0045] The position where the angle between adjacent broken line segments is greater than 15° or the cumulative segment length exceeds 5 times the control point step length is set as the segment point, and the length of each segment is not less than 5 times the control point step length.

[0046] In step S4, the Kalman filter of the next segment uses the state estimate and covariance at the end of the previous segment as initial values, and performs arc-length linear weighted fusion on the parameters of a set number of control points on both sides of the segment boundary.

[0047] The technical solution of the present invention has one or more of the following beneficial technical effects:

[0048] (1) This invention uses bidirectional Kalman filtering to perform forward and reverse recursive smoothing of the local correction parameter sequence along the arc length of the seam, thereby achieving noise-resistant smooth estimation of the correction parameters along the seam direction. This effectively suppresses the parameter jitter caused by random noise and abnormal ground features in the local statistical process, significantly reduces the risk of color patches and segment boundaries appearing along the seam in the corrected image, and improves the visual consistency and parameter stability of color balance at the seam.

[0049] (2) This invention automatically extracts the main direction of the joint through principal component analysis and processes the polyline joint into a one-dimensional topological sequence by combining a unified arc length parameterization mechanism. This decouples the topological order of the filter sequence from the two-dimensional geometric direction of each control point, solving the technical problem that existing local correction methods are difficult to uniformly parameterize at arbitrary directions and polyline joints with multiple turns, and that parameter mutations are easy to occur at segmentation points. This invention achieves universal adaptation to complex joint forms such as near-straight lines, curves, and even polylines, and smooth parameter transition at turning points.

[0050] (3) This invention adopts a two-stage architecture of global alignment followed by local smoothing correction. It combines the responsibility separation design of anisotropic spatial weight interpolation (wide diffusion along the seam direction and narrow diffusion perpendicular to the seam direction) and Alpha mixing coefficient. It can complete the pixel-level high-precision color correction of the entire buffer by performing calculations only at the seam control points. This greatly reduces the computational overhead and memory usage of traditional local methods for pixel-by-pixel window statistics. While meeting the timeliness requirements of large-format remote sensing image engineering production, it also takes into account the natural smoothness of color transition at the seam and the fidelity of the original radiometric features of the image.

[0051] Furthermore, other beneficial effects of the present invention will be mentioned in the specific embodiments. Attached Figure Description

[0052] Figure 1 This is a schematic diagram illustrating the spatial concept of seam buffer zone, control point, and anisotropic diffusion in one embodiment of the present invention;

[0053] Figure 2 This is an overall flowchart of one embodiment of the present invention. Detailed Implementation

[0054] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0055] To facilitate a clear description of the technical solutions in the embodiments of the present invention, the terms "first" and "second" are used to distinguish identical or similar items with essentially the same function and effect. Those skilled in the art will understand that the terms "first" and "second" do not limit the quantity or execution order.

[0056] This invention provides an adaptive color equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering. The method employs a two-stage architecture combining "global pre-alignment" and "local smoothing correction," transferring the Kalman filtering concept from time-series signal processing to the spatial seam correction scenario. The correction parameter sequence distributed along the seam arc length is equated to a one-dimensional time series. Bidirectional Kalman filtering suppresses local statistical noise while preserving the true gradual change trend of radiation differences. Simultaneously, a unified arc length parameterization mechanism adapts to arbitrary-direction, polygonal seams with inflections. Finally, anisotropic spatial interpolation extends the smoothing parameters to the entire seam buffer, achieving a natural and smooth color transition at the seam.

[0057] Figure 1This is a schematic diagram illustrating the spatial concept of seam buffer zones, control points, and anisotropic diffusion in one embodiment of the present invention. Figure 1 As shown, it includes the seam line, the buffer zone boundary formed by inner and outer dashed lines, various control points distributed along the seam, the tangential and normal directions obtained by uniform arc length parameterization at the seam zigzag turning point, and the tangential / normal diffusion radius required for subsequent interpolation.

[0058] The present invention will be further described in detail below with reference to specific embodiments.

[0059] Figure 2 This is an overall flowchart of one embodiment of the present invention. Figure 2 As shown, the reference image and the image to be corrected are first converted to the Lab color space (the Lab color space is a color space designed based on the human visual perception system, where the L component represents brightness, and the a and b components represent color contrast dimensions). Systematic overall radiation bias is eliminated by global statistical alignment. Then, Principal Component Analysis (PCA) is used to automatically extract the geometric direction of arbitrary seam lines, and a unified arc length parameterization mechanism is introduced to process the polyline seam into a one-dimensional topological sequence. Next, local statistics and sampling are performed at the seam control points to construct the observation parameter sequence to be smoothed. Then, bidirectional Kalman filtering is used to perform bidirectional recursive smoothing of the parameter sequence, obtaining the final correction parameters that smoothly transition along the seam direction under the premise of noise resistance. Finally, the correction parameters are interpolated to each pixel in the area to be corrected through anisotropic spatial weights, and combined with the mixing coefficient (e.g., alpha mixing coefficient) that decays with the distance from the pixel to the seam line normal, adaptive equalization correction of color radiation is completed, outputting a mosaic image with natural color transition and no local jumps.

[0060] Specifically, step S1: First, seam buffer generation and global pre-alignment are performed. This step uses Adaptive Instance Normalization (AdaIN) to achieve global linear pre-alignment. Using the seam line specified manually or output by a preceding algorithm as the center line, it extends W / 2 pixels to both sides along the normal direction, forming a strip-shaped buffer with a total width of W. The total width W of the seam buffer is related to the image ground sampling distance (GSD), and the formula for calculating the total width W of the seam buffer is:

[0061] ;

[0062] Where W is the total width of the seam buffer zone, in pixels; D is the desired width of the ground transition zone, ranging from 20 to 100 meters; and GSD is the ground sampling distance of the image. This indicates a rounding up operation; the typical value range of W at the pixel scale is 128 to 512 pixels, and in this embodiment, W is 256 pixels. After generating the initial buffer, the intersection of the buffer polygon and the geometric overlap polygon of the two digital orthophoto maps (DOM) is calculated, and the part exceeding the overlap area is clipped to ensure that the buffer falls completely within the geometric overlap area of ​​the two images; if the normal width of the overlap area is less than W, the buffer width is automatically narrowed to the normal width of the overlap area itself to avoid statistical failure due to the lack of corresponding reference pixels.

[0063] After completing the buffer construction, the buffer pixels of the reference image (hereinafter referred to as Scene A) and the image to be corrected (hereinafter referred to as Scene B) are converted to the Lab color space. The mean and standard deviation of the L, a, and b channels of the two images within the buffer are calculated respectively. Global pre-alignment of the image to be corrected is performed through a global linear transformation. The formula for the global linear transformation is:

[0064]

[0065] Where L', a', and b' are the pixel values ​​of the L, a, and b channels of the image to be corrected after global pre-alignment, respectively; L, a, and b are the original pixel values ​​of the L, a, and b channels of the image to be corrected before pre-alignment, respectively; the superscript forms ^L, ^a, and ^b correspond to the luminance channel, a color contrast channel, and b color contrast channel in the Lab color space, respectively; μ_B^L, μ_B^a, and μ_B^b are the global mean values ​​of the L, a, and b channels of the image to be corrected within the buffer, respectively; σ_B^L, σ_B^a, and σ_B^b are the global standard deviation values ​​of the L, a, and b channels of the image to be corrected within the buffer, respectively; σ_A^L, σ_A^a, and σ_A^b are the global standard deviation values ​​of the L, a, and b channels of the reference image within the buffer, respectively; and μ_A^L, μ_A^a, and μ_A^b are the global mean values ​​of the L, a, and b channels of the reference image within the buffer, respectively.

[0066] This step employs a two-stage design: global first, then local. This design pre-emptively eliminates the systematic overall radiometric bias between two images through a linear transformation of global statistics, converting residual radiometric differences into small-amplitude, localized, gradually varying signals along the seams. Its advantages are twofold: firstly, it ensures that subsequent local observations fall within the effective small-bias range of the linear correction model, preventing the scale parameter from deviating significantly from 1 due to excessive initial bias, thus avoiding noise amplification and numerical instability; secondly, this global alignment process ensures the successful validity of the random walk assumption in the subsequent Kalman filter state equation, facilitating stable engineering tuning of process noise and observation noise parameters. Conversely, if this global stage is skipped and local correction is performed directly, the large systematic overall bias will force significant fluctuations in the local correction parameters, causing the true radiometrically varying signal and local statistical noise to overlap, making it difficult for the filter to effectively distinguish between the two. Therefore, the organic combination of the two-stage design is a crucial prerequisite for the effective execution of subsequent parameter filtering and spatial interpolation processes.

[0067] Step S2: Extract the main direction of the seam and parameterize the arc length of the arbitrary-direction polyline. For seams with near-straight-line shapes, perform Principal Component Analysis (PCA) on the spatial coordinates of all pixels within the buffer zone to solve for the principal eigenvectors of the coordinate covariance matrix, obtaining the global tangent direction vector t and normal direction vector n of the seam line. The criteria for determining near-straight-line seams are: the ratio of the secondary eigenvalue to the principal eigenvalue of the covariance matrix is ​​less than 0.05, and the maximum turning angle of adjacent line segments of the seam polyline is less than 15°; when the criteria are met, it is treated as a single straight-line segment, and a set of tangent direction vectors and normal direction vectors are shared globally. This process does not require manual specification of the seam direction angle and can automatically adapt to near-straight-line seams with arbitrary directions.

[0068] For polygonal or curved joints that do not meet the near-straightness criteria, a unified arc length parameterization method is used as a preferred embodiment of the present invention: the cumulative arc length s is calculated along the joint polygonal line, and all control points are sorted according to their arc lengths to form a single one-dimensional parameter sequence. The subsequent Kalman filtering process is executed once on this unified arc length sequence. The sequence naturally spans all inflection points and remains continuous, with parameters automatically and smoothly transitioning at inflection points without the need for segmented connection processing. Specifically, the local tangent direction t_i of each control point is taken as the bisector of the angle between the directions of the preceding and following polygonal segments, and the tangent directions of adjacent control points are smoothed using a moving average with a window size of 3. The local normal direction n_i of the control point is obtained by rotating the tangent direction t_i counterclockwise by 90°, thus avoiding abrupt changes in the normal direction at inflection points. This design decouples the one-dimensional topological order of the filtering sequence from the two-dimensional geometric direction of each control point, which is the core mechanism of the present invention for adapting to joints with arbitrary orientations.

[0069] As an alternative implementation, if segmentation of the seam is required in some implementation scenarios, such as when the seam is extremely long and needs to be processed in blocks, the segmentation criteria are as follows: when the angle between adjacent broken line segments is greater than 15°, or when the cumulative segment length exceeds 5 times the control point step size, it is set as a candidate segmentation point, and the length of each segment is not less than 5 times the control point step size. If the length is insufficient, it is merged with the adjacent segment. The Kalman filter connection rule for each segment after segmentation is as follows: the state estimate X and covariance P at the end of the previous segment are used as the initial values ​​for the next segment filtering, instead of re-initializing, and the parameters of the three control points on both sides of the segment boundary are subjected to arc length weighted linear fusion to ensure the first-order smooth transition of cross-segment parameters.

[0070] Step S3: Perform seam control point sampling and local statistics. Project the buffer pixels along the tangent direction (for near-straight lines) or along the arc length of a polygonal line (for polygonal lines), and sort them. Select N control points {P1, P2, P3, P4, P5, P6, P7, …, PN} uniformly according to a fixed step size Δs. Figure 1 The circular marker on the joint line is shown. The neighborhood radius R and the step size Δs satisfy: This ensures that the neighborhoods of adjacent control points partially overlap, and that the observation parameter sequences are spatially continuous and seamless; simultaneously, it satisfies R≤W / 2, ensuring that the neighborhood range does not exceed the buffer boundary. The step size Δs is determined based on the desired ground interval, preferably ranging from 50 to 200 meters, converted to pixel units; at the pixel scale, the typical value range of Δs is 64 to 256 pixels. In this embodiment, Δs is set to 128 pixels, R to 128 pixels, and W to 256 pixels, satisfying the above two constraints.

[0071] Since the buffer is located within the geometric overlap area of ​​the two images, after geometric registration, the same pixel position in both images (A and B) contains the same pixel value. Pairing can be performed position by position. Therefore, local statistics do not require classifying pixels by side of the seam. For the same circular neighborhood centered on the control point with radius R, the mean and standard deviation of the Lab three channels are calculated for the same pixels in both the A and B layers, naturally forming statistical results based on the pairing of corresponding points in the overlap area. Only when the overlap area data is provided in a single-sided valid form in certain implementation scenarios does it degenerate into classifying sides by the sign of the signed normal distance d⊥ from the pixel to the seam line. The side with d⊥>0 is assigned to the reference image, and the side with d⊥<0 is assigned to the image to be corrected, and statistics are performed separately.

[0072] During the local statistical process, four levels of outlier removal are performed sequentially: invalid pixel removal, excluding completely transparent pixels, pixels outside the effective mosaic area, and pixels with cloud or cloud shadow masks if the mask is available; saturated pixel removal, removing bright, saturated, and extremely dark pixels with any RGB channel value greater than or equal to 254 or less than or equal to 1; and outlier removal, first calculating the initial mean μ and standard deviation σ of all effective pixels in the neighborhood, and then removing outliers. After re-performing the pixels, statistics are re-executed, and the influence of outliers such as small-area clouds, shadows, and water flares is suppressed through two rounds of 3σ clipping. Here, μ_L is the initial mean of the L-channel pixel values ​​in the neighborhood, and σ_L is the initial standard deviation of the L-channel pixel values ​​in the neighborhood. Observation confidence is recorded, and the number of effective pixels ultimately participating in the statistics within the neighborhood is denoted as n_i. When n_i is less than 10% of the total number of pixels in the neighborhood, the control point is marked as a low-confidence observation. In subsequent filtering steps, its observation noise is adaptively amplified instead of directly removing the control point, thus avoiding gaps in the observation parameter sequence. After completing the statistics and outlier processing, the local ideal correction parameter at the control point is calculated, i.e., the observation parameter vector Z_i at the i-th control point, and:

[0073] Z_i=[scale_L, shift_L, scale_a, shift_a, scale_b, shift_b]_i;

[0074] Where Z_i is the observation parameter vector (local ideal correction parameter) at the i-th control point; the subscript i is the control point number sorted along the joint arc length; scale_L, scale_a, and scale_b are the local scale correction factors for channels L, a, and b, respectively; shift_L, shift_a, and shift_b are the local offset correction amounts for channels L, a, and b, respectively.

[0075] The scale factor and offset are calculated using the following formula: scale = σ_A / σ_B. In the formula, σ_A and μ_A are the local standard deviation and local mean of the corresponding channel of the reference image in the neighborhood of the control point, respectively; σ_B and μ_B are the local standard deviation and local mean of the corresponding channel of the image to be corrected in the neighborhood of the control point, respectively; a very small positive constant is added when calculating σ_B to avoid numerical errors caused by the denominator being 0.

[0076] Step S4: Performing bidirectional Kalman filtering smoothing is the core step of this invention. The correction parameter sequence {Z1, Z2, ..., ZN} is considered as a one-dimensional dynamic system state evolving along the seam direction (i.e., arc length order), and a corresponding state-space model is established, with the following state equation:

[0077] ;

[0078] Where i is the control point number sorted along the joint arc length; X(i) is the system state vector at the i-th control point, i.e., the smoothed correction parameter vector; X(i-1) is the system state vector at the (i-1)-th control point; w(i) is the process noise vector at the i-th control point, which follows a normal distribution with a mean of 0 and a covariance matrix of Q; Q is the process noise covariance matrix, which characterizes the degree of random fluctuation in the evolution of the correction parameters along the joint direction.

[0079] In this invention, the observation equation is:

[0080] ;

[0081] Where Z(i) is the observation parameter vector at the i-th control point; X(i) is the system state vector at the i-th control point; v(i) is the observation noise vector at the i-th control point, which follows a normal distribution with a mean of 0 and a covariance matrix of R_i; R_i is the observation noise covariance matrix corresponding to the i-th control point, which characterizes the degree of random error in the local statistical results.

[0082] Before performing filtering, the relevant engineering parameters are tuned: First, the observation noise R is estimated using the robust residuals of the observation parameter sequence. For each parameter component, a median filter with a window size of 5 is applied to the observation parameter sequence {Z_i} to obtain the trend term. The median absolute deviation (MAD) of the difference between the observed value and the trend term is calculated, taking R = (1.4826 × MAD)², where R is the baseline observation noise covariance of the high-confidence control point; MAD is the median absolute deviation of the difference between the observed value and the median filtered trend term; 1 .4826 is the standard conversion coefficient from median absolute deviation to standard deviation under a normal distribution. For the aforementioned low-confidence control points, the observation noise is adaptively amplified according to the formula R_i = R × (n_ref / n_i), where R_i is the adaptively adjusted observation noise covariance of the i-th low-confidence control point; R is the baseline observation noise covariance; n_ref is the reference effective pixel count, which is half the total number of pixels in the neighborhood of the control point; and n_i is the number of effective pixels actually participating in the statistics in the neighborhood of the i-th control point. The fewer the effective pixels, the greater the observation noise, and the lower the weight of the corresponding observation.

[0083] Second, the process noise Q is tuned using the noise ratio q = Q / R, where q is the ratio of process noise to observation noise (noise ratio), used to control the smoothing strength of the filter; Q is the process noise covariance matrix; and R is the reference observation noise covariance matrix. The smaller the q value, the stronger the smoothing strength and the slower the parameter tracking speed. The preferred range for q is q∈[0.01, 0.5]. In this embodiment, q is set to 0.1, and its equivalent smoothing bandwidth is approximately... The distance between approximately 3 control points is linked to the control point step size Δs.

[0084] Third, the initial value is set to the mean of the first three observations, X(0), and the initial covariance is P(0) = R. A larger initial covariance is used to make the filtering converge quickly, and the estimation results at the beginning of the sequence are compensated by the subsequent inverse filtering. Fourth, the six parameter components each perform the filtering process independently, sharing the same noise ratio q value; the scale factor component can be filtered in the logarithmic domain, that is, the log(scale) is filtered and then the exponential is taken to restore it, ensuring that the restored scale factor is always positive.

[0085] After the parameters are tuned, perform forward Kalman filtering recursion and reverse Kalman filtering recursion on the parameter sequence respectively (hereinafter referred to as forward filtering and reverse filtering):

[0086] (a) Forward filtering performs prediction and update steps point by point from i=1 to i=N, calculating the prediction covariance during the prediction phase: Where P_pred is the one-step prediction covariance matrix; P(i-1) is the state covariance matrix after the (i-1)th control point update; and Q is the process noise covariance matrix.

[0087] (b). Calculate the Kalman gain again: K = P_pred / (P_pred + R_i); where K is the Kalman gain; R_i is the observation noise covariance matrix of the i-th control point;

[0088] (c). Calculate the state estimate during the update phase: X(i) = X(i-1) + K × (Z(i) - X(i-1)); where X(i) is the updated state vector of the i-th control point; Z(i) is the observation parameter vector of the i-th control point;

[0089] (d). Then update the covariance: P(i) = (1-K) × P_pred; where P(i) is the updated state covariance matrix of the i-th control point.

[0090] The reverse filtering process proceeds from i=N to i=1, using the final result of the forward filtering as the initial value, and repeats the prediction and update process (a)~(d). Finally, the arithmetic mean of the forward Kalman filter recursion result and the reverse Kalman filter recursion result is taken to obtain a bidirectionally smoothed correction parameter sequence {X1, …, XN} along the seam direction. Bidirectional filtering eliminates the inherent phase lag problem of unidirectional recursive filtering, ensuring good smoothness of parameters at both the start and end points of the seam. This filtering process is a single-pass recursive operation with a time complexity of O(N). The total complexity of the six parameter components is O(6N), where N is the number of control points, typically ranging from tens to hundreds. The computational load is extremely low, meeting the engineering application requirements of large-format remote sensing imagery.

[0091] Step S5: Perform anisotropic spatial weighted interpolation. After obtaining the smoothing correction parameters at the control points, the parameter interpolation needs to be extended to all pixels in the entire buffer. For any pixel p within the region to be corrected, calculate its tangential distance component d∥ and normal distance component d⊥ to each seam control point P_i. The calculation formulas are as follows:

[0092]

[0093] d∥ is the tangential distance component of pixel p relative to the i-th control point; d⊥ is the normal distance component of pixel p relative to the i-th control point; p is the spatial coordinate vector of the pixel to be corrected; P_i is the spatial coordinate vector of the i-th control point; t_i is the unit vector of the tangent direction at the i-th control point; n_i is the unit vector of the normal direction at the i-th control point; “⋅” indicates the vector dot product operation.

[0094] Furthermore, an anisotropic Gaussian kernel function is used to calculate the interpolation weights for each control point corresponding to the pixel. The weight formula is as follows:

[0095] ;

[0096] w_i(p) is the interpolation weight of the i-th control point corresponding to pixel p; The natural exponential function is used; d∥ is the tangential distance component of pixel p relative to the i-th control point; d⊥ is the normal distance component of pixel p relative to the i-th control point; σ∥ is the tangential diffusion radius, σ⊥ is the normal diffusion radius, in pixels, and its diffusion pattern is as follows: Figure 1 As shown by the dashed line of the ellipse.

[0097] In this invention, the tangential diffusion radius σ∥ is proportional to the control point step size Δs, with a value range of σ∥ = (1.5~3)×Δs. In this embodiment, σ∥ = 2Δs is used to ensure that the interpolation weight of any pixel covers 2 to 3 adjacent control points along the seam direction, achieving a continuous transition of parameters along the seam direction. The normal diffusion radius σ⊥ is proportional to the buffer width W, with a value range of σ⊥ = (0.25~0.5)×(W / 2). In this embodiment, σ⊥ = W / 4 is used to ensure that the correction effect is constrained within the buffer range in the normal direction. The ratio σ∥ / σ⊥ typically ranges from 4 to 16, naturally satisfying the characteristic that σ∥ is much larger than σ⊥, reflecting the physical constraint that parameters change slowly along the seam direction and decay rapidly perpendicular to the seam direction, without the need for additional direction judgment logic. The final correction parameter for pixel p is the normalized weighted sum of all control point parameters, calculated using the following formula:

[0098] ;

[0099] X(p) is the final correction parameter vector at pixel p; ∑ is the summation operation over all control points i; w_i(p) is the interpolation weight of the i-th control point corresponding to pixel p; X_i is the smoothing correction parameter vector of the i-th control point; the denominator is the sum of the interpolation weights of all control points, realizing the normalization of the weights.

[0100] In this step, the normalized interpolation weights are only responsible for the spatial distribution of the correction parameters between control points, that is, determining the correction parameter value used at pixel p; while the intensity envelope of the applied correction amount is independently undertaken by the Alpha mixing coefficients in the subsequent steps. The two work in series and their responsibilities do not overlap. First, the parameter X(p) at pixel p is obtained through interpolation in this step, and then the intensity of the applied correction amount is determined by the Alpha mixing coefficients.

[0101] Pixel-by-pixel color correction is applied. The Lab value of pixel p is linearly transformed using the interpolated local correction parameters. Taking the L channel as an example, the linearly corrected pixel value of the L channel at pixel p is:

[0102] ;

[0103] Where L_corr(p) is the pixel value of the L channel at pixel p after linear correction; scale_L(p) is the scale correction parameter for the L channel at pixel p; L(p) is the original pixel value of the L channel at pixel p before correction; shift_L(p) is the offset correction parameter for the L channel at pixel p; and the a and b channels undergo similar linear transformations. The correction intensity is then controlled by the alpha blending coefficients to obtain the final output pixel value of the L channel at pixel p:

[0104] ;

[0105] Where L_out(p) is the final output pixel value of the L channel at pixel p; α(p) is the correction intensity mixing coefficient at pixel p, with a value range of [0, 1]; L(p) is the original pixel value of the L channel at pixel p before correction; and L_corr(p) is the pixel value of the L channel at pixel p after linear correction.

[0106] The calculation logic for the Alpha mixing coefficient is as follows:

[0107] When |d⊥|≤W / 2, α(p)=cos²(π×|d⊥(p)| / W), where α(p) is the correction intensity mixing coefficient at pixel p, |d⊥(p)| is the absolute value of the normal distance from pixel p to the seam line, and W is the total width of the buffer; in other words, α(p) decays according to the cosine square law as the absolute value of the normal distance from pixel to the seam line |d⊥(p)| decreases.

[0108] When |d⊥|>W / 2, α=0, meaning the Alpha value on the seam line is 1, and full-strength correction is performed; the Alpha value at the outer edge of the buffer is 0, and the correction strength smoothly returns to zero. The role of the Alpha mixing coefficient is to control the correction strength envelope at the edge of the buffer, so that the correction effect smoothly disappears from the seam outwards. It is not used to balance global and local corrections, nor is it mixed with interpolation weights.

[0109] Finally, the corrected Lab values ​​are cropped to the legal range, where the L channel ranges from [0, 100] and the a and b channels range from [-128, 127]. The image is then converted back to the RGB color space using a standard color space conversion matrix and written back to the image to be corrected. The final output is a corrected image with natural seam transitions and no local tone jumps.

[0110] To objectively and quantitatively evaluate the actual correction effect of the adaptive color equalization method for arbitrary seams in remote sensing images based on spatiotemporal Kalman filtering provided in this embodiment, the following verification indicators are established: The residual color difference across the seam ΔE is specifically calculated by taking a pair of opposing windows along the seam at every step length Δs, with the seam line as the axis of symmetry and a normal offset of ±R / 2. The International Commission on Illumination Lab (CIELAB) color difference ΔE*ab corresponding to the Lab mean value in these windows is calculated, and the average and maximum values ​​are statistically analyzed along the seam direction. The accuracy of color correction is measured by comparing the values ​​before and after correction. The parameter smoothness is specifically evaluated by calculating the first-order difference variance of the correction parameter sequence. The ratio of the variance of the original observation value Z before filtering to the smoothed value X after filtering reflects the algorithm's noise resistance and smoothing ability. Algorithm time consumption and memory usage are also assessed. Under the same hardware configuration and test image conditions, the end-to-end processing time and peak memory usage of this method and the comparison method are statistically analyzed to measure the efficiency of engineering applications.

[0111] In specific comparative verification, the first verification scenario targeted long curved and zigzag seams. Artificial zigzag seam images with more than two turns and a length of no less than 20,000 pixels were selected for testing to examine the parameter continuity at the turns and the color consistency across the seam. Experimental results showed that before correction, the mean and maximum residual color difference across the seam were 5.89 and 7.45, respectively. After processing using the existing global linear color migration method and histogram matching method, both indicators decreased to 1.30 and 2.44. After processing using the local window adaptive method, i.e., Wallis filtering, the mean and maximum residual color difference could be reduced to 0.08 and 0.23, respectively. However, the first-order difference variance of the parameters was as high as 0.196, and the maximum jump in correction between adjacent windows reached 0.34. The correction result exhibited obvious visible patchy discontinuities, and the parameter jump at the seam turns also reached 0.14. In this embodiment, after processing the scene, the average and maximum residual color differences are 0.06 and 0.20, respectively, and the first-order difference variance of the parameters is only 0.011, which is 17.3 times lower than the Wallis method. The abrupt change between adjacent windows is reduced to 0.13, and the parameters at the transition points are smooth and continuous with no visible abrupt changes. This proves that this embodiment can simultaneously achieve the correction accuracy of the local method and the parameter smoothness and consistency of the global method at the joint of the polyline.

[0112] The second verification scenario addresses the condition where radiation differences gradually change along the seam direction. Images showing a continuous gradual change in radiation difference between adjacent flight zones due to differences in solar altitude angles were selected for testing to verify the algorithm's adaptive matching capability for spatial radiation gradations. The experiment statistically analyzed the average residual color difference for the first, middle, and last segments. The corresponding values ​​for the global linear color migration method were 2.36, 0.59, and 2.32, respectively, while the corresponding values ​​for the histogram matching method were 2.35, 0.59, and 2.32. Both global methods exhibited significant residual errors at both ends of the seam, demonstrating the inherent deficiency of global methods in lacking spatial awareness. In stark contrast, the average residual color difference in the first, middle, and last segments of this embodiment is 0.05, 0.05, and 0.06, respectively. The results along the seam are highly balanced, with the overall average and maximum residual color difference being only 0.05 and 0.14, respectively. Furthermore, the first-order difference variance of the parameter sequence is reduced by 30.8 times compared to the original observation, verifying the superior adaptability of this embodiment in complex spatial gradient scenarios.

[0113] The third verification scenario targets low-quality images containing noise and thin clouds. Images with grainy noise or thin cloud patches are selected, and the filtering noise robustness of this embodiment is measured using parameter smoothness and residual chromatic aberration as the main indicators. Without cloud masking and outlier removal, the average and maximum residual chromatic aberrations of the local window adaptive method are 0.95 and 5.20, respectively, with a first-order difference variance of 1.667. The maximum jump in adjacent window correction is 4.32, indicating that local statistical noise is directly introduced into the correction parameters, forming strong visible color artifacts on the image. Under the same condition of not enabling cloud masking, the average and maximum residual chromatic aberrations of this embodiment are 0.63 and 3.31, respectively, and the first-order difference variance of the parameter is reduced to 0.057. The bidirectional Kalman filtering reduces this indicator by 29.3 times compared to the original observation value, improving the noise robustness by more than an order of magnitude. Furthermore, after enabling the outlier removal and cloud masking processes described above in this embodiment, the average residual color difference of the correction results is reduced to 0.10, the maximum residual color difference is reduced to 0.49, and the effective pixel ratio of some cloud control windows is reduced to 0.55, triggering the confidence-driven adaptive amplification mechanism for observation noise R in this embodiment. The experimental results indirectly verify the effectiveness of this low-confidence observation weighting design.

[0114] To better illustrate the present invention, numerous specific details have been provided in the detailed embodiments described above. Those skilled in the art should understand that the present invention can be practiced even without certain specific details. In some instances, methods, means, elements, and circuits well known to those skilled in the art have not been described in detail in order to highlight the main points of the present invention.

[0115] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A color adaptive equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering, characterized in that, include: Step S1: Obtain pixel data of the reference image and the image to be corrected within the seam buffer, and perform global linear pre-alignment on the image to be corrected based on the global statistics within the seam buffer; Step S2: Set control points along the seam line at a fixed step size, map the control points into a one-dimensional sequence ordered along the seam arc length using a unified arc length parameterization method, and determine the tangent and normal at each control point; Step S3: Perform local statistics within the preset neighborhood of each control point, calculate the local correction parameters at each control point, and obtain the observation parameter sequence corresponding to the arc length order of the one-dimensional sequence; Step S4: Establish a Kalman filter state space model, take the observation parameter sequence as the observation input of the state space model, and perform forward Kalman filter recursion and reverse Kalman filter recursion in sequence along the arc length order of the one-dimensional sequence. The arithmetic mean of the forward Kalman filter recursion result and the reverse Kalman filter recursion result is used to obtain the smoothed correction parameter sequence at each control point. Based on the tangent and normal directions at each control point, the smoothing correction parameters of each control point are interpolated to each pixel in the seam buffer using anisotropic spatial weights, thereby obtaining the correction parameters of each pixel in the seam buffer. Step S5: Use the correction parameters of each pixel to perform a linear transformation on the pixels of the image to be corrected, and combine the mixing coefficient with the attenuation along the normal distance of the seam to control the correction intensity, and output the color-corrected image.

2. The color adaptive equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering according to claim 1, characterized in that, In step S1, the global linear pre-alignment specifically includes: The pixels of the reference image and the image to be corrected located in the seam buffer are converted to the Lab color space, and the global mean and global standard deviation of the L, a, and b channels are calculated respectively. A linear transformation is performed on the pixel values ​​of each channel of the image to be corrected, so that the mean and standard deviation of the transformed image in the L, a, and b channels are aligned with the mean and standard deviation of the L, a, and b channels of the reference image, respectively.

3. The color adaptive equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering according to claim 1, characterized in that, In step S2, the unified arc length parameterization method is specifically as follows: For joints with a broken or curved shape, the cumulative arc length is calculated along the broken line of the joint, and the control points are sorted in ascending order of arc length to form the one-dimensional sequence. The local tangent direction of each control point in the middle of the polyline is taken as the bisector of the angle between the directions of the polyline segments before and after the control point, and the direction of the first and last control points is taken as the direction of the polyline segment in which they are located; and the tangent direction of adjacent control points is smoothed by moving average, and the local normal direction of each control point is obtained by rotating the local tangent direction of the control point counterclockwise by 90°.

4. The color adaptive equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering according to claim 1, characterized in that, In step S2, for seams with near-straight-line shapes, the tangential and normal directions are determined by principal component analysis: Principal component analysis is performed on the spatial coordinates of pixels within the seam buffer to solve for the principal eigenvectors of the coordinate covariance matrix, thereby obtaining the overall tangent direction vector and normal direction vector of the seam line. The criteria for determining a near-straight line shape are: the ratio of the secondary eigenvalue to the principal eigenvalue of the covariance matrix is ​​less than 0.05, and the maximum turning angle of adjacent line segments of the seam polygonal line is less than 15°.

5. The color adaptive equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering according to claim 1, characterized in that, In step S3, the local statistics specifically include: Within a circular neighborhood centered on the control point and with a neighborhood radius proportional to the sampling step size, the mean and standard deviation of the Lab three channels are statistically analyzed for the corresponding pixels of the reference image and the image to be corrected, and the local correction parameters are calculated. During the local statistical process, the following processing is performed: invalid pixel removal, saturated pixel removal, and outlier removal are performed sequentially, and observation confidence scores are recorded simultaneously. The outlier removal method employs a multi-iterative removal approach based on the initial mean and standard deviation. After removing outliers whose pixel values ​​deviate from the initial mean by more than a preset multiple, the statistical analysis is performed again. The observation confidence record is to record the number of valid pixels that ultimately participate in the statistics within the neighborhood, and control points with a number of valid pixels lower than a preset threshold are marked as low-confidence observations.

6. The color adaptive equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering according to claim 1, characterized in that, Step S4 specifically includes: The forward Kalman filter recursion is performed sequentially from the first control point to the last control point: calculate the predicted covariance based on the state covariance and process noise of the previous control point; calculate the Kalman gain based on the predicted covariance and observation noise; update the current state estimate based on the Kalman gain, the previous state estimate, and the current observation; and update the current state covariance based on the Kalman gain and the predicted covariance. The reverse Kalman filter recursion, from the last control point to the first control point, uses the state estimate and covariance at the end of the forward Kalman filter recursion as initial values ​​and performs the same prediction and update steps.

7. The color adaptive equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering according to claim 6, characterized in that, Step S4 also includes an engineering parameter tuning step: The observation noise covariance is estimated based on the robust residuals of the observation sequence. First, median filtering is performed on the observation sequence to obtain the trend term, and then the observation noise covariance is calculated based on the deviation between the observations and the trend term. The process noise covariance is tuned by setting the ratio of process noise to observation noise, which is controlled within the range of 0.01 to 0.

5. For control points marked as low-confidence observations, the corresponding observation noise covariance is adaptively amplified based on the ratio of the actual effective pixel count to the reference effective pixel count.

8. The color adaptive equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering according to claim 1, characterized in that, In step S4, the specific method for performing interpolation using anisotropic spatial weights is as follows: Calculate the tangential and normal distances from any pixel within the seam buffer to each control point; The interpolation weights for each control point corresponding to the pixel are calculated using an anisotropic Gaussian kernel function. The interpolation weights are obtained by performing an exponential function operation based on the ratio of the tangential distance to the tangential diffusion radius and the ratio of the normal distance to the normal diffusion radius. Wherein, the tangential diffusion radius is greater than the normal diffusion radius; The final correction parameter at the pixel is equal to the normalized weighted sum of the smooth correction parameters at all control points.

9. The color adaptive equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering according to claim 1, characterized in that: In step S5, controlling the correction intensity using the mixing coefficient that attenuates along the seam normal distance specifically includes: obtaining the absolute value of the normal distance from the pixel to the seam line; when the absolute value of the normal distance is not greater than half the width of the seam buffer, the mixing coefficient attenuates with the normal distance according to the square cosine law; when the absolute value of the normal distance is greater than half the width of the seam buffer, the mixing coefficient is 0; and, The final output pixel value is obtained by weighted summing of the original pixel value and the linearly corrected pixel value using the mixing coefficient.

10. The color adaptive equalization method for arbitrary-direction seams based on spatiotemporal Kalman filtering according to claim 1, characterized in that, It also includes the steps for handling the joint segment connection: The position where the angle between adjacent broken line segments is greater than 15° or the cumulative segment length exceeds 5 times the control point step length is set as the segment point, and the length of each segment is not less than 5 times the control point step length. In step S4, the Kalman filter of the next segment uses the state estimate and covariance at the end of the previous segment as initial values, and performs arc-length linear weighted fusion on the parameters of a set number of control points on both sides of the segment boundary.