A method for measuring whole-layer atmospheric turbulence based on centroid drift and light intensity scintillation

Abstract

The present invention discloses a method for measuring whole-layer atmospheric turbulence based on centroid drift and light intensity flicker. It belongs to the technical field of satellite-to-ground laser communication. The atmospheric turbulence parameter measurement step includes: using a high-resolution imaging lens and a high-speed visible light camera as a receiving device to collect the starlight spot image after atmospheric turbulence; then, the edge detection of the star map is performed by Zernike moment to accurately analyze and identify the shape features in the star map, and realize sub-pixel level edge detection; on this basis, dynamic Lanczos interpolation is performed on the edge area of ​​the star point to realize resampling and retain the details of the star map, and the normalized amplitude enhancement weight method of the gray gradient is used to improve the accuracy of centroid positioning, improve the anti-interference ability and achieve higher positioning accuracy; finally, by calculating the light intensity flicker variance and the centroid drift variance, the whole-layer atmospheric turbulence intensity in the atmospheric channel transmission condition during satellite-to-ground laser communication is obtained.

Description

Technical Field

[0001] The invention relates to the field of satellite-to-ground laser communication technology, and in particular to a whole-layer atmospheric turbulence measurement method based on centroid drift and light intensity flicker. Background Art

[0002] Satellite-to-ground laser communication is a key link in the integration of space, air, and land. The free space channel in which it is located is unstable and random, which affects the phase and amplitude of the laser beam, causing problems such as beam wavefront distortion, light intensity flicker, center of mass drift, and arrival angle fluctuations, affecting the stability of the laser communication system. Atmospheric turbulence parameters are one of the important parameter bases for evaluating and monitoring the transmission performance of satellite-to-ground laser communication links. It is necessary to monitor the atmospheric turbulence channel conditions in real time to provide a reference for the scheduling of the communication system.

[0003] At present, most atmospheric turbulence parameter measurement methods use customized optical equipment, such as laser scintillometers or meteorological sensors, and rely on meteorological data including temperature, relative humidity, air pressure, wind speed, etc. for estimation. For example, temperature pulsators and large-aperture scintillators are used. Among them, scintillators use a two-end optical device consisting of a transmitter and a receiver at both ends of the turbulent medium to observe scintillation and realize the measurement of atmospheric turbulence. For the inclined paths and moving target platforms commonly found in satellite-to-ground laser communication links, it is not convenient to use the two-end measurement method. This requires further exploration of new atmospheric turbulence parameter measurement methods. In the field of astronomical imaging, the plane waves emitted by stars propagate in atmospheric turbulence, which will produce wavefront distortion, thereby changing the propagation direction and energy of the light waves. In order to study the impact of atmospheric turbulence on light wave propagation, researchers use differential image motion monitors (DIMMs) to measure seeing parameters, mainly at night. DIMMs play a key role in satellite-to-ground laser communications, making it possible to measure turbulence in the entire layer of the atmosphere.

[0004] The applicant, the Institute of Space Information Innovation of the Chinese Academy of Sciences, disclosed in its patent application document "A real-time calculation device and method for atmospheric optical turbulence based on ground measurement" (application date: June 21, 2023, application number: 202310740978.0, publication number: CN116484167B) that by measuring and counting the proportion of low-altitude atmospheric optical turbulence in the entire atmospheric turbulence, real-time atmospheric turbulence intensity measurement can be indirectly achieved. This atmospheric turbulence measurement method has the following shortcomings: the spatial coverage of the measurement is limited, and it cannot provide information on the turbulence parameters of the entire layer of the atmosphere.

[0005] The applicant, Yunnan Observatory of the Chinese Academy of Sciences, disclosed in its patent application document "Method for measuring atmospheric turbulence in a local area based on fluctuations in the laser wavefront arrival angle" (application date: May 20, 2024, application number: 202410623753.1, publication number: CN118244379B) that by measuring the offset of the center of the laser spot on different baselines, the covariance of the laser wavefront arrival angle of different baseline lengths is calculated to obtain the atmospheric turbulence profile and the distribution of atmospheric turbulence intensity with height. This atmospheric turbulence measurement method has the following shortcomings: This method can measure the distribution of turbulence intensity that changes with height, but is not suitable for multi-path measurement of the entire layer of atmospheric turbulence.

[0006] By capturing the faint changes in the light intensity and position fluctuations of stars, data on atmospheric channel transmission conditions in clear sky areas can be accurately collected. At the same time, the light spots of multiple stars can be extracted in the observation field, and the observation results of each star can reflect the intensity of atmospheric turbulence on the path. Summary of the invention

[0007] In order to address the challenges of limited spatial range of whole-layer atmospheric turbulence parameter measurement and insufficient positioning accuracy of star spot centroid in satellite-to-ground laser communications, the present invention proposes a whole-layer atmospheric turbulence measurement method based on centroid drift and light intensity scintillation, which aims to improve the robustness, accuracy and convenience of whole-layer atmospheric turbulence measurement.

[0008] The above purpose is achieved through the following technical solutions:

[0009] The present invention provides a method for measuring whole-layer atmospheric turbulence based on centroid drift and light intensity flicker, the steps comprising:

[0010] S1. Collect star images after atmospheric turbulence;

[0011] S2. Use Zernike moments to detect the edge of the star image to obtain the shape characteristics of the star image, separate the star spots from the background, and locate the edge information of the star points;

[0012] S3. Perform dynamic Lanczos interpolation to resample the edges of the star points;

[0013] S4. Use the normalized amplitude enhancement weight of the gray gradient to locate the centroid of the star spot on the interpolated star image;

[0014] S5. Calculate the intensity flicker based on the edge detection result of the star image using the Zernike moment in step S2, calculate the light spot centroid drift variance based on the centroid positioning result in step S4, and measure the turbulence intensity of the entire atmosphere based on the intensity flicker and the light spot centroid drift.

[0015] Furthermore, the specific implementation steps of collecting the star map after atmospheric turbulence in step S1 include:

[0016] The S1.1 uses a high-resolution imaging lens and a high-speed visible light camera to photograph stars;

[0017] S1.2 Continuously collect multiple frames of star images, and the acquired star spots occupy multiple pixels.

[0018] Furthermore, the step S2 of using the Zernike moment to perform edge detection on the star image to obtain the shape features of the star image includes the following specific implementation steps:

[0019] S2.1 After reading the star map, the Zernike moments of multiple pixels around each pixel are calculated, and the shape characteristics of the star point are obtained using the Zernike moments. The Zernike moments are complex numbers defined on the unit circle and decomposed into radial and angular functions, which are composed of a set of orthogonal functions to the unit circle. The internal complex polynomial is derived from the Nick moments of the two-dimensional continuous function with n rows and m columns. The calculation formula is as follows:

[0020] (1),

[0021] in, represents a continuous image function, is the coordinate point The polar coordinate form of From the origin to the coordinate point The vector length of is the angle of the vector in the counterclockwise direction, is a two-dimensional polynomial The conjugate polynomial of is expressed as:

[0022] (2),

[0023] where j is the mode number, is the Nick moment polynomial of n rows and m columns, which is defined as:

[0024] (3),

[0025] Among them, s is from 0 to , and Constraints should be met and ,and The value of is an even number. In digital images, the integral is replaced by the sum:

[0026] (4),

[0027] S2.2 Establish an ideal model for sub-pixel edge detection, there is a straight line contained in the unit circle , The grayscale values ​​on both sides are and , From the origin to the straight line The vertical distance for The angle between the direction and the x-axis; calculate the three Zernike moments ,in represents the average brightness of the star spot, Indicates the tilted shape of the star spot; Indicates the circularity of the star spot; then calculates four positioning parameters ,in Represents the background gray value, Indicates the step height, represents the vertical distance from the origin to the edge, Represents the angle between the vertical line and the horizontal line, thereby completing the edge positioning of the star spot, that is:

[0028] { d = Z 2 0 Z ′ 11 k = 3 Z ′ 11 2 ( 1 − d 2 ) 3 / 2 φ = t a n − 1 [ ( I m [ Z 1 1 ] ) / ( R e [ Z 1 1 ] ) ] h = [ Z 0 0 − k π 2 + k a r c s i n d + k d ( 1 − d 2 ) ] / π (5),

[0029] The Zernike moment of rotation is .

[0030] S2.3 calculates the Zernike moment by convolution of the template on the image. The larger the template size, the higher the solution accuracy. Assuming that the selected template size is N×N, N is the diameter of the circle, the sub-pixel position of the edge calculated by the sub-pixel edge detection algorithm of the Zernike moment is as follows:

[0031] (6),

[0032] in, is the sub-pixel coordinate of the edge, is the origin coordinate.

[0033] S2.4 uses the threshold segmentation method to separate the star spots from the background based on the grayscale difference between the star spots and the background in the star map. is the threshold value, when When , the pixel is an edge point, the boundary point is marked as 1, and the other points are marked as 0 to locate the edge information of the star point.

[0034] Furthermore, the step S3 of performing dynamic Lanczos interpolation on the edge area of ​​the star map to achieve resampling includes:

[0035] S3.1 takes the interpolation point as the center and obtains the edge area of ​​the star spot exceeding the grayscale threshold as a window;

[0036] S3.2 Use the Lanczos kernel function to calculate the weight of each pixel in the window for interpolation calculation:

[0037] Lanczos interpolation is performed by reconstructing the kernel of the filter definition:

[0038] (7),

[0039] in, is the normalized sinc function, is a positive integer, when =3 is Lanczos3 interpolation. By calculating the discrete convolution between the input image and the Lanczos3 kernel, the corresponding value of the expanded pixel can be obtained. The Lanczos3 interpolation formula is as follows:

[0040] (8),

[0041] in, Represents a variable, represents Lanczos3 interpolation, let is the image with coordinates Pixel points; Lanczos3 interpolation function in the edge area of ​​the star map It is expressed as follows:

[0042] (9),

[0043] S3.3 adds the values ​​of all pixels in the window multiplied by the weights to obtain the interpolation result of the star spot.

[0044] Furthermore, the centroid positioning is performed using the normalized amplitude enhancement weight of the grayscale gradient in step S4, and the specific implementation steps include:

[0045] S4.1 Use the grayscale centroid method to calculate the gradient value at the pixel and obtain the gradient matrix and grayscale matrix;

[0046] Assume the size of the star map is P×Q, and For each pixel, the horizontal and vertical coordinates are The output grayscale value of the CCD at each pixel, the theoretical centroid coordinates for:

[0047] (10),

[0048] After segmenting the interpolated star map, it is divided into two parts: inside and outside the circle. Then, based on the segmentation of the gradient and grayscale matrix, the new intensity value is calculated. for The gradient value at the pixel is the first-order derivative of the grayscale matrix. for:

[0049] (11),

[0050] for The partial derivatives of the components, for The partial derivatives of the components, is the amplitude of the gradient:

[0051] (12),

[0052] S4.2 Use the normalized amplitude value of the gradient to enhance the weight near the center of the star point and obtain a new gray value distribution;

[0053] Using the normalized distribution of gradients To further enhance the star point center The normalized amplitude matrices of the nearby weights, gradient matrices, and grayscale matrices follow the following formula:

[0054] (13),

[0055] This is the distribution of grayscale values ​​of the star map after Lanczos3 interpolation.

[0056] S4.3 The gray value distribution obtained , update the centroid coordinates after Lanczos3 interpolation, the updated centroid coordinates of the spot for:

[0057] (14).

[0058] Furthermore, the atmospheric turbulence intensity is measured based on the light intensity flicker and the light spot centroid drift in step S5, and the specific implementation steps include:

[0059] S5.1 Based on the centroid positioning result of step S4, the centroid drift variance of the star spot is calculated as , is the drift variance of the spot mass center in the horizontal direction, is the vertical drift variance, is the receiving aperture size, thus obtaining the spot centroid drift variance:

[0060] (15),

[0061] in, Indicates the transmission distance. represents the refractive index structure constant of the atmosphere;

[0062] S5.2 Calculate the light intensity flicker based on the edge detection result of the star image using the Zernike moment in step S2. The sum of the grayscale of all pixels on each frame of the image is used as the light intensity of the frame. The flicker index Determines the intensity fluctuation, that is, the normalized light intensity fluctuation variance, and inverts to derive the Rytov variance , also known as the logarithmic light intensity fluctuation variance, the logarithmic light intensity fluctuation variance is an evaluation index of the atmospheric turbulence intensity. It is expressed as:

[0063] (16),

[0064] Where < > represents the statistical average, I represents the total grayscale value of the image taken by the camera, represents the distance from the observation point to the optical axis, Indicates light intensity No. Moment, and The relationship between them is expressed as:

[0065] (17),

[0066] in, Indicates the variance of light intensity fluctuations, namely Rytov variance; represents the wave number of light waves, represents the wavelength of light, From the Kolmogorov refractive index fluctuation power spectrum, we get:

[0067] (18),

[0068] in, The light intensity without atmospheric turbulence is calculated by using the flicker factor of the light intensity, i.e. the normalized fluctuation variance, to evaluate the fluctuation of the light intensity.

[0069] S5.3 Based on the centroid drift and light intensity flicker of the spot, the turbulence intensity of the whole atmosphere is measured, that is, the atmospheric refractive index structure constant is calculated. :

[0070] (19).

[0071] Beneficial effects:

[0072] Traditional methods may not be able to evaluate the impact of the entire atmospheric turbulence on system performance from a wide area, have limitations in the measurement range, and may introduce distortion when processing data. The present invention uses imaging technology, image processing algorithms, and atmospheric optical principles to provide an effective solution for measuring atmospheric turbulence parameters in satellite-to-ground laser communications. Compared with the prior art, the present invention has the following beneficial effects:

[0073] The present invention is based on a measurement method of light intensity flicker and spot centroid drift, and aims to obtain the whole-layer atmospheric turbulence intensity of the satellite-to-ground laser communication system. This method can expand the space observation range, increase the convenience of wide-area measurement, and improve the robustness and accuracy of the whole-layer atmospheric turbulence measurement, providing an important basis and support for optimizing the satellite-to-ground laser communication system. By adopting Zernike moment edge detection, it is possible to effectively describe shape characteristics, suppress aliasing, retain details and reduce distortion caused by interpolation. At the same time, dynamic Lanczos3 interpolation and grayscale gradient normalized amplitude enhancement weights are used to improve anti-interference ability and positioning accuracy. The present invention further lays a foundation for research on multipath measurement of the whole-layer atmospheric turbulence intensity, anisotropy of atmospheric turbulence, and spatiotemporal correlation of atmospheric turbulence intensity. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 A flow chart is designed for a whole-layer atmospheric turbulence measurement method based on mass center drift and light intensity flickering of the present invention;

[0075] Figure 2 It is a schematic diagram of measuring the whole-layer atmospheric turbulence of the present invention;

[0076] Figure 3 It is a schematic diagram of the transmission of the star point through the atmospheric turbulence of the present invention;

[0077] Figure 4 In the figure, (a) is a schematic diagram of the gray value matrix of the present invention, and (b) is a schematic diagram of the sub-pixel distribution of Gaussian pixels;

[0078] Figure 5 This is a schematic diagram of the star spot pixel and centroid positioning of the present invention: Figure 4 In the figure, (a) is the star point pixel with the first brightness of the spot, (b) is the centroid location of the first brightness of the spot, (c) is the star point pixel with the second brightness of the spot, (d) is the centroid location of the second brightness of the spot, (e) is the star point pixel with the third brightness of the spot, (f) is the centroid location of the third brightness of the spot, (g) is the star point pixel with the fourth brightness of the spot, and (h) is the centroid location of the fourth brightness of the spot;

[0079] Figure 6 In the figure, (a) is a probability density distribution diagram of the flicker factor of the present invention, and (b) is a probability density distribution diagram of the average light intensity;

[0080] Figure 7 For the present invention based on Vega Measurement results. DETAILED DESCRIPTION

[0081] By combining the drawings in the embodiments of the present invention, the technical scheme is clearly and completely described. It should be pointed out that the embodiments only represent a part of the present invention, not all of it. According to these embodiments, other embodiments obtained by ordinary technicians without innovative work are all within the scope of protection of the present invention.

[0082] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0083] like Figure 1 As shown, it is a schematic diagram of the method flow of an embodiment of the present invention. In order to realize the measurement of the whole layer of atmospheric turbulence between the satellite and the ground, step 001 first uses a visible light camera and a high-resolution focal length lens to collect a star map. Step 002 uses the Zernike moment to detect the edge of the star spot to obtain the shape characteristics of the star point. Step 003 performs dynamic Lanczos3 interpolation on the inner edge area of ​​the star spot to achieve resampling and retain the details of the star map; step 004 uses the normalized amplitude enhancement weight method of the grayscale gradient to improve the accuracy of the center of mass positioning, and achieve better anti-interference ability and higher positioning accuracy; step 005 calculates the turbulence intensity of the whole layer of the atmosphere under the current communication link by the light intensity flicker variance and the center of mass drift variance. The specific steps include:

[0084] Step 1: Use a visible light camera and a high-resolution lens to capture star images that pass through atmospheric turbulence.

[0085] Identify the location of the satellite-to-ground laser communication test, determine the meteorological conditions at the time, as well as the wavelength range, image size and field of view of the camera observation, and set the optimal exposure time, sampling frequency and other parameters of the observation camera.

[0086] In this embodiment, the collection location is Changchun, and a high-resolution imaging lens with a focal length of 14mm, 50mm or 75mm and a high-speed visible light camera are used as the receiving device. The camera resolution can be set to 1920×1200, 2048×1536 pixels, the exposure time is 60ms, and the sampling frequency is 16.65Hz. This camera has a USB 3.0 interface, and the response at a wavelength of 550nm is more significant, which provides good conditions for capturing starlight spot images after atmospheric turbulence.

[0087] like Figure 2As shown, it is a schematic diagram of the whole-layer atmospheric turbulence measurement architecture provided in this embodiment. When light passes through the atmosphere, the position of the light spot will drift due to the change of the refractive index caused by atmospheric turbulence, and the intensity of the light will also fluctuate randomly. The light intensity flicker provides information about the change of turbulence intensity, and the light spot centroid drift provides information about the turbulence position. The method based on light intensity flicker and light spot centroid drift can make comprehensive use of these two observation methods. Multiple star spots can be extracted from the observed field of view, and the observation results of each star can reflect the atmospheric turbulence intensity on the path. By observing and analyzing the light spots of multiple stars, the atmospheric turbulence conditions in different directions and angles can be comprehensively considered, which reflects the atmospheric turbulence intensity in the wide-area transmission space to a certain extent.

[0088] Continuously collect multiple frames of star images at the communication time, and consider the impact of atmospheric turbulence on laser amplitude and phase. By observing and recording the intensity flickering and center of mass drift of star spots, we can have a more comprehensive understanding of the impact of atmospheric turbulence and increase the spatial coverage of the equipment to obtain star images. Figure 3 The star point transmission diagram shown in Figure 1 shows the process of the star point undergoing channel fading. In this process, channel fading introduces two main types of noise, namely multiplicative noise and additive noise. The attenuation and distortion effects of the star point signal during transmission will cause the amplitude and phase of the star point signal to change. Additive noise will be added to the star point signal during transmission, thus affecting the quality and reliability of the signal.

[0089] Step 2: Use Zernike moments to perform edge detection on the star image to obtain the shape features of the star image.

[0090] After reading the star map, the Zernike moments of multiple pixels around each pixel are calculated, and the shape characteristics of the star point are obtained using the Zernike moments. The Zernike moments are complex numbers defined on the unit circle and can be decomposed into radial and angular functions, which are composed of a set of orthogonal functions to the unit circle. The internal complex polynomial is derived from the Nick moments of the two-dimensional continuous function with n rows and m columns. The calculation formula is as follows:

[0091] (1),

[0092] in, represents a continuous image function, is the coordinate point The polar coordinate form of From the origin to the coordinate point The vector length of is the angle of the vector in the counterclockwise direction, is a two-dimensional polynomial The conjugate polynomial of is expressed as:

[0093] (2),

[0094] where j is the mode number, is the Nick moment polynomial of n rows and m columns, which is defined as:

[0095] (3),

[0096] Among them, s is from 0 to , and Constraints should be met and ,and The value of is an even number. In digital images, the integral is replaced by the sum:

[0097] (4).

[0098] Establish an ideal model for sub-pixel edge detection, such as Figure 4 The figure shows the gray value matrix diagram and Gaussian pixel sub-pixel distribution diagram, which describes the sub-pixel distribution inside the Gaussian pixel and the process of distributing the star points on multiple pixels through the Gaussian distribution function. Based on the pixel grayscale histogram of the star point, a local area of ​​the original image with a higher resolution is displayed. Represents a single pixel, , represent the rows and columns of the matrix respectively, , Representing the vertical and horizontal directions respectively, by dispersing the star points to multiple pixels, the positions of the star points can be determined more accurately, thereby improving the positioning accuracy.

[0099] There exists a straight line contained in the unit circle , The grayscale values ​​on both sides are and , From the origin to the straight line The vertical distance for The angle between the direction and the x-axis; calculate the three Zernike moments ,in represents the average brightness of the star spot, Indicates the tilted shape of the star spot; Indicates the circularity of the star spot; then calculates four positioning parameters ,in Represents the background gray value, Indicates the step height, represents the vertical distance from the origin to the edge, Represents the angle between the vertical line and the horizontal line, thereby completing the edge positioning of the star spot, that is:

[0100] { d = Z 2 0 Z ′ 11 k = 3 Z ′ 11 2 ( 1 − d 2 ) 3 / 2 φ = t a n − 1 [ ( I m [ Z 1 1 ] ) / ( R e [ Z 1 1 ] ) ] h = [ Z 0 0 − k π 2 + k a r c s i n d + k d ( 1 − d 2 ) ] / π (5),

[0101] Among them, the Zernike moment of rotation .

[0102] The Zernike moment is calculated by convolution of the template with the image. The larger the template size, the higher the solution accuracy. Assuming that the selected template size is N×N dimensions, the sub-pixel position of the edge calculated by the sub-pixel edge detection algorithm of the Zernike moment is as follows:

[0103] (6),

[0104] in, is the sub-pixel coordinate of the edge, The threshold segmentation method is used to separate the star spots from the background according to the grayscale difference between the star spots and the background in the star map. is the threshold value. The condition threshold value can be set between 0 and 255 as required, such as the grayscale threshold value is 210. When , the pixel is an edge point. Mark the boundary point as 1 and the other points as 0 to more accurately locate the edge information of the star point.

[0105] like Figure 5 The figure shows a schematic diagram of star point pixel and centroid positioning, showing the centroid positioning of four groups of different star spot intensities. The present invention uses Zernike moments to detect the edge of the star map to obtain the shape characteristics of the star map and provide a more stable and effective description. Dynamic Lanczos3 interpolation is applied in the edge area of ​​the star map to retain details and reduce distortion caused by interpolation. The normalized amplitude enhancement weight of the grayscale gradient is used to perform centroid positioning to enhance positioning accuracy.

[0106] Step 3: Perform dynamic Lanczos3 interpolation to resample the edge area of ​​the star point.

[0107] The present invention uses Lanczos-3 interpolation technology to dynamically interpolate the edge area for upsampling. Lanczos is a spatial domain interpolation technology. Dynamic Lanczos3 interpolation is performed on the edge area of ​​the star map to achieve resampling. With the interpolation point as the center, a window is selected in the image. For each pixel in the target image, interpolation calculation is performed based on the surrounding neighboring pixels and the Lanczos filter. The Lanczos interpolation is composed of the reconstruction kernel of the filter. definition:

[0108] (7),

[0109] in, is the normalized sinc function, is a positive integer, when =3 is Lanczos3 interpolation. By calculating the discrete convolution between the input image and the Lanczos3 kernel, the corresponding value of the expanded pixel can be obtained.

[0110] Lanczos3 interpolation function in the edge area of ​​the star map It is expressed as follows:

[0111] (8),

[0112] in, Represents a variable, represents Lanczos3 interpolation, let is the image with coordinates Pixel points; Lanczos3 interpolation function in the edge area of ​​the star map It is expressed as follows:

[0113] (9),

[0114] The pixel values ​​of all star spots are multiplied by the weights, and then the results are added to obtain the interpolation value of the target pixel. This can effectively suppress aliasing, preserve details, and reduce distortion caused by interpolation.

[0115] Step 4: Use the normalized amplitude enhancement weight of the grayscale gradient to locate the centroid.

[0116] Assume that the coordinates of the theoretical center of mass are , the star map size is P×Q. and For each pixel, the horizontal and vertical coordinates are The output grayscale value of the CCD at each pixel, the theoretical centroid coordinates for:

[0117] (10),

[0118] After segmenting the interpolated star map, it is divided into two parts: the inner part and the outer part. Then, based on the segmentation of the gradient and grayscale matrix, the new intensity value is calculated. Using a two-dimensional Gaussian function Theoretical physics model as a star point:

[0119] =

[0120] (20),

[0121] in, is the proportional coefficient representing the grayscale amplitude. , They are and The variance of the direction, is the correlation coefficient, for The gradient value at the pixel is the first-order derivative of the grayscale matrix. for:

[0122] (11),

[0123] for The partial derivatives of the components, for The partial derivatives of the components, is the amplitude of the gradient:

[0124] (12),

[0125] Using the normalized distribution of gradients To further enhance the star point center The normalized amplitude matrices of the nearby weights, gradient matrices, and grayscale matrices follow the following formula:

[0126] (13),

[0127] This is the distribution of grayscale values ​​of the star map after Lanczos3 interpolation.

[0128] The gray value distribution obtained by , update the centroid coordinates after Lanczos3 interpolation, the updated centroid coordinates of the spot for:

[0129] (14)

[0130] Step 5: Atmospheric turbulence measurement based on light intensity flicker and light spot centroid drift.

[0131] The current center of mass coordinates have been obtained from the previous steps. Due to the fluctuation of the arrival angle, the center of mass drifts. It is necessary to calculate the center of mass displacement variance of the spot according to the star images of the previous and next frames: , is the drift variance of the spot mass center in the horizontal direction, is the vertical drift variance, is the receiving aperture size, thus obtaining the spot centroid drift variance:

[0132] (15),

[0133] Then calculate the light intensity flicker for the star point area detected by Zermike moment edge detection. The sum of the grayscale of all pixels on each frame of the image can be used as the light intensity of the frame. Flicker Index Determines the intensity fluctuation, that is, the normalized light intensity fluctuation variance, and inverts to derive the Rytov variance , also known as the logarithmic light intensity fluctuation variance, the logarithmic light intensity fluctuation variance is an evaluation index of the atmospheric turbulence intensity. It is expressed as:

[0134] (16),

[0135] Where < > represents the statistical average, I represents the total grayscale value of the image taken by the camera, represents the distance from the observation point to the optical axis, Indicates light intensity No. Moment, and The relationship between them is expressed as:

[0136] (17),

[0137] in, Indicates the variance of light intensity fluctuations, namely Rytov variance; represents the wave number of light waves, represents the wavelength of light, From the Kolmogorov refractive index fluctuation power spectrum, we get:

[0138] (18),

[0139] in, is the light intensity when there is no atmospheric turbulence, such as Figure 6 The following is a probability density distribution diagram of the flicker factor and average light intensity. The value range of the flicker factor and average light intensity is divided into a certain number of equal intervals, and the number of items falling in each interval is calculated. The central values ​​of all intervals form a sequence. , the number of flicker factors and average light intensities in all intervals also forms a sequence , we can get the histogram of the sample by taking X and Y as the horizontal and vertical coordinates. It shows the characteristics of Gaussian distribution, which is in line with expectations.

[0140] Finally, the variance of light intensity scintillation and centroid drift can be obtained to obtain the turbulence intensity of the entire atmosphere under the current space link, that is, the atmospheric refractive index structure constant. :

[0141] (19)

[0142] like Figure 7 As shown, the measurement results of Vega are provided. Figure 7 The data collection in Changchun is shown in Figure 2. Value is 10 -17 ~10 -16 Fluctuates between 2.308×10 -17 , showing weak turbulence characteristics. In contrast, the data collected in Changchun on October 11, 2024, observed Value is 10 -16 ~10 -15 Fluctuates between 2.693×10 -16 , showing moderate turbulence characteristics, and the measurement results reveal the changes in atmospheric turbulence at different time points.

Claims

1. A method for measuring whole-layer atmospheric turbulence based on centroid drift and light intensity flicker, characterized in that the steps include: S1. Collect star images after atmospheric turbulence; S2. Use Zernike moments to detect the edge of the star image to obtain the shape characteristics of the star image, separate the star spots from the background, and locate the edge information of the star points; S3. Perform dynamic Lanczos interpolation to resample the edges of the star points; S4. Using the normalized amplitude enhancement weight of the gray gradient to locate the centroid of the star spot on the interpolated star image; S5. Calculate the intensity flicker according to the edge detection result of the star map using the Zernike moment in step S2, calculate the variance of the spot mass center drift according to the mass center positioning result in step S4, and measure the turbulence intensity of the entire atmosphere based on the intensity flicker and the spot mass center drift; The step S2 uses the Zernike moment to perform edge detection on the star image to obtain the shape features of the star image. The specific implementation steps include: S2.1 After reading the star map, the Zernike moments of multiple pixels around each pixel are calculated, and the shape characteristics of the star point are obtained using the Zernike moments. The Zernike moments are complex numbers defined on the unit circle and decomposed into radial and angular functions. They are composed of a set of orthogonal to the unit circle U:x 2 +y 2 ≤1 is derived from the complex polynomial inside, the Nick moment Z of the two-dimensional continuous function with n rows and m columns nm The calculation formula is as follows: Where f(x, y) represents the continuous image function, (ρ, θ) is the polar coordinate form of the coordinate point (x, y), ρ is the length of the vector from the origin to the coordinate point (x, y), and θ is the angle of the vector in the counterclockwise direction. is a two-dimensional polynomial V nm The conjugate polynomial of (ρ, θ) is expressed as: V nm (ρ,θ)=R nm (p)e jmθ (2) Where j is the mode number, R nm (ρ) is the Nick moment polynomial with n rows and m columns, which is defined as: Among them, s is from 0 to m and n should satisfy the constraints mod(nm, 2) = 0 and 0 ≤ m ≤ n, and the value of n - |m| is an even number. In digital images, the integral is replaced by the summation: S2.2 Establish an ideal model for sub-pixel edge detection. There is a straight line L contained in the unit circle. Circle , L Circle The grayscale values ​​on both sides are a and a+k respectively, and d is the distance from the origin to the straight line L Circle The vertical distance is the angle between the d direction and the x axis; calculate the three Zernike moments Z nm (Z 00 , Z 11 , Z 20 ), where Z 00 represents the average brightness of the star spot, Z 11 Represents the tilted shape of the star spot; Z 21 Indicates the circularity of the star spot; then calculate the four positioning parameters Where h represents the background gray value, k represents the step height, and d represents the vertical distance from the origin to the edge. Represents the angle between the vertical line and the horizontal line, thereby completing the edge positioning of the star spot, that is: The Zernike moment of rotation is S2.3 calculates the Zernike moment by convolution of the template on the image. The larger the template size, the higher the solution accuracy. Assuming that the selected template size is N×N, N is the diameter of the circle, the sub-pixel position of the edge calculated by the sub-pixel edge detection algorithm of the Zernike moment is as follows: Among them, (x s ,y s ) is the sub-pixel coordinate of the edge, (x0, y0) is the origin coordinate; S2.4 uses the threshold segmentation method to separate the star spots from the background based on the grayscale difference between the star spots and the background in the star map. t is the threshold value, when d≤d t When , the pixel is an edge point, the boundary point is marked as 1, and the other points are marked as 0, and the edge information of the star point is located; In step S3, dynamic Lanczos interpolation is performed in the edge area of ​​the star map to achieve resampling, and the specific implementation steps include: S3.1 takes the interpolation point as the center and obtains the edge area of ​​the star spot exceeding the grayscale threshold as a window; S3.2 uses the Lanczos kernel function to calculate the weight of each pixel in the window for interpolation calculation: The Lanczos interpolation is defined by the reconstruction kernel L(x) of the filter: Wherein, sin c(x) is a standardized sinc function, a is a positive integer, and when a=3, it is Lanczos3 interpolation. By calculating the discrete convolution between the input image and the Lanczos3 kernel, the corresponding value of the expanded pixel can be obtained; the Lanczos3 interpolation formula is as follows: Among them, x Lanczos represents a variable, L(x Lanczos ) represents Lanczos3 interpolation, let s uv is the pixel with coordinates (u, v) in the image; Lanczos3 interpolation function s in the edge area of ​​the star map Lanczos3 (x Lan ,y Lan ) is expressed as follows: S3.3 adds the values ​​of all pixels in the window and the value multiplied by the weight to obtain the interpolation result of the star spot; The centroid positioning is performed using the normalized amplitude enhancement weight of the grayscale gradient in step S4, and the specific implementation steps include: S4.1 uses the grayscale centroid method to calculate the gradient value at the pixel and obtain the gradient matrix and grayscale matrix; Assume that the size of the star map is P×Q, x i and j For each pixel, the horizontal and vertical coordinates are given by f(x i ,y j ) is the output grayscale value of the CCD at each pixel, and the theoretical centroid coordinate (x theory ,y theory )for: After segmenting the interpolated star map, it is divided into two parts: inside and outside the circle. Then, based on the segmentation of the gradient and grayscale matrix, the new intensity value is calculated. For Lanczos3 (x Lan ,y Lan )The gradient value at the pixel is the first-order derivative of the grayscale matrix, for: g x For x Lan Partial derivatives of components, g y for y Lan The partial derivative of the component, M(x GV ,y GV ) is the amplitude of the gradient: S4.2 uses the normalized amplitude value of the gradient to enhance the weight near the center of the star point to obtain a new gray value distribution; Using the normalized distribution M of the gradient normalized (x Lan ,y Lan ) to further enhance the star point center (x Center ,y Center ), the normalized amplitude matrices of the gradient matrix and the grayscale matrix follow the following formula: S(x Lan ,y Lan ) is the distribution of gray values ​​of the star map after Lanczos3 interpolation; S4.3 The gray value distribution S(x Lan ,y Lan ), update the centroid coordinates after Lanczos3 interpolation, the updated centroid coordinates of the spot (x Current ,y Current )for: 。 2. A method for measuring whole-layer atmospheric turbulence based on centroid drift and light intensity scintillation according to claim 1, characterized in that The specific implementation steps of collecting the star map after atmospheric turbulence in step S1 include: The S1.1 uses a high-resolution imaging lens and a high-speed visible light camera to photograph stars; S1.2 continuously collects multiple frames of star images, and the acquired star spots occupy multiple pixels.

3. The whole-layer atmospheric turbulence measurement method based on centroid drift and light intensity scintillation according to claim 1 is characterized in that The step S5 of measuring the atmospheric turbulence intensity based on the light intensity flicker and the light spot centroid drift includes the following specific implementation steps: S5.1 Calculate the centroid drift variance of the star spot using the centroid positioning result of step S4: is the drift variance of the spot mass center in the horizontal direction, is the drift variance in the vertical direction, and D is the receiving aperture size, thus obtaining the drift variance of the spot centroid: Where L represents the transmission distance, represents the refractive index structure constant of the atmosphere; S5.2 Calculate the light intensity flicker based on the edge detection result of the star image using the Zernike moment in step S2. The sum of the grayscale of all pixels on each frame of the image is taken as the light intensity of the frame. The flicker index Determines the intensity fluctuation, that is, the normalized light intensity fluctuation variance, and inverts to derive the Rytov variance Also known as logarithmic light intensity fluctuation variance, logarithmic light intensity fluctuation variance is an evaluation index of atmospheric turbulence intensity. It is expressed as: Among them, <> represents the statistical average, I represents the sum of the grayscale values ​​of the image taken by the camera, and r represents the distance from the observation point to the optical axis. 2 > indicates the second th Moment, and The relationship between them is expressed as:​ in, Indicates the variance of light intensity fluctuations, namely Rytov variance; represents the wave number of light waves, λ represents the wavelength of light, From the Kolmogorov refractive index fluctuation power spectrum, we get: Among them, I0 is the light intensity when there is no atmospheric turbulence, and the light intensity fluctuation is evaluated using the scintillation factor of the light intensity, that is, the normalized fluctuation variance; S5.3 measures the turbulence intensity of the entire atmosphere based on the centroid drift of the light spot and the scintillation of the light intensity, that is, calculates the atmospheric refractive index structure constant

Images