Terahertz cloud measurement radar calibration method based on channel perception

Abstract

The invention relates to the technical field of radar reflectivity factor calibration, and discloses a terahertz cloud measurement radar calibration method based on channel perception, which comprises the following steps: S1, acquiring two types of input data; s2, constructing a multi-scale channel model; s3, constructing a distance error estimation model; s4, constructing a dynamic time correction model; and S5, constructing a calibration module. According to the method, a high-precision physical channel model is constructed to independently predict path attenuation of signals, the concepts of gain deviation and path attenuation factors are put forward, the influence of thermal noise attenuated by radar on radar waves in the atmospheric environment is reflected respectively, and a Kalman filtering algorithm is used to estimate the path attenuation of the signals. And filtering errors generated in the time direction, and finally, calibrating the radar reflectivity factor by using the channel data by means of a calibration formula.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of weather forecasting, and particularly relates to a channel-aware terahertz cloud radar calibration method. BACKGROUND

[0002] With the intensification of global climate change trends, accurate weather forecasting and climate change research have become a major strategic demand for national development. Clouds, especially high-altitude ice clouds that have a significant impact on the Earth's radiation balance, are one of the most critical but least understood components in current climate models. Therefore, accurate detection of clouds is a scientific problem that needs to be broken through.

[0003] Terahertz cloud radar is considered as a key development direction of next-generation atmospheric remote sensing technology due to its high sensitivity to cloud microphysical parameters (such as particle concentration, phase state, etc.) because of its wavelength comparable to cloud particle size. However, this high sensitivity also poses a serious challenge, which is that the terahertz wave will be strongly attenuated by meteorological elements such as water vapor and cloud particles when penetrating the atmosphere, and this attenuation changes dramatically with meteorological conditions. This complex path attenuation introduces a huge, time-varying measurement error to the radar echo signal, causing the inverted radar reflectivity factor to deviate significantly from the true value, thereby fundamentally restricting the detection accuracy of the radar.

[0004] In order to ensure the accuracy of the measurement data, the radar system must be accurately calibrated. However, the existing radar calibration techniques mainly have the following defects:

[0005] (1) The traditional segmented calibration method tests and calibrates each subsystem of the radar transmitter, receiver, etc. independently in the laboratory environment. The fundamental defect is that it fails to achieve absolute calibration of the whole link "end to end", and cannot evaluate the overall performance of the system under real observation conditions. The errors of each subsystem will accumulate in actual work, resulting in significant systematic deviation in the final measurement results.

[0006] (2) The existing end-to-end static calibration method is currently more widely used, mainly including the use of known radiation sources "sun method", or in a specific weather conditions (such as stable light rain) "light rain method", and the use of standard targets with known radar cross section (RCS) (such as corner reflector, metal ball) for calibration. The common fatal defect of these methods is the "static" characteristics and strict application conditions. For example, the "sun method" requires clear sky, "light rain method" depends on the stability of precipitation conditions, and the standard target method has the problems of difficult deployment, complex process, and easy to be blocked by the environment. Therefore, these methods cannot meet the needs of all-weather, continuous operation of the cloud radar, and cannot deal with the problems of time drift of system parameters due to temperature changes, or nonlinear errors introduced by dynamic changes of the environment. They have the limitations of "limited application scenarios, poor dynamic performance, and complex calibration process".

[0007] (3) Although there have been some studies on weather radar calibration at lower frequencies (such as Ka band, W band), there is no public report on the calibration of terahertz cloud radar in existing patents and technologies. This shows that there is no effective calibration solution for the strong attenuation and high dynamic characteristics of the terahertz frequency band. The main difficulty is that the targets measured by the radar (cloud, water vapor) are the main source of measurement error (path attenuation). That is, in order to accurately measure the cloud, the attenuation caused by the cloud must be accurately known; but to accurately know the attenuation, the detailed characteristics of the cloud need to be known in advance. Existing technologies try to break this cycle by introducing an external reference target (sun, corner reflector) that is not related to the cloud, but as mentioned earlier, these external references are often unreliable or unavailable in actual applications. Details are as follows:

[0008] Prior art one

[0009] The core idea of traditional end-to-end calibration technology is to use an external standard target with known and stable radar cross section (RCS) as a reference, and by comparing the measured return power of the radar with the theoretical return power, the overall gain and error of the radar system are inversely calculated and calibrated. The specific implementation scheme usually includes: 1) Fix the standard target such as corner reflector or metal ball on the tower, top of the mountain or building far away. The radar periodically measures the target to realize the monitoring and calibration of the system performance. 2) Use a tethered balloon to hang a metal ball, or use a UAV to carry a calibration ball, and move it to a specific position in the radar detection airspace for calibration. This way is more flexible than fixed targets in deployment.

[0010] The principle of the method provided by the prior art one is direct, but in practical application, especially in a business cloud measurement radar system requiring all-weather and continuous operation, there are fundamental defects: 1) The cost of building and maintaining a tower or platform for placing a standard target is high. Using a balloon or an unmanned aerial vehicle as a mobile platform has a complex operation process and is limited by airspace control, wind and other weather conditions, making it difficult to achieve regular and automated calibration operations. 2) The fixed target is easily blocked by terrain and buildings, and can only be calibrated for a fixed path between the radar and the target, and cannot reflect the real attenuation on other detection paths. The mobile target is severely affected by weather (such as strong winds and precipitation) and cannot work in bad weather conditions. 3) Such methods are essentially static or quasi-static. They can only provide one calibration at a specific time point and cannot track and compensate for system gain drift caused by environmental temperature changes or path attenuation changes caused by rapid changes in atmospheric conditions in real time and continuously. This is fatal to cloud measurement tasks that require high time resolution and high stability.

[0011] Prior art two

[0012] In order to overcome the problem of difficult deployment of external standard targets, some researchers have developed methods for calibration using stable or known radiation sources in nature. Mainly include the sun method, light rain method, dry snow method.

[0013] The sun method uses the sun as a stable and known microwave source of radiant flux. By accurately pointing the radar antenna at the sun, the received solar radiation power is measured and compared with the known theoretical value of the solar radiation flux, thereby calibrating the gain of the radar receiving chain.

[0014] The light rain method / dry snow method is mainly used for differential reflectivity (ZDR) calibration of dual-polarization radars. The physical basis is that under certain weather conditions (such as drizzle or dry snow), the shape of the precipitation particles is close to spherical, and their theoretical ZDR value should tend to 0 dB. Therefore, by measuring the degree of deviation of ZDR from 0 dB at this time, the ZDR deviation of the system can be calibrated.

[0015] Although the method provided by the prior art two eliminates the trouble of deploying external targets, its inherent limitations make it also unable to meet the calibration needs of terahertz cloud measurement radars. Specifically:

[0016] 1) The "sun method" requires clear and cloudless conditions, and any cloud cover will make it invalid. The "light rain method" and "dry snow method" completely rely on the occurrence of specific types of precipitation weather, and the calibration opportunity is accidental and unpredictable. Therefore, these methods cannot achieve all-weather and continuous calibration.

[0017] 2) The "sun method" can only calibrate the receiving link of the radar, and cannot cover the transmitting link, so it cannot calibrate the ZDR error introduced by the transmitting channel imbalance. The "light rain method" can calibrate the total ZDR of the system, but the calibration result is easily affected by the water accumulation on the top of the antenna cover or the structure of the flight warning light, introducing new errors.

[0018] 3) The accuracy of the "light rain method" is highly dependent on the assumption that "raindrops are ideal spheres", but the actual precipitation has a very complex raindrop spectrum distribution and morphology, resulting in a large range of true values of ZDR, thereby affecting the calibration accuracy.

[0019] 4) Like the standard target method, these methods are essentially static and occasional, and cannot solve the problem of severe and dynamic changes in the atmospheric path attenuation caused by water vapor and cloud particles in the terahertz band.

[0020] Overall, the existing technology cannot solve the technical problem of accurately, real-time, dynamic and automatic end-to-end calibration of the terahertz cloud measuring radar in a complex and dynamically changing weather environment. The existing methods cannot effectively compensate for the dynamic drift of the complex and variable atmospheric path attenuation and system parameters in the terahertz band due to their static and non-adaptive nature, which seriously restricts the practical application performance of the terahertz cloud measuring radar. SUMMARY

[0021] In order to overcome or alleviate one or more of the above technical problems, the purpose of the present application is to provide a channel-aware terahertz cloud measuring radar calibration method, which aims to solve the technical problems of unstable echo intensity, system gain error and large calibration error of the existing terahertz cloud measuring radar caused by atmospheric attenuation and system drift in a complex weather environment. The core goal is to realize online estimation and dynamic correction of radar system errors, thereby completing end-to-end dynamic self-calibration of the radar reflectivity factor, significantly improving the accuracy, stability and environmental adaptability of the radar measurement.

[0022] The present application provides the following technical solutions:

[0023] A channel-aware terahertz cloud measuring radar calibration method, comprising the following steps:

[0024] S1, two types of input data are obtained, one is the original, uncalibrated echo data collected by the terahertz cloud measuring radar, and the echo data is the radar reflectivity factor of each range bin; the other is the related environmental parameters, such as temperature, air pressure, relative humidity, etc. measured by the ground meteorological station;

[0025] S2, a multi-scale channel model is constructed, the environmental parameters are input into the multi-scale channel model, the multi-scale channel model is constructed to contain micro-particle scattering, meso-meteorological disturbance and macro-terrain turbulence attenuation, and the comprehensive path attenuation experienced by the radar signal at each distance bin on the propagation path can be calculated;

[0026] S3, a distance error estimation model is constructed, the path attenuation predicted by the multi-scale channel model and the original radar data are input into the error estimation model, and the error estimation model calculates the system gain deviation and the path attenuation factor at the current moment through a specific algorithm;

[0027] S4, a dynamic time correction model is constructed, the original error parameter sequence is estimated and input into the dynamic time correction model, the dynamic time correction model adopts Kalman filtering algorithm to perform time smoothing processing on the original error parameter, and outputs the correction value which is stable in time and can reflect the slow drift of the system;

[0028] S5, a calibration module is constructed, the smoothed correction value is applied to the echo data obtained in step S1, the final accurately calibrated radar reflectivity factor is calculated according to the calibration formula, and the final accurately calibrated radar reflectivity factor is output.

[0029] According to some embodiments, step S2 of constructing a multi-scale channel model includes the following steps:

[0030] S2-1, micro-scale channel model

[0031] The size distribution of cloud droplets is approximated by a Gamma distribution, and by setting different size distribution parameters, the cloud droplet size distribution in typical meteorological scenarios of thin clouds and thick clouds is simulated respectively, and the specific expression is as follows:

[0032]

[0033] Wherein, n(r) represents the number of liquid droplets in unit volume with a radius of r~r+dr, r mod is the mode radius of each type of particle, N is the particle number density, α and γ are constants describing the slope of the distribution, and a is a standardization constant.

[0034] For cloud particles, the Mie scattering theory is used to calculate the scattering characteristics, the multiple scattering effect in the scattering field is ignored during the calculation, and the cloud particles are approximated as spherical particles, and the cloud particle extinction efficiency is calculated as follows:

[0035]

[0036] Wherein, Q e is the cloud particle extinction efficiency, a n and b nis a function of the size parameter and the complex refractive index of the particle, x = 2πr / λ represents the size parameter, where λ represents the wavelength, for cloud particles with a certain size distribution, the calculation of the signal attenuation caused by the unit distance is as follows:

[0037]

[0038] where Q e is the extinction efficiency of cloud particles, m is the refractive index of the particle relative to the surrounding medium, r is the radius of the particle, and v is the frequency.

[0039] S2-2, constructing a mesoscale channel model

[0040] In order to better describe the change of temperature and humidity in the region and its influence on the calculation of attenuation, the present application adopts a two-dimensional Gaussian random field method to model the relative humidity in the research area, generates a humidity distribution field with spatial correlation and natural fluctuation characteristics by superimposing multiple Gaussian kernel functions, thereby depicting the random fluctuation of humidity in geographical space, and the specific formula is as follows:

[0041]

[0042] where RH0 is the basic humidity value, A k is the amplitude of the kth Gaussian disturbance, (x k , y k ) is the position of the kth disturbance center, L corr is the correlation length, which determines the spatial scale of humidity change, and N is the number of Gaussian kernels, and the temperature T(h) at height h is represented as:

[0043] T(h) = T0 + Γ·h (5)

[0044] In the formula, T0 is the temperature at the reference height, Γ is the temperature decrement rate, which is-0.0065K / m in the standard atmosphere, h is the height, and the unit is m;

[0045] S2-3, constructing a macro-scale channel model, the steps are as follows:

[0046] First, atmospheric absorption

[0047] For the absorption attenuation of terahertz signals in the atmosphere, the calculation is based on the ITU-RP.676 recommended model, the total attenuation γ(f) of the signal in the atmosphere comes from including water vapor and oxygen, and the calculation formula is:

[0048] γ(f) = γ0(f) + γ w (f) = 0.1820fN''(f) (dB / km) (6)

[0049] where γ0(f) represents the characteristic attenuation of dry air to the signal, γ w (f) represents the characteristic attenuation under certain water vapor density, f represents the frequency, unit is GHz; N"(f) is the hypothetical part of the complex refractive index related to the frequency, as follows:

[0050]

[0051] where S i is the intensity of the i-th line, F i is the curve shape factor, and the summation is extended to all lines, N D "(f) is the dry continuous band of nitrogen absorption and Debye spectrum caused by atmospheric pressure;

[0052] Substitute the temperature field T(x, y) and humidity field RH(x, y) generated by the mesoscale channel model into the macroscopic absorption attenuation formula to obtain the attenuation coefficient of the spatial position (x, y, h):

[0053] γ(x, y, h) = γ0(f, T(h)) + γ w (f, RH(h)) (8)

[0054] In the formula, γ(x, y, h) represents the signal attenuation coefficient of the position, unit is dB / km, the temperature and humidity dependence is provided by the mesoscale model, and the macroscopic atmospheric absorption formula is used for correction calculation;

[0055] Second, atmospheric turbulence

[0056] The refractive index structure parameter is used to describe the change of atmospheric turbulence attenuation with altitude:

[0057]

[0058] where v represents the wind speed, unit is m / s, h represents the altitude, unit is m, represents the C n 2 Normal value on the horizontal ground;

[0059] The Rytov approximation method considers that the atmospheric turbulence attenuation can be represented as:

[0060]

[0061] where L represents the propagation distance, unit is m, k represents the signal wave number;

[0062] And the Andrew method represents the atmospheric turbulence attenuation as:

[0063]

[0064] where D represents the diameter of the receiver, σ I 2 represents the scintillation index;

[0065]

[0066] where C n 2 represents the refractive index structure parameter, k represents the signal wave number, and L represents the propagation distance, in meters;

[0067] Third, multipath effect

[0068] For the multipath effect, the Jakes model is used for simulation, and the complex envelope corresponding to the Rayleigh fading channel is described as follows:

[0069]

[0070] where w d represents the maximum Doppler shift, φ n is a random initial phase, θ n represents the incident angle of the nth path, and M represents the number of paths.

[0071] According to some embodiments, step S3 constructs a distance error estimation model, including the following steps:

[0072] S3-1, gain deviation The estimation formula is as follows:

[0073]

[0074] where α i represents the atmospheric channel attenuation of the ith distance library at the current temperature, humidity, pressure and signal frequency; Z meas,i is the reflectivity factor measured by the radar corresponding to the distance library, and N represents the cloud droplet particle concentration obtained by backstepping according to the reflectivity data, as shown in the following formula:

[0075] N = 10 5 · Z miscalib + 2 × 10 6 (15)

[0076] where Z miscalib represents the measured value of the radar reflectivity factor;

[0077] Attenuation factor The estimation formula is as follows:

[0078]

[0079] where k radark represents the slope of the uncalibrated radar reflectivity data versus normalized range common k represents the slope of the channel model attenuation versus normalized range

[0080] S3-2, calibration of the range-dependent radar reflectivity factor, the formula for calibration is:

[0081]

[0082] where Z cal is the calibrated value of the radar reflectivity factor, d is the detection distance, D max is the upper limit of the detection distance, is the gain deviation, is the attenuation factor.

[0083] According to some embodiments, step S4 constructs a dynamic time correction model, including the following steps:

[0084] The noisy original estimation sequence is regarded as a quasi-periodic signal, a first-order random walk model is established in the state space through the Kalman prediction equation, and the initial value is set to be the gain deviation 1.5 dB and the attenuation factor 0.75; then the gain matrix is driven by real-time observation residual to perform least mean square error correction; when the observation noise is large, the high-frequency correction amount is automatically compressed, and when the noise is small, the real drift is quickly tracked. Finally, the uncertainty is estimated through covariance update, the new gain deviation and attenuation factor are output every time a new radar sample comes, and the online dynamic correction of the radar reflectivity factor is completed; the core formula includes:

[0085] Prediction state equation:

[0086]

[0087] Prediction covariance equation:

[0088]

[0089] Kalman gain equation:

[0090]

[0091] Update state equation:

[0092]

[0093] Update covariance:

[0094] P k = (I-K k H)P k (22)

[0095] where, is the prior state estimate, is the prior state estimate. is the prior estimate covariance matrix, P k-1 is the posterior estimate covariance matrix. Q is the process noise covariance matrix. K k is the Kalman gain matrix, floating between 0 and 1. H is the observation matrix, and R is the observation noise covariance matrix, z k is the actual observation vector. I is the identity matrix.

[0096] According to some embodiments, step S5 constructs a calibration module for calibrating the radar reflectivity factor after obtaining the estimated system gain deviation and path attenuation factor, and the formula is:

[0097]

[0098] wherein, Z cal (t) is the calibrated value of the radar reflectivity factor at a specific distance-time direction, Z miscalib (t) is the uncalibrated value of the radar reflectivity factor at a specific distance-time direction, d is the calibration distance, t is the detection time, and D max is the upper limit of the detection distance.

[0099] Compared with the prior art, the present application has the following beneficial effects:

[0100] The method for calibrating a terahertz cloud measurement radar based on channel perception provided by the present application quantifies and compensates complex atmospheric attenuation effects more accurately than existing technical solutions by constructing a physically driven multi-scale channel model, thereby greatly improving the measurement accuracy of the radar reflectivity factor.

[0101] In addition, the present application can be run online and continuously, and can estimate and correct errors introduced by environmental changes and system drift in real time. The Kalman filtering mechanism ensures that the calibration parameters can follow the dynamic changes of the system state, which is not achieved by all static calibration methods.

[0102] Furthermore, the present application completely eliminates the dependence on external calibration targets such as the sun and corner reflectors, and realizes true end-to-end dynamic self-calibration. The entire calibration process is fully automated and does not require human intervention, greatly simplifying the operation process and reducing the difficulty of operation and maintenance.

[0103] Finally, since there is no need to deploy and maintain expensive external calibration equipment, the present application not only achieves higher performance, but also significantly reduces the construction and operation cost of high-precision terahertz cloud measurement radar systems, laying a foundation for the widespread application and networking observation of this technology. BRIEF DESCRIPTION OF DRAWINGS

[0104] Figure 1 is the overall technical solution flowchart provided by the embodiments of the present application.

[0105] Figure 2 The application provides a fading cause and a research method.

[0106] Figure 3 The application provides a change of fading of different scales with a propagation distance under multi-scale channel modeling.

[0107] Figure 4 The application provides a multi-scale channel modeling system.

[0108] Figure 5 The application provides a time direction estimation of gain deviation and fading factors.

[0109] Figure 6 The application provides a radar reflectivity factor calibration effect comparison under a thin cloud scene.

[0110] Figure 7 The application provides a radar reflectivity factor calibration effect comparison under a thick cloud / rain scene.

[0111] Figure 8 The application provides a radar reflectivity factor calibration comparison with time frames at 800 m, 2000 m and 3200 m.

[0112] Figure 9 The application provides a measured radar reflectivity factor calibration comparison. DETAILED DESCRIPTION

[0113] The application constructs a micro-meso-macro three-dimensional integrated terahertz channel modeling framework in terms of channel modeling, and comprehensively integrates multiple sources such as particle scattering, meteorological disturbance and terrain reflection. For the micro scale, the distribution of cloud droplets is modeled by using the Gamma distribution theory, and the Mie scattering theory is used to calculate the attenuation of cloud droplets. For the meso scale, a humidity distribution field with spatial correlation and natural fluctuation characteristics is generated by superimposing multiple Gaussian kernel functions, so as to depict the random fluctuation of humidity in geographical space. Meanwhile, the vertical gradient and spatial variation of temperature are also considered. The output temperature and humidity distribution will be used as input parameters for the calculation of atmospheric absorption attenuation at the macro scale, to correct the atmospheric absorption model. For the macro scale, the application calculates the absorption attenuation of the terahertz signal in the atmosphere based on the ITU-R P.676 recommendation model. The macro scale also introduces the attenuation of atmospheric turbulence and multipath effect. The application uses the Andrew method to quantify the atmospheric turbulence attenuation, and uses the Jakes model to simulate the multipath effect. Finally, the application realizes the quantitative calculation of micro, meso and macro three attenuations under certain conditions. The user can directly input the atmospheric environmental factors to intuitively obtain the total attenuation. The calculation of the total attenuation can effectively help to solve the problem that the existing calibration scheme is greatly affected by the environment.

[0114] The core technical principle of the application is that instead of relying on external reference targets, a high-precision physical channel model is constructed to independently predict the path attenuation of the signal, and a physical mapping relationship between the atmospheric channel attenuation characteristics and the radar measurement error is established based on this, thereby breaking the coupling between the measurement target and the measurement error.

[0115] The application will be described in detail below in conjunction with the embodiments and drawings, but it should be understood that the embodiments and drawings are only used to exemplarily describe the application, and cannot constitute any limitation on the protection scope of the application. All reasonable transformations and combinations within the scope of the inventive concept of the application fall within the protection scope of the application.

[0116] The application will be further described below in conjunction with the drawings.

[0117] Embodiment 1

[0118] As Figure 1 , the embodiment discloses a channel-aware terahertz cloud measurement radar calibration method, which mainly includes the following steps:

[0119] S1, two types of input data are obtained, one is the original and uncalibrated echo data collected by the terahertz cloud measurement radar, and the echo data is the radar reflectivity factor of each distance bin; the other is the related environmental parameters, such as the temperature, air pressure and relative humidity measured by the ground meteorological station;

[0120] S2, a multi-scale channel model is constructed, environmental parameters are input into the multi-scale channel model, the multi-scale channel model includes atmospheric molecular attenuation, water cloud particle scattering, multipath effect and atmospheric turbulence attenuation, and the comprehensive path attenuation experienced by the radar signal at each distance bin on the propagation path is calculated;

[0121] S3, a distance error estimation model is constructed, the path attenuation predicted by the multi-scale channel model and the original radar data are input into the error estimation model, and the model calculates the system gain deviation and the path attenuation factor at the current moment through a specific algorithm;

[0122] S4, a dynamic time correction model is constructed, the original error parameter sequence is estimated and input into the dynamic time correction model, the dynamic time correction model uses a Kalman filter algorithm to perform time smoothing processing on the original error parameters, and outputs a correction value that is stable in time and can reflect the slow drift of the system;

[0123] S5, a calibration module is constructed, the smoothed correction value is applied to the echo data obtained in step S1, the final accurately calibrated radar reflectivity factor is calculated according to a calibration formula, and is output.

[0124] Specifically:

[0125] According to Figure 2 , in step S2, a multi-scale channel model is constructed:

[0126] To more accurately describe the attenuation of the terahertz radar signal in the atmospheric channel during propagation, the embodiment constructs a channel model that covers particle scattering, meteorological disturbance, atmospheric absorption, atmospheric turbulence and multipath effect, and includes different scales from micro to macro:

[0127] S2-1, micro-scale channel model

[0128] The micro-scale mainly studies the change of channel parameters in a very small space. In this scale range, the formation and distribution of cloud droplets, water droplets and other scatterers are relatively uniform. At this time, the micro-scale channel model needs to be able to describe the uniform and stable scattering characteristics of cloud droplets, water droplets and other particles. Particle distribution function (PDF) can be used to describe the size distribution of these particles, for example, Gamma distribution can be used to approximate the size distribution of cloud droplets, and by setting different size distribution parameters, the cloud droplet size distribution in two typical meteorological scenarios, thin cloud and thick cloud / rain, can be simulated. The specific expression is as follows:

[0129]

[0130] Where n(r) represents the number of liquid droplets in unit volume with a radius of r~r+dr, r modis the mode radius of each type of particle, N is the particle number density, a and g are constants describing the slope of the distribution, and a is a normalization constant.

[0131] For cloud particles, the Mie scattering theory is used to calculate the scattering characteristics, and the multiple scattering effect in the scattering field is ignored in the calculation process. The cloud particles are approximated as spherical particles, and the cloud particle extinction efficiency is calculated as follows:

[0132]

[0133] where Q e n and b n are functions of the size parameter and the complex refractive index of the particle, and x = 2πr / λ represents the size parameter, where λ represents the wavelength. For cloud particles with a certain size distribution, the calculation of the signal attenuation caused by a unit distance is as follows:

[0134]

[0135] where Q e is the extinction efficiency of the cloud particle, m is the refractive index of the particle relative to the surrounding medium, r is the radius of the particle, and v is the frequency.

[0136] S2-2, Constructing mesoscale channel model

[0137] At the mesoscale, the distribution of channel parameters is influenced by more extensive but still local weather. By simulating the influence of local meteorological disturbances on terahertz signal propagation, a channel model is established that can reflect changes in meteorological factors such as humidity and temperature. Specifically, a two-dimensional Gaussian random field method is used to model the relative humidity in the study area. By superimposing multiple Gaussian kernel functions, a humidity distribution field with spatial correlation and natural fluctuation characteristics is generated, thereby depicting the random fluctuations of humidity in geographical space. This method can effectively reflect the non-uniformity of humidity distribution caused by mesoscale meteorological processes such as weather systems and local convection. The specific formula is as follows:

[0138]

[0139] where RH0 is the basic humidity value, A k is the amplitude of the kth Gaussian disturbance, (x k , y k ) is the position of the kth disturbance center, L corr is the correlation length, which determines the spatial scale of humidity change, and N is the number of Gaussian kernels. In addition, the model also considers the vertical gradient and spatial variation of temperature, further improving the description of atmospheric parameters. The temperature at height h is represented as

[0140] ​T(h) = To + Γ · h (5)

[0141] where To is the temperature at the reference height (usually ground), Γ is the temperature lapse rate, which is usually taken as -0.0065 K / m in standard atmosphere, and h is the height (unit: m). The temperature and humidity distribution output by the mesoscale model will be used as input parameters for the macro-scale atmospheric absorption attenuation calculation to correct the atmospheric absorption model, so as to reflect the influence of local disturbance on signal path attenuation.

[0142] S2-3, constructing a macro-scale channel model

[0143] (1) Atmospheric absorption

[0144] For the absorption attenuation of terahertz signals in the atmosphere, the embodiment is based on the ITU-R P.676 recommendation model for calculation. The total attenuation γ(f) of the signal in the atmosphere mainly comes from water vapor and oxygen, and its calculation formula is:

[0145] γ(f) = γ0(f) + γ w (f) = 0.1820 f N"(f) (dB / km) (6)

[0146] where γ0(f) (unit: dB / km) represents the characteristic attenuation of dry air (i.e. oxygen) to the signal, γ w (f) (unit: dB / km) represents the characteristic attenuation under a certain water vapor density, f represents the frequency, and the unit is GHz; N"(f) is the hypothetical part of the complex refractive index related to the frequency:

[0147]

[0148] where S i is the intensity of the i-th line, F i is the curve shape factor, and the sum is extended to all lines, N D "(f) is the dry continuous band of nitrogen absorption and Debye spectrum caused by atmospheric pressure.

[0149] In order to reflect the influence of mesoscale local disturbance on path attenuation, the embodiment substitutes the temperature field T(x, y) and humidity field RH(x, y) distribution generated by the mesoscale model into the macro-scale absorption attenuation formula to obtain the attenuation coefficient of the spatial position (x, y, h):

[0150] γ(x, y, h) = γ0(f, T(h)) + γ w (f, RH(h)) (8)

[0151] where γ(x, y, h) represents the signal attenuation coefficient (dB / km) at the location, the temperature and humidity dependence is provided by the mesoscale model, and the macroscopic atmospheric absorption formula is used for correction calculation.

[0152] (2) Atmospheric turbulence

[0153] The attenuation caused by atmospheric turbulence is mainly due to the heating effect of solar radiation on the ground surface under clear weather, which leads to the increase of air temperature near the ground and the decrease of air density, resulting in the temperature and density difference between the ground and the upper air. This spatial non-uniformity promotes the upward movement and disturbance of air, forming turbulent motion, thereby causing the turbulence attenuation of electromagnetic waves in the propagation process. In order to characterize the atmospheric turbulence, the refractive index structure parameter is used to describe the change of atmospheric turbulence attenuation with altitude:

[0154]

[0155] where v represents the wind speed in m / s, h represents the altitude in m, represents the C n 2 The normal value on the horizontal ground surface.

[0156] Due to the nonlinear behavior of atmospheric quantities such as temperature, pressure and wind speed, it is complex to evaluate turbulence. The Rytov approximation method considers that the atmospheric turbulence attenuation can be represented as:

[0157]

[0158] where L represents the propagation distance in m, and k represents the signal wave number.

[0159] However, this method does not consider the receiver aperture effect, which is very important for turbulence research.

[0160] The Andrew method represents the atmospheric turbulence attenuation as:

[0161]

[0162] where D represents the receiver diameter, σ I 2 represents the scintillation index.

[0163]

[0164] where C n 2 represents the refractive index structure parameter, k represents the signal wave number, and L represents the propagation distance in m;

[0165] (3) Multipath effect

[0166] The THz cloud radar signal inevitably scatters and reflects with rough surfaces such as mountains and cities, and then produces multipath effect. For the multipath effect, the Jakes model is used for simulation, and the complex envelope corresponding to the Rayleigh fading channel described by the Jakes model is:

[0167]

[0168] where w d represents the maximum Doppler shift, φ n is a random initial phase, θ n represents the incidence angle of the nth path, and M represents the number of paths.

[0169] As Figure 3 is the multi-scale channel model constructed in this embodiment, which includes the variation relationship of micro, meso, macro and total attenuation with the increase of propagation distance.

[0170] Through the comprehensive of the above-mentioned micro-scale attenuation, meso-scale attenuation and macro-scale attenuation, the signal attenuation in the thin cloud and thick cloud / rain two scenarios is obtained, as shown in Figure 4 .

[0171] In step S3, an error estimation model is constructed for predicting the gain deviation and the attenuation factor.

[0172] In order to more accurately calibrate the reflectivity factor of the THz cloud radar, the embodiment proposes an estimation method of gain deviation and attenuation factor based on the multi-scale channel model. The method is based on the estimation of the error caused by the multi-scale channel model attenuation data and the radar thermal noise, and realizes the calibration of the radar reflectivity factor. The method includes the following steps:

[0173] S3-1, estimation of gain deviation The estimation formula is as follows:

[0174]

[0175] where α i represents the atmospheric channel attenuation of the ith distance library under the current temperature, humidity, pressure and signal frequency, which is proposed by the multi-scale channel modeling system. Z meas,i is the radar measured reflectivity factor corresponding to the distance library, N represents the cloud droplet particle concentration obtained by inversely calculating the reflectivity data, and is usually calculated by:

[0176] N=10 5 ·Z miscalib +2×10 6 (15)

[0177] where Z miscalib represents the measured value of the radar reflectivity factor.

[0178] Attenuation factor The estimation formula is:

[0179]

[0180] Wherein, k radar represents the slope of the uncalibrated radar reflectivity data and the normalized distance, k common represents the slope of the channel model attenuation and the normalized distance.

[0181] S3-2, calibration of radar reflectivity factor, the embodiment proposes a new formula for calibrating radar reflectivity factor using channel data:

[0182]

[0183] Wherein, Z cal is the calibration value of the radar reflectivity factor, d is the detection distance, D max is the upper limit of the detection distance, and in the embodiment, the upper limit is 4000 meters, is the gain deviation, is the attenuation factor.

[0184] Based on the above formula, the embodiment can realize real-time and efficient dynamic calibration of the reflectivity factor of the terahertz cloud radar within a predetermined detection distance range, thereby significantly improving the accuracy and stability of the radar measurement.

[0185] In step S4, a dynamic time correction model is constructed for dynamic calibration of the system gain deviation and the path attenuation factor at different time frames.

[0186] The embodiment proposes a dynamic calibration method based on channel observation for the problem that the radar system gain deviation and the path attenuation factor change dynamically with time, resulting in inaccurate reflectivity factor. The core technical principle is: the noisy original estimation sequence is regarded as a quasi-periodic signal, a first-order random walk model is established in the state space through the Kalman prediction equation, and the initial value is set to be the gain deviation 1.5 dB and the attenuation factor 0.75; then the real-time observation residual is used to drive the gain matrix for least mean square error correction; when the observation noise is large, the high-frequency correction amount is automatically compressed, and when the noise is small, the real drift is quickly tracked. Finally, the uncertainty is estimated through covariance update, the new gain deviation and attenuation factor are output every time a new radar sample comes, and the online dynamic correction of the radar reflectivity factor is completed; the core formula includes:

[0187] Predicted state:

[0188]

[0189] Predicted covariance:

[0190]

[0191] Kalman gain:

[0192]

[0193] Updated state:

[0194]

[0195] Updated covariance:

[0196]

[0197] where, is the prior state estimate, is the previous posterior state estimate. is the prior estimate covariance matrix, P k-1 is the posterior estimate covariance matrix. Q is the process noise covariance matrix. K k is the Kalman gain matrix, floating between 0 and 1. H is the observation matrix, R is the observation noise covariance matrix, z k is the actual observation vector. I is the identity matrix.

[0198] Specifically, as shown in Figure 5 , the dynamic calibration of the system gain deviation and the path attenuation factor in different time frames is shown, and compared with the true value.

[0199] In step S5, a calibration module is constructed, which is used to calibrate the radar reflectivity factor after the estimated system gain deviation and path attenuation factor are obtained. The formula for calibrating the radar reflectivity factor is:

[0200]

[0201] where, Z cal (t) is the calibrated value of the radar reflectivity factor at a specific distance and time, Z miscalib (t) is the uncalibrated value of the radar reflectivity factor at a specific distance and time, d is the calibration distance, t is the detection time, D max is the upper limit of the detection distance.

[0202] Embodiment 2:

[0203] In this embodiment, the calibration method provided in Embodiment 1 is used to show the calibration effect of the radar reflectivity factor under different distance libraries (thin cloud and thick cloud scenarios).

[0204] In this case, two typical meteorological conditions, thin cloud and thick cloud, are selected. First, the calibration effect of the radar reflectivity factor under different distance libraries is shown as follows: Figure 4The multi-scale channel model is shown (temperature 15℃, humidity 60% and 90% respectively). Then the real value, uncalibrated data and calibrated data of the radar echo under the same distance library are respectively compared and shown. As shown in Figure 6 and Figure 7 As shown in the uncalibrated case, the radar reflectivity factor presents systematic deviation with the increase of distance, and the error is significant; after using the dynamic calibration method based on channel perception, the calibrated data is highly consistent with the real data in the full distance library range, effectively eliminating the systematic error accumulated with distance, and significantly improving the measurement accuracy and consistency. The calibration effect of this embodiment is shown in Table 1. It can be seen that the average error under thin cloud condition is reduced from 2.89dB before calibration to 0.41dB after calibration, and the average error is reduced by 85.96%. The root mean square error (RMSE) is reduced from 2.96dB before calibration to 0.52dB after calibration, and the RMSE is reduced by 82.51%. The average error under thick cloud / rain condition is reduced from 2.48dB before calibration to 0.58dB after calibration, and the average error is reduced by 76.71%. The RMSE is reduced from 2.71dB before calibration to 0.59dB after calibration, and the RMSE is reduced by 78.27%.

[0205] Table 1 Performance comparison before and after calibration (thin cloud, thick cloud / rain)

[0206]

[0207] Example 3

[0208] In this embodiment, the radar reflectivity factor calibration effect under different time sequences is shown.

[0209] As shown in Figure 8As shown, the case selects three representative distance points - 800 meters (short distance point), 2000 meters (medium distance point) and 3200 meters (long distance point), and analyzes the comparison of real data, uncalibrated data and calibrated data after using the method of the application on the time frame (continuous 150 sampling time points). The results show that the uncalibrated radar reflectivity factor has obvious systematic deviation and drift, and after using the online calibration method based on channel perception, the calibrated data can closely track the real value, significantly reducing the mean square error and average error under the time sequence (see Table 2). It can be seen that the average error of 800 meters distance is reduced from 1.93 dB to 0.07 dB (96.16% reduction), and the RMSE is reduced from 1.89 dB to 0.06 dB (96.83% reduction); the average error of 2000 meters distance is reduced from 2.16 dB to 0.14 dB (93.38% reduction), and the RMSE is reduced from 2.12 dB to 0.12 dB (94.39% reduction); the average error of 3200 meters distance is reduced from 1.95 dB to 0.12 dB (93.80% reduction), and the RMSE is reduced from 1.91 dB to 0.09 dB (95.09% reduction). Reflecting the improvement of the calibration method provided by embodiment 1 on the long-term stability of radar measurement. This case further verifies the real-time adaptability and robustness of the scheme provided by embodiment 1 in dynamic weather environment.

[0210] Table 2 Comparison of radar reflectivity factor time series calibration performance

[0211]

[0212] Embodiment 4

[0213] As Figure 9 The radar reflectivity factor before and after calibration is obtained from the measured data. The measured scene temperature is 23℃, and the humidity is 95%. In the process of processing the measured data, first, the uncalibrated radar reflectivity factor is obtained by means of the radar formula and the input parameters of the radar system itself and the transmitted waveform, such as: transmitted peak power, horizontal and vertical beam width, transmitted pulse width, antenna gain, etc. Then, the estimated values of the attenuation factor and the gain deviation are obtained by means of the weather conditions at that time and the radar detection distance data. The final calibrated radar reflectivity factor is obtained by combining the formula.

[0214] The above embodiments are only preferred embodiments of the present application, and the protection scope of the present application is not limited to the above embodiments. Any technical solutions falling within the concept of the present application shall fall within the protection scope of the present application. It should be pointed out that, for ordinary skilled persons in the art, improvements and refinements without departing from the principles of the present application shall also be considered as the protection scope of the present application.

Claims

1. A terahertz cloud-measuring radar calibration method based on channel awareness, characterized in that: Includes the following steps: S1. Acquire two types of input data: one is the raw, uncalibrated echo data collected by the terahertz cloud radar, which is the radar reflectivity factor of each range library; the other is the relevant environmental parameters, such as temperature, air pressure, and relative humidity measured by ground weather stations. S2. Construct a multi-scale channel model, input the environmental parameters into the multi-scale channel model, the multi-scale channel model includes atmospheric molecule attenuation, water cloud particle scattering, multipath effect and atmospheric turbulence attenuation, and calculate the comprehensive path attenuation experienced by the radar signal at each range along the propagation path. S3. Construct a range error estimation model. Input the path attenuation predicted by the multi-scale channel model and the original radar data into the error estimation model. The error estimation model calculates the system gain deviation and path attenuation factor at the current moment using a specific algorithm. S4. Construct a dynamic time correction model, estimate the original error parameter sequence and feed it into the dynamic time correction model. The dynamic time correction model uses the Kalman filter algorithm to perform time smoothing on the original error parameters and outputs a time-stable correction value that can reflect the slow drift of the system. S5. Construct a calibration module, apply the smoothed correction value to the echo data obtained in step S1, calculate the final precisely calibrated radar reflectivity factor according to the calibration formula, and output it.

2. The terahertz cloud measurement radar calibration method based on channel awareness according to claim 1, characterized in that: Step S2 involves constructing a multi-scale channel model, including the following steps: S2-1, Microscale Channel Model The Gamma distribution is used to approximate the size distribution of cloud droplets. By setting different size distribution parameters, the droplet size distribution under typical meteorological scenarios of thin and dense clouds is simulated. The specific expressions are as follows: Where n(r) represents the number of droplets with radii between r and r+dr per unit volume, r mod α is the model radius of various particles, N is the particle number density, α and γ are constants describing the distribution slope, and a is the normalization constant. For cloud particles, the Mie scattering theory is used to calculate the scattering characteristics. Multiple scattering effects in the scattering field are ignored during the calculation, and cloud particles are approximated as spherical particles. The extinction efficiency of cloud particles is calculated as follows: Among them, Q e It is the extinction efficiency of cloud particles, a n and b n It is a function of the size parameter and the complex refractive index of the particle, where x = 2πr / λ represents the size parameter, and λ represents the wavelength. For cloud particles with a certain size distribution, the signal attenuation is calculated as follows: Among them, Q e is the extinction efficiency of cloud particles, m is the refractive index of the particle relative to the surrounding medium, r is the radius of the particle, and v is the frequency. S2-2, Constructing a mesoscale channel model A two-dimensional Gaussian random field method is used to model the relative humidity in the study area. By superimposing multiple Gaussian kernel functions, a humidity distribution field with spatial correlation and natural fluctuation characteristics is generated, thereby characterizing the random fluctuation of humidity in geographic space. The specific formula is as follows: Among them, RH0 is the basic humidity value, A k Let x be the amplitude of the k-th Gaussian perturbation. k ,y k Let L be the location of the k-th disturbance center. corr The relevant length determines the spatial scale of humidity change, N is the number of Gaussian kernels, and the temperature T(h) at height h is expressed as: T(h)=T0+Γ·h (5) In the formula, T0 is the temperature at the reference altitude, Γ is the temperature lapse rate, which is taken as -0.0065 K / m in standard atmosphere, and h is the altitude in meters. S2-3. Construct a macroscopic-scale channel model, the steps of which are as follows: First, atmospheric absorption The absorption and attenuation of terahertz signals in the atmosphere is calculated based on the ITU-RP.676 recommended model. The total attenuation γ(f) of the signal in the atmosphere comes from water vapor and oxygen, and its calculation formula is as follows: γ(f)=γ0(f)+γ w (f)=0.1820fN″(f) (dB / km) (6) Where γ0(f) represents the characteristic attenuation of the signal by dry air, γ w (f) represents the characteristic attenuation under a certain water vapor density condition, where f represents the frequency in GHz; N""f" is the assumed part of the composite refractive index related to this frequency, as shown in the following formula: Among them, S i F is the intensity of the i-th line. i It is the curve shape factor and the sum extended to all lines, N D (f) is the dry continuous band of nitrogen absorption and Debye spectrum caused by atmospheric pressure; Substituting the temperature field T(x,y) and humidity field RH(x,y) generated by the mesoscale channel model into the macroscale absorption attenuation formula, we obtain the attenuation coefficient at the spatial location (x,y,h): γ(x,y,h)=γ0(f,T(h))+γ w (f,RH(h)) (8) In the formula, γ(x,y,h) represents the signal attenuation coefficient at that location, in dB / km. The temperature and humidity dependence is provided by the mesoscopic model, and the macroscopic atmospheric absorption formula is used for correction calculation. Second, atmospheric turbulence Using refractive index structure parameters to describe the variation of atmospheric turbulence attenuation with altitude: Where v represents wind speed in m / s and h represents altitude in m. This represents C under the condition of a wind speed of 21 m / s. n 2 Normal values ​​on a level surface; The Retoff approximation method states that atmospheric turbulence attenuation is expressed as: Where L represents the propagation distance in meters (m) and k represents the signal wavenumber; Atmospheric turbulence attenuation can be expressed using the Andrew method as follows: Where D represents the receiver diameter, σ I 2 Indicates the flicker index; Among them, C n 2 The refractive index structure parameter represents the signal wavenumber, and L represents the propagation distance, with units of meters (m). Third, multipath effect For the multipath effect, the Jakes model is used for simulation, and the corresponding complex envelope describing the Rayleigh fading channel is: Among them, w d Represents the maximum Doppler frequency shift, φ n It is a random initial phase, θ n Let represent the angle of incidence for the nth path, and M represent the number of paths.

3. The terahertz cloud measurement radar calibration method based on channel awareness according to claim 2, characterized in that: Step S3 involves constructing an error estimation model, including the following steps: S3-1, Gain Deviation The estimation formula is as follows: Where, α i Z represents the atmospheric channel attenuation of the i-th distance library under the current temperature, humidity, pressure, and signal frequency; meas,i The reflectivity factor is obtained from radar measurements corresponding to the range library, where N represents the cloud droplet particle concentration derived from the reflectivity data, as shown in the following formula: N=10 5 ·Z miscalib +2×10 6 (15) Among them, Z miscalib The measured value representing the radar reflectivity factor; Attenuation factor The estimation formula is: Where, k radar k represents the slope of the uncalibrated radar reflectivity data relative to the normalized distance. common The slope representing the channel model attenuation and normalized distance; S3-2. Calibration of radar reflectivity factor. The calibration formula is as follows: Among them, Z cal Here, d represents the calibration value of the radar reflectivity factor, and d represents the detection range. max This is the upper limit of the detection distance. For gain bias, This is the attenuation factor.

4. The terahertz cloud measurement radar calibration method based on channel awareness according to claim 3, characterized in that: Step S4 involves constructing a dynamic time correction model, including the following steps: The noisy original estimated sequence is treated as a quasi-periodic signal. A first-order random walk model is established in the state space using the Kalman prediction equation, with initial values ​​set as a gain bias of 1.5 dB and an attenuation factor of 0.

75. The gain matrix is ​​then driven by real-time observation residuals for minimum mean square error correction. When observation noise is high, the high-frequency correction is automatically compressed; when noise is low, the actual drift is quickly tracked. Finally, the uncertainty is updated using covariance, enabling the output of the latest gain bias and attenuation factor for each new radar sample, thus achieving online dynamic correction of the radar reflectivity factor. The core formulas include: Predicted status: Predicting covariance: Kalman gain: Update status: Update covariance: in, For prior state estimation, This is the posterior state estimate for the previous step. To estimate the covariance matrix a priori, P k-1 Let K be the posterior estimated covariance matrix. Q is the process noise covariance matrix. k Let be the Kalman gain matrix, floating between 0 and 1. H is the observation matrix, R is the observation noise covariance matrix, and z is the Kalman gain matrix. k is the actual observed vector. I is the identity matrix.

5. The terahertz cloud measurement radar calibration method based on channel awareness according to claim 4, characterized in that: Step S5 constructs a calibration module, which is used to calibrate the radar reflectivity factor after obtaining the estimated system gain deviation and path attenuation factor. The formula is as follows: Among them, Z cal (t) represents the calibrated value of the radar reflectivity factor over time at a specific distance, where d is the detection distance, and D max This represents the upper limit of the detection range.

Images