A precise method for comparing ceilometer heights.

Abstract

This invention provides an accurate cloud height comparison method, belonging to the field of cloud height technology. The invention establishes an aerosol scattering feature discrimination matrix to eliminate the aerosol layer and retain the true cloud signal. It employs a maximum likelihood estimation-based Bayesian inference algorithm for cloud boundaries to obtain the posterior probability distribution of cloud base height. A two-layer game optimization model is constructed to improve cloud identification accuracy and minimize data registration errors. A cloud motion vector estimation algorithm is used to calculate cloud movement speed and perform spatiotemporal interpolation correction. Correlation coefficients and linear regression slopes are calculated to determine the pass / failability of the comparison. Incremental deviation analysis is performed on all maximum cloud height values ​​to determine full-range performance. This invention solves the technical problem of decreased comparison accuracy caused by spatiotemporal mismatch between the observation data of the reference instrument and the tested instrument due to dynamic cloud movement.

Description

Technical Field

[0001] This invention belongs to the field of cloud height measurement technology, and more specifically, relates to a precise cloud height measurement method. Background Technology

[0002] As a key piece of equipment for meteorological observation, the cloud base height meter measures cloud height using lidar technology and is widely used in aviation meteorological support and weather forecasting. Traditional cloud base height meter comparison methods involve installing a reference instrument and the tested instrument at a fixed location for simultaneous observation, and then evaluating the performance of the tested instrument by comparing the cloud base height data measured by the two instruments. In current cloud base height meter comparison practices, due to the continuous movement of clouds and the spatial distance between the reference and tested instruments, when the cloud movement speed is high, the two instruments may actually be observing cloud regions in different spatial locations at the same time, leading to spatiotemporal registration discrepancies in the observation data. In other words, existing technology suffers from a technical problem where the dynamic movement of clouds causes a spatiotemporal mismatch in the observation data between the reference and tested instruments, resulting in a decrease in comparison accuracy. Summary of the Invention

[0003] In view of this, the present invention provides a precise cloud height comparison method, which can solve the technical problem in the prior art where the dynamic movement of clouds during the cloud height comparison process causes a spatiotemporal mismatch between the observation data of the reference instrument and the test instrument, resulting in a decrease in the accuracy of the comparison.

[0004] This invention is implemented as follows: It provides a precise cloud height measurement method. A reference instrument and a test instrument are installed in a standard location, and nanosecond-level time synchronization is achieved through the BeiDou time synchronization system. The accuracy of the test instrument in measuring distances to fixed targets is verified. The reference instrument and the test instrument are activated to continuously observe and synchronously record the raw data of the laser echo signal. Wavelet multi-scale decomposition is performed on the raw laser echo signal data to extract candidate cloud positions, and depolarization ratio and color ratio parameters are calculated. An aerosol scattering feature discrimination matrix is ​​established to remove the aerosol layer and retain the true cloud signal. A Bayesian inference algorithm based on maximum likelihood estimation is used to obtain the posterior probability distribution of cloud base height. A two-layer game optimization model is constructed to improve cloud identification accuracy and minimize data registration error. A cloud motion vector estimation algorithm is used to calculate the cloud movement speed and perform spatiotemporal interpolation correction. The correlation coefficient and linear regression slope are calculated to determine the pass / failability of the comparison test.

[0005] Specifically, the distance measurement accuracy verification involves setting up 5 fixed-height target boards within the range of the astrolabe, measuring each fixed-height target board 10 times, calculating the standard deviation of the deviation, and proceeding to the next step when the standard deviation of the deviation is less than 50% of the instrument's nominal accuracy.

[0006] The fixed-height target plate is an aluminum alloy plate with a diffuse reflection coating, measuring 1 meter by 1 meter, and is installed on support rods at different heights, with the plate surface perpendicular to the laser emission direction.

[0007] The wavelet multi-scale decomposition uses Daubechies wavelet basis functions for multi-level decomposition, and a peak detection algorithm is applied to extract candidate cloud locations in the wavelet coefficients of the 3rd to 5th layers.

[0008] The peak detection algorithm searches for local maxima in each wavelet coefficient sequence. When the value of a point is more than twice the average of the 20 adjacent sampling points, the point is determined to correspond to a candidate cloud location.

[0009] Wherein, the depolarization ratio is the ratio of the intensity of the received vertically polarized light to the intensity of the parallelly polarized light, and the color ratio parameter is the ratio of the backscattering coefficients of the lidar at a wavelength of 532 nanometers and a wavelength of 1064 nanometers.

[0010] The aerosol scattering feature discrimination matrix is ​​a two-dimensional lookup table. Candidate cloud locations with a depolarization ratio of less than 0.15 and a color ratio parameter greater than 1.8 are marked as aerosol layers and discarded.

[0011] The cloud boundary Bayesian inference algorithm based on maximum likelihood estimation models the laser echo photon count as a Poisson distribution, constructs the prior distribution of cloud base height using historical climate statistics, and calculates the posterior probability distribution using the Bayesian formula.

[0012] The two-layer game optimization model includes an upper-layer model that aims to maximize cloud layer recognition accuracy and a lower-layer model that aims to minimize data registration error, and is solved iteratively through Stackelberg game.

[0013] The cloud motion vector estimation algorithm uses optical flow to calculate the horizontal movement speed of the cloud layer for three consecutive frames of echo images. When the cloud layer movement speed exceeds 5 meters per second, spatiotemporal interpolation correction is performed.

[0014] The spatiotemporal interpolation correction extrapolates the cloud position at the time of observation by the reference instrument to the expected position at the time of observation by the test instrument based on the cloud motion vector, thereby eliminating the registration error caused by cloud movement.

[0015] The correlation coefficient is calculated using the least squares method by fitting the cloud base height data of the reference instrument as the independent variable and the cloud base height data of the test instrument as the dependent variable, which is the Pearson correlation coefficient.

[0016] The linear regression slope is obtained by fitting a univariate linear equation with the reference instrument cloud base height data as the independent variable and the test instrument cloud base height data as the dependent variable using the least squares method.

[0017] Among them, the test instrument is considered qualified when the correlation coefficient is greater than 0.95 and the linear regression slope is in the range of 0.95 to 1.05.

[0018] Among them, all the maximum cloud height values ​​identified during the comparative test are sorted, and the incremental deviation of each point is calculated based on the minimum value. When the absolute value of the incremental deviation of more than 90% of the data points is less than 5% of the instrument's range, the test instrument is judged to be qualified for full range performance.

[0019] The BeiDou time synchronization system provides a unified time reference for the reference instrument and the instrument under test by receiving the Coordinated Universal Time signal broadcast by the BeiDou satellite, with a time synchronization accuracy better than 50 nanoseconds.

[0020] This invention achieves nanosecond-level time synchronization through the BeiDou time synchronization system to eliminate time deviations. It constructs an aerosol scattering feature discrimination matrix by combining wavelet multi-scale decomposition with depolarization ratio and color ratio parameters to eliminate non-cloud interference signals. It obtains the posterior probability distribution of cloud base height using a maximum likelihood estimation-based Bayesian inference algorithm to quantify measurement uncertainty, thus solving the problem of insufficient cloud signal recognition accuracy. Furthermore, this invention calculates cloud movement speed using a cloud motion vector estimation algorithm. When the cloud movement speed exceeds a threshold, it employs spatiotemporal interpolation correction technology to extrapolate cloud positions at different observation times to a unified reference spatiotemporal based on the cloud motion vector, eliminating registration errors caused by cloud movement. Finally, this invention couples and optimizes maximizing cloud recognition accuracy with minimizing data registration errors through a two-layer game optimization model. In summary, this invention solves the technical problem mentioned in the background art where dynamic cloud movement leads to spatiotemporal mismatch between the observation data of the reference instrument and the tested instrument, resulting in decreased comparative measurement accuracy. Attached Figure Description

[0021] Figure 1 This is a flowchart of the method of the present invention.

[0022] Figure 2 This is a distribution diagram of wavelet coefficients at different levels of wavelet multi-scale decomposition.

[0023] Figure 3 This is a distribution map of candidate cloud locations in the aerosol scattering feature discrimination matrix.

[0024] Figure 4 A three-dimensional surface plot showing the peak value of the a priori probability and the change in the signal-to-noise ratio for cloud identification confidence.

[0025] Figure 5 Scatter plot of cloud base height data for the reference instrument and the tested instrument before and after spatiotemporal interpolation correction. Detailed Implementation

[0026] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0027] like Figure 1 As shown, Figure 1 This is a flowchart of a precise ceilometer comparison method provided by the present invention. The method includes the following steps:

[0028] A precise method for comparing ceilometer heights includes the following steps:

[0029] S01. Install the reference instrument and the instrument under test in a standard venue, keeping the distance between the two instruments at more than 3 meters. Make the optical axis parallel by adjusting the base bolts, and use the Beidou time synchronization system to achieve nanosecond-level time synchronization.

[0030] S02. Verify the accuracy of the instrument under test for measuring distance at fixed targets. Set up 5 fixed-height target boards within the range of the astrolabe. Repeat the measurement 10 times for each fixed-height target board and record the deviation sequence between the measured distance and the actual distance.

[0031] S03. Calculate the mean deviation and standard deviation of all deviation values ​​in the test of the fixed target object. If the standard deviation is less than 50% of the instrument's nominal accuracy, proceed to the next step. Otherwise, calibrate the instrument under test and repeat step S02.

[0032] S04. Start the reference instrument and the test instrument to observe continuously for 12 hours, with the sampling frequency set to once per minute, and simultaneously record the raw data of the laser echo signal and the manually observed cloud-like data.

[0033] S05. Perform wavelet multi-scale decomposition on the original laser echo signal data, apply peak detection algorithm to extract cloud candidate positions in the wavelet coefficients of the 3rd to 5th layers, and calculate the depolarization ratio and color ratio parameters of each cloud candidate position.

[0034] S06. Establish an aerosol scattering feature discrimination matrix, mark the candidate cloud locations with a depolarization ratio of less than 0.15 and a color ratio parameter greater than 1.8 as aerosol layers and remove them, while retaining the real cloud signals;

[0035] S07. The real cloud signal is processed by the Bayesian inference algorithm for cloud boundary based on maximum likelihood estimation to obtain the posterior probability distribution and posterior uncertainty interval of the cloud base height of each layer.

[0036] S08. Register the simultaneous multi-layer cloud data identified by the reference instrument and the test instrument. When the cloud movement speed exceeds 5 meters per second, use the cloud motion vector estimation algorithm to perform spatiotemporal interpolation correction on the simultaneous multi-layer cloud data.

[0037] S09. Calculate the correlation coefficient between the cloud base height data of the reference instrument and the cloud base height data of the tested instrument after spatiotemporal interpolation correction. When the correlation coefficient is greater than 0.95 and the linear regression slope is in the range of 0.95 to 1.05, the tested instrument is deemed to have passed the comparative test.

[0038] S10. Sort all the maximum cloud height values ​​identified during the comparison test, and calculate the incremental deviation of each point based on the minimum value. When the absolute value of the incremental deviation of more than 90% of the data points is less than 5% of the instrument's range, the instrument under test is deemed to have qualified performance across the entire range.

[0039] The BeiDou timing system, by receiving Coordinated Universal Time (UTC) signals broadcast by the BeiDou satellites, provides a unified time reference for both the reference instrument and the instrument under test, with a timing accuracy better than 50 nanoseconds, ensuring that the time deviation of the synchronized observation data is negligible. The fixed-height target plate is an aluminum alloy flat plate with a diffuse reflection coating, measuring 1 meter by 1 meter, installed on support rods at different heights, with the plate surface perpendicular to the laser emission direction.

[0040] The distance measured is the distance data obtained by the instrument under test from a target board at a fixed height. The true distance is the actual distance between the target board and the instrument under test measured using a distance measuring device that has passed legal metrological verification. The deviation sequence is a series of all deviation values ​​obtained by subtracting the true distance from the measured distance, arranged in the order of measurement. The deviation mean is the arithmetic mean of all deviation values ​​in the deviation sequence. The deviation standard deviation is the square root of the sum of the squares of the differences between each deviation value and the deviation mean, divided by the length of the deviation sequence minus 1. The instrument's nominal accuracy refers to the maximum permissible measurement error value given in the instrument's manufacturer's technical specifications. For cloud heights less than 1000 meters, it is typically 10% of the measured value; for cloud heights greater than 1000 meters, it is 20%.

[0041] The raw laser echo signal data consists of a time-varying sequence of atmospheric backscattered photon counts received after the reference and tested instruments emitted laser pulses. Manually observed cloud records are textual information recorded by observers regarding the presence, morphology, and observation time of clouds. Wavelet multi-scale decomposition involves multi-level decomposition of the raw laser echo signal data using Daubechies wavelet basis functions. The third level corresponds to a spatial scale of 200 to 400 meters, and the fifth level corresponds to a spatial scale of 800 to 1600 meters. The wavelet coefficients at different levels reflect the signal characteristics of cloud layers of varying thicknesses.

[0042] The peak detection algorithm searches for local maxima in each wavelet coefficient sequence. When the value of a point is more than twice the average of the 20 adjacent sampling points, the point is determined to correspond to a candidate cloud location. The depolarization ratio is the ratio of the received vertically polarized light intensity to the parallel-polarized light intensity. The depolarization ratio of cloud droplets is typically less than 0.05, while that of aerosol particles is between 0.20 and 0.35. The color ratio parameter is the ratio of the backscattering coefficients of the lidar at wavelengths of 532 nm and 1064 nm. The color ratio parameter of water clouds is close to 1, while that of aerosol layers is typically greater than 2.

[0043] The aerosol scattering feature discrimination matrix is ​​a two-dimensional lookup table. The horizontal axis represents the depolarization ratio, divided into 50 intervals from 0 to 0.5, and the vertical axis represents the chromaticity ratio parameter, divided into 50 intervals from 0.5 to 3.0. Each element in the aerosol scattering feature discrimination matrix stores the target type label corresponding to the combination of the depolarization ratio and chromaticity ratio parameter. A target type label value of 1 indicates a real cloud layer, and a target type label value of 0 indicates an aerosol layer. The real cloud layer signal is the portion of the original laser echo signal data corresponding to the real cloud layer after removing the aerosol layer.

[0044] The maximum likelihood estimation-based Bayesian inference algorithm for cloud boundaries models laser echo photon counts as a Poisson distribution with parameter λ, where λ represents the true backscattering intensity. A prior distribution of cloud base height is constructed using historical climate statistics. The posterior probability distribution is calculated using the Bayesian formula combined with real cloud signals. The mean of the posterior probability distribution is used as the cloud base height estimate, and the standard deviation of the posterior probability distribution is used as the posterior uncertainty interval. Simultaneously, multi-layer cloud data is the set of data containing cloud base height and cloud quantity identified by the reference instrument and the tested instrument at the same observation time.

[0045] The cloud movement speed was calculated using optical flow on three consecutive echo images. Optical flow selects cloud feature points in the previous frame and searches for the closest corresponding grayscale value in the next frame; the difference between these two points is divided by the time interval to obtain the cloud movement speed. The cloud motion vector estimation algorithm calculates the horizontal and vertical movement speeds of the clouds based on the cloud movement speed, forming a cloud motion vector. Spatiotemporal interpolation correction extrapolates the cloud position from the reference instrument's observation time to the expected position on the test instrument based on the cloud motion vector, thus eliminating registration errors caused by cloud movement. The extrapolation formula is: the corrected height equals the observed height plus the vertical component of the cloud motion vector multiplied by the time difference.

[0046] The reference instrument's cloud base height data consisted of a sequence of estimated cloud base heights calculated by the reference instrument using a Bayesian inference algorithm based on maximum likelihood estimation. The test instrument's cloud base height data consisted of a sequence of estimated cloud base heights calculated by the test instrument using the same algorithm. The correlation coefficient was calculated using the least squares method to fit the reference instrument's cloud base height data as the independent variable and the test instrument's cloud base height data as the dependent variable, yielding a Pearson correlation coefficient. The linear regression slope was obtained by fitting a univariate linear equation using the least squares method to fit the reference instrument's cloud base height data as the independent variable and the test instrument's cloud base height data as the dependent variable. The linear regression slope reflects the proportional relationship between the measurement results of the reference instrument and the test instrument.

[0047] The total maximum cloud height is the set of maximum cloud base height data at all observation times during the comparison test. The incremental deviation is calculated as the increment of the reference instrument for the i-th maximum cloud height minus the increment of the tested instrument, where the increment is defined as the minimum value among the i-th maximum cloud height minus all maximum cloud height values. The instrument range is the upper limit of the cloud height measurement range specified in the technical specifications of the tested instrument minus the lower limit of the cloud height measurement range, typically between 30 meters and 7500 meters.

[0048] The specific implementation methods of the above steps are described in detail below.

[0049] The specific implementation of step S01 involves first installing the reference instrument and the tested instrument in a standard site that conforms to the construction specifications of marine observation stations. The site should be an open area, ensuring no obstructions of 10 meters in height within a 20-meter radius. The distance between the two instruments should be set to 3 to 5 meters to avoid mutual interference of the laser beams. During installation, the optical axis parallelism is achieved by adjusting three adjusting bolts on the base. The calibration method involves setting a common aiming point at a distance and adjusting the bolts until the laser beams of the two instruments converge at the same point, with the deviation angle controlled within 0.1 degrees. Time synchronization utilizes the timing function of the BeiDou satellite navigation system, receiving the Coordinated Universal Time (UTC) signal broadcast by the BeiDou satellite. The time reference is written into the instrument's internal clock through a hardware clock synchronization module, achieving a timing accuracy of 50 nanoseconds, ensuring that the time deviation between the sampling times of the two instruments is less than 100 nanoseconds. Verification of the fixed-target ranging accuracy was conducted immediately after instrument installation. Five fixed-height target plates at different heights (10m, 500m, 1000m, 1500m, and 2000m) were set up along the laser emission direction in the test area, covering 10% to 90% of the instrument's measurement range. The fixed-height target plates were 1m x 1m aluminum alloy flat plates coated with a Lambertian diffuse reflectance coating, with a reflectivity of 85% to 90%. The plate surface was installed perpendicular to the laser emission direction. The true distance to each fixed-height target plate was measured using a laser rangefinder that had passed national metrological verification, with a measurement accuracy better than 1mm. The measurement was repeated three times, and the average value was used as the true distance benchmark. The instrument measured each fixed-height target plate 10 times consecutively, with a 30-second interval between each measurement. The 10 distance measurements were recorded to form a deviation sequence. Calculate the mean and standard deviation of the deviation sequence. The standard deviation reflects the dispersion of the measurement and the magnitude of random error. When the standard deviation is less than 50% of the instrument's nominal accuracy, the distance measuring system of the tested instrument is considered to be working normally and can proceed to the subsequent comparison test process. Otherwise, the tested instrument needs to have its optical system cleaned, laser energy calibrated, or receiver sensitivity adjusted before re-measuring.

[0050] The specific implementation of step S02 involves simultaneously activating both the reference instrument and the tested instrument into continuous observation mode. The observation duration is set to 12 hours, and the sampling frequency is uniformly set to 1 minute per instance, i.e., triggering laser emission and data recording every 60 seconds to ensure that the observation period covers the complete weather change process from 1 hour after sunrise to 1 hour before sunset. Raw laser echo signal data is recorded using a high-speed data acquisition card with a sampling rate of 100 MHz, recording the backscattered photon count time series within 15 microseconds after each laser pulse emission. The raw data is stored in binary format. Manual cloud observation records are kept by the observer visually observing the sky every 30 minutes, recording the presence of clouds, cloud morphology, cloud thickness estimation, and weather phenomena at the current moment. These manual observation records serve as an independent reference for verifying the instrument's operational status. The raw laser echo signal data preprocessing employed a wavelet multi-scale decomposition method, using the Daubechies wavelet basis function for a 5-level decomposition. The first level corresponds to a spatial scale of 50-100 meters, the second to 100-200 meters, the third to 200-400 meters, the fourth to 400-800 meters, and the fifth to 800-1600 meters. The purpose of wavelet decomposition is to transform the laser echo signal from the time domain to the time-frequency domain. Wavelet coefficients at different levels correspond to signal characteristics of cloud layers with varying thicknesses; thin cloud signals are primarily reflected in the wavelet coefficients of the second and third levels, while thick cloud signals are primarily reflected in the wavelet coefficients of the fourth and fifth levels. The peak detection algorithm is executed in the wavelet coefficient sequence from layer 3 to layer 5. The algorithm principle is to find the local maxima of the wavelet coefficients. The judgment condition is that when the wavelet coefficient value of a certain point is greater than twice the average of the 10 sampling points before and after it, the point is identified as a candidate cloud position. The threshold of twice is determined based on the statistical characteristics of atmospheric turbulence noise, which can effectively distinguish the real cloud signal from the background noise. For each candidate cloud position, the depolarization ratio and color ratio parameters are calculated. The depolarization ratio is obtained by analyzing the intensity ratio of the vertical polarization component to the parallel polarization component in the received optical path. The color ratio parameter is calculated by the ratio of the backscattering coefficients at wavelengths of 532 nm and 1064 nm by a dual-wavelength lidar. The depolarization ratio and color ratio parameters are used for subsequent aerosol and cloud differentiation.

[0051] The specific implementation of step S03 involves constructing an aerosol scattering feature discrimination matrix to automatically identify aerosol layers and real cloud layers. The horizontal axis of the matrix represents the depolarization ratio parameter, with a value ranging from 0 to 0.5 divided into 50 equally spaced intervals, each interval having a width of 0.01. The vertical axis represents the color ratio parameter, with a value ranging from 0.5 to 3.0 divided into 50 equally spaced intervals, each interval having a width of 0.05. Each grid cell in the matrix stores a target type label. The label value is pre-assigned based on historical aerosol and cloud layer observation statistics. When the depolarization ratio is less than 0.15 and the color ratio parameter is greater than 1.8, the corresponding grid label value is 0, representing an aerosol layer; the remaining grid labels are 1, representing a real cloud layer. The thresholds of 0.15 and 1.8 are determined based on Mie scattering theory and cloud droplet particle scattering characteristics. The depolarization ratio and color ratio parameters of the candidate cloud locations extracted in step S02 are substituted into the aerosol scattering feature discrimination matrix to query the corresponding label values. Candidate cloud locations with a label value of 0 are identified as aerosol layers and removed from the dataset, while candidate cloud locations with a label value of 1 are retained as real cloud signals for subsequent processing. A Bayesian inference algorithm for cloud boundaries based on maximum likelihood estimation is used to accurately estimate the cloud base height of the real cloud signals. The algorithm models the laser echo photon count as a Poisson distribution process, with the Poisson distribution parameter representing the real backscattering intensity. The prior distribution of cloud base height for the observation area is extracted from a historical climate statistics database, and the prior distribution is represented by a Gaussian mixture model to indicate the probability distribution of cloud base height under different weather types. The Bayesian inference process combines the prior distribution with the currently observed real cloud signals and uses the Bayesian formula to calculate the posterior probability distribution of the cloud base height. The mean of the posterior probability distribution is used as the cloud base height estimate, and the standard deviation is used as the posterior uncertainty interval, which reflects the reliability of the estimation result. The two-layer game-theoretic optimization model is executed after obtaining a preliminary cloud base height estimate. The upper-layer model aims to maximize the cloud identification confidence, while the lower-layer model aims to minimize the registration error index. The upper-layer model's input parameters include the peak value of the posterior probability distribution, the signal-to-noise ratio (SNR), and the number of identified clouds. Its output parameter is the cloud identification confidence. The objective function is calculated by multiplying the posterior probability peak value normalization term, the logarithm of the SNR term, and the square root of the cloud number term. The lower-layer model's input parameters include cloud movement speed, time synchronization error, and spatial spacing. Its output parameter is the registration error index. Its objective function is calculated by summing the product of cloud movement speed and time synchronization error, and the square of the spatial spacing term. The two objective functions are linked by a coupling term, which is the product of the SNR and cloud movement speed. The value of this coupling term affects the identification confidence of the upper-layer model and the registration error of the lower-layer model, forming a two-way feedback mechanism.The Stackelberg game algorithm is used for iterative optimization. The upper-level model, as the leader, first calculates the recognition confidence based on the current signal-to-noise ratio and the number of clouds. The lower-level model, as the follower, adjusts the spatiotemporal interpolation parameters to optimize the registration error based on the recognition confidence output by the upper-level model. The iteration loop continues until the relative change of the two objective functions is less than 1%, at which point convergence is determined. The cloud recognition confidence and registration error indices after convergence are used for subsequent data quality assessment.

[0052] The specific implementation of step S04 involves calculating the cloud movement speed using the optical flow method. The principle of the optical flow method is to estimate the motion vector by analyzing the positional changes of identical feature points in consecutive frame images. Three consecutive frames of laser echo intensity images are selected, with a 60-second time interval between frames. Feature points at the cloud boundary are identified in the first frame image. The selection criterion for feature points is pixels with an echo intensity gradient greater than a threshold of 100. In the second and third frames, a template matching algorithm is used to search for the corresponding positions of the feature points. The template size is 11 pixels multiplied by 11 pixels, and the matching criterion is maximizing the normalized cross-correlation coefficient. The positional offset of the feature points between adjacent frames is divided by the time interval to obtain the cloud movement speed. The median of the movement speeds of all feature points is taken as the estimated cloud movement speed at that moment. The median filtering method can effectively suppress the influence of individual abnormal feature points. When the cloud movement speed exceeds 5 meters per second, the cloud is determined to be in a rapid movement state. At this time, the cloud observed by the reference instrument and the tested instrument may not be the same cloud body, requiring spatiotemporal interpolation correction to eliminate registration errors. The cloud motion vector estimation algorithm further calculates the vertical motion component based on the cloud movement speed. This is done by analyzing the rate of change of cloud echo intensity over time in the vertical direction. A positive rate of change indicates cloud base rise, while a negative rate indicates cloud base descent. The magnitude of the vertical motion speed is equal to the ratio of the vertical gradient of the echo intensity to the rate of change over time. Spatiotemporal interpolation correction extrapolates the cloud position from the reference instrument's observation time to the expected position at the test instrument's observation time. The extrapolation method adds the cloud base height measured by the reference instrument to the product of the vertical component of the cloud motion vector and the time difference. The time difference is the test instrument's observation time minus the reference instrument's observation time, typically within the range of -60 seconds to +60 seconds. The spatiotemporally interpolated reference instrument cloud base height data and the test instrument's cloud base height data are paired. The pairing principle is that data with the same timestamp or a time difference of less than 10 seconds form data pairs, creating paired datasets. The Pearson correlation coefficient of the paired datasets is calculated. The correlation coefficient reflects the degree of linear correlation between the measurements of the two instruments and is calculated by dividing the covariance of the paired data by the product of the standard deviations of the two datasets. The slope of the linear regression was calculated using the least squares method, with the cloud base height data of the reference instrument as the independent variable and the cloud base height data of the test instrument as the dependent variable. A univariate linear equation was fitted, and the slope of the linear regression reflected the consistency of the measurement scales of the two instruments. The judgment criterion was that if the correlation coefficient was greater than 0.95 and the linear regression slope was within the range of 0.95 to 1.05, the measurement results of the test instrument and the reference instrument were considered to have good consistency, and the test instrument was deemed to have passed the comparison test. Otherwise, the comparison test was deemed unqualified, and instrument calibration or re-comparison testing was required.

[0053] The specific implementation of step S05 involves a full-range performance evaluation of the maximum cloud base height data identified at all observation times throughout the entire comparative test. The evaluation aims to verify the linearity and long-term stability of the tested instrument across the entire measurement range. The maximum cloud base height data of the reference instrument and the tested instrument at 720 observation times during the comparative test are arranged in ascending order of numerical value, forming a maximum cloud base height sequence for the reference instrument and a maximum cloud base height sequence for the tested instrument. Using the minimum value of the two sequences as the reference point, the incremental value for each observation point is calculated. The incremental value is defined as the cloud base height at the current observation point minus the minimum value, reflecting the change in cloud height relative to the reference point. The incremental deviation for each observation point is calculated. The incremental deviation is defined as the incremental value of the reference instrument minus the incremental value of the tested instrument. A positive incremental deviation indicates that the tested instrument's measurement is too low, and a negative incremental deviation indicates that the tested instrument's measurement is too high. The number of data points with an absolute incremental deviation less than 5% of the instrument's range was counted. The instrument's range is 7500 meters minus 30 meters equals 7470 meters; 5% of the instrument's range equals 373.5 meters. When the proportion of data points with an absolute incremental deviation less than 373.5 meters exceeds 90% of the total data points, the tested instrument is deemed to have passed the full-range performance test. This 90% criterion is based on the requirements for cloud height measurement accuracy in marine meteorological observation operational specifications. The full-range performance evaluation results reflect the systematic and linear errors of the tested instrument under different cloud height conditions. By analyzing the incremental deviation distribution characteristics, the existence of fixed deviations or range-related deviations in the instrument is determined, providing a quantitative basis for subsequent instrument calibration and maintenance.

[0054] It should be noted that the key technical ideas of this invention include an aerosol cloud separation method based on wavelet multi-scale decomposition and a two-parameter discriminant matrix, a cloud boundary uncertainty quantification method based on Bayesian inference, and a spatiotemporal registration and correction method driven by a two-layer game optimization model. The aerosol cloud separation method based on wavelet multi-scale decomposition and a two-parameter discriminant matrix decomposes the laser echo signal at different spatial scales to extract feature information of clouds of different thicknesses. It then constructs a two-dimensional discriminant space by combining depolarization ratio and color ratio parameters. Compared with traditional single-threshold discriminant methods, this method can simultaneously utilize particle polarization characteristics and spectral scattering characteristics to achieve accurate differentiation between aerosols and clouds, effectively reducing the false cloud detection rate under high-concentration aerosol conditions and improving the accuracy of cloud identification in multi-layer thin cloud scenarios. The Bayesian inference-based cloud boundary uncertainty quantification method models the randomness of laser echo photon counting as a Poisson process, integrating historical climate statistics prior knowledge with real-time observation data. Through posterior probability distribution, it not only provides an estimate of cloud base height but also quantifies the uncertainty range of the estimation result. Compared to traditional fixed threshold methods that only output a single value, the Bayesian method provides measurement confidence information, enabling the comparison results to have credibility assessment capabilities. It is particularly suitable for robust identification of weak signal clouds and complex multi-layered cloud scenarios. The two-layer game-theoretic optimization model-driven spatiotemporal registration correction method models the optimization of cloud identification accuracy and the minimization of data registration error as a game relationship between upper and lower layers. It achieves synergy between the two optimization objectives through the coupling term of the signal-to-noise ratio and cloud movement speed. Compared to traditional methods that independently perform identification and registration, the game-theoretic model can dynamically adjust registration parameters according to the current signal quality, maintaining a balance between identification accuracy and registration precision under rapid cloud movement conditions, significantly improving the reliability of comparison measurements under dynamic weather conditions. The synergistic effect of the three technical approaches lies in the wavelet multi-scale decomposition providing high-quality cloud candidate signals, the Bayesian inference method achieving accurate boundary estimation and quantifying uncertainty based on this, and the two-layer game model optimizing the spatiotemporal registration strategy based on uncertainty information, forming a complete technical chain from signal preprocessing, feature extraction, boundary recognition to registration correction. Compared with the existing technology that processes each link independently, this invention achieves intelligent and accurate comparison of the entire process of cloud height measurement through multi-level information fusion and feedback optimization mechanism, and has significant advantages in comparison accuracy and stability under complex atmospheric conditions.

[0055] It should be noted that this invention also solves the following technical problem: During the comparison measurement of cloud height instruments, the laser backscattering signal of the aerosol layer is easily misidentified as a cloud signal, leading to the inclusion of false cloud data in the comparison results and affecting the accuracy assessment. This invention calculates the depolarization ratio and color ratio parameters of the laser echo signal, and utilizes the physical differences between cloud droplets and aerosol particles in polarization and wavelength scattering characteristics to establish an aerosol scattering feature discrimination matrix to achieve automatic classification of target types. Candidate positions with a depolarization ratio less than 0.15 and a color ratio parameter greater than 1.8 are marked as aerosol layers and eliminated, thus retaining the true cloud signal for subsequent comparison analysis. In addition, traditional cloud height instrument comparison methods only provide point estimates of cloud base height, which cannot quantify measurement uncertainty and make it difficult to assess the confidence level of the comparison results. This invention adopts a Bayesian inference algorithm for cloud boundaries based on maximum likelihood estimation, constructs a prior distribution of cloud base height using historical climate statistics, and calculates a posterior probability distribution by combining it with the true cloud signal. The standard deviation of the posterior probability distribution is used as the posterior uncertainty interval, providing a reliable basis for confidence assessment of the comparison results.

[0056] Specifically, the principle of this invention is as follows: The solution to this technical problem lies in establishing a complete technical chain for cloud dynamic characteristic perception and spatiotemporal registration correction. First, the nanosecond-level time reference provided by the BeiDou time service system ensures strict synchronization of the observation times of the two instruments, eliminating sources of time-dimensional deviation. Second, wavelet multi-scale decomposition can extract cloud features at different spatial scales. Depolarization ratio and color ratio parameters distinguish cloud droplets from aerosol particles from a physical scattering mechanism. The aerosol scattering feature discrimination matrix achieves rapid and accurate target classification through a two-dimensional lookup table, ensuring reliable identification of cloud signals. Third, the cloud motion vector estimation algorithm calculates the moving speed by tracking cloud feature points using optical flow. Spatiotemporal interpolation correction extrapolates the cloud position at the reference instrument's observation time to the subject instrument's observation time based on the cloud motion vector, achieving spatial-dimensional registration correction. Fourth, the two-layer game optimization model couples cloud identification and registration error optimization through a Stackelberg game framework. The upper-layer model optimizes the identification confidence, while the lower-layer model responds and optimizes the registration error. Iterative solution ensures that both objectives reach their optimal values ​​simultaneously.

[0057] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0058] The specific implementation of step S01 involves installing a reference instrument and the tested instrument in a standard location, maintaining a distance of at least 3 meters between them. The optical axes are paralleled using base adjustment bolts, and nanosecond-level time synchronization is achieved using the BeiDou time synchronization system. The BeiDou time synchronization system receives Coordinated Universal Time (UTC) signals broadcast by BeiDou satellites, providing a unified time reference for both the reference instrument and the tested instrument. The timing accuracy is better than 50 nanoseconds, ensuring negligible time deviations in the synchronized observation data. The ranging accuracy of the tested instrument against fixed targets is verified. Five fixed-height target plates are set up within a range of 10 to 2000 meters. Measurements are repeated 10 times for each fixed-height target plate, and the deviation sequence between the measured distance and the actual distance is recorded. The fixed-height target plates are aluminum alloy flat plates with a diffuse reflection coating, measuring 1 meter by 1 meter, installed on support rods at different heights, with the plate surface perpendicular to the laser emission direction. The formula for calculating the average deviation is as follows:

[0059] ;

[0060] In the formula, This is the mean deviation, in meters. The empirical value for the number of measurements is 10. For the first The deviation value of this measurement, in meters. This is the index for the number of measurements, with values ​​ranging from 1 to... . The calculation formula is: ,in For the first The distance measurement value is in meters. The distance is the actual distance, in meters. The formula for calculating the standard deviation is as follows:

[0061] ;

[0062] In the formula, represents the standard deviation of the deviation, in meters. When Proceed to the next step, in which The instrument's nominal accuracy is measured in meters. For cloud heights less than 1000 meters, the accuracy is typically 10% of the measured value, and for cloud heights greater than 1000 meters, it is 20%. Otherwise, the instrument under test should be calibrated and the measurement repeated.

[0063] The specific implementation of step S02 involves starting both the reference instrument and the tested instrument for continuous observation for 12 hours, with a sampling frequency of 1 minute per observation, simultaneously recording the raw laser echo signal data and manually observed cloud-like data. Wavelet multi-scale decomposition is performed on the raw laser echo signal data, using the Daubechies wavelet basis function for multi-level decomposition. The formula for calculating the spatial scale corresponding to layer wavelet decomposition is as follows:

[0064] ;

[0065] In the formula, For the first The spatial scale corresponding to each layer, in meters. This serves as the baseline spatial scale, typically taken as 100 meters. The wavelet decomposition level is denoted as . This represents the range resolution of the lidar, measured in meters, typically 15 meters. Layer 3 corresponds to a spatial scale of 200 to 400 meters, and layer 5 corresponds to a spatial scale of 800 to 1600 meters. A peak detection algorithm is applied to the wavelet coefficients from layers 3 to 5 to extract candidate cloud locations. The peak detection algorithm's criterion is that when the value at a certain point... When the value is greater than twice the average of the 20 adjacent sampling points, the point is determined to be a candidate cloud location. The determination formula is as follows:

[0066] ;

[0067] In the formula, For the first The first layer of wavelet coefficient sequence The values ​​at each sampling point are dimensionless. For sampling point index, The baseline wavelet coefficient value is usually set to 1. Indices are used for adjacent sampling points. The depolarization ratio and color ratio parameters are calculated for each candidate cloud location. The formula for calculating the depolarization ratio is as follows:

[0068] ;

[0069] In the formula, This refers to the depolarization ratio, in dimensionless units. The intensity of vertically polarized light is expressed in watts. The intensity of parallel polarized light is expressed in watts. The depolarization ratio of cloud droplets is typically less than 0.05, while that of aerosol particles ranges from 0.20 to 0.35. The formula for calculating the color ratio parameter is as follows:

[0070] ;

[0071] In the formula, These are color ratio parameters, in dimensionless units. The backscattering coefficient at a wavelength of 532 nanometers is expressed in meters per steradian. This is the backscattering coefficient at a wavelength of 1064 nanometers, expressed in meters per steradian. The color ratio parameter of water clouds is close to 1, while that of aerosol layers is typically greater than 2.

[0072] The specific implementation of step S03 involves establishing an aerosol scattering feature discrimination matrix. This matrix is ​​a two-dimensional lookup table, with the horizontal axis representing the depolarization ratio divided into 50 intervals from 0 to 0.5, and the vertical axis representing the color ratio parameter divided into 50 intervals from 0.5 to 3.0. (Matrix row index) The calculation formula is expressed as follows:

[0073] ;

[0074] In the formula, This represents the matrix row index corresponding to the depolarization ratio, with values ​​ranging from 1 to 50. This represents the minimum depolarization ratio; the default value is 0. This represents the maximum depolarization ratio; the default value is 0.5. This is a floor function. Matrix column index. The calculation formula is expressed as follows:

[0075] ;

[0076] In the formula, This represents the matrix column index corresponding to the color ratio parameter, with values ​​ranging from 1 to 50. This is the minimum value for the color ratio parameter; the default value is 0.5. This represents the maximum color ratio parameter; the default value is 3.0. Matrix elements. Store the target type label corresponding to the combination of depolarization ratio and color ratio parameters, when and hour, Indicates an aerosol layer, otherwise Representing the actual cloud layer. Candidate cloud locations marked as aerosol layers are removed, retaining only the actual cloud layer signal. A Bayesian inference algorithm for cloud boundaries based on maximum likelihood estimation is used to process the actual cloud layer signal, modeling the laser echo photon count as a Poisson distribution. The probability density function of the Poisson distribution is expressed as follows:

[0077] ;

[0078] In the formula, For given parameters Observed at time The probability of a photon, in dimensionless form. This is a count of laser echo photons, expressed in units of individual photons. The value represents the actual backscattering intensity, expressed in units of [number of particles]. The base of the natural logarithm is approximately 2.71828. A priori distribution of cloud base height is constructed using historical climate statistics. The prior distribution adopts a Gaussian distribution form, and is described as follows:

[0079] ;

[0080] In the formula, This represents the prior probability density of cloud base height, expressed in meters. This refers to the cloud base height, in meters. This is the a priori distribution mean, expressed in meters, and determined based on historical climate statistics. The standard deviation of the prior distribution, in meters, is determined based on historical climate statistics. Pi, approximately 3.14159. The posterior probability distribution is calculated using Bayes' theorem combined with real cloud data. The Bayes' theorem is expressed as follows:

[0081] ;

[0082] In the formula, Count a given number of photons The posterior probability density of cloud base height, in meters. Let be the likelihood function for a given cloud base height, in dimensionless units. The mean of the posterior probability distribution is used as the estimate of the cloud base height. The unit is meters, and the standard deviation of the posterior probability distribution is used as the posterior uncertainty interval. The unit is meters. A two-layer game-theoretic optimization model is constructed to improve cloud layer identification accuracy and minimize data registration error. The objective function of the upper-layer model is expressed as follows:

[0083] ;

[0084] In the formula, The confidence level for cloud identification is dimensionless. This is the baseline confidence level, typically set to 1. The peak value of the posterior probability distribution is given, and the unit is dimensionless. This is the baseline posterior probability value, typically taken as 0.8. Signal-to-noise ratio, expressed in decibels. The reference signal-to-noise ratio is typically set to 20 dB. The number of cloud layers identified, in layers. The baseline number of cloud layers is typically set to 3. The objective function of the lower-level model is expressed as follows:

[0085] ;

[0086] In the formula, The registration error index is expressed in meters. The reference registration error is typically taken as 10 meters. The speed of cloud movement is measured in meters per second. The baseline cloud movement speed is typically taken as 5 meters per second. This represents the time synchronization error, expressed in seconds. The reference time synchronization error is typically taken as a value of [value missing]. Second, This refers to the spatial spacing, in meters, with a default value of 3 meters. The baseline spatial spacing is typically taken as 5 meters. The two objective functions are related through a coupling term, which is described as follows:

[0087] ;

[0088] In the formula, The values ​​for the coupling terms are expressed in decibels per second. This is the baseline coupling term value, typically taken as 100 dBm / s. The upper-level model constraints are: and The constraints of the lower-level model are: Miqie Seconds. Through iterative solution using Stackelberg game, the upper-level model first determines the confidence level for cloud identification, and the lower-level model optimizes the registration error index based on the cloud identification confidence level. The iteration converges when the changes in both objective functions are less than 1%.

[0089] The specific implementation of step S04 involves using a cloud motion vector estimation algorithm to calculate the cloud movement speed, and employing optical flow to calculate the horizontal cloud movement speed using three consecutive echo images. The optical flow method selects cloud feature points in the previous frame and searches for the corresponding point with the closest grayscale value in the next frame. The formula for calculating the grayscale matching error is as follows:

[0090] ;

[0091] In the formula, This represents the grayscale matching error, expressed in dimensionless units. For a moment The image grayscale value function, with dimensionless units. For a moment The image grayscale value function, with dimensionless units. These are the coordinates of the feature points in the previous frame, in pixels. The coordinates of the search point in the next frame, in pixels. and The offset within the window. This is the radius of the search window, typically set to 5 pixels. Selecting this will... smallest As a corresponding point, the formula for calculating cloud movement speed is as follows:

[0092] ;

[0093] In the formula, The difference in horizontal position of the feature points is calculated using the following formula: The unit is meters. The vertical position difference of the feature points is calculated using the following formula: The unit is meters. This is a time interval, in seconds. Spatiotemporal interpolation correction was performed on the multi-layer cloud data identified by the reference instrument and the tested instrument at the same time, using meters per second. The corrected height calculation formula is as follows:

[0094] ;

[0095] In the formula, This is the corrected height, in meters. The observation height is in meters. This represents the vertical component of the cloud motion vector, expressed in meters per second. The time difference is expressed in seconds. Calculate the correlation coefficient and linear regression slope between the spatiotemporally interpolated cloud base height data of the reference instrument and the cloud base height data of the tested instrument. The mean of the cloud base height data of the reference instrument is also included. The calculation formula is expressed as follows:

[0096] ;

[0097] In the formula, The reference instrument's cloud base height data is the average value, in meters. For the number of observations, For reference instrument number Sub-cloud base height data, in meters. This is the observation count index, with values ​​ranging from 1 to... Mean cloud base height data of the tested instruments The calculation formula is expressed as follows:

[0098] ;

[0099] In the formula, The mean cloud base height data of the tested instrument is expressed in meters. For the tested instrument Sub-cloud base height data, in meters. The correlation coefficient calculation formula is as follows:

[0100] ;

[0101] In the formula, The correlation coefficient is dimensionless. The slope of the linear regression is obtained by least squares fitting, and the calculation formula is expressed as follows:

[0102] ;

[0103] In the formula, The slope is the linear regression slope, in dimensionless units. and The instrument under test is deemed to be qualified upon comparison.

[0104] The specific implementation of step S05 involves sorting all the maximum cloud height values ​​identified during the comparative measurement process, and calculating the incremental deviation of each point based on the minimum value. The formula for calculating the incremental deviation is as follows:

[0105] ;

[0106] In the formula, For the first The incremental deviation of the maximum cloud height, in meters. For the maximum cloud height index, For reference instrument number The maximum cloud height, in meters. This is the minimum of all maximum cloud height values ​​from the reference instrument, in meters. For the tested instrument The maximum cloud height, in meters. This is the minimum value among all maximum cloud height values ​​of the tested instrument, in meters. When the condition is met... The instrument under test was deemed to have qualified performance across its entire range at that time. To meet Number of data points The total number of all data points. The instrument's measuring range is in meters, and the calculation formula is: ,in This is the upper limit of the cloud height measurement range, typically 7500 meters. This is the lower limit of the cloud height measurement range, usually 30 meters.

[0107] The formula for calculating the signal-to-noise ratio is as follows:

[0108] ;

[0109] In the formula, Effective signal power, measured in watts. The noise power is expressed in watts. The Stackelberg game is a dynamic game model where the upper-level model, acting as the leader, makes the first decision, and the lower-level models, acting as followers, respond by observing the upper-level model's decision. The equilibrium point is determined through back induction.

[0110] To better understand and implement this invention, the following is a specific application scenario example 2: A marine meteorological observation team conducted a cloud height measurement comparison at a marine observation station in the South China Sea. The goal was to verify the performance of three newly purchased test instruments to ensure they met the requirements for operational marine meteorological observation. The observation station is equipped with one reference instrument, model CL51, which has passed national metrological verification. Its cloud height measurement range is 30 meters to 7500 meters, and its measurement accuracy is no more than 10% error when the cloud height is less than 1000 meters and no more than 20% error when the cloud height is greater than 1000 meters. The three test instruments, models CYY-1, CYY-2, and CYY-3, are all domestically produced laser cloud height instruments with the same technical specifications as the reference instrument.

[0111] The technical team first selected a suitable comparison testing site at the observation station. The site was located in an open area to the north of the station, with no buildings or trees taller than 10 meters obstructing the view within 30 meters. A reference instrument was installed at the center of the site. Test instrument CYY-1 was installed 3.5 meters east of the reference instrument, CYY-2 4 meters south, and CYY-3 3.8 meters west, forming a cross-shaped arrangement. Each instrument underwent optical axis parallel calibration using the adjusting bolts on its tripod base. A laser aiming target was set up 200 meters away from the instruments, and the optical axes of each instrument were adjusted until the laser beam converged at the center of the target, with the optical axis parallelism deviation controlled within 0.08 degrees. All instruments received satellite time signals via the BeiDou time synchronization module, achieving a time synchronization accuracy of 45 nanoseconds, ensuring consistent time reference across the four instruments.

[0112] In the fixed-target ranging accuracy verification phase, the technical team set up five fixed-height target plates along a straight line on the north side of the test site. Each target plate was a 1-meter by 1-meter aluminum alloy flat plate coated with a diffuse reflection coating, achieving a reflectivity of 88%. The heights of the five target plates were 10 meters, 480 meters, 980 meters, 1520 meters, and 1980 meters, respectively. A calibrated laser rangefinder with an accuracy of 0.5 millimeters was used to measure the actual distances to each of the five target plates. The tested instrument, CYY-1, measured each target plate 10 times, and the ranging results are shown in Table 1.

[0113] Table 1. Distance measurement results of the CYY-1 instrument on a fixed target.

[0114]

[0115] Based on the data in Table 1, the deviation sequence was calculated. The standard deviations of the deviations at the five test points were 0.08 m, 0.42 m, 0.95 m, 1.52 m, and 2.18 m, respectively, corresponding to the nominal accuracies of the instruments: 1 m, 48 m, 98 m, 304 m, and 396 m. Since the standard deviations were all less than 50% of the nominal accuracies, the CYY-1 ranging system under test was determined to be functioning normally. The test results for the CYY-2 and CYY-3 instruments were similar, with standard deviations ranging from 0.09 m to 2.35 m and from 0.07 m to 2.12 m, respectively, both meeting the conditions for entering the comparative testing process.

[0116] The continuous observation phase lasted 12 hours, from 8:00 AM to 8:00 PM, with a sampling frequency of once per minute. A total of 720 sets of observation data were collected by the four instruments. Weather conditions during the observation period were cloudy to overcast, with visibility ranging from 15 to 20 kilometers, wind speeds from 3 to 7 meters per second, temperatures from 26°C to 31°C, and relative humidity from 62% to 78%. Human observers recorded cloud cover every 30 minutes, for a total of 24 records, including 18 records with clouds and 6 records without clouds. The raw laser echo signal data was recorded using a high-speed data acquisition card at a sampling rate of 100 MHz. A photon count time series was recorded within 15 microseconds after each laser pulse emission, resulting in a raw data file size of approximately 86 gigabytes.

[0117] The technical team used wavelet multi-scale decomposition to preprocess the laser echo signal, selecting the Daubechies4 wavelet basis function for a 5-level decomposition. For example... Figure 2 As shown, the third-layer wavelet coefficients mainly reflect cloud signals with a thickness of 200 to 400 meters, the fourth-layer wavelet coefficients reflect cloud signals with a thickness of 400 to 800 meters, and the fifth-layer wavelet coefficients reflect cloud signals with a thickness of 800 to 1600 meters. A peak detection algorithm was applied to the third to fifth-layer wavelet coefficients, extracting a total of 1856 candidate cloud locations, averaging 2.58 candidate cloud locations per set of observation data. The depolarization ratio and color ratio parameters were calculated for each candidate cloud location; the depolarization ratio ranged from 0.03 to 0.42, and the color ratio parameter ranged from 0.78 to 2.65.

[0118] The aerosol scattering feature discrimination matrix divides the depolarization ratio from 0 to 0.5 into 50 intervals and the chromaticity ratio parameter from 0.5 to 3.0 into 50 intervals, forming 2500 grid cells. According to the discrimination criteria, 327 candidate cloud locations with a depolarization ratio less than 0.15 and a chromaticity ratio greater than 1.8 were identified as aerosol layers and discarded, retaining 1529 true cloud signals. Figure 3 As shown, the candidate cloud locations after removing the aerosol layer are mainly distributed in the region with a depolarization ratio of 0.03 to 0.12 and a color ratio parameter of 0.85 to 1.65, which is consistent with the scattering characteristics of water clouds.

[0119] A Bayesian inference algorithm for cloud boundaries based on maximum likelihood estimation was used to process 1529 real cloud signals. A historical climate statistics database provided the prior distribution of cloud base heights for the observation area. The prior distribution was a mixture model of three Gaussian distributions, with peak values ​​at 580 meters, 1250 meters, and 2100 meters, corresponding to typical heights of low, middle, and high clouds, respectively. Bayesian inference calculated the posterior probability distribution of the cloud base height for each layer. The mean of the posterior probability distribution was used as the cloud base height estimate, and the standard deviation was used as the posterior uncertainty interval. The reference instrument identified 142 sets of observation data with three cloud layers, 385 sets with two cloud layers, and 193 sets with one cloud layer. The tested instrument CYY-1 identified 94.3% of the cloud layers with a difference of no more than one layer from the reference instrument, and 87.6% of the cloud base heights differed by less than 50 meters.

[0120] A two-layer game-theoretic optimization model collaboratively optimizes cloud identification and data registration. The upper-layer model calculates the confidence level of cloud identification, with input parameters including the peak value of the posterior probability distribution, the signal-to-noise ratio (SNR), and the number of clouds. During the observation period, the SNR ranged from 18.5 dB to 32.7 dB, with an average of 24.3 dB. Figure 4 As shown, the confidence level for cloud identification increases with the increase of the posterior probability peak value and the signal-to-noise ratio. When the posterior probability peak value is greater than 0.75 and the signal-to-noise ratio is greater than 22 dB, the confidence level for cloud identification exceeds 0.90, indicating that the identification result is reliable. The lower-level model calculates the registration error index, with input parameters including cloud movement speed, time synchronization error, and spatial spacing. The cloud movement speed calculated by the optical flow method ranges from 1.2 m / s to 8.7 m / s, with 89 instances exceeding 5 m / s, accounting for 12.4% of the total observation times. The time synchronization error is controlled within a certain range using the BeiDou time synchronization system. The spatial spacing is 3.5 to 4 meters. The registration error range is 2.3 to 15.8 meters. When the cloud movement speed is less than 5 meters per second, the registration error is less than 8 meters. When the cloud movement speed is greater than 5 meters per second, the registration error increases to 10 to 16 meters, requiring spatiotemporal interpolation correction.

[0121] In the iterative solution of the Stackelberg game, the upper-level model first determines the cloud identification confidence level, and the lower-level model adjusts the spatiotemporal interpolation parameters based on the identification confidence level. Initially, the registration error index was 12.5 meters. After three iterations, the registration error index decreased to 5.8 meters; after four iterations, it was 5.6 meters, a change of 3.4%; after five iterations, it was 5.5 meters, a change of 1.8%; and after six iterations, it was 5.4 meters, a change of 0.9%, satisfying the convergence condition. The coupling term value increased from the initial 96.5 dBm / s to 112.3 dBm / s at convergence, reflecting the synergistic optimization effect of the signal-to-noise ratio and cloud movement speed.

[0122] Spatiotemporal interpolation correction processed 89 moments where cloud movement speeds exceeded 5 meters per second. The cloud motion vector estimation algorithm calculated the vertical motion speed range from -0.8 meters per second to +1.2 meters per second, with negative values ​​indicating cloud base descent and positive values ​​indicating cloud base elevation. For example... Figure 5 As shown, after spatiotemporal interpolation correction, the correlation coefficient between the cloud base height data of the reference instrument and the tested instrument CYY-1 increased from 0.923 to 0.967, and the linear regression slope was adjusted from 1.082 to 1.012. This meets the criteria of a correlation coefficient greater than 0.95 and a linear regression slope within the range of 0.95 to 1.05, thus the tested instrument CYY-1 is deemed qualified for the comparison test. The correlation coefficients of the tested instruments CYY-2 and CYY-3 are 0.971 and 0.965, respectively, and the linear regression slopes are 1.008 and 0.997, respectively, both meeting the qualification criteria for the comparison test.

[0123] During the full-range performance evaluation phase, the technical team extracted maximum cloud base height data from 720 observation times. The maximum cloud base height range for the reference instrument was 285 meters to 3420 meters, while the maximum cloud base height range for the tested instrument CYY-1 was 278 meters to 3385 meters. Incremental values ​​were calculated based on the minimum value. The increment range for the reference instrument was 0 meters to 3135 meters, and the increment range for the tested instrument CYY-1 was 0 meters to 3107 meters. A total of 672 data points had an absolute increment deviation of less than 373.5 meters, accounting for 93.3% of the total data points, exceeding the 90% threshold. Therefore, the tested instrument CYY-1 was deemed to have passed the full-range performance evaluation. The passing data point percentages for the tested instruments CYY-2 and CYY-3 were 91.8% and 92.6%, respectively, both meeting the full-range performance requirements.

[0124] The advancement of this invention over traditional cloud height measurement methods lies in addressing the core shortcomings of traditional methods through multi-level information fusion and dynamic optimization mechanisms. Traditional methods use fixed thresholds to determine cloud boundaries, which easily leads to false cloud detections when aerosol concentrations are high. This invention, however, utilizes a two-parameter discrimination matrix combining wavelet multi-scale decomposition with depolarization ratio and color ratio parameters. By leveraging multi-dimensional information from particle scattering characteristics, it achieves accurate differentiation between aerosols and clouds, solving the misjudgment problem from a physical mechanism perspective. Traditional methods only output a single cloud base height value, failing to assess the reliability of the measurement results. This invention employs a Bayesian inference framework to integrate historical climate statistics prior knowledge with real-time observation data. Through posterior probability distribution, it not only provides an estimated cloud base height but also quantifies the uncertainty interval, providing a basis for confidence assessment of the comparison results. Traditional methods employ fixed spatiotemporal registration parameters under conditions of rapid cloud movement, leading to a significant increase in registration error as cloud movement speed increases. This invention models the minimization of cloud identification accuracy and registration error as a collaborative optimization problem through a two-layer game optimization model. It utilizes the coupling mechanism between signal-to-noise ratio and cloud movement speed to achieve dynamic adaptive registration, minimizing registration error while maintaining identification accuracy, thus fundamentally improving the reliability of comparative measurements under dynamic weather conditions.

[0125] It should be noted that the variables involved in this invention are explained in detail in Tables 2 and 3.

[0126] Table 2. Variable Explanation Table (Part 1)

[0127]

[0128] Table 3. Variable Explanation Table (Part Two)

[0129]

[0130] 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 changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A precise method for comparing ceilometer heights, characterized in that, Reference instruments and test instruments were installed in a standard field and synchronized at the nanosecond level using the BeiDou time synchronization system. The accuracy of the test instrument in measuring fixed targets was verified. Specifically, five fixed-height target boards were set up within the measurement range, and each fixed-height target board was measured 10 times. The standard deviation of the deviation was calculated. When the standard deviation of the deviation was less than 50% of the instrument's nominal accuracy, the reference instrument and the test instrument were started to continuously observe and synchronously record the raw data of the laser echo signal. Wavelet multi-scale decomposition was performed on the raw data of the laser echo signal to extract candidate cloud positions and calculate the depolarization ratio and color ratio parameters. An aerosol scattering feature discrimination matrix was established to remove the aerosol layer and retain the real cloud signal. A Bayesian inference algorithm based on maximum likelihood estimation was used to obtain the posterior probability distribution of cloud base height. A two-layer game optimization model was constructed to improve the accuracy of cloud identification and minimize the data registration error. A cloud motion vector estimation algorithm was used to calculate the cloud movement speed and perform spatiotemporal interpolation correction. The correlation coefficient and linear regression slope were calculated to determine the pass rate of the comparison test. Among them, the cloud boundary Bayesian inference algorithm based on maximum likelihood estimation models the laser echo photon count as a Poisson distribution, constructs the prior distribution of cloud base height through historical climate statistics, and calculates the posterior probability distribution using the Bayesian formula; The aerosol scattering feature discrimination matrix is ​​a two-dimensional lookup table, with the horizontal axis representing the depolarization ratio and the vertical axis representing the color ratio parameter. Candidate cloud locations with a depolarization ratio less than 0.15 and a color ratio parameter greater than 1.8 are marked as aerosol layers and removed. The two-layer game optimization model includes an upper-layer model that aims to maximize cloud identification accuracy and a lower-layer model that aims to minimize data registration error. It is solved iteratively through Stackelberg game. The upper-layer model first determines the cloud identification confidence level, and the lower-layer model optimizes the registration error index based on the cloud identification confidence level. The model converges when the change in both objective functions is less than 1%. The objective function of the upper-level model is expressed as follows: ； In the formula, The confidence level for cloud identification is dimensionless. As the baseline confidence level, The peak value of the posterior probability distribution is given, and the unit is dimensionless. The baseline posterior probability value, Signal-to-noise ratio, expressed in decibels. As a reference signal-to-noise ratio, The number of cloud layers identified, in layers. The baseline number of cloud layers; The objective function of the lower-level model is expressed as follows: ； In the formula, The registration error index is expressed in meters. For the reference registration error, The speed of cloud movement is measured in meters per second. As a baseline cloud movement speed, This represents the time synchronization error, expressed in seconds. For the reference time synchronization error, Spatial distance, in meters. The reference spatial spacing.

2. The accurate astronomical height comparison method according to claim 1, characterized in that, The fixed-height target plate is an aluminum alloy plate with a diffuse reflection coating, measuring 1 meter by 1 meter, installed on support rods at different heights, with the plate surface perpendicular to the laser emission direction.

3. The accurate astronomical height comparison method according to claim 2, characterized in that, The wavelet multi-scale decomposition uses Daubechies wavelet basis functions for multi-level decomposition, and a peak detection algorithm is applied to extract candidate cloud locations in the wavelet coefficients of the 3rd to 5th levels.

4. The accurate astronomical height comparison method according to claim 3, characterized in that, The peak detection algorithm searches for local maxima in each wavelet coefficient sequence. When the value of a point is more than twice the average of the 20 adjacent sampling points, the point is determined to correspond to a candidate cloud location.

5. The accurate astrolabe height comparison method according to claim 4, characterized in that, The depolarization ratio is the ratio of the intensity of the received vertically polarized light to the intensity of the parallelly polarized light, and the color ratio parameter is the ratio of the backscattering coefficients of the lidar at a wavelength of 532 nm and 1064 nm.

6. The accurate astronomical height comparison method according to claim 5, characterized in that, The cloud motion vector estimation algorithm uses optical flow to calculate the horizontal movement speed of the cloud layer for three consecutive frames of echo images. When the cloud layer movement speed exceeds 5 meters per second, spatiotemporal interpolation correction is performed.

Images