Spaceborne IPDA Echo Noise Suppression Method for Timing Signal Image Processing

Abstract

This invention discloses a method for suppressing echo noise in spaceborne IPDA lidar, relating to the field of remote sensing technology, through image-based processing of time-series signals. The invention converts online and offline echo time-series signals into grayscale images. Multi-level denoising is performed on the offline denoised grayscale image to obtain a high-quality reference image. This reference image guides the local linear model filtering of the online original grayscale image to achieve structural information transfer. A weighted coefficient matrix is ​​then dynamically generated based on the surface reflectance matching degree of the two original images. Finally, the online structure-transferring image and the offline denoised image are adaptively weighted and fused. This invention, by guiding online filtering and adaptive fusion with the offline denoised grayscale image, effectively solves the inherent contradiction between denoising and fidelity preservation in existing single-signal denoising techniques. While significantly suppressing noise, it accurately preserves the true gas absorption structure information and surface detail features, providing a high-quality data foundation for greenhouse gas concentration inversion using spaceborne IPDA lidar.

Description

Technical Field

[0001] This invention relates to the field of remote sensing technology, specifically to a method for suppressing echo noise in spaceborne IPDA for time-series signal image processing. Background Technology

[0002] Spaceborne IPDA lidar is an active atmospheric remote sensing device that inverts the column concentration of greenhouse gases such as carbon dioxide in the atmosphere by emitting laser pulses of specific wavelengths towards the Earth's surface and receiving the echo signals. Spaceborne IPDA lidar typically emits two laser pulses with extremely similar wavelengths simultaneously: one located at the center of the gas absorption spectral line, and the other at the edge of the absorption spectral line. The two laser pulses travel along identical atmospheric transmission paths and are reflected back to the satellite from the surface. By comparing the difference in energy between the two echoes, common-mode interference such as atmospheric molecular scattering, aerosol attenuation, and changes in surface reflectivity can be eliminated, thereby accurately calculating the amount of gas absorbed along the path. However, the echo signals received by the spaceborne platform are extremely weak and severely contaminated by various noise sources such as solar background noise, detector shot noise, and thermal noise, resulting in a very low signal-to-noise ratio for single-pulse echoes, making them unsuitable for direct concentration inversion. Therefore, efficient noise suppression of spaceborne IPDA echo signals is a crucial prerequisite for improving the accuracy of gas concentration detection.

[0003] Currently, noise suppression techniques for lidar echo signals mainly fall into two categories: time-domain filtering methods and frequency-domain analysis methods. In time-domain filtering, commonly used denoising algorithms include classic statistical filtering methods such as moving average filtering, Kalman filtering, Wiener filtering, and median filtering. These methods suppress random noise by analyzing the time-domain statistical characteristics of the signal. In frequency-domain analysis, time-frequency analysis tools such as wavelet transform, empirical mode decomposition, Fourier transform, and Hilbert-Huang transform are widely used for denoising non-stationary signals, achieving noise suppression by separating noise and signal components in different frequency domains. Furthermore, machine learning-based denoising methods are gradually emerging, such as dimensionality reduction techniques like principal component analysis and independent component analysis, as well as intelligent algorithms like support vector machines and neural networks, which achieve adaptive denoising by learning the characteristic differences between noise and signal. At the hardware level, improving detector sensitivity, optimizing optical system design, and employing multi-pulse averaging techniques have also improved the signal-to-noise ratio of the echo signal to some extent.

[0004] However, the aforementioned existing technologies are all based on a single-signal processing framework, meaning they only perform noise suppression on the online echo signal itself. This single-signal processing approach has a fundamental flaw: without prior knowledge, the denoising algorithm cannot accurately distinguish between noise components and true structural information in the echo signal, leading to a dilemma of over-smoothing or noise residue during the denoising process. Specifically, while a strong denoising strategy can effectively suppress noise, it may also erase important structural details and edge features in the signal; conversely, while a weak denoising strategy can better preserve the signal structure, the noise suppression effect is often unsatisfactory. This inherent contradiction between denoising and fidelity preservation makes it difficult for existing single-signal denoising methods to achieve ideal noise suppression results in spaceborne IPDA lidar echo signal processing, severely restricting further improvements in lidar detection accuracy. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method for suppressing echo noise in spaceborne IPDA for time-series signal image processing.

[0006] To achieve the above objectives, the technical solution of the present invention is as follows: a spaceborne IPDA echo noise suppression method for time-series signal image processing, comprising:

[0007] The online and offline echo timing signals of a preset number of pulses from the spaceborne IPDA lidar are acquired and arranged according to the pulse sequence number to construct the online and offline original grayscale images, respectively.

[0008] The offline original grayscale image is denoised to obtain an offline denoised grayscale image. Using the offline denoised grayscale image as a reference image, the online original grayscale image is subjected to local linear model filtering. For each pixel, within a local window centered on it, the linear transformation coefficient that minimizes the weighted mean square error between the filtered output and the offline denoised grayscale image is solved. The filtered output value of the pixel is calculated using the linear transformation coefficient. The online structure-transfer grayscale image is obtained by processing pixel by pixel.

[0009] Based on the online and offline original grayscale images, the surface reflectance matching degree value is calculated pixel by pixel, and online and offline weighting coefficients are generated for each pixel location to obtain a weighting coefficient matrix;

[0010] Based on the weighting coefficient matrix, the online structure-transmitted grayscale image and the offline denoised grayscale image are weighted and fused pixel by pixel to generate an online denoised grayscale image.

[0011] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0012] 1. This invention constructs a dual-channel collaborative processing framework, which obtains reliable structural prior information through high-quality denoising of offline original grayscale images, and uses this information to guide the local linear model filtering of online original grayscale images, thus realizing the effective transmission of structural information. This technical solution based on the guided filtering idea enables the algorithm to accurately preserve the real surface reflectivity structure and edge features while suppressing noise, thus solving the inherent contradiction between denoising and fidelity preservation in a single signal processing framework.

[0013] 2. This invention arranges a continuous pulse time-series signal into a two-dimensional grayscale image according to the pulse sequence number, so that the originally discrete noise points present a regular texture in the image spatial domain. This allows for efficient noise reduction using mature digital image processing algorithms. The rich spatial context information in the two-dimensional neighborhood enables the local linear model to more accurately estimate the expected output of each pixel. By image processing, the time-series signal noise problem is transformed into a structured image processing problem, which significantly improves the signal-to-noise ratio gain of noise suppression compared to one-dimensional point-by-point processing.

[0014] 3. This invention proposes a dynamic weighted fusion mechanism based on surface reflectance matching degree. By calculating the normalized grayscale difference between the online and offline original images and converting it into a matching degree value, the consistency between the two channels at each pixel location can be accurately quantified. Based on this, a spatially continuous weighted coefficient matrix is ​​generated through smoothing processing and sigmoid function mapping. This achieves an intelligent fusion strategy that trusts the online signal to preserve gas absorption characteristics when the matching degree is high, while relying on the offline signal to avoid noise contamination when the matching degree is low. This spatially adaptive fusion mechanism ensures that the final output image achieves optimal processing results in different regions, effectively balancing noise suppression and signal fidelity.

[0015] 4. This invention employs a three-stage cascaded offline denoising strategy, sequentially removing periodic noise, Gaussian noise, and restoring edges to provide a high-fidelity reference benchmark for subsequent processing. The local linear model filtering utilizes adaptive local windows, regularization terms, adaptive weights, and multi-scale fusion, significantly improving adaptability to complex terrains and numerical stability. The matching degree calculation incorporates normalized grayscale differences, exponential decay, and S-shaped mapping, ensuring the fusion decision remains highly sensitive to uncertain regions while saturation at both ends guarantees stable output under extreme conditions. The synergistic effect of these multiple stages enables this invention to operate stably under the harsh conditions of low signal-to-noise ratio and strong interference on spaceborne platforms, and its multi-stage refined design enhances noise suppression and structural integrity under complex terrain conditions. Attached Figure Description

[0016] The disclosure of this invention is illustrated with reference to the accompanying drawings. It should be understood that the drawings are for illustrative purposes only and are not intended to limit the scope of protection of this invention. In the drawings, the same reference numerals are used to refer to the same parts. Wherein:

[0017] Figure 1 This is a flowchart of the method of the present invention;

[0018] Figure 2 This is a flowchart of the image processing of the present invention. Detailed Implementation

[0019] It is readily understood that, based on the technical solution of this invention, those skilled in the art can propose various interchangeable structural methods and implementations without altering the essential spirit of the invention. Therefore, the following detailed embodiments and accompanying drawings are merely illustrative examples of the technical solution of this invention and should not be considered as the entirety of the invention or as limitations or restrictions on the technical solution of this invention.

[0020] like Figure 1 and Figure 2 As shown, the specific implementation steps of the present invention include the following steps:

[0021] The online and offline echo timing signals of a preset number of pulses from the spaceborne IPDA lidar are acquired and arranged according to the pulse sequence number to construct the online and offline original grayscale images.

[0022] This step converts the one-dimensional echo timing signals of the two wavelength channels of the spaceborne IPDA lidar—namely, the online wavelength sensitive to gas absorption and the offline wavelength whose influence from target gas absorption is negligible—into a two-dimensional grayscale image. The online echo timing signal refers to the laser pulse echo signal corresponding to the center wavelength of the target gas absorption spectral line; the offline echo timing signal refers to the laser pulse echo signal corresponding to the flank of the target gas absorption spectral line and far from the center wavelength. This signal is not significantly attenuated by the target gas absorption along the transmission path, and therefore can serve as a pure background reference characterizing the surface reflection characteristics and atmospheric molecule / aerosol scattering characteristics.

[0023] Specifically, after each laser pulse is emitted, the sequence of amplitude changes over time corresponding to that pulse is received, which is the echo timing signal. Within a preset time window, a preset number of echo timing signals are continuously acquired. Then, the timing signal of each pulse is used as a row of the image, and each sampling point is used as a column of the image. They are arranged from top to bottom according to the order of pulse emission to obtain the online original grayscale image and the offline original grayscale image, respectively.

[0024] For example, if a pulse contains 1024 time sampling points, then each pulse corresponds to one row of the image, containing 1024 pixels; if there are 500 pulses, then a grayscale image of 500×1024 is formed.

[0025] After determining the preset quantity, it is necessary to further determine the corresponding time window to accurately control the start and end times of data reception during actual acquisition. The time window is the preset number of pulses divided by the pulse repetition frequency. The pulse repetition frequency is a fixed parameter of the spaceborne IPDA lidar system, representing the number of laser pulses emitted per second. Those skilled in the art will understand that this time window is an engineering approximation, and the actual error is within the system's tolerance range.

[0026] The time window determines the length of the ground along the orbit covered by the image processing: the window duration multiplied by the satellite speed equals the ground coverage distance. Therefore, through the adaptive pulse number, the time window is automatically adjusted, thus using a shorter window in flat areas and a longer window in complex areas.

[0027] In practical systems, to prevent excessively long time windows from causing non-negligible geometric errors due to Earth's rotation or satellite attitude changes, or excessively short windows from resulting in insufficient pulse acquisition for effective noise reduction, a minimum window duration and a maximum window duration can be preset. The minimum window duration is no less than the time interval between two pulses to ensure that at least enough pulses are acquired for interpolation completion. The maximum window duration is determined based on the preset maximum allowable ground position removal and the satellite's ground velocity. If the calculated window duration exceeds the range, it is limited to the boundary value, and the actual number of pulses acquired is adjusted accordingly to ensure that data acquisition meets both noise suppression requirements and the physical constraints of the spaceborne platform.

[0028] The number of preset pulses is determined specifically as follows:

[0029] Based on the terrain features of the target area scanned by the satellite-borne IPDA, the terrain complexity index of the target area is calculated. The terrain complexity index is obtained by calculating a weighted combination of the slope change rate and the standard deviation of elevation within the target area. The preset quantity is adjusted based on the terrain complexity index.

[0030] The method for adjusting the preset number based on the terrain complexity index is as follows: when the terrain complexity index is lower than the first threshold, the preset maximum number of pulses is used; when the terrain complexity index is higher than the second threshold, the preset minimum number of pulses is used; when the terrain complexity index is between the first threshold and the second threshold, the number of pulses and the terrain complexity index are inversely proportional in a piecewise linear relationship.

[0031] This step defines how the number of pulses is automatically adjusted according to the complexity of the terrain. Its design intent is to collect more pulses in flat areas to improve the signal-to-noise ratio, and to reduce the number of pulses in complex terrain areas to avoid smoothing out important spatial details.

[0032] The slope change rate is calculated based on the digital elevation model of the target area. The slope value is obtained by calculating the ratio of the elevation difference between adjacent grid points to the horizontal distance. Then, the degree of change of all slope values ​​in the entire target area is statistically analyzed. The average absolute deviation of the slope values ​​is used to quantify the severity of slope change, reflecting the spatial variation rate of local surface tilt: the larger the value, the more frequent and severe the topographic undulation.

[0033] The elevation standard deviation is calculated as the dispersion of the elevation values ​​of all grid points within the target area relative to the average elevation. It is calculated as the square root of the average of the sum of squares of the differences between the elevation of each grid point and the average elevation, reflecting the overall fluctuation range in the vertical direction.

[0034] The terrain complexity index is calculated by normalizing the slope change rate and elevation standard deviation to ensure they fall within the same numerical range, and then weighting and summing them according to a preset weight ratio.

[0035] Since the slope change rate directly affects the laser incident angle, and thus has a more direct effect on the change in surface reflectivity, a higher weight is given to the slope change rate. For example, the weight of the slope change rate can be set in the range of 0.6 to 0.8, while the weight of the elevation standard deviation is set in the range of 0.2 to 0.4. In actual implementation, the weights can be pre-adjusted within the range according to the geological characteristics of the target area, and the sum of the weights is 1.

[0036] The first threshold, the second threshold, the maximum number of pulses, and the minimum number of pulses are all fixed values ​​preset according to the system design parameters of the spaceborne IPDA lidar, and the second threshold is greater than the first threshold.

[0037] When the terrain complexity index is below the first threshold, it indicates that the terrain is relatively flat. At this time, the preset maximum number of pulses is used to improve the signal-to-noise ratio through multi-pulse coherent accumulation, making full use of the characteristics of flat areas with gentle spatial changes for efficient noise suppression. When the terrain complexity index is above the second threshold, it indicates that the terrain is extremely complex. At this time, the preset minimum number of pulses is used to avoid spatial detail blurring caused by multi-pulse averaging, ensuring that the edge and structural information of complex terrain is preserved. When the terrain complexity index is between the first and second thresholds, the number of pulses and the terrain complexity index have a piecewise linear inverse relationship. That is, the number of pulses decreases linearly as the terrain complexity index increases. This linear relationship ensures a smooth transition of the number of pulses under different terrain complexities.

[0038] For example, in plains areas with a terrain complexity index of 0.15, the system will automatically select the maximum number of pulses of 1000 for processing, thereby obtaining a higher signal-to-noise ratio and more stable signal quality by increasing the number of pulses; in mountainous areas with a terrain complexity index of 0.85, the minimum number of pulses of 200 will be selected, thereby prioritizing the preservation of terrain edges and detailed features by reducing the number of pulses, while reducing computational complexity and storage overhead.

[0039] Before arranging by pulse number, the process also includes pulse number continuity detection and interpolation completion, specifically:

[0040] Obtain the timing tag for each pulse, where the timing tag is the sequential number of the pulse emission time;

[0041] Scan the timing tags in chronological order to determine whether the timing tags of adjacent pulses are consecutive;

[0042] If a timing tag is detected to be missing, an interpolated echo timing signal for the missing pulse is generated by using an interpolation method based on the echo timing signals of a preset number of valid pulses before and after the missing position, and the interpolated echo timing signal is inserted into the corresponding missing position.

[0043] The completed pulse sequence is arranged in ascending order of the time tag, so that the constructed online original grayscale image and the offline original grayscale image have continuous pulse numbers in the row direction.

[0044] Pulse sequence continuity detection is used to address pulse loss issues caused by satellite platform jitter and data transmission packet loss. The preset number is usually set to 3 to 7, preferably 5, which means using 3 to 7 valid pulses before and after the missing position for interpolation.

[0045] The interpolation method can be linear interpolation, spline interpolation, or correlation-based interpolation. Taking linear interpolation as an example, if the i-th pulse is missing, but the (i-1)-th and (i+1)-th pulses exist, the interpolated echo signal of the i-th pulse can be calculated as the average of the gray values ​​of the corresponding time sampling points of the two valid pulses before and after it. For cases with multiple consecutive missing points, higher-order interpolation methods can be used. For example, if three pulses are consecutively missing—the i-th, (i+1)-th, and (i+2)-th—a polynomial fitting model can be constructed using five valid pulses before and after it for interpolation. This continuity detection and interpolation completion ensures that subsequent image processing will not produce artifacts or errors due to data discontinuity, guaranteeing the spatial continuity of the image.

[0046] In practice, the system first reads the time sequence labels of all pulses, then checks the continuity of the label sequence in chronological order. If a label jump is detected, such as jumping directly from 100 to 102, the missing label 101 is identified, and the interpolation process is initiated. After interpolation, all the completed pulses are arranged from top to bottom in ascending order of their time sequence labels, so that each row of the grayscale image corresponds to a consecutive pulse number, and adjacent rows spatially correspond to consecutive ground positions, providing a reliable data foundation for subsequent spatial filtering processing.

[0047] This step constitutes the front-end data preparation stage of the entire noise suppression method, converting the raw pulse echo sequence received by the spaceborne IPDA lidar into a standardized, continuous grayscale image format with adaptive spatial sampling density. Specifically, online and offline raw grayscale images are constructed by arranging the temporal signals of a certain number of pulses within a preset time window according to their sequence numbers, providing a unified data structure for subsequent image processing. Based on this, the number of pulses is adaptively determined according to the terrain complexity index of the target area, thereby controlling the spatial sampling density of the grayscale image. Furthermore, pulse sequence continuity detection and interpolation completion are performed before arrangement to ensure the correctness and integrity of the spatial adjacency relationships between image rows. This step provides high-quality, physically meaningful, and resource-efficient input data for subsequent core operations such as offline grayscale image denoising, online structure transfer filtering, and weighted fusion.

[0048] This step significantly improves the targeting and efficiency of data acquisition through an adaptive adjustment mechanism of terrain complexity index and pulse number. This allows the system to acquire more pulse data in flat terrain areas to improve signal-to-noise ratio and noise suppression, while appropriately reducing the pulse number in complex terrain areas to preserve fine spatial structure and edge details. Pulse sequence continuity detection and interpolation completion techniques effectively solve the problem of pulse loss or data interruption that may occur during actual observation. By generating echo timing signals of missing pulses through interpolation, the continuity of the constructed image in the row direction is guaranteed, thus avoiding artifacts or distortion caused by data discontinuity during image processing. By constructing two sets of original grayscale images, one online and one offline, a comparative benchmark is provided for subsequent noise suppression and quality assessment, enabling the system to more accurately identify and separate noise components. After the original one-dimensional signal is visualized, the previously difficult-to-process noise exhibits regular texture features in the spatial domain, facilitating the full utilization of mature digital image processing algorithms to achieve better noise suppression results than traditional one-dimensional signal processing methods.

[0049] In the implementation process, regarding the calculation of the terrain complexity index, in addition to using a weighted combination of slope change rate and elevation standard deviation, other terrain feature parameters such as curvature change rate, terrain undulation, and surface roughness can be introduced, or machine learning methods can be used to directly extract complexity features from the digital elevation model. Regarding the pulse number adjustment strategy, in addition to piecewise linear inverse proportional relationships, more complex adaptive algorithms such as nonlinear function relationships, fuzzy logic control, and neural network adaptive adjustment can be used. In the pulse sequence continuity detection and interpolation completion stage, in addition to interpolation methods based on a preset number of effective pulses, alternative techniques such as time series prediction models, signal reconstruction methods based on physical models, or spatial interpolation using data from adjacent scan lines can be used. Regarding image construction methods, in addition to constructing grayscale images by pulse sequence, different organizational methods can be used, such as gridded arrangement based on geographical location, dynamic arrangement based on timestamps, or fusion arrangement combining multi-source data.

[0050] The offline original grayscale image is denoised to obtain an offline denoised grayscale image. Using the offline denoised grayscale image as a reference image, the online original grayscale image is subjected to local linear model filtering. For each pixel, within a local window centered on it, the linear transformation coefficient that minimizes the weighted mean square error between the filtered output and the offline denoised grayscale image is solved. The filtered output value of the pixel is calculated using the linear transformation coefficient. The online structure-transfer grayscale image is obtained by processing pixel by pixel.

[0051] This step uses the high-quality denoising result of the offline denoised grayscale image as a reference benchmark to guide the filtering process of the online original grayscale image, thereby achieving effective noise suppression and accurate transmission of structural information. Specifically, the logic of the local linear model filtering process is as follows: within a local spatial window, there exists a local linear mapping relationship between the surface reflectivity structure in the noise-contaminated online echo signal and the surface reflectivity structure in the high-quality offline echo signal. By solving this mapping relationship, the reflectivity structure information that has been effectively recovered in the offline signal can be transmitted to the online signal, thereby achieving noise suppression and structure preservation. The local linear model filtering process borrows the core principle of guided filtering, that is, using the structural information of the guiding image to guide the filtering process of the target image.

[0052] In the actual calculation process, for each pixel in the image, a local window is first defined around it. Then, the grayscale values ​​of all pixels within this window are collected as the input feature vector. Simultaneously, the pixel value at the corresponding position in the offline denoised grayscale image is obtained as the reference target vector. Subsequently, an optimization problem is solved to determine the optimal linear transformation coefficients. The goal of this optimization problem is to minimize the weighted mean square error between the filtered output and the offline denoised grayscale image. Finally, the obtained linear transformation coefficients are used to filter the center pixel to obtain the final output value of that pixel. By performing this processing flow pixel by pixel, the entire online structure-transfer grayscale image is finally obtained. This image retains both the real-time characteristics of the original online image and the high-quality characteristics of the offline denoised grayscale image.

[0053] The noise reduction process specifically includes:

[0054] The offline original grayscale image is initially denoised using a nonlocal mean denoising algorithm to obtain a preliminary denoised grayscale image.

[0055] The initial denoised grayscale image is subjected to wavelet threshold denoising to obtain a wavelet denoised grayscale image;

[0056] Anisotropic diffusion filtering is applied to the wavelet-denoised grayscale image to obtain the offline denoised grayscale image;

[0057] The diffusion coefficient of the anisotropic diffusion filtering process is negatively correlated with the local gradient magnitude of the wavelet-denoised grayscale image.

[0058] This denoising process employs a three-stage cascade approach. The nonlocal mean denoising algorithm primarily targets structural and periodic noise in the image; wavelet threshold denoising mainly targets Gaussian white noise and shot noise; and anisotropic diffusion filtering further smooths flat areas while preserving edge structures. The three stages complement each other, progressively improving image quality.

[0059] First, a nonlocal means algorithm is used for initial denoising. This algorithm utilizes the repetitive texture structures present in the image, searching for similar blocks within a large area around each pixel, and suppressing noise through weighted averaging. This is used to remove periodic or quasi-periodic noise patterns. Typical parameters are a search window size of 21×21 pixels, a similar block size of 7×7 pixels, and a filter intensity parameter set based on an estimate of the image noise standard deviation, ranging from 0.6 to 0.9 times the noise standard deviation. After nonlocal means processing, large blocks of noise in the image are effectively suppressed, but edges may be slightly blurred.

[0060] Then, wavelet thresholding is applied to the preliminary denoising results: the sym8 wavelet basis is selected, and the decomposition is performed in 3 to 5 layers. A soft thresholding function is used to process the high-frequency coefficients of each layer. The threshold value is adaptively calculated based on the noise level, and is the noise standard deviation multiplied by the square root of twice the logarithm of the number of pixels. Wavelet thresholding can effectively remove Gaussian white noise and shot noise, while preserving high-frequency information at the edges relatively well.

[0061] Finally, anisotropic diffusion filtering is performed: this filter simulates the thermal diffusion equation, but the diffusion coefficient varies with the local gradient magnitude. Specifically, the diffusion coefficient is:

[0062] ;

[0063] in The diffusion coefficient is... This represents the local gradient magnitude of the image. The standard deviation of the gradient magnitude within the neighborhood window of the current pixel. To adjust the parameters, the range of values ​​is: When the local gradient magnitude is large, the diffusion coefficient is small, the diffusion process is suppressed, and the edges are preserved; when the local gradient magnitude is small, the diffusion coefficient is large, the diffusion process is accelerated, and the noise is effectively removed.

[0064] Through three-stage denoising, noise in the offline original grayscale image is suppressed step by step, ultimately obtaining a clear, flat, and smooth offline denoised grayscale image, providing a high-quality reference benchmark for subsequent online original grayscale image filtering processing.

[0065] The local window is an adaptive local window, and the construction of the adaptive local window is specifically as follows:

[0066] Calculate the local gradient magnitude of the offline denoised grayscale image at the center pixel location;

[0067] The shape and size of the local window are determined based on the local gradient magnitude, wherein the larger the local gradient magnitude, the smaller the size of the local window and the larger the aspect ratio, and the major axis of the local window is perpendicular to the gradient direction.

[0068] The adaptive local window is constructed based on the shape and size.

[0069] Since the texture and edge features of the original grayscale image vary greatly in different regions, a fixed-size square window can easily lead to over-smoothing of edge regions and insufficient denoising of flat regions. Therefore, an adaptive local window is introduced to adapt to different image region features by dynamically adjusting the shape and size of the window.

[0070] In practice, the local gradient magnitude of the center pixel in the offline denoised grayscale image is first calculated. The horizontal and vertical gradient components are calculated in the local 3×3 neighborhood using the Sobel or Prewitt operator. The gradient magnitude is obtained by calculating the square root of the sum of the squares of the horizontal and vertical gradient components.

[0071] Then, the parameters of the local window are determined based on the magnitude of the gradient, specifically:

[0072] When the gradient magnitude is less than the first gradient threshold, the first gradient threshold is set to 0.3 to 0.5 times the average gradient of the image, indicating that the region is a smooth region. At this time, the local window adopts a larger size and an isotropic shape, such as a 9×9 or 11×11 square window.

[0073] When the gradient magnitude is greater than the second gradient threshold, which is set to 1.5 to 2.0 times the average gradient of the image, it indicates that the region is a strong edge region. At this time, the local window adopts a shape with a smaller size and a larger aspect ratio, such as a 3×3 or 5×5 window and an elliptical window with a major axis to minor axis ratio of 3:1 or 4:1. The major axis is perpendicular to the gradient direction so that the major axis of the window is parallel to the edge direction in the image, thereby taking more points along the edge direction during filtering and avoiding detail blurring caused by crossing the edge.

[0074] When the gradient magnitude is between the first threshold and the second threshold, the size and shape of the local window are determined by linear interpolation.

[0075] For example, assuming the average gradient of the image is 50, the first threshold is set to 20, and the second threshold is set to 80, when the gradient magnitude of a pixel is 15, an 11×11 square window is used; when the gradient magnitude is 50, a 7×7 window with an ellipse having a major axis to minor axis ratio of 2:1 is used; and when the gradient magnitude is 100, a 3×3 window with an ellipse having a major axis to minor axis ratio of 4:1 and the major axis is perpendicular to the gradient direction is used.

[0076] The specific process for solving the linear transformation coefficients is as follows:

[0077] Construct the grayscale value vector of the pixels within the local window and the reference grayscale value vector of the offline denoised grayscale image at the corresponding position;

[0078] A regularization coefficient to prevent overfitting is introduced into the weighted mean square error to obtain a regularized weighted mean square error.

[0079] Based on the gray value vector, the reference gray value vector, and the regularized weighted mean square error, the regularized linear regression equation is solved to obtain the linear transformation coefficients.

[0080] When solving for linear transformation coefficients within a local window, the conventional least-squares solution may encounter numerical instability, especially when the pixel grayscale variation within the window is extremely small. In this case, the coefficient matrix of the normal equation, i.e., the product matrix of the transpose of the design matrix, approaches singularity, leading to abnormally large or numerically unstable coefficient estimates. To address this issue, a regularization coefficient is introduced into the weighted mean square error. This regularization coefficient is the square of the Euclidean norm of the weight coefficient vector in the linear transformation coefficients multiplied by a small positive number. The ridge parameter, also known as the regularization coefficient, is used to prevent overfitting by constraining the magnitude of the weight coefficients and to improve the condition number of the coefficient matrix of the normal equation, thereby enhancing numerical stability.

[0081] The regularization coefficient The range of values ​​is Preferred When the pixel grayscale within the window changes significantly, Take a smaller value to reduce the constraints on coefficient estimation; when the pixel grayscale change within the window is small, A larger value is chosen to enhance numerical stability. In practical applications, It can adaptively adjust based on the grayscale standard deviation of a local window.

[0082] In practice, the grayscale information of all pixels is first collected within the local window, and two vectors are constructed: one is an input vector containing the grayscale values ​​of all pixels within the window, and the other is a reference vector containing the grayscale values ​​of the corresponding pixels in the offline denoised grayscale image. The lengths of these two vectors are the same, and both are equal to the total number of pixels within the local window.

[0083] Subsequently, a linear transformation model is established. The goal of this model is to find a set of optimal transformation parameters that minimize the difference between the output result obtained after linearly transforming the input vector with these parameters and the reference vector. This set of transformation parameters includes two parts: a weight coefficient vector and a bias term.

[0084] To find the optimal transformation parameters, an optimization objective function needs to be constructed, which consists of two parts: the first part is a weighted mean squared error term, which measures the difference between the transformation output and the reference value, where the error term of each pixel is weighted according to its adaptive weight; the second part is a regularization coefficient, which is used to constrain the magnitude of the weight coefficients to prevent the model from becoming too complex and causing overfitting.

[0085] In the actual solution process, a step-by-step solution strategy is adopted. First, the bias term is calculated. This bias term is obtained by weighted averaging, that is, by multiplying the adaptive weight of each pixel by its corresponding reference gray value, then summing all weighted results and dividing by the total weight to obtain the bias term. This bias term is then substituted into the optimization objective function to further solve for the weight coefficient vector.

[0086] The solution for the weight coefficient vector is specifically obtained by constructing a matrix equation containing the input vector, weight matrix, and regularization coefficients, and then solving the inverse matrix of the equation to obtain the optimal weight coefficients. This process ensures that the weight coefficient vector that minimizes the objective function is found under all constraints.

[0087] The obtained weight coefficient vector and bias term are the linear transformation coefficients, which are used to calculate the filtered output value of each pixel. The weight coefficient vector is responsible for weighted combination of the input gray value vector, and the bias term is responsible for providing overall offset adjustment.

[0088] The weights in the weighted mean square error are adaptive weights, and the calculation of the adaptive weights includes:

[0089] Calculate the spatial distance and grayscale difference between each pixel within the local window and the center pixel;

[0090] The adaptive weight is calculated based on the spatial distance and the grayscale difference;

[0091] The principle for allocating the adaptive weights is as follows: pixels that are closer to the center pixel in terms of spatial location and whose grayscale value is closer to the center pixel's grayscale value are assigned higher weights; conversely, pixels that are farther away in terms of spatial distance or have a greater difference in grayscale value are assigned lower weights.

[0092] The adaptive weights are used to calculate the weighted mean square error.

[0093] This adaptive weighting mechanism takes into account two key factors: the spatial relationship between pixels and the degree of grayscale similarity, thereby accurately evaluating the contribution of each pixel within the local window to the filtering process of the central pixel.

[0094] In the specific implementation process, it is first necessary to calculate the spatial distance between each pixel and the center pixel. The Euclidean distance metric is used, which is to calculate the straight-line distance between two pixels in a two-dimensional coordinate system. The closer the pixel is, the stronger its spatial correlation with the center pixel.

[0095] Secondly, the grayscale difference between each pixel and the center pixel is calculated by taking the absolute difference between the grayscale values ​​of the two pixels. The smaller the grayscale difference, the more similar the two pixels are in terms of visual features.

[0096] After obtaining the spatial distance and grayscale difference, Gaussian weighting was applied to both factors. For spatial distance, a Gaussian function based on the window radius was used for weighting; for grayscale difference, a Gaussian function based on the local window grayscale standard deviation was used. The Gaussian parameter for spatial distance was set to a value equal to the window radius. The Gaussian parameter for the grayscale difference is set to the local window grayscale standard deviation. times.

[0097] The final adaptive weight is the product of the spatial weight and the grayscale weight. This ensures that only pixels that are spatially close and have similar grayscale features receive a higher weight value. Near edges, due to the large grayscale difference, even spatially close pixels are assigned a lower weight, thus protecting the edge structure from damage. In smooth regions, due to the small grayscale difference, pixels in the neighborhood receive a higher weight, fully utilizing neighborhood information for noise suppression.

[0098] For example, in a 7x7 local window, if a pixel is close to the center pixel and has a small difference in grayscale value, then that pixel will play a greater role in the filtering process; conversely, if a pixel is close in space but has a large difference in grayscale value, or if the grayscale values ​​are similar but the spatial distance is far, its weight will be reduced accordingly.

[0099] The local linear model filtering process is a multi-scale local linear model filtering process, which includes:

[0100] The local windows are constructed at multiple preset scales to obtain local windows of multiple scales;

[0101] The corresponding linear transformation coefficients are solved independently within a local window at each scale, thus obtaining the linear transformation coefficients for each of the multiple scales.

[0102] For each scale, the center pixel is predicted using its corresponding linear transformation coefficients, and the difference between the prediction result and the offline denoised grayscale image at the corresponding position is calculated as the prediction error for that scale.

[0103] Based on the magnitude of the prediction error at each scale, a corresponding scale weight is assigned to each scale, with the scale having a larger prediction error receiving a larger weight.

[0104] The linear transformation coefficients of each scale are weighted according to their corresponding scale weights to obtain the fused linear transformation coefficients, and the final filtered output value is calculated using the fused coefficients.

[0105] A single-scale local window is insufficient to capture both large-scale structures and small-scale details simultaneously. Multi-scale local linear model filtering significantly improves the robustness of the filtering process and its adaptability to different image features by performing local linear modeling at multiple different scales and adaptively fusing the results based on the prediction performance at each scale.

[0106] In the specific implementation process, multiple scales are pre-set for processing, with a scale range of [missing information]. indivual.

[0107] For example, four scales can be selected: the first scale corresponds to a small 3×3 window, the second scale corresponds to a medium 5×5 window, the third scale corresponds to a large 7×7 window, and the fourth scale corresponds to a large 9×9 window. Each scale represents a different receptive field size and can capture local features of different ranges.

[0108] For each preset scale, the aforementioned adaptive local window construction method is independently applied. That is, the shape and size of the window at that scale are determined based on the local gradient magnitude of the offline denoised grayscale image at that pixel. The linear transformation coefficients at that scale are solved, including the weight coefficient vector and the bias term. The center pixel is predicted based on the linear transformation coefficients to obtain the predicted value at that scale. This predicted value is calculated based on the local window information of the current scale.

[0109] To evaluate the prediction quality at each scale, the prediction error for that scale needs to be calculated. The prediction error is obtained by comparing the difference between the predicted value at that scale and the reference gray value at the corresponding position in the offline denoised grayscale image. The absolute error is used as the metric. The smaller the prediction error, the better the prediction effect of that scale in the current local area. Therefore, a higher weight is given to that scale.

[0110] The weights are allocated using an exponential decay function, which establishes a non-linear mapping relationship to convert prediction errors into scale weights. When the prediction error is 0, the corresponding weight reaches its theoretical maximum value; as the prediction error increases, the weights decay exponentially.

[0111] The exponent in the exponential decay function is the negative of the ratio of prediction error to decay factor. The calculation process specifically includes: first, determining the decay factor; then, normalizing the prediction error for each scale and calculating the ratio of prediction error to decay factor; next, using this ratio as the negative exponent of the exponential function to calculate the initial weights; finally, normalizing the initial weights for all scales by dividing the weight of each scale by the sum of all scale weights to ensure that the final weight sum is 1. This normalization process ensures that the sum of weights for all scales is 1, guaranteeing that the fused coefficients are still effective linear transformation coefficients, ensuring that the final filtered output value is within a reasonable range, and allowing the contributions of different scales to be accurately reflected in the final result proportionally.

[0112] The decay factor controls the rate of exponential decay, and its value ranges from the maximum value among all scale prediction errors at the current pixel. The range of values ​​for 'a' is 1. The range of values ​​for b is When the attenuation factor is set to the maximum prediction error When the exponential decay curve is steeper, the weight differences corresponding to different prediction errors are more significant; when the decay factor is set to the maximum prediction error... At this time, the exponential decay curve is relatively flat, and the weight differences corresponding to different prediction errors are relatively smaller. If the decay factor is lower than... This leads to an overly extreme weight allocation, causing the scale with the smallest prediction error to receive almost all the weights; if the decay factor is higher than... This leads to a tendency for weight allocation to become averaged, making it impossible to effectively distinguish the differences in prediction quality at different scales.

[0113] After obtaining the weights for each scale, the linear transformation coefficients of each scale are weighted and fused according to their corresponding weights. Specifically, the weight coefficient vector of each scale is multiplied by its corresponding scale weight, and then the results of all scales are summed to obtain the fused weight coefficient vector. Similarly, the bias term of each scale is multiplied by its corresponding scale weight, and then the results of all scales are summed to obtain the fused bias term. The final filtered output value is calculated using the fused weight coefficient vector and the bias term. This filtered output value integrates information from multiple scales and can adaptively reflect the scale combination most suitable for the local features of the current image.

[0114] This step first performs multi-level denoising on the offline original grayscale image to obtain a high-quality offline denoised grayscale image, providing a reliable reference benchmark for subsequent online filtering of the original grayscale image. Then, using this offline denoised grayscale image as a reference, local linear model filtering is performed on the online original grayscale image, realizing the transfer of structural information from the offline high-quality image to the online real-time image. This lays a solid foundation for the noise suppression and image quality improvement of the entire scheme. Specifically, this step, through a series of refined designs such as constructing an adaptive local window, introducing regularization coefficients, adopting an adaptive weighting mechanism, and implementing a multi-scale fusion strategy, ensures that the filtering process can adaptively adjust the processing parameters according to the local features of different image regions. This effectively suppresses noise while preserving the structural details and edge information of the image to the greatest extent. Simultaneously, the output of this step, the online structural transfer grayscale image, will serve as an important input for subsequent processing steps, providing crucial data support for achieving the final goal of the entire scheme.

[0115] This step employs a three-stage cascaded denoising process, progressively suppressing shot noise, thermal noise, and background noise in the offline original grayscale image. This provides a high-fidelity structural prior for subsequent guided filtering, significantly improving the accuracy of online original grayscale image denoising. The adaptive local window adjusts its shape and size based on the gradient. Small, flat elliptical windows are used in edge regions and elongated along the edge direction to effectively prevent cross-edge smoothing, while large square windows are used in flat regions to enhance denoising capabilities, thus achieving a balance between detail preservation and noise suppression. The introduction of regularization coefficients eliminates the ill-conditioned problem of coefficient calculation in flat regions, preventing output value runaway and ensuring the numerical stability of the algorithm. Adaptive weights, by simultaneously considering spatial distance and grayscale similarity, make the filtering focus more on neighborhoods homogeneous with the center pixel, further strengthening edge preservation capabilities. The multi-scale fusion mechanism automatically allocates weights based on prediction errors, enabling the filtering results to balance large-scale structural consistency with small-scale texture details, exhibiting stronger robustness and adaptability in online original grayscale image processing on complex terrains.

[0116] In the implementation process, for denoising, in addition to the three-stage processing flow of nonlocal mean denoising, wavelet threshold denoising, and anisotropic diffusion filtering, other mature denoising algorithms such as total variation denoising, bilateral filtering, guided filtering, and BM3D denoising can also be used, or deep learning methods such as convolutional neural networks can be used for end-to-end denoising. For local window construction, in addition to adaptive adjustment based on gradient magnitude, window parameters can be determined based on texture features, curvature features, or a fixed window combined with edge detection. For solving linear transform coefficients, in addition to regularized linear regression, ridge regression, LASSO regression, support vector regression, or nonlinear regression methods such as multinomial regression and neural network regression can be used as alternative techniques. For weight calculation, in addition to Gaussian weighting based on spatial distance and gray-level differences, other distance metrics such as Mahalanobis distance can be used, or multidimensional features such as color space and texture features can be introduced for comprehensive weight calculation. In terms of multi-scale processing, in addition to constructing windows of different sizes in the spatial domain, multi-scale analysis and fusion can also be achieved in the frequency domain, wavelet domain, or through image pyramid structures.

[0117] Based on the online and offline original grayscale images, the surface reflectance matching degree value is calculated pixel by pixel, and online and offline weighting coefficients are generated for each pixel location to obtain a weighting coefficient matrix.

[0118] After obtaining the online structure-transfer grayscale image and the offline denoised grayscale image, it is necessary to determine how to fuse the two into the final online denoising result. In order to obtain the best fusion effect, the weights need to be dynamically allocated according to the consistency between the online and offline original images at each pixel position.

[0119] If the grayscale values ​​of the two images at a certain pixel are very close, it indicates that the area is less affected by gas absorption or has a low noise level, and the online structure can be trusted to transmit the image more readily. Conversely, if the two images differ significantly, it indicates that the online channel may be subject to strong interference or abrupt changes in surface reflectance, and the offline denoised grayscale image should be relied upon more. Therefore, this step calculates the surface reflectance matching degree value pixel by pixel based on the online and offline original grayscale images, quantifying the degree of grayscale similarity between the two channels at the same location. Based on this matching degree value, an online weighting coefficient and an offline weighting coefficient are generated for each pixel, and the sum of the two is 1. The weighting coefficients of all pixels are arranged according to the original image size to obtain the weighting coefficient matrix.

[0120] The specific steps for generating online and offline weighting coefficients for each pixel location are as follows:

[0121] A smooth surface reflectance matching degree map is generated based on the surface reflectance matching degree value;

[0122] Each pixel value of the smooth surface reflectance matching degree map is mapped to an online weighting coefficient through a sigmoid function. The output value of the sigmoid function increases monotonically with the input value, and has the maximum slope in the middle range of the input value, and tends to saturate at both ends.

[0123] The offline weighting coefficient is calculated based on the online weighting coefficient, such that the sum of the online weighting coefficient and the offline weighting coefficient is 1.

[0124] This step further defines how to derive the final online weighting coefficients from the original pixel-by-pixel matching degree values. Since the original matching degree values ​​may contain isolated noise points, directly using the original matching degree values ​​may cause drastic and discontinuous weight jumps due to image noise, which is not conducive to generating visually smooth and structurally coherent fusion results. Therefore, it is necessary to first generate a matching degree map of the same size as the original image based on the matching degree values, and then smooth the map to obtain a smooth surface reflectance matching degree map, which suppresses local high-frequency noise while preserving the overall spatial distribution characteristics of the matching degree.

[0125] The value of each pixel in the smooth matching degree map is mapped to an online weighted coefficient through a sigmoid function. Specifically, the sigmoid function is:

[0126] ;

[0127] Where S is the online weighting coefficient, with a value range of [0, 1]; x is the matching degree value of the current pixel in the smooth surface reflectance matching degree map, with a value range of [0, 1]. This is the slope control parameter, and its value range is... This parameter controls the steepness of the S-curve in the middle region. When the value of this parameter is large, the curve is steeper, which means that small changes in the middle matching degree range will cause significant adjustments to the weights; when the value of this parameter is small, the curve is relatively flat, and the changes in weights are more moderate. This is the center point of the sigmoid function, which is the input matching degree value corresponding to the function output of 0.5.

[0128] The boundary of the intermediate interval is determined by two threshold parameters: a low matching threshold and a high matching threshold. This is the midpoint between the low matching threshold and the high matching threshold.

[0129] The low matching degree threshold ranges from 0.2 to 0.4. When the matching degree value is lower than this threshold, the system determines that it is in a mismatch state. At this time, the S-shaped function output tends to saturate, and the corresponding online weighting coefficient is close to 0, indicating that the system almost completely distrusts online data and mainly relies on offline benchmark data.

[0130] The high matching degree threshold ranges from 0.6 to 0.8. When the matching degree value is higher than this threshold, the system determines that it is in a matching state. At this time, the function output also tends to saturate, and the corresponding online weighting coefficient is close to 1, indicating that the system highly trusts the online data.

[0131] The range of the intermediate interval is defined by the two thresholds mentioned above, from the low matching degree threshold to the high matching degree threshold. According to the value range of the thresholds, the intermediate interval is the narrowest when the low matching degree threshold is the maximum value of 0.4 and the high matching degree threshold is the minimum value of 0.6, with a width of 0.2; the intermediate interval is the widest when the low matching degree threshold is the minimum value of 0.2 and the high matching degree threshold is the maximum value of 0.8, with a width of 0.6. Within this intermediate interval, small changes in the matching degree value will cause significant adjustments to the online weighting coefficient, making the system highly sensitive to pixels in an uncertain state.

[0132] In practical implementation, in scenarios with high noise levels or frequent surface changes, the intermediate range can be appropriately expanded to enable the system to maintain sensitivity to uncertain situations for a longer period of time; in scenarios with high data quality or relatively stable surfaces, the intermediate range can be appropriately narrowed to enable the system to make clear matching or non-matching judgments more quickly. Through this flexible threshold range setting, this solution can adapt to the needs of different application scenarios.

[0133] The generation of the smooth surface reflectance matching degree map is specifically as follows:

[0134] For each pixel location, calculate the normalized grayscale difference between the grayscale value of the online original grayscale image at that pixel location and the grayscale value of the offline original grayscale image at that pixel location;

[0135] Based on the normalized grayscale difference, the surface reflectance matching degree value is calculated using an exponential decay function;

[0136] A surface reflectance matching degree map is generated based on the surface reflectance matching degree value;

[0137] The surface reflectance matching degree map is subjected to local smoothing filtering to obtain a smoothed surface reflectance matching degree map;

[0138] The local smoothing filtering process is as follows: taking each pixel as the center, within its neighborhood window, a weighted average is performed based on the matching degree values ​​of each pixel within the window to obtain the smoothed matching degree value of that pixel; wherein, the weight of each pixel within the window is determined based on its spatial distance from the center pixel, and the closer the spatial distance, the greater the weight of the pixel.

[0139] This step further refines the calculation process and smoothing method for the surface reflectance matching degree value. The calculation of normalized grayscale difference is the foundation of the entire process. The core purpose is to eliminate the influence of the overall brightness difference that may exist between the online original grayscale image and the offline denoised grayscale image, so that the system can focus on the relative change between the two images. The specific calculation method is as follows: For each pixel position in the image, first calculate the absolute value of the difference between the grayscale value of the online original grayscale image and the offline denoised grayscale image at that position. Then divide this difference by the sum of the grayscale values ​​of the two images at that position, while adding a very small constant to the denominator to prevent division by zero error. The value range is... Through this normalization process, the resulting difference value is constrained to the range of 0 to 1, where 0 indicates that the two images are completely identical at that pixel position, and 1 indicates that the difference is the largest.

[0140] This scheme uses an exponential decay function to convert the normalized grayscale difference values ​​calculated above into matching scores. The smaller the difference between two images, the higher the matching score should be. Therefore, the matching score is calculated using an exponential decay method, that is, when the difference value is 0, the matching score is 1. As the difference value increases, the matching score decays rapidly according to an exponential law. This exponential decay characteristic makes the system more sensitive to small changes in difference, and is used to distinguish different degrees of matching.

[0141] To control the sensitivity of the matching degree to changes in the difference value, this scheme introduces an adjustable control parameter. The range of this control parameter is determined based on the statistical characteristics of the image itself. Specifically, the range of this parameter is set to the maximum normalized grayscale difference value of all pixels in the entire image. The maximum normalized grayscale difference value is calculated by traversing the entire online and offline original grayscale images, calculating the normalized grayscale difference value at all pixel locations, and extracting the maximum value among them.

[0142] The parameter setting method based on the characteristics of the image itself is adaptive: if there are large difference regions in the image, the range of values ​​for the control parameters will be increased accordingly, so that the matching degree calculation is less sensitive to differences and avoids misjudging normal surface changes as anomalies; if the differences in the image are generally small, the range of values ​​for the control parameters will be reduced accordingly, so that the matching degree calculation is more sensitive to small differences and can detect potential anomalies in a timely manner.

[0143] The local smoothing filtering process employs a spatial weighted averaging technique, such as Gaussian filtering. Specifically, a local neighborhood window, such as a 3×3 or 5×5 square window, is defined around each pixel. The weighted average of all pixels within this window is then calculated as the smoothed matching degree value of the center pixel. The weight allocation follows a Gaussian distribution, meaning that pixels closer to the center pixel receive higher weights, while pixels farther away receive lower weights. By utilizing the spatial continuity of ground features, the matching degree calculation error caused by random noise or outliers in individual pixels can be effectively eliminated, making the final weight allocation more spatially reasonable and stable, and avoiding unnatural abrupt changes or noise spots in the weight map.

[0144] This step dynamically calculates the fusion weight of each pixel based on the grayscale similarity between the online and offline original images, thus providing a reliable weight matrix for the subsequent weighted fusion of the online structure-transfer image and the offline denoised grayscale image. Specifically, firstly, the surface reflectance matching degree value is obtained by calculating the normalized grayscale difference pixel by pixel and then performing exponential decay. Next, the matching degree map is locally smoothed to eliminate isolated noise. Then, the smoothed matching degree is mapped to online weighting coefficients using a sigmoid function, and offline weighting coefficients are automatically determined. This series of operations ensures that the weight of each pixel can truly reflect the consistency of the two channels, while taking into account spatial continuity and nonlinear enhancement.

[0145] This step converts grayscale differences into matching scores using an exponential decay function, with the decay coefficient positively correlated with the noise level of the offline denoised grayscale image. This results in a more severe penalty for matching scores in noisier regions, thus avoiding misrepresentation of the online channel under low signal-to-noise ratio conditions. Local smoothing filtering effectively suppresses weight abrupt changes caused by isolated noise points, ensuring a continuous spatial transition of weighting coefficients and preventing artifacts or blocky discontinuities in the fused image. The S-shaped mapping imparts a large slope in the intermediate matching score interval, enhancing the distinguishability of moderately consistent regions and making the fusion decision clearer. Simultaneously, saturation at both ends ensures weight stability under extreme conditions. The final weighting coefficient matrix adaptively balances the contributions of the online structure transfer image and the offline denoised grayscale image, preserving gas absorption features in the online signal when the matching score is high, and relying on the clean background of the offline denoised grayscale image when the matching score is low, thereby achieving the optimal trade-off between noise suppression and signal fidelity.

[0146] In the implementation process, for matching degree calculation, a ratio can be used with a negative exponent after logarithmic transformation; alternatively, absolute difference can be used directly without normalization, but the effect is affected by brightness. The attenuation coefficient in the exponential decay function can be fixed and not adaptive with noise level; linear attenuation can also be used instead. For smoothing, median filtering or bilateral filtering can be used instead of weighted average filtering, and the window shape can be square, circular, or cross-shaped; alternatively, a more robust mapping function can be used without smoothing. For mapping to weighted coefficients, linear mapping or piecewise linear functions can be used instead of sigmoid functions. In addition, the smoothing step can be omitted, and nonlinear mapping can be performed directly on the original matching degree map; or machine learning methods can be used to learn the mapping function from the data. The generation of weighted coefficients can also be independent of online and offline original images, and instead use an online structure to convey the difference between the image and the offline denoised grayscale image.

[0147] Based on the weighting coefficient matrix, the online structure-transmitted grayscale image and the offline denoised grayscale image are weighted and fused pixel by pixel to generate an online denoised grayscale image.

[0148] This step is the final fusion stage of the entire image processing workflow. Its core task is to fuse the online structured grayscale image and the offline denoised grayscale image obtained from the previous processing according to the weighting coefficient matrix dynamically calculated pixel by pixel, thereby generating the final high-quality online denoised grayscale image.

[0149] Specifically, for each pixel in the image, there are three values: first, the grayscale value transmitted by the online structure of the pixel, which preserves the gas absorption structure of the online channel and the noise has been initially suppressed; second, the grayscale value of the pixel after offline denoising, which is a clean surface reflection background with almost no gas absorption; and third, the online weighting coefficient corresponding to the pixel, which ranges from 0 to 1 and reflects the degree of trust in the online channel.

[0150] In practice, the weighted fusion operation employs a pixel-by-pixel linear weighting method. First, the online and offline weighting coefficients for each pixel are read from the weighting coefficient matrix. For each pixel location in the image, the fused grayscale value equals the grayscale value at that location in the online structure-transferred grayscale image multiplied by the online weighting coefficient, plus the grayscale value at that location in the offline denoised grayscale image multiplied by the offline weighting coefficient. Since the sum of the online and offline weighting coefficients is 1, this fusion process is essentially a convex combination of the two images.

[0151] The numerical distribution of the weighting coefficient matrix directly determines the spatial adaptability of the fusion strategy. Specifically, in areas with high surface reflectance matching, i.e., areas where the online original grayscale image and the offline denoised grayscale image are highly similar, the high-quality data of the offline denoised grayscale image is preferred because the surface conditions in these areas are relatively stable, and the denoising effect of the offline denoised grayscale image is more reliable. In areas with low surface reflectance matching, i.e., areas where the online original grayscale image and the offline denoised grayscale image differ significantly, the structural information of the online structure-transmitting grayscale image is preferred because these areas may have real surface changes or the online data has higher credibility.

[0152] For example, if the online structure-transmitted grayscale value of a pixel is 145, the offline denoising grayscale value is 138, and the online weighting coefficient for that pixel in the weighting coefficient matrix is ​​0.65, then the final online denoising grayscale value is... The value is approximately 143, rounded down. This value falls between the two and is biased towards the online value. At another pixel, the online structure transfer grayscale value is 200, the offline denoising value is 80, and the online weighting coefficient is only 0.2, so the final value is... It is closer to the offline value, avoiding the outlier value of 200 that may be caused by noise.

[0153] This step combines two optimized images from different sources: an online structure-transfer grayscale image (preserving the gas absorption structure but potentially containing residual noise) and an offline denoised grayscale image (clean but lacking gas absorption information). These are then convexly combined using a dynamically calculated weighting coefficient matrix, pixel by pixel, to generate a final online denoised grayscale image that combines the advantages of both while minimizing their disadvantages. Specifically, it utilizes the previously generated weighting coefficients to precisely quantify and apply the confidence level of each pixel. This ensures that the online structure-transfer image is prioritized in areas of high matching to preserve the gas absorption signal, while the offline denoised grayscale image is prioritized in areas of low matching to avoid noise contamination. After this step, the output image can be directly used for subsequent greenhouse gas column concentration inversion or other quantitative analyses.

[0154] This step effectively preserves the true structural information and detailed features of the original online grayscale image, avoiding structural blurring or detail loss problems that may occur with traditional denoising methods. It can fully utilize the high-quality denoising results of the offline denoised grayscale image, providing clearer and lower-noise image quality in areas with stable terrain. Through an adaptive weighting mechanism, it can dynamically adjust the fusion strategy according to the specific matching situation of each pixel position, achieving optimal information selection in space.

[0155] In the implementation process, nonlinear weighted fusion methods can also be used, such as weighted fusion based on the sigmoid function, to make the weight transition smoother; multi-scale weighted fusion methods can also be used, performing weighted fusion at different scale levels and then integrating the results of each scale; image patch-based weighted fusion methods can also be used, performing weighted calculations on an image patch basis to better maintain the integrity of the local structure; in addition, deep learning methods can be introduced, using convolutional neural networks to automatically learn the optimal weighted fusion strategy to further improve the fusion quality.

[0156] The technical scope of this invention is not limited to the content described above. Those skilled in the art can make various modifications and variations to the above embodiments without departing from the technical concept of this invention, and all such modifications and variations should fall within the protection scope of this invention.

Claims

1. A method for suppressing echo noise in a spaceborne IPDA for time-series signal image processing, characterized in that, Specifically, it includes: The online and offline echo timing signals of a preset number of pulses from the spaceborne IPDA lidar are acquired and arranged according to the pulse sequence number to construct the online and offline original grayscale images, respectively. The offline original grayscale image is denoised to obtain an offline denoised grayscale image; Using the offline denoised grayscale image as a reference image, a local linear model filtering process is performed on the online original grayscale image. For each pixel, within a local window centered on it, the linear transformation coefficient that minimizes the weighted mean square error between the filtered output and the offline denoised grayscale image is solved. The filtered output value of the pixel is calculated using the linear transformation coefficient. The online structure-transfer grayscale image is obtained by processing pixel by pixel. Based on the online and offline original grayscale images, the surface reflectance matching degree value is calculated pixel by pixel, and online and offline weighting coefficients are generated for each pixel location to obtain a weighting coefficient matrix; Based on the weighting coefficient matrix, the online structure-transmitted grayscale image and the offline denoised grayscale image are weighted and fused pixel by pixel to generate an online denoised grayscale image.

2. The spaceborne IPDA echo noise suppression method for time-series signal image processing according to claim 1, characterized in that: The number of preset pulses is determined specifically as follows: Based on the terrain features of the target area scanned by the satellite-borne IPDA, the terrain complexity index of the target area is calculated. The terrain complexity index is obtained by calculating a weighted combination of the slope change rate and the standard deviation of elevation within the target area. The preset quantity is adjusted based on the terrain complexity index. The method for adjusting the preset number based on the terrain complexity index is as follows: when the terrain complexity index is lower than the first threshold, the preset maximum number of pulses is used; when the terrain complexity index is higher than the second threshold, the preset minimum number of pulses is used; when the terrain complexity index is between the first threshold and the second threshold, the number of pulses and the terrain complexity index are inversely proportional in a piecewise linear relationship.

3. The spaceborne IPDA echo noise suppression method for time-series signal image processing according to claim 1, characterized in that: Before arranging by pulse number, the process also includes pulse number continuity detection and interpolation completion, specifically: Obtain the timing tag for each pulse, where the timing tag is the sequential number of the pulse emission time; Scan the timing tags in chronological order to determine whether the timing tags of adjacent pulses are consecutive; If a timing tag is detected to be missing, an interpolated echo timing signal for the missing pulse is generated by using an interpolation method based on the echo timing signals of a preset number of valid pulses before and after the missing position, and the interpolated echo timing signal is inserted into the corresponding missing position. The completed pulse sequence is arranged in ascending order of the time tag, so that the constructed online original grayscale image and the offline original grayscale image have continuous pulse numbers in the row direction.

4. The spaceborne IPDA echo noise suppression method for time-series signal image processing according to claim 1, characterized in that: The noise reduction process specifically includes: The offline original grayscale image is initially denoised using a nonlocal mean denoising algorithm to obtain a preliminary denoised grayscale image. The initial denoised grayscale image is subjected to wavelet threshold denoising to obtain a wavelet denoised grayscale image; Anisotropic diffusion filtering is applied to the wavelet-denoised grayscale image to obtain the offline denoised grayscale image; The diffusion coefficient of the anisotropic diffusion filtering process is negatively correlated with the local gradient magnitude of the wavelet-denoised grayscale image.

5. The spaceborne IPDA echo noise suppression method for time-series signal image processing according to claim 1, characterized in that: The local window is an adaptive local window, and the construction of the adaptive local window is specifically as follows: Calculate the local gradient magnitude of the offline denoised grayscale image at the center pixel location; The shape and size of the local window are determined based on the local gradient magnitude, wherein the larger the local gradient magnitude, the smaller the size of the local window and the larger the aspect ratio, and the major axis of the local window is perpendicular to the gradient direction. The adaptive local window is constructed based on the shape and size.

6. The spaceborne IPDA echo noise suppression method for time-series signal image processing according to claim 1, characterized in that: The specific process for solving the linear transformation coefficients is as follows: Construct the grayscale value vector of the pixels within the local window and the reference grayscale value vector of the offline denoised grayscale image at the corresponding position; A regularization coefficient to prevent overfitting is introduced into the weighted mean square error to obtain a regularized weighted mean square error. Based on the gray value vector, the reference gray value vector, and the regularized weighted mean square error, the regularized linear regression equation is solved to obtain the linear transformation coefficients.

7. The spaceborne IPDA echo noise suppression method for time-series signal image processing according to claim 1, characterized in that: The weights in the weighted mean square error are adaptive weights, and the calculation of the adaptive weights includes: Calculate the spatial distance and grayscale difference between each pixel within the local window and the center pixel; The adaptive weight is calculated based on the spatial distance and the grayscale difference; The principle for allocating the adaptive weights is as follows: pixels that are closer to the center pixel in terms of spatial location and whose grayscale value is closer to the center pixel's grayscale value are assigned higher weights; conversely, pixels that are farther away in terms of spatial distance or have a greater difference in grayscale value are assigned lower weights. The adaptive weights are used to calculate the weighted mean square error.

8. The spaceborne IPDA echo noise suppression method for time-series signal image processing according to claim 1, characterized in that: The local linear model filtering process is a multi-scale local linear model filtering process, which includes: The local windows are constructed at multiple preset scales to obtain local windows of multiple scales; The corresponding linear transformation coefficients are solved independently within a local window at each scale, thus obtaining the linear transformation coefficients for each of the multiple scales. For each scale, the center pixel is predicted using its corresponding linear transformation coefficients, and the difference between the prediction result and the offline denoised grayscale image at the corresponding position is calculated as the prediction error for that scale. Based on the magnitude of the prediction error at each scale, a corresponding scale weight is assigned to each scale, with the scale having a larger prediction error receiving a larger weight. The linear transformation coefficients of each scale are weighted according to their corresponding scale weights to obtain the fused linear transformation coefficients, and the final filtered output value is calculated using the fused coefficients.

9. The spaceborne IPDA echo noise suppression method for time-series signal image processing according to claim 1, characterized in that: The specific steps for generating online and offline weighting coefficients for each pixel location are as follows: A smooth surface reflectance matching degree map is generated based on the surface reflectance matching degree value; Each pixel value of the smooth surface reflectance matching degree map is mapped to an online weighting coefficient through a sigmoid function. The output value of the sigmoid function increases monotonically with the input value, and has the maximum slope in the middle range of the input value, and tends to saturate at both ends. The offline weighting coefficient is calculated based on the online weighting coefficient, such that the sum of the online weighting coefficient and the offline weighting coefficient is 1.

10. The spaceborne IPDA echo noise suppression method for time-series signal image processing according to claim 9, characterized in that: The generation of the smooth surface reflectance matching degree map is specifically as follows: For each pixel location, calculate the normalized grayscale difference between the grayscale value of the online original grayscale image at that pixel location and the grayscale value of the offline original grayscale image at that pixel location; Based on the normalized grayscale difference, the surface reflectance matching degree value is calculated using an exponential decay function; A surface reflectance matching degree map is generated based on the surface reflectance matching degree value; The surface reflectance matching degree map is subjected to local smoothing filtering to obtain a smoothed surface reflectance matching degree map; The local smoothing filtering process is as follows: taking each pixel as the center, within its neighborhood window, a weighted average is performed based on the matching degree values ​​of each pixel within the window to obtain the smoothed matching degree value of that pixel; wherein, the weight of each pixel within the window is determined based on its spatial distance from the center pixel, and the closer the spatial distance, the greater the weight of the pixel.

Images