A method for retrieving auroral brightness based on satellite remote sensing data

Abstract

The application discloses an aurora brightness inversion method based on satellite remote sensing data, and comprises the following steps: constructing an atmospheric transmittance lookup table based on VIIRS night light data, inverting DNB band surface reflectivity and calculating moonlight radiation brightness; establishing a radiation transfer equation containing aurora uplink radiation, surface reflected aurora radiation and surface reflected moonlight radiation, and solving aurora real brightness pixel by pixel; and removing abnormal high brightness signals and urban artificial light interference, and obtaining spatially continuous and stable annual aurora brightness distribution through annual scale synthesis and spatial interpolation smoothing. The application can realize wide-range, quantitative and long-time sequence aurora brightness monitoring, effectively overcome the shortcomings of ground-based observation and numerical simulation, and provide reliable technical means for spatial weather research and polar environment monitoring.

Description

Technical Field

[0001] This invention relates to the field of Earth observation remote sensing data processing and spatial information technology, specifically to a method for quantitatively retrieving nighttime aurora brightness using multi-source satellite remote sensing data. This method is used to separate and quantitatively retrieve the true aurora radiance from VIIRS nighttime light remote sensing data for aurora signals that are self-luminous in the upper atmosphere. Background Technology

[0002] Auroras are a typical luminous phenomenon occurring in the upper atmosphere of high-latitude regions. They are formed when high-energy charged particles from the solar wind enter the atmosphere along geomagnetic field lines and collide with atmospheric components, producing radiant radiation. They are an important direct indicator of space weather activity. The spatial distribution and radiation intensity of auroras can reflect magnetosphere-ionosphere coupling processes, the intensity of high-energy particle deposition, and the degree of geomagnetic disturbance. They have significant scientific and applied value in space weather monitoring, geomagnetic storm process analysis, polar ionospheric environment research, and spacecraft safety operation assessment. Furthermore, as a unique natural spectacle, auroras have extremely high ornamental and tourism value; understanding their brightness distribution is also important for promoting the development of specialty tourism. Therefore, conducting large-scale, quantitative, and long-term continuous auroral brightness monitoring is a crucial foundational task in space environment and space weather research.

[0003] Currently, aurora brightness monitoring mainly relies on ground-based optical observation systems and numerical simulation methods based on physical mechanisms. Ground-based observations typically utilize equipment such as all-space imagers and spectrometers to acquire local aurora images or radiation information, offering the advantage of high observation accuracy. However, limitations such as sparse station deployment, limited observation fields of view, and cloud cover make it difficult to obtain large-scale, continuous aurora brightness distribution information and meet the long-term monitoring needs at regional or even global scales. On the other hand, aurora numerical simulation methods are usually based on magnetosphere-ionosphere coupling models and high-energy particle transport models, which are highly dependent on solar wind parameters, magnetic field structure, and ionospheric background fields. The model parameters are complex and have significant uncertainties. Furthermore, it is difficult to directly establish a pixel-scale correspondence between the simulation results and actual observation data, and it is also difficult to generate aurora brightness products with clear spatial resolution.

[0004] Auroras, as a weak, spontaneous radiation signal from the upper atmosphere, differ significantly from conventional nighttime low-light sources such as city lights, firelight, and moonlight reflection. For example, auroras are emitted by the atmosphere itself, not by surface reflection or artificial light sources, resulting in fundamentally different radiation transmission paths. Auroral radiation sources are located at high altitudes, making correction methods for surface-oriented light sources unsuitable. Auroral signals are weak and easily masked by strong interference from moonlight, surface reflection, town lights, and heat sources. Furthermore, auroras are transient and spatially discontinuous, making it impossible to preserve the true radiation intensity using conventional nighttime light synthesis and noise reduction methods. Therefore, conventional methods for nighttime low-light inversion, moonlight correction, and nighttime light product processing cannot accurately invert their true intensity and cannot be directly used for quantitative inversion of auroral brightness. Chinese invention patent CN108959347B discloses a method for determining auroral observation candidate areas based on topographic and nighttime light data, but it only filters observation locations and does not involve quantitative inversion of auroral radiation brightness, thus failing to obtain large-scale, continuously distributed auroral brightness data. Another researcher proposed a nighttime low-light remote sensing radiometric correction technique, in which low-light includes auroras. Specifically, it constructs a lunar spectral irradiance model based on VIIRS / DNB data to achieve accurate calculation of lunar radiant flux. However, this method is only used for lunar correction and cannot quantitatively separate the true brightness of auroras from nighttime light signals, thus failing to form a complete method for auroral brightness inversion and interference removal.

[0005] Overall, existing methods for monitoring aurora brightness are significantly inadequate in terms of spatial coverage, continuous monitoring capabilities, and observation consistency. For a long time, there has been a lack of a universal method that can quantitatively retrieve aurora brightness using operational satellite remote sensing data. Summary of the Invention

[0006] This invention aims to address the problems existing in the process of aurora brightness inversion based on nighttime light remote sensing data, such as surface reflection, moonlight illumination, atmospheric attenuation, nighttime combustion, and human activity interference. It proposes an aurora brightness inversion method based on satellite remote sensing data to achieve separation, interference suppression, and spatial continuous reconstruction of nighttime aurora radiation information, thereby improving the accuracy and stability of aurora brightness inversion.

[0007] To achieve the above-mentioned technical objectives, the technical solution adopted by the present invention is as follows:

[0008] A method for inverting aurora brightness based on satellite remote sensing data, the method comprising the following steps:

[0009] Step A: Based on aerosol optical thickness and observation geometric parameters, construct an atmospheric transmittance lookup table to retrieve the atmospheric transmittance τ(θ) in the satellite observation direction. S Atmospheric transmittance τ(θ) in the direction of moonlight incidence M Atmospheric transmittance τ(θ) in the direction of aurora radiation A );

[0010] Step B involves constructing a DNB-band surface reflectance inversion model based on multispectral reflectance to obtain time-continuous surface reflectance data. ;

[0011] Step C: Calculate the lunar radiance at the top of the atmosphere based on the lunar phase angle, Earth-Moon distance, and lunar optical parameters. ;

[0012] Step D involves constructing a nighttime radiative transfer equation that includes direct auroral radiation, auroral surface reflected radiation, and moonlight surface reflected radiation. For each pixel in the nighttime remote sensing image, based on the nighttime remote sensing sensor's entrance pupil radiance, the moonlight radiance at the top of the atmosphere, the surface reflectance, and the atmospheric transmittance in all directions, the equation is transformed and solved pixel by pixel to obtain the auroral radiance of each pixel. Generate an image of aurora radiation brightness:

[0013] ;

[0014] In the formula, L is the entrance pupil radiance of the remote sensing sensor at night. M The brightness of lunar radiation at the top of the atmosphere;

[0015] Step E involves sequentially removing non-aurora interference, synthesizing time series data, and optimizing spatial smoothness in the aurora radiance image to obtain a spatially continuous and stable annual-scale aurora radiance spatial distribution result.

[0016] Step A further includes:

[0017] Step A1: Preset atmospheric state parameters and observation geometric parameters. Use aerosol optical thickness as the atmospheric state input variable and the observation zenith angle as the geometric observation condition variable. Select multiple parameter combinations within their respective preset ranges according to a set step size. Calculate the atmospheric transmittance corresponding to each parameter combination using the 6S radiative transfer model. Based on the calculation results and corresponding parameter values, construct a two-dimensional atmospheric transmittance lookup table of aerosol optical thickness and observation zenith angle. The aerosol optical thickness range covers all operating conditions from clean to heavily polluted atmospheres, and the observation zenith angle range covers the entire viewing angle variation range of the large-angle tilt observation of the image edge observed by the VIIRS sensor.

[0018] Step A2: Based on MERRA-2 data, aerosol optical thickness is extracted as the atmospheric state input. Using VIIRS nighttime light remote sensing data to observe the zenith angle, lunar zenith angle, and auroral zenith angle as geometric observation conditions, the atmospheric transmittance τ(θ) in the satellite observation direction is calculated using bilinear interpolation based on the atmospheric transmittance lookup table constructed in Step A1. S Atmospheric transmittance τ(θ) in the direction of moonlight incidenceM Atmospheric transmittance τ(θ) in the aurora radiation direction A ).

[0019] Step B further includes:

[0020] Typical land cover spectra of water bodies, snow, vegetation, and soil were selected from the ECOSTRESS land cover spectral library. The typical reflectance of each type of typical land cover in each band was calculated by combining the spectral response functions of the DNB band and the VIIRS multispectral data in the visible light band. A conversion model from multispectral reflectance to DNB band surface reflectance was established with the typical reflectance of the DNB band as the dependent variable and the typical reflectance of the VIIRS multispectral visible light band as the independent variable.

[0021] The actual observed multispectral visible light band reflectance is read from the VIIRS surface reflectance data, the actual DNB band surface reflectance is calculated using the conversion model, and clear-sky pixels are selected using the quality control band to generate a DNB band surface reflectance time series.

[0022] For missing values ​​in the DNB band surface reflectance time series, an initial time window is selected centered on the time phase corresponding to the missing value for interpolation. If the number of valid observations within the window does not reach a preset threshold, the time window is gradually expanded by a preset step size until the number of valid observations meets the requirement. The median of all valid values ​​within the window is used as the interpolated value to obtain temporally continuous surface reflectance data. .

[0023] Furthermore, in step C, the lunar radiance at the top of the atmosphere... The calculation formula is:

[0024] ;

[0025] In the formula, This represents the brightness of lunar radiation at the top of the atmosphere. The solar constant, The average albedo of the moon, Lunar phase angle, The radius of the moon, This is the distance between the Earth and the Moon.

[0026] Furthermore, in step E, the process of sequentially removing non-auroral interference from the auroral radiance image includes the following steps:

[0027] Based on the VIIRS thermal anomaly product, thermal anomaly flare pixels generated by open flame, natural gas combustion or industrial heat source are extracted from aurora radiation brightness images. A first radius buffer is constructed with the thermal anomaly flare pixels as the center, and the pixels in the aurora radiation brightness image located in the first radius buffer are masked to remove abnormal high brightness signal interference.

[0028] Based on the land cover product, the impermeable surface coverage at a preset resolution is calculated. Pixels with impermeable surface coverage greater than a preset ratio are identified as urban pixels. A second radius buffer is constructed with the urban pixels as the center, and the pixels in the aurora radiation brightness image located in the second radius buffer are masked to remove interference from artificial light signals in urban areas.

[0029] Furthermore, in step E, the process of time-series synthesis of the aurora radiance image includes the following steps:

[0030] For auroral radiance images that have eliminated non-auroral interference, the mean and standard deviation of auroral radiance over time are calculated pixel by pixel. Outliers in the time series are eliminated based on the three-standard-deviation criterion, and only observations that deviate from the mean by no more than three standard deviations are retained. The mean of the effective observations after outlier elimination is synthesized to obtain the annual-scale auroral radiance composite result.

[0031] Furthermore, in step E, for spatially missing pixels formed after removing non-aurora interference, an adaptive neighborhood interpolation method based on Gaussian weights is used:

[0032] Centered on the missing pixel, an initial 3×3 neighborhood window is selected. If the number of valid values ​​within the window is less than a preset threshold, the neighborhood window is gradually expanded with a fixed compensation until the number of valid values ​​within the window meets the requirement. Within the neighborhood window, Gaussian weights are assigned based on the spatial distance between each valid pixel and the center pixel.

[0033] ;

[0034] In the formula, Let be the Gaussian weight of the i-th valid pixel within the neighborhood window. The distance from the i-th valid pixel to the center pixel. This is the Gaussian kernel scaling parameter, which takes the value of 1 / 2 of the neighborhood window radius;

[0035] By utilizing the spatial correlation of neighboring pixels, a locally weighted estimation of the missing values ​​is performed to obtain the aurora brightness interpolation value of the central missing pixel:

[0036] ;

[0037] In the formula, This is the aurora brightness interpolation value. The aurora brightness of the i-th valid pixel within the neighborhood window;

[0038] A 21×21 sliding neighborhood window is constructed with the target pixel as the center. All valid pixel values ​​are extracted within this window. After removing high and low value outliers, the intermediate samples are weighted and averaged using Gaussian weights. The calculation result is used as a smoothed estimate of the target pixel, and finally, the spatially continuous and stable annual-scale aurora brightness spatial distribution result is obtained.

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

[0040] The aurora brightness inversion method based on satellite remote sensing data of this invention can realize large-scale, continuous spatial distribution monitoring of aurora brightness using operational nighttime light remote sensing data. Compared with traditional aurora monitoring methods that rely on ground-based optical observations or numerical simulations, this invention effectively overcomes the limitations of sparse distribution of ground-based observation stations, limited spatial representativeness, and cloud obstruction. It also avoids the problems of high dependence on complex physical parameters and external driving data, and the uncertainty of numerical simulation methods. This invention can more accurately and stably obtain the continuous spatial distribution of large-scale aurora brightness, providing a reliable data foundation for polar space environment research. Its main innovations are: 1) Under a unified radiative transfer framework, the apparent radiance of nighttime light is decomposed into auroral radiation, surface reflection of moonlight, and atmospheric transmission, constructing an auroral radiance inversion model and realizing the effective extraction of auroral radiation signals from nighttime light remote sensing data; 2) A comprehensive background correction method based on nighttime light remote sensing is proposed. By modeling the surface reflectance in the DNB band, combining the moonlight radiance calculated based on lunar phase and observation geometry, and the atmospheric transmittance calculated by the lookup table method, the influence of moonlight radiation, atmospheric effects, and surface reflection characteristics is eliminated. Furthermore, by combining thermal anomaly and land cover data, interference factors such as abnormally bright areas and urban artificial light are eliminated, achieving the elimination of multi-source background influence; 3) By combining annual-scale synthesis and spatial interpolation and smoothing processing, the stability and spatial continuity of auroral brightness results under complex observation conditions are improved, enhancing the usability of the results. The method of this invention has a complete process and stable data sources, with good versatility and scalability, providing a practical technical means for conducting long-term, continuous auroral brightness remote sensing monitoring. Attached Figure Description

[0041] Figure 1 This is a flowchart of the aurora brightness inversion method based on satellite remote sensing data of the present invention;

[0042] Figure 2NPP / VIIRS DNB radiance maps (a. February 26, 2026, b. March 1, 2020, c. March 17, 2020, d. March 22, 2020);

[0043] Figure 3 Aurora radiance maps obtained after removing the effects of atmosphere, moonlight, and surface reflectivity (a. February 26, 2026, b. March 1, 2020, c. March 17, 2020, d. March 22, 2020).

[0044] Figure 4 The aurora radiation brightness map for 2020 was obtained by masking nighttime light anomalies and urban areas and then performing annual composite processing.

[0045] Figure 5 This is the final composite aurora radiation brightness map for 2020. Detailed Implementation

[0046] The embodiments of the present invention will be described in further detail below with reference to the accompanying drawings.

[0047] This invention discloses a method for inverting aurora brightness based on satellite remote sensing data, the method comprising the following steps:

[0048] Step A: Based on aerosol optical thickness and observation geometric parameters, construct an atmospheric transmittance lookup table to retrieve the atmospheric transmittance τ(θ) in the satellite observation direction. S Atmospheric transmittance τ(θ) in the direction of moonlight incidence M Atmospheric transmittance τ(θ) in the direction of aurora radiation A );

[0049] Step B involves constructing a DNB-band surface reflectance inversion model based on multispectral reflectance to obtain time-continuous surface reflectance data. ;

[0050] Step C: Calculate the lunar radiance at the top of the atmosphere based on the lunar phase angle, Earth-Moon distance, and lunar optical parameters. ;

[0051] Step D involves constructing a nighttime radiative transfer equation that includes direct auroral radiation, auroral surface reflected radiation, and moonlight surface reflected radiation. For each pixel in the nighttime remote sensing image, based on the nighttime remote sensing sensor's entrance pupil radiance, the moonlight radiance at the top of the atmosphere, the surface reflectance, and the atmospheric transmittance in all directions, the equation is transformed and solved pixel by pixel to obtain the auroral radiance of each pixel. Generate an image of aurora radiation brightness:

[0052] ;

[0053] In the formula, L is the entrance pupil radiance of the remote sensing sensor at night. M The brightness of lunar radiation at the top of the atmosphere;

[0054] Step E involves sequentially removing non-aurora interference, synthesizing time series data, and optimizing spatial smoothness in the aurora radiance image to obtain a spatially continuous and stable annual-scale aurora radiance spatial distribution result.

[0055] 1) Calculation of atmospheric transmittance

[0056] Nighttime light remote sensing signals are affected by atmospheric scattering and absorption during propagation, leading to attenuation of radiative energy and enhancement of path radiation. Therefore, it is necessary to calculate atmospheric transmittance from different observation directions to effectively correct for atmospheric effects. Using aerosol optical thickness and observation zenith angle as input variables, the atmospheric transmittance under different combinations of conditions is calculated using the 6S radiative transfer model, and an atmospheric transmittance lookup table is constructed. For parameter settings, the aerosol optical thickness ranges from 0 to 2 with a step size of 0.2, covering conditions from clean to heavily polluted atmospheres; the observation zenith angle ranges from 0° to 70° with a step size of 10°, reflecting the angular change from the VIIRS sensor zenith angle to the large-angle tilt observation at the image edge. For each parameter combination, the corresponding atmospheric transmittance is calculated, and a two-dimensional atmospheric transmittance lookup table is constructed based on the calculation results and corresponding parameter values.

[0057] Using aerosol optical thickness extracted from MERRA-2 data as atmospheric state input, and combining zenith angles observed from VIIRS nighttime light remote sensing data (lunar zenith angle and auroral zenith angle approximated as 0° considering that auroral radiation mainly comes from above) as observation geometric parameters, the atmospheric transmittance τ(θ) in the satellite observation direction was calculated using bilinear interpolation based on the constructed atmospheric transmittance lookup table. S Atmospheric transmittance τ(θ) in the direction of moonlight incidence M Atmospheric transmittance τ(θ) in the direction of aurora radiation A ).

[0058] 2) Calculation of surface reflectance

[0059] The reflection of moonlight and auroras by the Earth's surface is a crucial factor affecting nighttime light observation; therefore, it is necessary to calculate the surface reflectance of the VIIRS nighttime light remote sensing band (DNB band). Typical spectra of water bodies, snow, vegetation, and soil were selected from the ECOSTRESS ground object spectral library. The reflectance of each band was calculated by combining the spectral response functions of the DNB band and the visible light bands (M4, M5, and M7 bands) from VIIRS multispectral data. Using the DNB band reflectance as the dependent variable and the M4, M5, and M7 band reflectances as independent variables, a multispectral reflectance to DNB reflectance conversion model was established.

[0060] ;

[0061] In the formula, ρ DNB ρ represents the surface reflectance in the DNB band. M4 ρ M5 ρ M7 These are the surface reflectance values ​​for the M4, M5, and M7 bands, respectively.

[0062] The reflectance of the M4, M5 and M7 bands was read from the VIIRS surface reflectance data. The corresponding DNB band reflectance was calculated using the above conversion model. Quality control bands were used for quality screening, and only clear sky pixels were retained.

[0063] To address missing values ​​in the reflectance time series, interpolation was performed using 5 time phases before and after the current time phase (a total of 11 time phases). When the number of valid observations within a window was less than 3, the time window range was gradually expanded to ±9, ±13, ±17, etc., until the number of valid observations within the window met the requirements. The median of all valid values ​​within the time window was calculated as the interpolated value, thus obtaining temporally continuous surface reflectance data.

[0064] 3) Calculation of lunar radiance at the top of the atmosphere

[0065] Moonlight is the primary natural light source at night, and its intensity is significantly affected by factors such as lunar phase, observation angle, and Earth-Moon distance. Under near-full moon conditions, the moonlight radiation signal is very strong, while under crescent moon conditions, the signal is very weak. These variations in moonlight radiation interfere with the quantitative inversion of aurora radiation. Therefore, quantitative calculations of moonlight radiation are necessary to accurately eliminate its influence during aurora inversion. The formula for calculating the brightness of moonlight radiation in the upper atmosphere is:

[0066] ;

[0067] In the formula, L M E represents the brightness of lunar radiation at the top of the atmosphere. sun A is the solar constant. M R is the average lunar albedo, α is the lunar phase angle, and R is the lunar phase angle. M Let d be the radius of the moon. EM This is the distance between the Earth and the Moon.

[0068] 4) Aurora Radiance Inversion

[0069] The radiation signals received by remote sensing sensors at night mainly consist of aurora up-reaching radiation, aurora down-reaching radiation reflected from the Earth's surface, and moonlight reflected from the Earth's surface. Their functional relationship can be expressed as:

[0070] ;

[0071] In the formula, L is the entrance pupil radiance of the remote sensing sensor at night, i.e., the apparent radiance; τ(θ) A ), τ(θ) S ), τ(θ) M ρ represents the atmospheric transmittance at the corresponding aurora zenith angle, satellite zenith angle, and lunar zenith angle, respectively. DNB L represents the surface reflectance in the DNB band. A L represents the aurora radiation brightness. M This represents the brightness of lunar radiation at the top of the atmosphere.

[0072] By transforming the radiative transfer equation above, we can obtain the aurora radiance L. A The expression function:

[0073] ;

[0074] The apparent radiance of the nighttime light was read at 500m resolution from the VNP46A1 nighttime light radiance product, and combined with the previously calculated surface reflectance ρ. DNB Lunar radiance at the top of the atmosphere (L) M Atmospheric transmittance τ(θ) in the direction of satellite observation S Atmospheric transmittance τ(θ) in the direction of moonlight incidence M Atmospheric transmittance τ(θ) in the direction of aurora radiation A ), calculate the aurora radiance at 500m resolution pixel by pixel.

[0075] 5) Removal of abnormally bright signals

[0076] Nighttime forest fires, natural gas combustion, and industrial heat sources generate abnormally bright signals, whose radiation characteristics differ significantly from aurora signals, necessitating identification and removal. Based on the VIIRS thermal anomaly product VNP14A1, flare pixels caused by fires and natural gas combustion were extracted. Considering the spillover effect of high-brightness light, a buffer zone with a radius of 5 km was constructed around these flare pixels. Subsequently, pixel values ​​within the buffer zone in the aurora radiance calculation results were masked to eliminate interference from abnormally bright signals.

[0077] 6) Removal of artificial light signals in urban areas

[0078] Artificial lighting in urban areas (such as road lighting and building lighting) has a stable spatial distribution and high brightness, which can easily have a significant impact on weak radiation signals such as auroras, so it is necessary to eliminate artificial light signals in urban areas.

[0079] Impermeable surface cover at 500 m resolution was calculated using Dynamic World land cover products, and urban pixels were extracted based on an impermeable surface cover greater than 10%. Considering the spillover effect of urban light, a buffer zone with a radius of 10 km was constructed for urban pixels. Then, pixel values ​​within the buffer zone in the aurora radiance calculation results were masked to eliminate interference from artificial light signals from urban areas.

[0080] 7) Annual-scale synthesis

[0081] Due to the significant temporal instability of auroral radiation signals, and the influence of various factors such as cloud cover, moonlight interference, atmospheric conditions, and sensor noise, daily inversion results often exhibit large random fluctuations and uncertainties, making it difficult to stably characterize the true distribution features of regional auroral radiation. Therefore, annual-scale composite processing is performed based on the daily auroral radiance obtained from the previous calculations and processing.

[0082] The mean and standard deviation of aurora brightness over time are calculated pixel by pixel. Outliers in the time series are removed using the three-standard-deviation (3σ) criterion, meaning only observations that deviate from the mean by no more than three standard deviations are retained. After removing outliers, the effective values ​​of each pixel are statistically analyzed and averaged to synthesize the annual aurora radiance.

[0083] 8) Spatial interpolation and smoothing

[0084] After obtaining the annual aurora radiance composite results, missing value imputation and spatial smoothing were performed to further improve the spatial continuity and stability of the data.

[0085] To address spatially missing pixels in urban areas and abnormally bright regions after masking, an adaptive neighborhood interpolation method based on Gaussian weights is employed. Centered on the missing pixel, valid pixels are initially selected within a 3×3 neighborhood for interpolation. If fewer than 5 valid values ​​are found within the window, the neighborhood window is progressively expanded to 5×5, 7×7, 9×9, and so on, until at least 5 valid values ​​are found. Within the neighborhood window, Gaussian weights are assigned based on the spatial distance between each valid pixel and the center pixel. The spatial correlation of neighboring pixels is then used to locally weight and estimate the missing value.

[0086] ;

[0087] In the formula, w is the aurora brightness interpolation value. i L is the Gaussian weight of the i-th valid pixel within the neighborhood window. M,i Let be the aurora brightness of the i-th valid pixel within the neighborhood window.

[0088] Gaussian weight wi The calculation formula is:

[0089] ;

[0090] In the formula, wi is the Gaussian weight of the i-th valid pixel within the neighborhood window, and d i σ is the distance from the i-th valid pixel to the center pixel, and σ is the Gaussian kernel scaling parameter, which is half of the window radius to ensure that the weights decay smoothly with distance.

[0091] Finally, to reduce spatial anomaly fluctuations caused by noise, local anomalous radiation, and inversion errors, spatial smoothing was performed on the interpolated aurora radiation data. A 21×21 sliding neighborhood window was constructed with the target pixel as the center. All valid pixel values ​​were extracted within this window, and the middle 50% of the samples were selected according to the numerical size (i.e., high and low value anomalies were removed). Based on this, a weighted average was calculated using Gaussian weights to serve as the smoothed estimate of the target pixel, ultimately yielding a spatially continuous and stable annual-scale aurora brightness spatial distribution result.

[0092] Example

[0093] This invention retrieves aurora brightness in a high-latitude region (longitude -120°~-110°, latitude 55°~65°) as the study area. The following are the specific implementation steps of an example. See the technical flow... Figure 1 .

[0094] 1) Obtain the 2020 VIIRS daily apparent radiance product VNP46A1, VIIRS daily surface reflectance product VNP09GA, VIIRS daily thermal anomaly product VNP14A1, MERRA-2 daily reanalysis aerosol optical thickness product M2T1NXAER, and Dynamic World annual land cover product for the study area.

[0095] 2) Using aerosol optical thickness and observed zenith angle as input variables, atmospheric transmittance was calculated using the 6S radiative transfer model, and a two-dimensional lookup table of atmospheric transmittance was constructed. Aerosol optical thickness was extracted from MERRA-2 data, and satellite-observed zenith angle, lunar zenith angle, and auroral zenith angle (approximately 0°) were extracted from VNP46A1 data. Based on the constructed atmospheric transmittance lookup table, bilinear interpolation was used to calculate the atmospheric transmittance τ(θ) at the corresponding satellite-observed zenith angle, lunar zenith angle, and auroral zenith angle. S ), τ(θ) M ), τ(θ) A ).

[0096] 3) Extract the surface reflectance of the M4, M5, and M7 bands from the VNP09GA data, substitute them into Formula 1 to calculate the surface reflectance of the DNB band, and then perform masking processing based on the cloud mask extracted from the QF1 band in the VNP09GA product. For missing reflectance values ​​after masking, a time neighborhood adaptive interpolation method is used. Centered on the current time phase, a time window is first constructed by selecting 5 time phases before and after the current time phase; when there are fewer than 3 valid observations in the window, the time window range is gradually expanded until there are no fewer than 3 valid observations. The median of the valid values ​​in the window is calculated as the interpolation value to generate surface reflectance data without missing values.

[0097] 4) Extract the lunar phase angle from the VNP46A1 data, and calculate the Earth-Moon distance d based on the astronomical ephemeris corresponding to the observation date. EM The lunar radiance L at the top of the atmosphere was calculated using Formula 2. M .

[0098] 5) Extract the apparent radiance of the DNB band at a resolution of 500m from the VNP46A1 data, and perform masking processing using the cloud mask extracted from the QF_Cloud_Mask band to obtain the apparent radiance of the nighttime under clear skies. Figure 2 A spatial distribution map of DNB radiance in the study area on a given date is provided.

[0099] 6) Combining the atmospheric transmittance corresponding to the previously calculated apparent radiance of nighttime light, surface reflectance, lunar radiance at the top of the atmosphere, satellite zenith angle, lunar zenith angle, and aurora zenith angle, the aurora radiance L is calculated pixel by pixel using Formula 4. A . Figure 3 A spatial distribution map of aurora radiance in the study area on *month*day, after removing atmospheric effects, moonlight, and surface reflection, is presented.

[0100] 7) Extract thermal anomaly pixels from VNP14A1 data, construct a buffer with a radius of 5km based on these pixels, and perform masking on the aurora radiance pixel values ​​within the buffer.

[0101] 8) Extract impermeable surface pixels at 10m resolution from Dynamic World land cover data, calculate the impermeable surface coverage at 500m resolution, and use 10% as a threshold to identify urban pixels. Construct a 10km buffer zone based on urban pixels, and perform masking on the aurora radiance within the buffer zone.

[0102] 9) Calculate the mean and standard deviation of the aurora brightness time series pixel by pixel, and remove outliers based on the 3σ criterion. Calculate the average value of the remaining valid observations as a composite value to obtain the annual aurora radiance. Figure 4 The aurora radiance obtained by synthesizing data on *month*day*year in the study area is given.

[0103] 10) To address the missing values ​​in the annual synthetic aurora radiation brightness caused by the removal of abnormally bright areas, towns, and water bodies, an adaptive neighborhood interpolation method based on Gaussian weights is adopted. The weights are calculated according to Formula 6 based on the distance between the effective value pixels in the neighborhood window and the center pixel of the window. Then, Formula 5 is used to locally weight the missing values ​​to obtain the interpolated aurora radiation brightness.

[0104] 11) Construct a 21×21 sliding window centered on each pixel, extract the effective value pixels within the window and sort them by value, retain the middle 50% of the samples, and combine them with Gaussian weights to perform a weighted average as the smoothed estimate of the current pixel, and smooth the aurora radiance. Figure 5 The spatial distribution map of the final annual synthetic aurora radiance is given.

[0105] Although preferred embodiments of this application have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments as well as all changes and modifications falling within the scope of this application.

[0106] Obviously, those skilled in the art can make various modifications and variations to this application without departing from the spirit and scope of this application. Therefore, if such modifications and variations fall within the scope of the claims of this application and their equivalents, this application also intends to include such modifications and variations.

Claims

1. A method for inverting aurora brightness based on satellite remote sensing data, characterized in that, The method includes the following steps: Step A: Based on aerosol optical thickness and observation geometric parameters, construct an atmospheric transmittance lookup table to retrieve the atmospheric transmittance τ(θ) in the satellite observation direction. S Atmospheric transmittance τ(θ) at the direction of moonlight incidence M Atmospheric transmittance τ(θ) in the aurora radiation direction A ); Step B involves constructing a DNB-band surface reflectance inversion model based on multispectral reflectance to obtain time-continuous surface reflectance data. ; Step C: Calculate the lunar radiance at the top of the atmosphere based on the lunar phase angle, Earth-Moon distance, and lunar optical parameters. ; Step D involves constructing a nighttime radiative transfer equation that includes direct auroral radiation, auroral surface reflected radiation, and moonlight surface reflected radiation. For each pixel in the nighttime remote sensing image, based on the nighttime remote sensing sensor's entrance pupil radiance, the moonlight radiance at the top of the atmosphere, the surface reflectance, and the atmospheric transmittance in all directions, the equation is transformed and solved pixel by pixel to obtain the auroral radiance of each pixel. Generate an image of aurora radiation brightness: ； In the formula, L is the entrance pupil radiance of the remote sensing sensor at night. M The brightness of lunar radiation at the top of the atmosphere; Step E involves sequentially removing non-aurora interference, synthesizing time series data, and optimizing spatial smoothness in the aurora radiance image to obtain a spatially continuous and stable annual-scale aurora radiance spatial distribution result.

2. The aurora brightness inversion method based on satellite remote sensing data according to claim 1, characterized in that, Step A further includes: Step A1: Preset atmospheric state parameters and observation geometric parameters. Use aerosol optical thickness as the atmospheric state input variable and the observation zenith angle as the geometric observation condition variable. Select multiple parameter combinations within their respective preset ranges according to a set step size. Calculate the atmospheric transmittance corresponding to each parameter combination using the 6S radiative transfer model. Based on the calculation results and corresponding parameter values, construct a two-dimensional atmospheric transmittance lookup table of aerosol optical thickness and observation zenith angle. The aerosol optical thickness range covers all operating conditions from clean to heavily polluted atmospheres, and the observation zenith angle range covers the entire viewing angle variation range of the large-angle tilt observation of the image edge observed by the VIIRS sensor. Step A2: Based on MERRA-2 data, aerosol optical thickness is extracted as the atmospheric state input. Using VIIRS nighttime light remote sensing data to observe the zenith angle, lunar zenith angle, and auroral zenith angle as geometric observation conditions, the atmospheric transmittance τ(θ) in the satellite observation direction is calculated using bilinear interpolation based on the atmospheric transmittance lookup table constructed in Step A1. S Atmospheric transmittance τ(θ) at the direction of moonlight incidence M Atmospheric transmittance τ(θ) in the aurora radiation direction A ).

3. The aurora brightness inversion method based on satellite remote sensing data according to claim 1, characterized in that, Step B further includes: Typical land cover spectra of water bodies, snow, vegetation, and soil were selected from the ECOSTRESS land cover spectral library. The typical reflectance of each type of typical land cover in each band was calculated by combining the spectral response functions of the DNB band and the VIIRS multispectral data in the visible light band. A conversion model from multispectral reflectance to DNB band surface reflectance was established with the typical reflectance of the DNB band as the dependent variable and the typical reflectance of the VIIRS multispectral visible light band as the independent variable. The actual observed multispectral visible light band reflectance is read from the VIIRS surface reflectance data, the actual DNB band surface reflectance is calculated using the conversion model, and clear-sky pixels are selected using the quality control band to generate a DNB band surface reflectance time series. For missing values ​​in the DNB band surface reflectance time series, an initial time window is selected centered on the time phase corresponding to the missing value for interpolation. If the number of valid observations within the window does not reach a preset threshold, the time window is gradually expanded by a preset step size until the number of valid observations meets the requirement. The median of all valid values ​​within the window is used as the interpolated value to obtain temporally continuous surface reflectance data. .

4. The aurora brightness inversion method based on satellite remote sensing data according to claim 1, characterized in that, In step C, the lunar radiance at the top of the atmosphere... The calculation formula is: ； In the formula, This represents the brightness of lunar radiation at the top of the atmosphere. The solar constant, The average albedo of the moon, Lunar phase angle, The radius of the moon, This is the distance between the Earth and the Moon.

5. The aurora brightness inversion method based on satellite remote sensing data according to claim 1, characterized in that, Step E, the process of sequentially removing non-aurora interference from the aurora radiance image, includes the following steps: Based on the VIIRS thermal anomaly product, thermal anomaly flare pixels generated by open flame, natural gas combustion or industrial heat source are extracted from aurora radiation brightness images. A first radius buffer is constructed with the thermal anomaly flare pixels as the center, and the pixels in the aurora radiation brightness image located in the first radius buffer are masked to remove abnormal high brightness signal interference. Based on the land cover product, the impermeable surface coverage at a preset resolution is calculated. Pixels with impermeable surface coverage greater than a preset ratio are identified as urban pixels. A second radius buffer is constructed with the urban pixels as the center, and the pixels in the aurora radiation brightness image located in the second radius buffer are masked to remove interference from artificial light signals in urban areas.

6. The aurora brightness inversion method based on satellite remote sensing data according to claim 1, characterized in that, Step E, the process of time-series synthesis of aurora radiance images, includes the following steps: For auroral radiance images that have eliminated non-auroral interference, the mean and standard deviation of auroral radiance over time are calculated pixel by pixel. Outliers in the time series are eliminated based on the three-standard-deviation criterion, and only observations that deviate from the mean by no more than three standard deviations are retained. The mean of the effective observations after outlier elimination is synthesized to obtain the annual-scale auroral radiance composite result.

7. The aurora brightness inversion method based on satellite remote sensing data according to claim 1, characterized in that, In step E, for spatially missing pixels formed after removing non-aurora interference, an adaptive neighborhood interpolation method based on Gaussian weights is used: Centered on the missing pixel, an initial 3×3 neighborhood window is selected. If the number of valid values ​​within the window is less than a preset threshold, the neighborhood window is gradually expanded with a fixed compensation until the number of valid values ​​within the window meets the requirement. Within the neighborhood window, Gaussian weights are assigned based on the spatial distance between each valid pixel and the center pixel. ； In the formula, Let be the Gaussian weight of the i-th valid pixel within the neighborhood window. The distance from the i-th valid pixel to the center pixel. This is the Gaussian kernel scaling parameter, which takes the value of 1 / 2 of the neighborhood window radius; By utilizing the spatial correlation of neighboring pixels, a locally weighted estimation of the missing values ​​is performed to obtain the aurora brightness interpolation value of the central missing pixel: ； In the formula, This is the aurora brightness interpolation value. The aurora brightness of the i-th valid pixel within the neighborhood window; A 21×21 sliding neighborhood window is constructed with the target pixel as the center. All valid pixel values ​​are extracted within this window. After removing high and low value outliers, the intermediate samples are weighted and averaged using Gaussian weights. The calculation result is used as a smoothed estimate of the target pixel, and finally, the spatially continuous and stable annual-scale aurora brightness spatial distribution result is obtained.

Images