Cutting height offset calculation method

Abstract

The invention relates to the technical field of atmospheric remote sensing detection, in particular to a cut height offset calculation method, which comprises the following steps of: acquiring limb scattered radiation corresponding to a plurality of predetermined engineering control cut heights; randomly initializing a plurality of cut-up particles, wherein the position of each cut-up particle represents a guess value of cut-up offset; defining a fitness function; multiple rounds of iteration are executed, and in each round of iteration, for each cut-up particle, the position of the cut-up particle is calculated according to the individual extreme value and the global extreme value of the cut-up particle; based on a fitness function and the updated position of each cut-up particle, calculating a fitness value corresponding to each cut-up particle in the current iteration, and updating an individual extreme value of each cut-up particle and a global extreme value of all cut-up particles; and after a preset iteration ending condition is met, taking the global extreme value of the last iteration as the finally obtained cut-up offset. According to the method, the precision and the range of an inversion result are effectively improved, complex model calculation is avoided, and meanwhile, the accuracy of an optimization process is improved.

Description

Technical Field

[0001] The invention relates to the technical field of atmospheric remote sensing detection, and in particular to a method for calculating a cutting height offset. Background Art

[0002] In satellite remote sensing, the tangent height refers to the shortest distance between the tangent line of the earth's surface in satellite remote sensing limb detection technology and the tangent line of the atmospheric height observed by satellite instruments, that is, the vertical distance from the observation point to the atmosphere. Figure 1 As shown, Figure 1 SCIAMACHY stands for Scanning Imaging Absorption Spectrometer for Atmospheric Composition, and Envisat stands for a satellite launched by the European Space Agency for earth observation. The measurement of tangent height is an important remote sensing observation method, which is widely used in atmospheric and climate research, such as atmospheric composition analysis and climate change research. Due to various complex atmospheric effects and instrument factors, such as atmospheric refraction and instrument parameter deviation, the observed tangent height often deviates from the actual tangent height. This deviation will affect the accuracy of atmospheric detection, especially for hyperspectral atmospheric detection instruments. The accuracy of tangent height has an important impact on the quality of observation data and subsequent atmospheric property inversion results. Therefore, calculating the tangent height offset and correcting the tangent height are of vital importance to improving the accuracy and reliability of atmospheric detection data.

[0003] The existing methods for calculating the cut height offset are mainly based on theoretical calculations of atmospheric models. However, theoretical calculations based on atmospheric models are often difficult to accurately reflect the actual atmospheric state due to the limitations of the accuracy and complexity of the models. Summary of the invention

[0004] In order to solve the above problem, the present invention provides a method for calculating cutting height offset.

[0005] The method for calculating the cutting height offset provided by the present invention comprises:

[0006] Get multiple predetermined engineering control cut heights h i The corresponding limb scattered radiation;

[0007] Randomly initialize multiple height-cutting particles, and the position of each height-cutting particle represents the guessed value of the height-cutting offset Δh;

[0008] The fitness function f is defined according to all engineering control cut heights, the edge scattered radiation corresponding to each engineering control cut height, and the cut height offset Δh;

[0009] Perform multiple iterations, and in each iteration:

[0010] For each cutting height particle, the speed and position of the cutting height particle are calculated according to the individual extreme value of the cutting height particle and the global extreme value of all cutting height particles;

[0011] Based on the fitness function f and the updated position of each height-cutting particle, the fitness value corresponding to each height-cutting particle in the current iteration is calculated;

[0012] According to the fitness value corresponding to each height-cutting particle in the current iteration, the individual extreme value of each height-cutting particle and the global extreme value of all height-cutting particles are updated;

[0013] After the preset iteration end condition is met, the global extreme value of the last iteration is used as the final cutting height offset.

[0014] Furthermore, according to all engineering control cut heights, the edge scattered radiation corresponding to each engineering control cut height, and the cut height offset Δh, the step of defining the fitness function f includes:

[0015] Define the fitness function f as:

[0016]

[0017] in, is the normalized radiation corresponding to the i-th engineering control cut height, is the normalized inverted radiation corresponding to the ith engineering control cut height, and I represents the number of engineering control cut heights.

[0018] Furthermore, the normalized radiation corresponding to the i-th engineering control cut height is for:

[0019]

[0020] Where λ is the wavelength, h i represents the height of the i-th engineering control, R(λ, h i ) represents the edge scattered radiance with wavelength λ at the i-th engineering control cut height, R(λ, h ref ) represents the wavelength λ at the predefined reference cut height h ref The limb scattered radiance at .

[0021] Furthermore, the normalized inversion radiation corresponding to the i-th engineering control cutting height is The way to obtain is:

[0022] The inversion cutting height corresponding to the i-th engineering control cutting height is calculated according to the following formula

[0023]

[0024] Among them, h i Indicates the control height of the i-th engineering project;

[0025] Using the preset radiative transfer model, according to Calculate the inversion radiation corresponding to the i-th engineering control cut height

[0026] Calculate the normalized inversion radiation corresponding to the i-th engineering control cut height

[0027]

[0028] in, Represents the inverted radiation corresponding to the predefined reference cut height.

[0029] Furthermore, the step of calculating the speed and position of the height-cutting particle according to the individual extreme value of the height-cutting particle and the global extreme value of all height-cutting particles includes:

[0030] The velocity of the j-th cutting particle at the k+1th iteration for:

[0031]

[0032] Where c1 and c2 are the preset first learning rate and second learning rate, rand1 and rand2 represent different random numbers, and p j represents the individual extreme value of the jth cutting height particle, represents the global extreme value at the kth iteration, ω k represents the inertia weight factor of the kth iteration, represents the position of the j-th cutting height particle at the k-th iteration;

[0033] The position of the j-th cutting particle at the k+1th iteration for:

[0034] Furthermore, the inertia weight factor ω of the kth iteration is k for:

[0035]

[0036] Among them, ω max Represents the maximum value of the preset inertia weight factor, ω min Indicates the minimum value of the preset inertia weight factor, K max Indicates the preset maximum number of iterations.

[0037] Furthermore, based on the fitness function f and the updated position of each height-cutting particle, the step of calculating the fitness value corresponding to each height-cutting particle in the current iteration includes:

[0038] The updated position of each cutting height particle is used as the cutting height offset Δh;

[0039] Control the height h for each project i , according to the cutting height offset Δh, h is calculated i Corresponding inversion cut height

[0040] Inverse cutting height based on all engineering control cutting heights And the fitness function f, calculate the fitness value of each cutting height particle corresponding to the current iteration.

[0041] Furthermore, according to the fitness value corresponding to each height-cutting particle in the current iteration, the step of updating the individual extreme value of each height-cutting particle includes:

[0042] For each cut height particle:

[0043] Determine whether the fitness value corresponding to the current iteration is greater than the fitness value corresponding to the current individual extreme value.

[0044] Small value;

[0045] If yes, then the individual extreme value of the height-cut particle is updated to the updated position of the height-cut particle;

[0046] If not, the individual extreme value of the cut height particle is not changed.

[0047] Furthermore, according to the fitness value corresponding to each height-cutting particle in the current iteration, the step of updating the global extreme value of all height-cutting particles includes:

[0048] For each cut height particle:

[0049] Determine whether the fitness value corresponding to the current iteration is smaller than the fitness value corresponding to the current global extreme value;

[0050] If yes, then the global extreme value is updated to the updated position of the height-cut particle;

[0051] If not, the global extremum is not changed.

[0052] Furthermore, the preset iteration end condition refers to:

[0053] The preset maximum number of iterations is reached, or the change in the global extremum in consecutive z iterations is less than a preset threshold, where z is a preset positive integer constant.

[0054] One or more technical solutions provided in the embodiments of the present application have at least the following technical effects or advantages:

[0055] The method for calculating the shear height offset proposed in the present invention transforms the shear height inversion problem into an optimization problem, utilizes the limb radiation characteristics of the UV-B band in the upper stratosphere, and adopts a particle swarm optimization algorithm to perform shear height correction, thereby effectively improving the accuracy and range of the inversion result. By designing a fitness function and using the sum of squared errors as an individual quality metric, the optimization process is guided, thus avoiding complex model calculations and improving the accuracy of the optimization process. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 A schematic diagram of the edge observation geometry disclosed in the background technology of the present invention;

[0057] Figure 2 The SCIAMACHY 32980-32992 track cutting height diagram calculated based on the method of the present invention is provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0058] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. Before describing in detail the technical solutions of each embodiment of the present invention, the nouns and terms involved are explained. In this specification, components with the same name or the same number represent similar or identical structures and are only for illustrative purposes.

[0059] The present invention converts the shear height inversion problem into an optimization problem, utilizes the edge radiation characteristics of the UV-B band in the upper stratosphere, and adopts a particle swarm optimization algorithm to perform shear height correction. In addition, the normalized radiation samples at multiple engineering tangent heights are composed into an inversion vector as the target to be fitted in the optimization process. Then, a fitness function is designed, and the sum of squared errors is used as an individual quality metric to guide the optimization process. Finally, the shear height offset is updated by the particle swarm optimization algorithm. During the whole process, the shear height offset is continuously iteratively updated until the iterative convergence condition is met, thereby obtaining the optimized shear height offset and a more accurate inversion shear height result.

[0060] The present invention is explained below in conjunction with specific embodiments.

[0061] 1. Data Acquisition

[0062] A plurality of cut heights located in the upper stratosphere corresponding to the remote sensing instrument are obtained as the engineering control cut heights, and remote sensing data observed by the remote sensing instrument at each engineering control cut height. The present invention only uses data in the upper stratosphere to reduce interference from aerosols.

[0063] Remote sensing instruments are devices for limb observation, which can measure atmospheric limb scattered radiation at different tangent heights. These different tangent heights are called the engineering control tangent heights of the remote sensing instrument.

[0064] In this embodiment, the remote sensing instrument corresponds to five engineering control cut heights, namely 37km, 40km, 43km, 46km and 49km. Under each engineering control cut height, the remote sensing instrument will observe the edge scattered radiation of the UV-B band in the upper stratosphere as remote sensing data, including the wavelength of the ultraviolet band and the edge scattered radiance corresponding to the wavelength.

[0065] 2. Inversion vector definition

[0066] In order to optimize the calculation of the cut height offset based on the particle swarm optimization algorithm, the remote sensing data needs to be normalized to the reference cut height to reduce the influence of the lower atmosphere on the radiance.

[0067] The reference cutting height is any one of a plurality of engineering control cutting heights, and the selection method may be the engineering control cutting height with the smallest difference from the mean of all engineering control cutting heights. The reference cutting height selected in this embodiment is the third one of the five engineering control cutting heights: 43 km.

[0068] Normalize the limb scattered radiation corresponding to the i-th engineering control cut height to the reference cut height h ref On the top, we get the normalized radiation corresponding to the i-th engineering control cut height

[0069]

[0070] Where λ is the wavelength, h i represents the height of the i-th engineering control, R(λ, h i ) represents the edge scattered radiance observed by the remote sensing instrument at the wavelength λ at the i-th engineering control height, R(λ, h ref ) indicates that the remote sensing instrument is at a wavelength of λ and a reference height h ref The limb scattered radiance observed at .

[0071] Combine the normalized radiation corresponding to all engineering control cut heights into the observation vector Y obs :

[0072]

[0073] Wherein, I represents the number of engineering control cut heights. In this embodiment, the wavelength is selected to be 300 nm.

[0074] 3. Model Construction

[0075] The present invention transforms the cutting height inversion problem into an optimization problem. The goal of this optimization problem is to find an optimal cutting height offset Δh so that the simulation vector Y calculated based on the cutting height offset Δh is mod As close as possible to the observation vector Y obsIn this optimization problem, the height-cutting offset Δh is encoded as the position of the height-cutting particle in the particle swarm optimization algorithm, and the simulation vector Y mod and the observation vector Y obs The sum of squared errors is used to construct the fitness function f. In this way, the particle swarm optimization algorithm can find the optimal solution by iteratively updating the position of the particles, that is, the height offset Δh.

[0076] The cutting height offset Δh is:

[0077]

[0078] in, It represents the inverted cut height obtained by inverting the cut height of the i-th engineering control according to the present invention.

[0079] The fitness function f is:

[0080]

[0081] in, is the normalized radiation corresponding to the i-th engineering control cut height, is the normalized inverted radiation corresponding to the i-th engineering control cut height.

[0082] The present invention uses a radiation transfer model to invert the cutting height corresponding to the i-th engineering control cutting height. Calculate the corresponding inverted radiation The radiation transfer models include MODTRAN (MODerate resolution atmospheric TRANsmission, medium resolution atmospheric transfer model), SCIATRAN (radiative TRANsfer model for SCIAMACHY, SCIAMACHY radiation transfer model), etc. How to calculate the radiation corresponding to the cut height based on the radiation transfer model is common knowledge in the field and will not be elaborated in the present invention.

[0083] After that, the normalized inversion radiation corresponding to the i-th engineering control cut height is calculated

[0084]

[0085] in, Represents the inverted radiation corresponding to the predefined reference cut height.

[0086] 4. Find the optimal solution

[0087] The particle swarm optimization algorithm of the present invention adopts a full model, in which the position of each height-cutting particle represents a guess value of the height-cutting offset Δh, and the position of the jth height-cutting particle at the first iteration is for:

[0088]

[0089] Among them, x min represents the minimum position in the search space, x max represents the maximum position of the search space, and rand represents a random number in the range of [0, 1]. The minimum and maximum positions of the search space are both pre-set values ​​based on domain knowledge and historical data.

[0090] The velocity of the j-th cutting particle at the k+1th iteration for:

[0091]

[0092] Among them, c1 and c2 are the preset first learning rate and second learning rate, respectively, which control whether the cutting height particle is more inclined to track its individual optimal position or the global optimal position, rand1 and rand2 represent random numbers generated in the range of [0,1], and p j represents the individual extreme value of the jth cutting height particle, represents the global extreme value at the kth iteration, ω k represents the inertia weight factor of the kth iteration, Indicates the position of the jth cutting height particle at the kth iteration. The individual extreme value refers to the optimal position reached by a cutting height particle during the iteration process, that is, the position with the minimum corresponding fitness value, which is initialized to the initial position of the corresponding cutting height particle. The global extreme value refers to the optimal position reached by all cutting height particles during the iteration process, that is, the position of the cutting height particle with the minimum corresponding fitness value among the individual extreme values ​​of all cutting height particles. The fitness value refers to the value calculated according to the fitness function f.

[0093] Among them, the inertia weight factor ω k for:

[0094]

[0095] Among them, ω max Represents the maximum value of the preset inertia weight factor, ω min Indicates the minimum value of the preset inertia weight factor, K max Indicates the preset maximum number of iterations.

[0096] The updated position of the j-th cutting particle at the k+1th iteration for:

[0097]

[0098] In each iteration, each particle updates its speed and position according to the above formula. The speed update reflects the "social influence" (attraction of global extreme values) and "individual influence" (attraction of individual extreme values) on the particle, while taking into account the historical search direction of the particle (i.e. the inertia of the particle). The position update is based on the current position and speed of the particle.

[0099] After updating the position of the cutting height particle, the updated position of each cutting height particle is used as the cutting height offset Δh, and the cutting height is calculated inversely in combination with the corresponding engineering control cutting height. Based on the inversion cutting height Calculate the fitness function f corresponding to the updated position and obtain the fitness value corresponding to the height-cutting particle in this round of iteration.

[0100] Determine whether the fitness value corresponding to the updated position of a height-cutting particle is smaller than the fitness value corresponding to the current individual extreme value of the height-cutting particle. If so, update the individual extreme value of the height-cutting particle to the updated position of the height-cutting particle; if not, do not change the individual extreme value of the height-cutting particle.

[0101] Determine whether the fitness value corresponding to the updated position of a height-cutting particle is smaller than the fitness value corresponding to the current global extreme value. If so, update the global extreme value to the updated position of the height-cutting particle; if not, do not change the global extreme value.

[0102] This iterative process will continue until a certain convergence condition is met: the preset maximum number of iterations is reached or the change in the global extreme value in consecutive z iterations is less than a preset threshold, where z is a preset positive integer constant.

[0103] When the iteration is finished, the global extreme value corresponding to the last iteration is used as the height cutting offset finally obtained by the present invention.

[0104] 5. Effect comparison

[0105] Please refer to Table 1, which shows the average value of the SCIAMACHY (Scanning Imaging Absorption Spectrometer for Atmospheric CHartographY) 32980-32992 track cut height offset obtained by the method provided by the present invention and the cut height comparison result with the SCIAMACHY data provided by ESA (European Space Agency). The calculation method of the cut height correction accuracy Acc is:

[0106] Acc=|TH ret -TH ESA |;

[0107] Among them, TH ret is the cutting height calculated by the cutting height offset obtained by the method provided by the present invention, TH ESA The cut height of Level 1B (L1B) data provided by ESA.

[0108] Table 1 Average values ​​of cutting height offset and precision obtained by the present invention

[0109]

[0110]

[0111] It can be seen from Table 1 that the average offset of the height cut inversion of each track is highly consistent, indicating that the method of the present invention has good stability; the average accuracy of all tracks is 0.23km, which is less than 250m, meeting the requirements of height cut correction.

[0112] Figure 2 The SCIAMACHY 32980-32992 track cut height calculated according to the method of the present invention is shown, and each line segment represents a track. Figure 2 It can be seen that the method of the present invention is applicable in the orbital initial stage and in the latitude zone south of 30° in the southern hemisphere. The orbital initial stage generally refers to SZA>90°, SZA (Solar Zenith Angle).

[0113] In summary, the present invention converts the shear height offset calculation problem into an optimization problem, and only uses the edge scattering radiation characteristics at a specific wavelength (300nm in this embodiment) to improve the speed of shear height inversion; the shear height edge scattering radiation characteristics in the upper stratosphere are used to reduce the interference of aerosols; the radiation normalization calculation reduces the interference of instrument effects and stray light. The square sum of the errors between the normalized radiation and the normalized inverted radiation is used as the fitness function to evaluate the quality of particles, and has a high inversion accuracy.

[0114] The above-described embodiments are merely descriptions of preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A method for calculating cutting height offset, characterized in that: include: Get multiple predetermined engineering control cut heights h i The corresponding limb scattered radiation; Randomly initialize multiple height-cutting particles, and the position of each height-cutting particle represents the guessed value of the height-cutting offset Δh; The fitness function f is defined according to all engineering control cut heights, the edge scattered radiation corresponding to each engineering control cut height, and the cut height offset Δh; Perform multiple iterations, and in each iteration: For each cutting height particle, the speed and position of the cutting height particle are calculated according to the individual extreme value of the cutting height particle and the global extreme value of all cutting height particles; Based on the fitness function f and the updated position of each height-cutting particle, the fitness value corresponding to each height-cutting particle in the current iteration is calculated; According to the fitness value corresponding to each height-cutting particle in the current iteration, the individual extreme value of each height-cutting particle and the global extreme value of all height-cutting particles are updated; After the preset iteration end condition is met, the global extreme value of the last iteration is used as the final cutting height offset.

2. The method for calculating the cutting height offset according to claim 1, characterized in that: According to all engineering control cut heights, the edge scattered radiation corresponding to each engineering control cut height, and the cut height offset Δh, the steps of defining the fitness function f include: Define the fitness function f as: in, is the normalized radiation corresponding to the i-th engineering control cut height, is the normalized inverted radiation corresponding to the ith engineering control cut height, and I represents the number of engineering control cut heights.

3. The method for calculating the cutting height offset according to claim 2, characterized in that: Normalized radiation corresponding to the i-th engineering control cut height for: Where λ is the wavelength, h i represents the height of the i-th engineering control, R(λ, h i ) represents the edge scattered radiance with wavelength λ at the i-th engineering control cut height, R(λ, h ref ) represents the wavelength λ at the predefined reference cut height h ref The limb scattered radiance at .

4. The method for calculating the cutting height offset according to claim 2, characterized in that: Normalized inversion radiation corresponding to the i-th engineering control cutting height The way to obtain is: The inversion cutting height corresponding to the i-th engineering control cutting height is calculated according to the following formula Among them, h i Indicates the control height of the i-th engineering project; Using the preset radiative transfer model, according to Calculate the inversion radiation corresponding to the i-th engineering control cut height Calculate the normalized inversion radiation corresponding to the i-th engineering control cut height in, Represents the inverted radiation corresponding to the predefined reference cut height.

5. The method for calculating the cutting height offset according to claim 1, characterized in that: The step of calculating the speed and position of the cutting particle according to the individual extreme value of the cutting particle and the global extreme value of all cutting particles comprises: The velocity of the j-th cutting particle at the k+1th iteration for: Where c1 and c2 are the preset first learning rate and second learning rate, rand1 and rand2 represent different random numbers, and p j represents the individual extreme value of the jth cutting height particle, represents the global extreme value at the kth iteration, ω k represents the inertia weight factor of the kth iteration, represents the position of the j-th cutting height particle at the k-th iteration; The position of the j-th cutting particle at the k+1th iteration for:

6. The method for calculating the cutting height offset according to claim 5, characterized in that: The inertia weight factor ω of the kth iteration k for: Among them, ω max Represents the maximum value of the preset inertia weight factor, ω min Indicates the minimum value of the preset inertia weight factor, K max Indicates the preset maximum number of iterations.

7. The method for calculating the cutting height offset according to claim 1, characterized in that: Based on the fitness function f and the updated position of each height-cutting particle, the step of calculating the fitness value corresponding to each height-cutting particle in the current iteration includes: The updated position of each cutting height particle is used as the cutting height offset Δh; Control the height h for each project i , according to the cutting height offset Δj, h is calculated i Corresponding inversion cut height Inverse cutting height based on all engineering control cutting heights And the fitness function f, calculate the fitness value of each cutting height particle corresponding to the current iteration.

8. The method for calculating the cutting height offset according to claim 1, characterized in that: According to the fitness value corresponding to each height-cutting particle in the current iteration, the steps of updating the individual extreme value of each height-cutting particle include: For each cut height particle: Determine whether the fitness value corresponding to the current iteration is smaller than the fitness value corresponding to the current individual extreme value; If yes, then the individual extreme value of the height-cut particle is updated to the updated position of the height-cut particle; If not, the individual extreme value of the cut height particle is not changed.

9. The method for calculating the cutting height offset according to claim 1, characterized in that: According to the fitness value corresponding to each height-cutting particle in the current iteration, the steps of updating the global extreme value of all height-cutting particles include: For each cut height particle: Determine whether the fitness value corresponding to the current iteration is smaller than the fitness value corresponding to the current global extreme value; If yes, then the global extreme value is updated to the updated position of the height-cut particle; If not, the global extremum is not changed.

10. The method for calculating the cutting height offset according to claim 1, characterized in that: The preset iteration end conditions are: The preset maximum number of iterations is reached, or the change in the global extremum in consecutive z iterations is less than a preset threshold, where z is a preset positive integer constant.

Images