A method and system for constructing an atmospheric model based on virtual temperature correction

Abstract

The application relates to the technical field of atmospheric parameter inversion model, and provides a method and system for constructing an atmospheric model based on virtual temperature correction, which solves the problem of insufficient adaptability of a dynamic scene. The application obtains a radiation brightness observation value of a satellite thermal infrared band, synchronously records a satellite observation zenith angle and an azimuth angle, and collects surface coverage type data, a dew point temperature and corresponding height information. Upward and downward radiation components are eliminated. A normalized vegetation coverage index is calculated based on satellite red light and near-infrared band reflectivity data, a surface coverage type is combined to construct a specific radiation rate mapping relationship, and a surface specific radiation rate spatial distribution is inverted. Vertical layered atmospheric data are generated based on the dew point temperature and the height, and a radiation transmission path correction parameter is constructed in combination with the observation angle. Finally, an atmospheric radiation transmission model is established, and a continuous spatialized virtual temperature field is output. The application solves the problem of insufficient adaptability of a dynamic scene in the prior art, and realizes high-precision inversion of a virtual temperature vertical profile of the atmosphere.

Description

Technical Field

[0001] The present application relates to the technical field of atmospheric parameter inversion models, and in particular to a method and system for constructing an atmospheric model based on virtual temperature correction. Background Art

[0002] The coordinated inversion of the large-scale surface thermal environment and the atmospheric virtual temperature field requires addressing the core issue of wide-area surface temperature monitoring and simultaneous analysis of vertical atmospheric thermal states. This technology must cover hundreds of square kilometers with a spatial resolution of 100 meters, while simultaneously modeling the surface thermal radiation characteristics and the distribution of virtual temperatures at different atmospheric levels.

[0003] Currently, the mainstream solution is the multi-channel split-window algorithm. This approach relies on radiance data from two adjacent bands of a satellite's thermal infrared sensor, using preset empirical coefficients to separate atmospheric water vapor interference to directly estimate surface temperature. The core process utilizes the statistical characteristics of atmospheric absorption differences between the two bands to construct a linear combination equation to convert radiance values ​​into temperature values.

[0004] However, existing methods have some significant shortcomings. First, the calibration of empirical coefficients relies on historical atmospheric conditions or standard models, making it difficult to dynamically respond to sudden atmospheric disturbances, resulting in systematic deviations in temperature inversion results. Second, this solution only outputs surface temperature at the surface level and cannot invert the virtual temperature distribution at different atmospheric altitudes. This results in a missing vertical dimension for thermal environment analysis and fails to meet the requirements of meteorological models for thermal field structures. Summary of the Invention

[0005] The present application provides a method and system for constructing an atmospheric model based on virtual temperature correction, which is used to solve the problem of insufficient adaptability to dynamic scenes in the existing technology, and realizes the inversion of the atmospheric virtual temperature vertical profile with high precision.

[0006] In a first aspect, the present application provides a method for constructing an atmospheric model based on virtual temperature correction, comprising:

[0007] Obtain the satellite's thermal infrared band radiation brightness observation values, and simultaneously record the satellite's observation zenith angle and azimuth angle data, while obtaining the surface cover type classification data, dew point temperature and corresponding height information;

[0008] Performing atmospheric path radiation distortion correction on the radiance observation value using a radiation calibration coefficient, and generating a calibrated radiance observation value by eliminating the atmospheric up-going radiation and down-going reflected radiation components;

[0009] Based on the reflectance data of the satellite red light band and the near-infrared band, a normalized vegetation cover index is calculated. In combination with the surface cover type classification data, a mapping relationship between the surface cover type and the emissivity is constructed. The nonlinear relationship between the vegetation cover index and the emissivity is fitted through the mapping relationship to invert the spatial distribution data of the surface emissivity.

[0010] Based on the dew point temperature and the corresponding height information, vertically layered atmospheric profile data is generated, the spatial distribution data is integrated with the atmospheric profile data, and correction parameters of the atmospheric radiation transmission path are constructed in combination with the angle data;

[0011] An atmospheric radiation transfer model is established based on the correction parameters, and the calibrated radiation brightness observation value is input into the atmospheric radiation transfer model. The vertical profile of the atmospheric virtual temperature is inverted by solving the radiation transfer equation, and the spatial distribution data of the atmospheric virtual temperature field is output.

[0012] Optionally, based on the reflectance data of the satellite red light band and the near-infrared band, a normalized vegetation cover index is calculated, and combined with the surface cover type classification data, a mapping relationship between the surface cover type and the emissivity is constructed, and the nonlinear relationship between the vegetation cover index and the emissivity is fitted through the mapping relationship to invert the spatial distribution data of the surface emissivity, including:

[0013] Obtaining the reflectance data of the satellite red light band and the near infrared band, for each pixel, substituting the values ​​of the reflectance data of the red light band and the near infrared band into the normalized difference calculation formula to generate the normalized vegetation cover index of each pixel;

[0014] Combined with the surface cover type classification data, the surface cover types are divided into four categories: water bodies, vegetation, bare soil, and artificial buildings. A mapping function between the vegetation cover index and the emissivity is established for each surface cover type, where the vegetation type is fitted with an exponential function and the non-vegetation type is fitted with a piecewise linear function.

[0015] Based on the mapping function, the vegetation cover index is input into the fitting function of the corresponding surface cover type, the emissivity value of each pixel is output, and the emissivity values ​​of all pixels are integrated to generate continuous spatial distribution data of surface emissivity.

[0016] Optionally, based on the dew point temperature and the corresponding height information, vertically layered atmospheric profile data is generated, the spatial distribution data is fused with the atmospheric profile data, and the correction parameters of the atmospheric radiation transmission path are constructed in combination with the angle data, including:

[0017] Based on the dew point temperature and the corresponding height information, generate atmospheric profile data by vertical height layer, wherein each layer includes temperature, water vapor density, and air pressure parameters;

[0018] Matching the spatial distribution data of the surface emissivity with the atmospheric profile data according to geographical location and fusing them into a three-dimensional gridded data set;

[0019] Calculating the slant path length of electromagnetic waves passing through each atmospheric layer based on the atmospheric profile parameters in the three-dimensional gridded dataset and the angular data of the zenith angle and azimuth angle observed by the satellite;

[0020] Calculating the atmospheric attenuation coefficient and scattering weight layer by layer based on the oblique path length and the atmospheric parameters of each layer in the three-dimensional gridded data set;

[0021] The attenuation coefficients and scattering weights of all layers are integrated to generate correction parameters of the atmospheric radiation transmission path that vary with altitude and angle.

[0022] Optionally, the vertical profile of the atmospheric virtual temperature is inverted by solving the radiation transfer equation, and the spatial distribution data of the atmospheric virtual temperature field is output, including:

[0023] According to the calibrated radiance observation value, the dry delay component data and the wet delay component data are solved by solving the radiative transfer equation using the atmospheric radiative transfer model;

[0024] Based on the dry delay component and wet delay component data, combined with the atmospheric correction parameter table, calculate the contribution degree data of each atmospheric layer to the dry delay component and the wet delay component;

[0025] Based on the contribution data, the virtual temperature value of each atmospheric layer is adjusted by iteratively solving the radiation transfer equation until the error between the calculated radiance value output by the atmospheric radiation transfer model and the calibrated radiance observation value is minimized, thereby generating an optimized vertical virtual temperature sequence;

[0026] According to the optimized vertical virtual temperature sequence, atmospheric virtual temperature vertical profile data of a single geographical point is formed, and the atmospheric virtual temperature vertical profile data of all geographical points are integrated to obtain a continuous atmospheric virtual temperature field spatial distribution map.

[0027] Optionally, according to the calibrated radiance observation value, solving the radiative transfer equation using the atmospheric radiative transfer model to obtain dry delay component data and wet delay component data includes:

[0028] Establishing a radiation transfer relationship to express the calibrated radiance observation value as a linear superposition of a dry delay component and a wet delay component, wherein the dry delay component corresponds to the absorption effect of dry atmospheric components and the wet delay component corresponds to the absorption effect of water vapor components;

[0029] Performing spectral decomposition on the calibrated radiance observation value, separating independent spectral response segments for different absorption characteristic bands, and obtaining radiance decomposition values ​​corresponding to each spectral response segment;

[0030] Based on the radiation transfer relationship, a spectral response segment dominated by dry atmospheric components is selected, and the corresponding radiation brightness decomposition value is substituted into the independent expression of the dry delay component to generate dry delay component data;

[0031] Based on the radiation transfer relationship, a spectral response segment dominated by water vapor components is selected, and the corresponding radiation brightness decomposition value is substituted into the independent expression of the wet delay component to generate wet delay component data.

[0032] Optionally, establishing an atmospheric radiation transfer model according to the correction parameters and inputting the calibrated radiance observation value into the atmospheric radiation transfer model includes:

[0033] Constructing a vertically layered atmospheric radiation transfer model framework, dividing the atmosphere into multiple continuous layers according to altitude, with each layer corresponding to the attenuation coefficient and scattering weight in the correction parameters;

[0034] In the model framework, the angle data is converted into the inclined path length of the electromagnetic wave propagation path in each layer, and the cumulative impact factor of the path is generated by combining the attenuation coefficient and scattering weight of each layer;

[0035] The calibrated radiation brightness observation value is converted into an equivalent radiation energy input value, and the equivalent radiation energy input value is input into an atmospheric radiation transfer model containing the path cumulative impact factor.

[0036] Optionally, the radiometric brightness observation value is corrected for atmospheric path radiation distortion using a radiometric calibration coefficient, and a calibrated radiometric brightness observation value is generated by eliminating the atmospheric upgoing radiation and downgoing reflected radiation components, including:

[0037] The radiance observation value is converted into an equivalent radiant energy value based on a radiation calibration coefficient, and the equivalent radiant energy value is separated into an atmospheric upward radiation component and a downward radiation component reflected by the surface, wherein the upward radiation component is calculated by a top-of-atmosphere radiation model, and the downward radiation component is determined by multiplying the surface reflectivity and the atmospheric downward radiation;

[0038] The superposition value of the upgoing radiation component and the downgoing radiation component is subtracted from the equivalent radiation energy value to generate a calibrated radiation brightness observation value representing the target surface radiation.

[0039] In a second aspect, the present application provides an atmospheric model construction system based on virtual temperature correction, comprising:

[0040] The acquisition module is used to obtain the radiance observation value of the satellite thermal infrared band, and synchronously record the angle data of the satellite observation zenith angle and azimuth angle, and obtain the surface cover type classification data, dew point temperature and corresponding height information;

[0041] a calibration module for performing atmospheric path radiation distortion correction on the radiance observation value using a radiation calibration coefficient, and generating a calibrated radiance observation value by eliminating the atmospheric up-going radiation and down-going reflected radiation components;

[0042] An inversion module is used to calculate the normalized vegetation cover index based on the reflectance data of the satellite red light band and the near-infrared band, and to construct a mapping relationship between the surface cover type and the emissivity in combination with the surface cover type classification data. The nonlinear relationship between the vegetation cover index and the emissivity is fitted through the mapping relationship to invert the spatial distribution data of the surface emissivity.

[0043] A construction module is used to generate vertically layered atmospheric profile data based on the dew point temperature and the corresponding height information, fuse the spatial distribution data with the atmospheric profile data, and construct correction parameters of the atmospheric radiation transmission path in combination with the angle data;

[0044] The output module is used to establish an atmospheric radiation transfer model based on the correction parameters, input the calibrated radiation brightness observation value into the atmospheric radiation transfer model, invert the atmospheric virtual temperature vertical profile by solving the radiation transfer equation, and output the atmospheric virtual temperature field spatial distribution data.

[0045] In a third aspect, the present application provides a computing device comprising a processing component and a storage component; the storage component stores one or more computer instructions; the one or more computer instructions are used to be called and executed by the processing component to implement a method for constructing an atmospheric model based on virtual temperature correction as described in the first aspect above.

[0046] In a fourth aspect, the present application provides a computer storage medium storing a computer program. When the computer program is executed by a computer, it implements a method for constructing an atmospheric model based on virtual temperature correction as described in the first aspect.

[0047] This application obtains the radiation brightness observation values ​​and angle data of the satellite thermal infrared band, and simultaneously collects the surface type, dew point temperature and altitude information, laying a multi-source data foundation for collaborative inversion; uses the radiation calibration coefficient to eliminate the upward and downward radiation distortion of the atmosphere to generate calibrated radiation brightness data, thereby improving the physical credibility of the observation value; calculates the normalized vegetation index based on the reflectivity of the red and near-infrared bands, constructs the emissivity mapping relationship in combination with the surface type and fits the nonlinear function, thereby realizing the accurate quantification of complex surface thermal radiation characteristics; generates layered atmospheric profile data based on dew point temperature and altitude, integrates the spatial distribution of surface emissivity and constructs the radiation transfer path correction parameters in combination with the observation angle, thereby significantly improving the accuracy of atmospheric thermal radiation path modeling; establishes a radiation transfer model based on the correction parameters and inputs the calibration data, and inverts the vertical profile of the atmospheric virtual temperature by solving the radiation equation. Finally, the three-dimensional spatialized virtual temperature field distribution is output, realizing the full-dimensional characterization of the vertical thermal structure from the surface to the atmosphere.

[0048] Furthermore, for each pixel, the red and near-infrared band reflectances are substituted into the normalized difference formula to generate a vegetation cover index, and the vegetation density distribution is quantified based on spectral differences; the surface types are divided into four categories: water bodies, vegetation, bare soil, and artificial buildings. The vegetation type is fitted with an exponential function, and the non-vegetation type is designed with a piecewise linear function to design differentiated mapping relationships, so as to solve the adaptation problem of the thermal radiation mechanism of different land types in a targeted manner; the vegetation index is input into the fitting function of the corresponding surface type to output the single pixel emissivity, and all pixel data are integrated to generate a continuous spatial distribution map, which ensures the spatial consistency of thermal radiation parameters on complex underlying surfaces.

[0049] These and other aspects of the present application will become more readily apparent from the description of the following embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, a brief introduction will be given below to the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without any creative work.

[0051] Figure 1 A flowchart of a method for constructing an atmospheric model based on virtual temperature correction provided by the present application is shown;

[0052] Figure 2 A scene diagram showing a method for constructing an atmospheric model based on virtual temperature correction provided by the present application is shown;

[0053] Figure 3 A schematic diagram of the structure of an atmospheric model construction system based on virtual temperature correction provided by the present application is shown;

[0054] Figure 4 A schematic structural diagram of a computing device provided by the present application is shown. DETAILED DESCRIPTION

[0055] In order to enable those skilled in the art to better understand the solution of the present application, the technical solution in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application.

[0056] In some of the processes described in the specification and claims of this application and the above-mentioned figures, multiple operations that appear in a specific order are included, but it should be clearly understood that these operations may not be executed in the order in which they appear in this document or may be executed in parallel. The serial numbers of the operations, such as 101, 102, etc., are only used to distinguish between different operations, and the serial numbers themselves do not represent any order of execution. In addition, these processes may include more or fewer operations, and these operations may be executed in sequence or in parallel. It should be noted that the descriptions of "first", "second", etc. in this document are used to distinguish different messages, devices, modules, etc., and do not represent a sequential order, nor do they limit "first" and "second" to being different types.

[0057] Research has found that traditional atmospheric temperature detection methods have significant flaws: First, their empirical calculation parameters rely heavily on historical weather or theoretical models, making them difficult to dynamically adjust in the face of sudden meteorological changes such as rainstorms and dust storms, resulting in large deviations between satellite-measured surface temperatures and actual temperatures. Second, existing solutions can only output surface temperature results for the surface layer, failing to capture the three-dimensional temperature distribution from the ground to high altitudes. This results in a lack of critical thermal field structure data for weather forecasts and disaster warnings. Therefore, a new modeling method that can accurately invert the vertical temperature distribution of the atmosphere is urgently needed.

[0058] In response to the above problems, the present invention proposes a virtual temperature correction modeling method that dynamically integrates multi-source data. The core of the method is to use the collaborative analysis of satellite observation data and meteorological parameters to construct an atmospheric thermal radiation propagation path with self-corrected errors. Specifically, the infrared heat signal received by the satellite is first combined with the solar angle information to eliminate the measurement interference caused by atmospheric scattering; then, the emissivity characteristics are automatically matched according to the vegetation coverage status and surface type; then, the water vapor temperature information at different altitudes in the atmosphere is integrated to generate a "vertical slice of the atmosphere" data chain; finally, based on the corrected thermal radiation propagation formula, a layered temperature distribution map from the surface to high altitude is inverted. This method successfully breaks through traditional limitations: by replacing empirical coefficients with real-time data drive, the temperature inversion accuracy under sudden meteorological disturbances is improved; at the same time, the innovative output of the atmospheric vertical virtual temperature field fills the gap in the demand for three-dimensional thermal structure in meteorological models.

[0059] This application introduces a virtual temperature correction method to improve the simulated atmospheric environment's temperature, pressure, atmospheric density, and other environmental parameters. This aims to address the inability of standard atmospheric models to meet the requirements of a realistic training environment under actual meteorological conditions (such as temperature deviation, humidity, wind speed, and terrain effects). This also addresses the issue of incomplete aircraft performance and aerodynamic effects during simulator training under complex meteorological conditions. By introducing virtual temperature (wet air temperature) and humidity to correct the standard atmospheric model, simulator training scenarios become more realistic, providing civil aviation pilots with a more realistic flight training environment.

[0060] The following will be combined with the drawings in the embodiments of this application to clearly and completely describe the technical solutions in the embodiments of this application. Obviously, the embodiments described are only part of the embodiments of this application, not all of the embodiments. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without making creative efforts are within the scope of protection of this application.

[0061] Figure 1 A flowchart of a method for constructing an atmospheric model based on virtual temperature correction is provided for the embodiment of the present application. Figure 1 As shown, the method includes:

[0062] 101. Obtain the satellite's thermal infrared band radiation brightness observation value, and simultaneously record the satellite's observation zenith angle and azimuth angle data, and simultaneously obtain the surface cover type classification data, dew point temperature and corresponding height information;

[0063] In the above steps, the satellite thermal infrared band radiance observation value refers to the surface radiation energy intensity received by the satellite sensor in the thermal infrared band, usually 8 to 14 microns. The original digital value needs to be converted into a physical radiation value through radiation calibration, with the unit being watts per square meter per steradian per micron ( ); the satellite observation zenith angle and azimuth describe the angle between the satellite line of sight and the surface normal respectively to and the horizontal azimuth of the satellite relative to true north to , used to correct for the effects of atmospheric path radiation; surface cover classification data refers to information on the classification of features such as vegetation, water bodies, and buildings. Vegetation indices such as NDVI are calculated using the visible-near-infrared band and directly influence surface emissivity estimates. Dew point temperature and its corresponding altitude refer to the temperature at which water vapor condenses in the air and its altitude. These temperatures are obtained through meteorological sounding data or atmospheric profile models and are used to calculate atmospheric transmittance and path radiation.

[0064] In the embodiment of the present application, first, the original satellite thermal infrared image, such as Band 10 of Landsat-8, is processed by radiometric calibration, and the digital value is converted into the top-of-atmosphere radiance value using the gain value and offset value in the sensor parameters. The formula is: ,in, is the top-of-atmosphere radiance value, gain is the gain value, offset is the offset value, and DN is a digital value. The gain and offset values ​​are extracted from the pre-set image metadata file. For example, if the DN value of a pixel is 5000, the gain is 0.001, and the offset is -0.1, the radiance value is calculated as .

[0065] Synchronously calculate the observation zenith angle and azimuth, and calculate the pixel coordinates based on the satellite geometric imaging model Take two surfaces of different elevations, such as , the spatial point is obtained by the direct calculation model and The latitude, longitude and elevation observations are made, and the zenith angle is calculated by calculating the geocentric and Connect the same vector The angle is obtained by the formula:

[0066]

[0067] in, is the zenith angle, For spatial points and The square of the distance, for Elevation surface at ,for The elevation surface at The radius of the Earth is 6371 km, for example, When the distance is 150 kilometers, the calculation is ; Azimuth Then through the vector The angle between the projection on the horizontal plane and the north direction is determined.

[0068] Secondly, the normalized vegetation index is calculated using the visible-near infrared bands of the same-view satellite, such as the Band 4 red band and Band 5 near infrared band of Landsat-8. The formula is NDVI, where NDVI is the vegetation index, Band 4 is the red band reflectance value, and Band 5 is the near infrared band reflectance value. Divide vegetation cover :like ,but , combined with Calculating surface emissivity , such as vegetation ; Then, input the image into an atmospheric transport model such as NASA's ATMCO tool to obtain the time and center longitude and latitude, such as , generate atmospheric profile parameters, dew point temperature Query the temperature of the water vapor saturation layer, such as 2000 meters , corresponding to the height Direct output, atmospheric transmittance and atmospheric radiance Calculated by model integration as ; Finally, the radiance , zenith angle , azimuth Write the raster file by pixel coordinates, surface cover type and emissivity Generate classified raster, dew point temperature and height Attached to metadata as a global parameter to complete data synchronization storage.

[0069] In actual application, during a satellite observation mission in region A, researchers obtained raw observation values ​​using the thermal infrared sensor of satellite B. When the sensor recorded a raw value of 500 for a pixel, the thermal infrared radiation brightness value was calculated using the gain coefficient 0.000335 and the offset value 0.12 provided in the metadata file accompanying the satellite data: first multiply the gain coefficient by the raw value 0.000335×500=0.1675, and then add the offset value 0.1675+0.12=0.287, and the thermal infrared radiation brightness value of the pixel was 0.2875 watts per square meter per steradian per micron ( ). At the same time, the satellite platform automatically records the observation geometry data of the pixel: According to the satellite orbit parameters and pixel coordinates, the observation zenith angle is calculated to be 28 degrees through spatial geometry. The specific calculation process is: take the ground points P1 and P2 corresponding to the pixel at an altitude of 0 meters and 1000 meters, combined with the earth's radius of 6371 kilometers, and use the formula The angle value is obtained, and the azimuth is corrected to 152 degrees based on the angle between the satellite's flight direction and true north. Furthermore, the vegetation index is calculated using the visible light band data from the C satellite: when the red light band reflectance value is 0.08 and the near-infrared band is 0.33, the vegetation index = (0.33-0.08) / (0.33+0.08) = 0.61. Referring to the typical soil value of 0.05 and the typical vegetation value of 0.75, the vegetation coverage is calculated: , and then the surface emissivity is calculated to be 0.004×0.84+0.986=0.989. The dew point temperature data is obtained from the D meteorological data platform. The dew point temperature is recorded at the corresponding altitude of 3000 meters and the pressure layer of about 700 hPa. All data are ultimately integrated into atmospheric correction models to support surface temperature inversion and environmental monitoring applications.

[0070] In the overall approach described in step 101, thermal infrared band radiance observations are combined with multi-source parameters to achieve precise quantification of surface thermal radiation and temperature inversion. Zenith and azimuth data eliminate atmospheric scattering and surface anisotropy interference. Land cover classification significantly improves the accuracy of emissivity calculations. Dew point temperature and altitude information drive atmospheric water vapor attenuation correction, removing atmospheric radiation contamination. Ultimately, high-precision surface temperature spatial distribution data is obtained, effectively supporting diverse applications such as analyzing urban heat island thermal patterns, identifying hidden forest fires, monitoring geothermal anomalies, and dynamically assessing ecological and environmental thermal characteristics.

[0071] 102. Performing atmospheric path radiation distortion correction on the radiance observation value using a radiation calibration coefficient, and generating a calibrated radiance observation value by eliminating the atmospheric up-going radiation and down-going reflected radiation components;

[0072] Optionally, step 102 may specifically include the following steps:

[0073] 1021. Convert the radiation brightness observation value into an equivalent radiation energy value based on a radiation calibration coefficient, and separate the equivalent radiation energy value into an atmospheric upward radiation component and a surface reflected downward radiation component.

[0074] The upward radiation component is calculated by a top-of-atmosphere radiation model, and the downward radiation component is determined by multiplying the surface reflectivity by the downward radiation of the atmosphere.

[0075] 1022. Subtract the superposition value of the uplink radiation component and the downlink radiation component from the equivalent radiation energy value to generate a calibrated radiation brightness observation value representing the target surface radiation.

[0076] In the above steps, the radiation calibration coefficient refers to the radiation calibration coefficient generated by the standard radiation source of the ground blackbody radiation calibration field through the least squares fitting algorithm, which is used to establish a quantitative relationship between the digital value and the actual radiation energy; the atmospheric upward radiation component refers to the radiation energy scattered or emitted upward by the atmosphere itself, which needs to be calculated by the top of the atmosphere radiation model such as MODTRAN, and is mainly affected by aerosols and water vapor absorption; the downward radiation component reflected by the surface refers to the energy reflected by the surface after the solar radiation is scattered downward by the atmosphere, and is determined by the product of the surface reflectivity and the atmospheric downward radiation; the calibrated radiation brightness observation value refers to the physical quantity that characterizes the actual surface radiation after eliminating atmospheric interference, and the unit is watts per square meter per steradian per micron ( ).

[0077] In the embodiment of the present application, first, in step 1021, the radiation brightness observation value of the satellite thermal infrared band is converted into an equivalent radiation energy value using the radiation calibration coefficient. The formula is: , where gain and offset are extracted from satellite metadata files such as Landsat-8 MTL.txt. For example, if the DN value of a pixel is 6500, the gain is 0.0008, and the offset is -0.05, the equivalent radiation energy value is calculated as

[0078] ; Then, the equivalent radiation energy value is separated into the atmospheric upward radiation component Downward radiation component reflected from the ground , where the upward radiation component is calculated by inputting parameters such as dew point temperature and atmospheric transmittance into the top of atmosphere radiation model such as MODTRAN. For example, the input dew point temperature Atmospheric transmittance , model output ; The downward radiation component is transmitted through the surface reflectivity Downward radiation from the atmosphere The product of is determined by the formula , where the surface reflectivity Deduced from vegetation cover type data such as water bodies , bare soil , atmospheric downward radiation Calculated by the solar altitude angle and the atmospheric profile model, for example, a vegetation pixel ,but .

[0079] Secondly, the equivalent radiation energy value obtained in step 1022 Subtract the uplink radiation component Downlink radiation component The superposition value of generates the calibrated radiation brightness observation value, the formula is ,in, is the calibrated radiance observation value, is the equivalent radiation energy value, is the uplink radiation component, is the downward radiation component. For example, the above pixel is calculated as , which directly represents the true radiation characteristics of the surface; the specific calculation process needs to be performed pixel by pixel: first read As in 5.15, add the model output As 2.3 and based on reflectivity calculation For example, 1.55, the final difference is used as the calibration result, for example, a farmland pixel , then the calibration value , which can be used for subsequent surface temperature inversion or environmental monitoring applications.

[0080] In a practical application, researchers used the SDGSAT-1 satellite's thermal infrared spectrometer (TIS) to acquire thermal infrared data with a spatial resolution of 30 meters. They then used radiometric calibration to convert the raw observations (DN values) into radiance values. For example, when a pixel has a DN value of 850, the radiance is calculated using a gain factor of 0.00012 and an offset value of 0.05: Then, the atmospheric interference was eliminated based on the split window algorithm: first, the brightness temperature difference of the three thermal infrared channels of TIS, TIS2: 10.3–11.3μm and TIS3: 11.5–12.5μm, such as the brightness temperature of TIS2 301.5K and TIS3 299.8K, was used in combination with the atmospheric water vapor content parameter. Calculate the empirical coefficients A=2.1, B=-0.8, and finally invert the surface temperature: At the same time, vegetation coverage was calculated using the Landsat-8 visible light band: when the red light band reflectance is 0.08 and the near-infrared band reflectance is 0.33, the NDVI value is (0.33-0.08) / (0.33+0.08)=0.61; combined with the typical soil NDVI value of 0.05 and the vegetation NDVI value of 0.75, the vegetation coverage was calculated: , and then the surface emissivity is calculated as: ε=0.004×0.84+0.986=0.989. To improve the accuracy, the atmospheric parameters provided by ERA5 reanalysis data, such as the dew point temperature at 3000 meters, are integrated. , and remove the atmospheric upward radiation based on the radiation transfer equation Downward radiation reflected from the surface ((1-ε) ). For example, if 、 =0.075, ε=0.989, then the reflected component is: (1-0.989)×0.075=0.000825, and the surface radiation brightness after calibration is: Finally, the generated surface temperature and atmospheric virtual temperature field collaborative analysis showed that the temperature in urban built-up areas is higher than that in suburban areas. , and high temperature area The spatial overlap with areas with low vegetation coverage <0.3 is 89%, providing a quantitative basis for alleviating the urban heat island effect.

[0081] In the overall solution of step 102 above, the system eliminates the two-way superposition noise pollution of the upward radiation of the atmosphere and the downward radiation component reflected by the surface, thereby achieving deep purification of the original thermal infrared radiation brightness observation value. The generated calibrated radiation brightness observation value completely retains the true thermal radiation background characteristics of the target surface, and its spatial consistency and spectral fidelity are significantly optimized, especially in areas with drastic fluctuations in atmospheric water vapor content or complex terrain areas, which can effectively avoid the annihilation of thermal anomaly signals and temperature inversion deviations. This process directly constructs a standardized surface outgoing radiation dataset without atmospheric pollution interference, providing millimeter-level thermal radiation intensity benchmark data support for key links such as thermal infrared emissivity separation, temperature-emissivity collaborative inversion, and surface energy balance modeling, greatly improving the scientificity and reliability of dynamic monitoring data of surface thermal characteristics.

[0082] 103. Based on the reflectance data of the satellite red light band and the near-infrared band, calculate the normalized vegetation cover index, combine the surface cover type classification data, construct a mapping relationship between the surface cover type and the emissivity, and fit the nonlinear relationship between the vegetation cover index and the emissivity through the mapping relationship to invert the spatial distribution data of the surface emissivity;

[0083] Optionally, step 103 may specifically include the following steps:

[0084] 1031. Obtaining reflectance data of the satellite red light band and the near infrared band. For each pixel, substituting the reflectance data of the red light band and the near infrared band into a normalized difference calculation formula to generate a normalized vegetation cover index for each pixel.

[0085] 1032. Based on the surface cover type classification data, the surface cover types are divided into four categories: water bodies, vegetation, bare soil, and artificial buildings. A mapping function between the vegetation cover index and the emissivity is established for each surface cover type, wherein the vegetation type is fitted using an exponential function, and the non-vegetation type is fitted using a piecewise linear function.

[0086] 1033. Based on the mapping function, the vegetation cover index is input into the fitting function of the corresponding surface cover type, the emissivity value of each pixel is output, and the emissivity values ​​of all pixels are integrated to generate continuous spatial distribution data of the surface emissivity.

[0087] In the above steps, the normalized vegetation cover index is a quantitative indicator calculated by dividing the reflectance difference between the satellite red light band and the near-infrared band by their sum, which is used to reflect the vegetation cover density, and the value range is from -1 to 1; the normalized difference operation formula is the mathematical expression for calculating the index; the surface cover type classification data is the information of the four pre-classified land features: water bodies, vegetation, bare soil, and artificial buildings; the emissivity is the ratio of the actual thermal radiation capacity of the surface to the radiation capacity of an ideal black body at the same temperature, and the value range is 0 to 1; the mapping function is a mathematical model that establishes the relationship between NDVI and emissivity; the exponential function is a nonlinear fitting function for vegetation classes; the piecewise linear function is a piecewise mathematical fitting function for non-vegetation classes, such as water bodies, bare soil, and artificial buildings; the fitting function refers to the specific form of the optimized mapping function; the emissivity value is the radiation capacity value calculated for each pixel; the spatial distribution data is a continuous raster image generated by integrating the emissivity values ​​of all pixels according to their geographical location, which is used for surface temperature inversion.

[0088] In the embodiment of the present application, first, the reflectance data of the satellite red light band Band 4 and the near infrared band Band 5 are obtained through step 1031, and the normalized difference formula is substituted into each pixel in the image to calculate its NDVI value. The formula is NDVI Among them, the red light band reflectance Band4 represents the absorption intensity of red light by vegetation. The lower the value, the stronger the absorption. The near-infrared band reflectance Band5 represents the scattering ability of vegetation cells. The higher the value, the denser the vegetation. For example, the red light reflectance of a farmland pixel is 0.18 and the near-infrared reflectance is 0.55. The calculation process is: first calculate the numerator difference as 0.55-0.18=0.37, then calculate the denominator sum as 0.55+0.18=0.73, and the final NDVI value is , which indicates a medium vegetation cover level. After all pixels are calculated, an NDVI raster layer is formed and directly input into the subsequent mapping function construction step.

[0089] Secondly, in step 1032, the mapping function between NDVI and emissivity is constructed for the four types of ground objects by combining the surface cover type data. The vegetation type is fitted with an exponential function, and the formula is: , where a, b, and c are optimization parameters. For example, set a = 0.96, b = -0.05, c = -3.0 to capture the non-linear saturation effect caused by leaf stacking. When the NDVI increases from 0.2 to 0.6, the growth rate of the emissivity slows down. For non-vegetation types, a piecewise linear function is used: water is fixed at ε = 0.995, close to a black body; for bare soil, ε is set to 0.97 when NDVI ≤ 0, and linearly decreases to 0.96 when 0 < NDVI < 0.2 to reflect the influence of soil moisture. The materials of artificial buildings are uniform, and ε is constant at 0.96. For example, for a sparse grassland pixel with a surface type of vegetation and NDVI = 0.15, substituting into the exponential function, the emissivity is calculated as 0.96 + (-0.05) × exp(-3.0 × 0.15) = 0.96 - 0.05 × 0.637 ≈ 0.928. For a sandy pixel with a surface type of bare soil and NDVI = -0.05, ε = 0.97 is directly assigned.

[0090] Finally, through step 1033, the emissivity is calculated pixel by pixel based on the mapping function. The specific process is to read the surface type and NDVI value of the pixel. For vegetation types, the exponential function is called, and for non-vegetation types, values are assigned according to preset rules. All results are integrated into spatial distribution data. For example, for a wetland pixel classified as water, ε = 0.995 is directly output; for an urban building pixel, ε = 0.96 is output; for a forest pixel classified as vegetation and NDVI = 0.7, substituting into the formula gives ε = 0.96 - 0.05 × exp(-3.0 × 0.7) = 0.96 - 0.05 × 0.122 ≈ 0.954. The processing process covers all pixels, and after generating a raster image, the spatial variation of the emissivity is presented. For example, the high value in the water area is 0.995, and in the vegetation area, it gradually increases with the increase of NDVI, providing input parameters for temperature inversion. For example, for a bare soil pixel with NDVI = 0.18, since 0 < 0.18 < 0.2, the emissivity is linearly calculated as 0.97 - 0.1 × 0.20.18 = 0.973.

[0091] In practical applications, in a certain urban ecological environment monitoring project, researchers used Landsat-8 satellite data to conduct an inversion of the spatial distribution of the surface emissivity. First, the reflectance data of the satellite's red light band and near-infrared band were obtained, and the normalized difference vegetation index (NDVI) was calculated for each pixel: when the red light reflectance is 0.08 and the near-infrared reflectance is 0.33, substituting into the formula (near-infrared reflectance - red light reflectance) / (near-infrared reflectance + red light reflectance), the NDVI value is (0.33 - 0.08) / (0.33 + 0.08) = 0.61. Subsequently, combined with the surface cover type classification data, it was pre-divided into four categories: water, vegetation, bare soil, and artificial buildings, and the mapping relationship between the vegetation cover index and the emissivity was constructed. For the calculation of the vegetation cover degree, the pure vegetation NDVI threshold of 0.70 and the bare soil NDVI threshold of 0.05 were set. When the NDVI of a certain pixel is 0.61, the vegetation cover degree is calculated Based on the classification results, the vegetation class uses an exponential function to fit the emissivity: A piecewise linear function was used for non-vegetation classes. For example, the bare soil class had an emissivity of ε = 0.96 + 0.02 × 0.05 = 0.961. When NDVI = 0.05, the water class was directly assigned a fixed value of 0.995. Finally, a mapping function was used to traverse all pixels and integrate them to generate continuous spatial distribution data of surface emissivity. The results showed that emissivity in densely vegetated areas was generally higher than 0.98, while urban built-up areas were concentrated between 0.90 and 0.95. Water bodies maintained a stable value of 0.995. This data significantly improved the accuracy of subsequent surface temperature inversion parameters.

[0092] In the overall solution of step 103 above, a dynamic nonlinear mapping model of vegetation cover index and emissivity is constructed by collaboratively analyzing the spectral characteristics of red and near-infrared reflectance and the classification data of surface cover types, thereby achieving refined inversion of the emissivity of mixed pixels. Specifically, the normalized difference vegetation index is used to quantify the spatial heterogeneity of vegetation cover structure. Based on the four types of ground cover patterns of water bodies, vegetation, bare soil, and artificial buildings, functional response mechanisms of emissivity and vegetation index are established respectively: an exponential function is used in the vegetation cover area to accurately capture the saturation effect of vegetation canopy emissivity as coverage increases, and a piecewise linear function is used in the non-vegetated area to accurately characterize the nonlinear effects of bare soil porosity, building materials, and water surface fluctuation on emissivity. The generated continuous spatial emissivity distribution field effectively solves the scale confusion problem in the traditional single pixel emissivity assignment, significantly improves the accuracy of thermal radiation characteristics characterization in complex underlying surface areas such as high-density vegetation urban areas and bare soil-vegetation transition zones, and provides millimeter-level precision emissivity parameter field support for key links such as surface temperature inversion, thermal inertia model construction, and surface energy closure analysis.

[0093] 104. Generate vertically layered atmospheric profile data based on the dew point temperature and the corresponding height information, fuse the spatial distribution data with the atmospheric profile data, and combine with the angle data to construct correction parameters for the atmospheric radiation transmission path;

[0094] Optionally, step 104 may specifically include the following steps:

[0095] 1041. Based on the dew point temperature and the corresponding altitude information, generate atmospheric profile data by vertical altitude layer, wherein each layer includes temperature, water vapor density, and air pressure parameters;

[0096] 1042. Match the spatial distribution data of the surface emissivity with the atmospheric profile data according to geographic location, and fuse them into a three-dimensional gridded data set;

[0097] 1043. Calculate the slant path length of the electromagnetic wave through each atmospheric layer based on the atmospheric profile parameters in the three-dimensional gridded dataset and the angle data of the satellite observation zenith angle and azimuth angle;

[0098] 1044. Calculate the atmospheric attenuation coefficient and scattering weight layer by layer based on the oblique path length and the atmospheric parameters of each layer in the three-dimensional gridded dataset;

[0099] 1045. Integrate the attenuation coefficients and scattering weights of all layers to generate correction parameters of the atmospheric radiation transmission path that vary with altitude and angle.

[0100] In the above steps, vertically layered atmospheric profile data refers to a multi-layer atmospheric parameter set divided by altitude, each layer containing temperature, water vapor density, and air pressure, which is generated by inferring dew point temperature and corresponding altitude information; the three-dimensional gridded dataset is a three-dimensional parameter matrix formed by matching and fusing the two-dimensional plane of surface emissivity spatial distribution data with the vertically layered atmospheric profile data according to geographical location; the slant path length describes the slant distance that electromagnetic waves pass through each atmospheric layer, which is calculated by the zenith angle and layer height of satellite observations; the atmospheric attenuation coefficient represents the energy loss rate of electromagnetic waves on the slant path, including absorption and scattering; the scattering weight describes the redistribution ratio of radiation energy by atmospheric particles; the correction parameter is the path radiation correction value after integrating the attenuation and scattering effects of all layers, which is used to eliminate the interference of the atmosphere on the surface thermal radiation observation.

[0101] In the embodiment of the present application, first, in step 1041, vertically layered atmospheric profile data is generated based on the dew point temperature and the corresponding altitude information at a fixed altitude interval, such as every 100 meters. The specific process is: N layers are divided from the surface to the top of the atmosphere, such as 10km, and the center height of each layer is , calculate the temperature of each layer using the atmospheric state equation ,in, is the dew point temperature, is the temperature lapse rate, take , is the center height of each floor, is the corresponding height information. Water vapor density Among them, water vapor pressure Dew point temperature formula is reversed, air pressure in, is the sea level pressure, is the center height of each layer. For example, Michu , then the temperature of the 2000-meter layer is directly taken as Meter layer temperature .

[0102] Next, the surface emissivity is calculated by step 1042. The spatial distribution data of the two-dimensional shed, with a resolution of 30 meters, is fused with the atmospheric profile data: the operation process is to calculate the center coordinates (longitude and latitude) of each surface pixel. ,latitude ), vertically superimpose the atmospheric profile layer data of the corresponding position, such as 200-meter intervals, to form a three-dimensional grid unit , cell properties include surface emissivity , the height of this floor ,temperature , water vapor density , air pressure For example, a farmland pixel of , at altitude Associated atmospheric parameters .

[0103] Next, step 1043 is combined with the satellite observation zenith angle Calculate the oblique path length of each layer: the path length of the electromagnetic wave passing through the kth layer Degrees, among which is the layer thickness, fixed value such as 200 meters, unit m, Function angles need to be converted to radians. For example, when hour, .

[0104] Then, the atmospheric attenuation coefficient is calculated layer by layer in step 1044. unit With scattering weight Dimensionless: The attenuation coefficient formula is ,in is the standard water vapor density, is the standard atmospheric pressure, is the water vapor density, is the attenuation coefficient, For example, when the air pressure hour, The formula for the scattering weight is ,in The unit is m for the current layer height. For example hour, .

[0105] Finally, step 1045 integrates all layer parameters to generate path correction parameters, and the total correction amount ,unit The contribution of each layer is superimposed, and the formula is ,in is the attenuation coefficient of the kth layer, is the oblique path length, is the scattering weight. For example, a pixel is in 5 layers of atmosphere: the first layer contributes , Layer 2 contribution Total correction amount accumulated , which will be input into the surface temperature inversion model to eliminate atmospheric interference.

[0106] In practical applications, in a geothermal environment and atmosphere collaborative monitoring project, researchers first used ground meteorological stations and sounding data to obtain dew point temperature and corresponding height information. For example, when the ground dew point temperature at a certain observation point is When combined with the measured data of different altitude layers, such as the dew point temperature at 2000 meters , dew point temperature at 5000 meters , generate atmospheric profile data by vertical height layer, where each layer contains temperature, water vapor density, and air pressure parameters, such as the temperature at 1000 meters altitude , water vapor density , air pressure 900hPa. Subsequently, the spatial distribution data of surface emissivity, such as the 0.85~0.99 spatial grid data obtained by vegetation index inversion, are matched with the above atmospheric profile data according to geographical location and fused into a three-dimensional gridded data set. For example, the surface emissivity grid with a resolution of 1km is aligned with the 47-layer atmospheric parameters in the vertical direction from the surface to the 100hPa isobaric surface through longitude and latitude coordinates to form an integrated data model from the surface to the high altitude. Based on the zenith angle of satellite observation (such as ) and azimuth (e.g. ) angle data, combined with the atmospheric layer height parameters in the three-dimensional grid, calculate the oblique path length of the electromagnetic wave passing through each atmospheric layer. For example, when the satellite zenith angle is When the kth layer of atmosphere is at Slant path length in meters , the calculation formula is: Based on the oblique path length and the atmospheric parameters of each layer, temperature, water vapor density, and air pressure, the atmospheric attenuation coefficient and scattering weight are calculated layer by layer. For example, the calculation of the water vapor attenuation coefficient is combined with the water vapor density of the layer. , using the empirical formula:

[0107] At the same time, the scattering weight is calculated based on the aerosol optical thickness (such as 0.1) and the scattering phase function (such as 0.85). The attenuation coefficients of all layers are integrated with the scattering weights to generate correction parameters that vary with height and angle. For example, the attenuation coefficients of 47 layers of the atmosphere ( ) and the scattering weight ( ), construct the atmospheric transmittance correction function: Finally, the output parameters are used to correct the radiation signal received by the satellite and improve the accuracy of surface temperature inversion.

[0108] In the overall solution of step 104 above, vertically layered atmospheric profile parameters are inverted using dew-point temperature and altitude information, and a three-dimensional gridded atmospheric radiation transmission path correction model is constructed in conjunction with the spatial distribution field of surface emissivity and satellite observation geometric angle data. The oblique path length and height-dependent attenuation coefficient are calculated based on the layered temperature, humidity, and pressure profiles to analyze the cumulative optical path effect of electromagnetic waves passing through each atmospheric layer. The horizontal distribution field of emissivity is integrated to achieve spatial coupling of surface-atmosphere parameters and quantify the attenuation weights of path scattering, absorption, and refraction on thermal radiation propagation. Ultimately, an angle-adaptive radiation transmission correction parameter set is generated, significantly improving the accuracy of eliminating atmospheric thermal radiation interference in high-humidity areas and complex terrain areas, providing a millimeter-level precision atmospheric path correction benchmark for thermal infrared surface temperature inversion.

[0109] 105. Establish an atmospheric radiation transfer model based on the correction parameters, input the calibrated radiation brightness observation value into the atmospheric radiation transfer model, invert the atmospheric virtual temperature vertical profile by solving the radiation transfer equation, and output the atmospheric virtual temperature field spatial distribution data.

[0110] Optionally, step 105 may specifically include the following steps:

[0111] 1051. Construct a vertically layered atmospheric radiation transfer model framework, dividing the atmosphere into multiple continuous layers according to altitude, with each layer corresponding to the attenuation coefficient and scattering weight in the correction parameters;

[0112] 1052. In the model framework, the angle data is converted into the inclined path length of the electromagnetic wave propagation path in each layer, and the path cumulative impact factor is generated by combining the attenuation coefficient and scattering weight of each layer;

[0113] 1053. Convert the calibrated radiation brightness observation value into an equivalent radiation energy input value, and input the equivalent radiation energy input value into an atmospheric radiation transfer model containing the path cumulative impact factor.

[0114] 1054. Calculate dry delay component data and wet delay component data by solving a radiation transfer equation using the atmospheric radiation transfer model according to the calibrated radiation brightness observation value.

[0115] Among them, step 1054 may specifically include the following processes: establishing a radiation transfer relationship, expressing the calibrated radiation brightness observation value as a linear superposition of a dry delay component and a wet delay component, wherein the dry delay component corresponds to the absorption effect of the dry atmospheric component, and the wet delay component corresponds to the absorption effect of the water vapor component; performing spectral decomposition on the calibrated radiation brightness observation value, separating independent spectral response segments for different absorption characteristic bands, and obtaining the radiation brightness decomposition value corresponding to each spectral response segment; based on the radiation transfer relationship, selecting the spectral response segment dominated by the dry atmospheric component, substituting the corresponding radiation brightness decomposition value into the independent expression of the dry delay component, and generating dry delay component data; based on the radiation transfer relationship, selecting the spectral response segment dominated by the water vapor component, substituting the corresponding radiation brightness decomposition value into the independent expression of the wet delay component, and generating wet delay component data.

[0116] 1055. Based on the dry delay component and the wet delay component data and in combination with the atmospheric correction parameter table, calculate the contribution degree data of each atmospheric layer to the dry delay component and the wet delay component;

[0117] 1056. Based on the contribution data, adjust the virtual temperature value of each atmospheric layer by iteratively solving the radiation transfer equation until the error between the calculated radiance value output by the atmospheric radiation transfer model and the calibrated radiance observation value is minimized, thereby generating an optimized vertical virtual temperature sequence;

[0118] 1057. Based on the optimized vertical virtual temperature sequence, atmospheric virtual temperature vertical profile data of a single geographic point is formed, and the atmospheric virtual temperature vertical profile data of all geographic points are integrated to obtain a continuous atmospheric virtual temperature field spatial distribution map.

[0119] In the above steps, the vertically layered atmospheric radiation transmission model framework refers to a multi-layer calculation structure that divides the atmosphere into fixed height intervals; the inclined path length describes the geometric distance that the electromagnetic wave obliquely passes through each layer of the atmosphere; the path cumulative impact factor characterizes the total attenuation effect of the radiation energy on the entire transmission path; the equivalent radiation energy input value refers to the actual radiation intensity of the surface after eliminating atmospheric interference; the dry delay component reflects the absorption effect of dry atmospheric components (nitrogen, oxygen, etc.) on radiation; the wet delay component describes the absorption effect of water vapor on radiation; the radiation transmission relationship is the mathematical expression of the observed radiation formed by the superposition of dry and wet delay components; Spectral decomposition refers to the technology of separating and processing radiation signals according to different absorption bands; spectral response segment refers to the specific wavelength range that is strongly absorbed by water vapor or dry atmosphere; the atmospheric correction parameter table contains the response coefficients of each layer of the atmosphere to different delay components; the contribution degree data quantifies the proportion of the influence of each layer of the atmosphere on the total delay component; the virtual temperature value is a temperature parameter that represents the thermal state of moist air equivalent to the properties of dry air; the vertical virtual temperature sequence refers to the set of virtual temperature values ​​at different altitudes; the atmospheric virtual temperature vertical profile data is the virtual temperature distribution at different altitudes of a single geographical location; the atmospheric virtual temperature field spatial distribution map is the three-dimensional temperature field formed by the vertical profiles of all geographical points.

[0120] In the embodiment of the present application, first, step 1051 is used to construct a radiation transfer model framework including a multi-layer structure, and the atmosphere is divided into continuous layers from the surface to the top layer at fixed height intervals such as 200 meters. Each layer is associated with the attenuation coefficient αk (unit: ) and the scattering weight wk (dimensionless). The model framework presets the initial virtual temperature reference value of each layer (e.g. 300K near the ground and 200K at the top layer) as the starting point for subsequent iterative optimization.

[0121] Secondly, step 1052 combines the satellite observation zenith angle Calculate the inclined path length of each layer using the formula:

[0122]

[0123] in, is the inclined path length of each layer, is the layer thickness, is the zenith angle. Then the path cumulative impact factor is generated ,in, is the attenuation coefficient of the i-th layer of atmosphere, which represents the rate of loss of radiation energy per unit path length (units ); is the slant path length of the electromagnetic wave through the i-th layer of the atmosphere, calculated from the layer thickness and the zenith angle of the satellite observation (unit: km); is the scattering weight of the i-th layer of atmosphere, which characterizes the scattering enhancement effect of particles such as aerosols on radiation (dimensionless); is the basic attenuation caused by the i-th layer, describing the radiation energy absorption loss. For example, when When single layer , the cumulative impact factor of the first three layers .

[0124] Then, the calibrated radiance observation value is converted to unit , directly loaded into the model as the equivalent radiation energy input value, where The atmospheric up-radiation and surface reflection components have been eliminated, such as the input value of a certain pixel .

[0125] Then, the radiation signal is decomposed in step 1054. First, the radiation transmission relationship is established.

[0126] is the dry delay component, is the wet delay component; then, the input radiation value is spectrally decomposed and characteristic bands are selected, such as the dry atmosphere dominant band , water vapor dominant band ;Radiance decomposition values ​​separated by dominant bands, such as dry bands , respectively, into the independent expressions: dry delay component ,in, is the correction parameter. Example calculation: hour, Wet delay component ,in, Example calculation for water vapor density: hour, .

[0127] Next, the contribution of each layer is calculated based on the atmospheric correction parameter table in step 1055. ,in, is the dry delay response coefficient, is the dry delay component, such as the ground layer For example, the ground layer contribution Wet delay contribution: contribution of the kth layer

[0128] ,in is the wet delay response coefficient, is the wet delay component, such as the ground layer , calculated .in and The value changes with the height, and the response coefficient of the high-rise layer decreases, such as the 3km layer .

[0129] Then, the virtual temperature value is iteratively optimized through step 1056. First, the virtual temperature of each layer is initialized. , such as 298K for the surface and 290K for the 1km layer, and then calculate the model output radiation value , if the error , such as the initial

[0130] , adjust the virtual temperature according to the contribution ratio: ,in, is the step size coefficient,

[0131] For thermal radiation sensitivity, calculate the adjustment amount Then update the virtual temperature, such as the surface layer is adjusted to Finally, repeat the above steps until the error , output optimized virtual temperature sequence, such as 298.4K for the surface and 290.2K for the 1km layer.

[0132] Finally, a three-dimensional virtual temperature field is generated through step 1057. Single point vertical profile: For each geographic coordinate, such as 118.5° east longitude, 118.5° north latitude , output height sequence virtual temperature values, such as 0km: 298.4K, 1km: 290.2K, 3km: 265.7K on the surface; then perform spatial integration and interpolate all pixel points to form a continuous three-dimensional grid (longitude × latitude × altitude); finally, visualize the output to generate a spatial distribution map, which can show the thermal structure characteristics.

[0133] In practical applications, in a certain city thermal environment monitoring project, researchers used drones equipped with thermal infrared sensors with a resolution of 0.5 meters and multispectral sensors to collect data, and inverted the surface temperature and atmospheric virtual temperature field through a deep learning model. First, the thermal infrared images obtained by the drone were radiometrically calibrated to generate a brightness temperature orthophoto map with a temperature range of 283.6-315.2 K. The normalized vegetation index was calculated based on the multispectral data, with an average NDVI of 0.61. The spatial distribution of the surface emissivity was inverted using an empirical model, with a range of 0.96-0.99 in vegetation areas and 0.89-0.93 in built-up areas. The low-altitude atmospheric parameters were simultaneously integrated, with an atmospheric transmittance of 0.88 and an upward radiation of 4.82. , downlink radiation 6.37 Subsequently, a three-dimensional radiation transfer model was constructed to eliminate the effects of building shading. The resulting spatial distribution of surface temperature was 293.5 K ± 1.2 K in vegetation-covered areas, 305.7 K ± 2.8 K in urban built-up areas, and 290.1 ​​K ± 0.6 K in water bodies. Simultaneously, a vertical profile of atmospheric virtual temperature was inverted using brightness temperature data from the Fengyun-3E microwave thermometer (MWTS-3) satellite using a pre-trained deep neural network with a hidden layer structure of 700-700-700-700. The results showed that the virtual temperature field at the lower troposphere (1000 hPa) was consistent with the ERA5 reanalysis data, with a root mean square error of 0.6 K. In the inversion of the typhoon system's warm core structure, the virtual temperature deviation at 300 hPa was less than 0.8 K. The resulting spatial distribution of the atmospheric virtual temperature field revealed that the core area of ​​the urban heat island had a virtual temperature 3.5 K higher than that of the suburbs, while the typhoon's warm core had a virtual temperature of -15.2 K at 500 hPa. Verification shows that the absolute error between the surface temperature inversion results and the ground measured values ​​is less than 1.5 K, and the correlation coefficient between the atmospheric virtual temperature profile and the sounding data exceeds 0.98.

[0134] In the overall solution of step 105, a vertically layered atmospheric radiation transfer model is constructed, deeply coupling the calibrated thermal radiation brightness observations with the three-dimensional path correction parameters to achieve accurate modeling and spatial reconstruction of the atmospheric virtual temperature field. This ultimately generates highly spatially coherent three-dimensional distribution data of the atmospheric virtual temperature field, which clearly demonstrates the thermal structure of the boundary layer, the vertical gradient of the inversion layer, and the thermodynamic fluctuations at the tropopause. This provides key atmospheric thermodynamic parameter support for analyzing the triggering mechanism of severe convective weather, tracking the thermal dynamics of volcanic ash clouds, and modeling the global energy balance.

[0135] The following is a complete embodiment of steps 101 to 105:

[0136] like Figure 2 As shown in the figure, in a geothermal resource exploration project in a certain area, researchers used the Resource-1 02E star thermal infrared sensor to carry out surface temperature inversion applications. Figure 2 The thermal infrared sensor of Ziyuan-1 02E in the upper center (shown as 21 in the figure) is The zenith angle (shown as 22 in the figure) is used to observe the ground.

[0137] First, obtain the passing satellite thermal infrared image, with a band range of 8.0-10.0μm and a spatial resolution of 16 meters, and simultaneously record the satellite observation zenith angle , and collect the dew point temperature measured by the ground weather station on that day The original thermal infrared imagery was radiometrically calibrated, converting the DN values ​​into radiance values. The pixel brightness temperature range was 283.6-317.4K. A split-window algorithm was then used to eliminate interference from upwelling atmospheric radiation and downwelling reflected radiation, generating a calibrated radiance dataset.

[0138] Subsequently, based on the red light band and near-infrared band data of the Tongjing satellite, the normalized difference vegetation index NDVI was calculated with an average value of 0.38. Combined with the surface cover classification map, the four categories of farmland, forest land, bare land, and water body were constructed to construct a nonlinear inversion model of emissivity, and the spatial distribution data of emissivity was output. Among them, the surface cover type was distinguished by different texture patterns: the farmland area (as shown in Figure 27) was marked with a grid texture with an emissivity of 0.92-0.95; the forest area (as shown in Figure 28) was represented by a circular texture with an NDVI average of 0.38; the bare land area (as shown in Figure 29) was distinguished by a dot texture with an emissivity of 0.87-0.89; the water area (as shown in Figure 210) was marked with a wavy texture, corresponding to a dew point temperature of -2.1℃.

[0139] At the same time, 20 layers of atmospheric profiles are generated based on the vertical distribution of dew point temperature (this application is a simplified drawing, Figure 2 Only four layers are divided, but actual applications may require the division of multiple atmospheric profiles). The layered structure of the atmospheric profile is as follows from top to bottom: the 3000m altitude layer (as shown in Figure 23) is the upper atmospheric layer, with a temperature of 268.4K±0.9K and an attenuation coefficient of 0.32; the 1500m altitude layer (as shown in Figure 24) is the key layer of the middle atmosphere, contributing 12.7% to the wet delay; the 1000hPa layer (as shown in Figure 25) is the standard atmospheric pressure reference layer, with a temperature of 5.2℃ and a water vapor density of 4.3g / m³; the surface layer (as shown in Figure 26) is the temperature inversion target layer, with a temperature of 291.3K±1.7K. Each layer contains temperature, water vapor density, and air pressure parameters, such as the 1000hPa layer temperature , water vapor density , fuse the emissivity data with the atmospheric profile, and calculate the slant path attenuation coefficient based on the satellite observation angle. The attenuation coefficient at the 3000-meter altitude layer is 0.32, and the atmospheric correction parameter set is generated.

[0140] The calibrated radiance data was input into a vertically layered atmospheric radiation transfer model. Spectral decomposition was used to separate the dry delay component (the 10-12μm band dominated by oxygen absorption) from the wet delay component (the 8-9μm band dominated by water vapor absorption). The contribution weights of each layer were quantified using an atmospheric correction parameter table. The 1500-meter layer contributed 12.7% to the wet delay. An iterative optimization algorithm was used to adjust the initial virtual temperature values ​​of each layer. When the root mean square error between the model output radiance and the observed value dropped to 1.85 K, an optimized vertical virtual temperature sequence was generated, such as a surface temperature of 291.3 K ± 1.7 K and a temperature of 268.4 K ± 0.9 K at the 3000-meter altitude layer. Finally, a spatial distribution map of surface temperature was generated with a 2-kilometer grid resolution. Verification showed that the inversion results, when compared with the measured data from 30 ground temperature measurement points, had an average absolute error of 1.2 K. The map clearly reveals that the surface temperatures of the three geothermal anomaly areas were 3.5 K higher than those of the surrounding areas. Drilling confirmed that the geothermal water temperatures at two of the three areas, at 1500 meters underground, reached 2.5 K. The three geothermal anomalies in the underground layer, marked by dashed circles, are: Anomaly Area 1 (shown as 211 in the figure), located above farmland, with a temperature 3.5K higher than the surrounding area; Anomaly Area 2 (shown as 212 in the figure), located above bare land, with a temperature 3.6K higher than the surrounding area; and Anomaly Area 3 (shown as 213 in the figure), located near a body of water, with a temperature 3.5K higher than the surrounding area. The two vertical black lines represent drilling verification points: Verification Point 1 (shown as 214 in the figure) and Verification Point 2 (shown as 215 in the figure) both discovered 62°C geothermal water 1500m below the surface, validating the accuracy of the remote sensing identification results and confirming the effective support of the inversion results for geothermal resource exploration.

[0141] Figure 3 The present invention provides a schematic diagram of a system for constructing an atmospheric model based on virtual temperature correction. Figure 3 As shown, the system includes:

[0142] The acquisition module 31 is used to obtain the radiance observation value of the satellite thermal infrared band, and simultaneously record the angle data of the satellite observation zenith angle and azimuth angle, and obtain the surface cover type classification data, dew point temperature and corresponding height information;

[0143] A calibration module 32 is configured to perform atmospheric path radiation distortion correction on the radiance observation value using a radiation calibration coefficient, and generate a calibrated radiance observation value by eliminating the atmospheric up-going radiation and down-going reflected radiation components;

[0144] An inversion module 33 is configured to calculate a normalized vegetation cover index based on the reflectance data of the satellite red light band and the near-infrared band, construct a mapping relationship between the surface cover type and the emissivity in combination with the surface cover type classification data, and fit the nonlinear relationship between the vegetation cover index and the emissivity through the mapping relationship to invert the spatial distribution data of the surface emissivity;

[0145] A construction module 34 is configured to generate vertically layered atmospheric profile data based on the dew point temperature and the corresponding height information, fuse the spatial distribution data with the atmospheric profile data, and construct correction parameters for the atmospheric radiation transmission path in combination with the angle data;

[0146] The output module 35 is used to establish an atmospheric radiation transfer model based on the correction parameters, input the calibrated radiation brightness observation value into the atmospheric radiation transfer model, invert the atmospheric virtual temperature vertical profile by solving the radiation transfer equation, and output the atmospheric virtual temperature field spatial distribution data.

[0147] Figure 3 The atmospheric model construction system based on virtual temperature correction can be executed Figure 1 The implementation principle and technical effects of the atmospheric model construction method based on virtual temperature correction described in the illustrated embodiment will not be elaborated on here. The specific manner in which each module and unit performs operations in the atmospheric model construction system based on virtual temperature correction in the above embodiment has been described in detail in the embodiments of the method and will not be elaborated on here.

[0148] In one possible design, Figure 3 The atmospheric model construction system based on virtual temperature correction of the embodiment shown can be implemented as a computing device, such as Figure 4 As shown, the computing device may include a storage component 41 and a processing component 42;

[0149] The storage component 41 stores one or more computer instructions, wherein the one or more computer instructions are called and executed by the processing component 42 .

[0150] The processing component 42 is used for the above Figure 1 The embodiment provides a method for constructing an atmospheric model based on virtual temperature correction.

[0151] The processing component 42 may include one or more processors to execute computer instructions to perform all or part of the steps in the above method. Of course, the processing component may also be implemented as one or more application-specific integrated circuits (ASICs), digital signal processors (DSPs), digital signal processing devices (DSPDs), programmable logic devices (PLDs), field programmable gate arrays (FPGAs), controllers, microcontrollers, microprocessors, or other electronic components to perform the above method.

[0152] The storage component 41 is configured to store various types of data to support operations at the terminal. The storage component can be implemented by any type of volatile or non-volatile memory device, or a combination thereof, such as static random access memory (SRAM), electrically erasable programmable read-only memory (EEPROM), erasable programmable read-only memory (EPROM), programmable read-only memory (PROM), read-only memory (ROM), magnetic memory, flash memory, magnetic disk, or optical disk.

[0153] Of course, a computing device may also include other components, such as input / output interfaces, display components, communication components, etc.

[0154] The input / output interface provides an interface between the processing component and the peripheral interface module, which can be an output device, an input device, etc.

[0155] The communication component is configured to facilitate, among other things, wired or wireless communications between the computing device and other devices.

[0156] Among them, the computing device can be a physical device or an elastic computing host provided by a cloud computing platform, etc. In this case, the computing device can refer to a cloud server, and the above-mentioned processing components, storage components, etc. can be basic server resources rented or purchased from the cloud computing platform.

[0157] The present application also provides a computer storage medium storing a computer program, wherein the computer program can achieve the above-mentioned Figure 1 The illustrated embodiment is a method for constructing an atmospheric model based on virtual temperature correction.

[0158] Those skilled in the art will clearly understand that, for the convenience and brevity of description, the specific working processes of the systems, devices and units described above can refer to the corresponding processes in the aforementioned method embodiments and will not be repeated here.

[0159] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units, i.e., they may be located in one location or distributed across multiple network units. Some or all of the modules may be selected based on actual needs to achieve the objectives of the present embodiment. Persons of ordinary skill in the art will be able to understand and implement the present invention without inventive effort.

[0160] Through the above description of the embodiments, those skilled in the art will clearly understand that each embodiment can be implemented using software plus a necessary general-purpose hardware platform, or of course, hardware. Based on this understanding, the essence of the above technical solution, or the portion that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, a magnetic disk, or an optical disk, and includes a number of instructions for causing a computer device (such as a personal computer or server) to execute the methods described in each embodiment or certain portions of the embodiments.

[0161] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present application, rather than to limit them. Although the present application has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for constructing an atmospheric model based on virtual temperature correction, characterized in that: include: Obtain the satellite's thermal infrared band radiation brightness observation values, and simultaneously record the satellite's observation zenith angle and azimuth angle data, while obtaining the surface cover type classification data, dew point temperature and corresponding height information; Performing atmospheric path radiation distortion correction on the radiance observation value using a radiation calibration coefficient, and generating a calibrated radiance observation value by eliminating the atmospheric up-going radiation and down-going reflected radiation components; Based on the reflectance data of the satellite red light band and the near-infrared band, a normalized vegetation cover index is calculated. In combination with the surface cover type classification data, a mapping relationship between the surface cover type and the emissivity is constructed. The nonlinear relationship between the vegetation cover index and the emissivity is fitted through the mapping relationship to invert the spatial distribution data of the surface emissivity. Based on the dew point temperature and the corresponding height information, vertically layered atmospheric profile data is generated, the spatial distribution data is integrated with the atmospheric profile data, and correction parameters of the atmospheric radiation transmission path are constructed in combination with the angle data; An atmospheric radiation transfer model is established based on the correction parameters, and the calibrated radiation brightness observation value is input into the atmospheric radiation transfer model. The vertical profile of the atmospheric virtual temperature is inverted by solving the radiation transfer equation, and the spatial distribution data of the atmospheric virtual temperature field is output.

2. The method according to claim 1, characterized in that Based on the reflectance data of the satellite red light band and the near-infrared band, the normalized vegetation cover index is calculated. Combined with the surface cover type classification data, a mapping relationship between the surface cover type and the emissivity is constructed. The nonlinear relationship between the vegetation cover index and the emissivity is fitted through the mapping relationship to invert the spatial distribution data of the surface emissivity, including: Obtaining the reflectance data of the satellite red light band and the near infrared band, for each pixel, substituting the values ​​of the reflectance data of the red light band and the near infrared band into the normalized difference calculation formula to generate the normalized vegetation cover index of each pixel; Combined with the surface cover type classification data, the surface cover types are divided into four categories: water bodies, vegetation, bare soil, and artificial buildings. A mapping function between the vegetation cover index and the emissivity is established for each surface cover type, where the vegetation type is fitted with an exponential function and the non-vegetation type is fitted with a piecewise linear function. Based on the mapping function, the vegetation cover index is input into the fitting function of the corresponding surface cover type, the emissivity value of each pixel is output, and the emissivity values ​​of all pixels are integrated to generate continuous spatial distribution data of surface emissivity.

3. The method according to claim 1, characterized in that Based on the dew point temperature and the corresponding height information, vertically layered atmospheric profile data is generated, the spatial distribution data is fused with the atmospheric profile data, and the correction parameters of the atmospheric radiation transmission path are constructed in combination with the angle data, including: Based on the dew point temperature and the corresponding height information, generate atmospheric profile data by vertical height layer, wherein each layer includes temperature, water vapor density, and air pressure parameters; Matching the spatial distribution data of the surface emissivity with the atmospheric profile data according to geographical location and fusing them into a three-dimensional gridded data set; Calculating the slant path length of electromagnetic waves passing through each atmospheric layer based on the atmospheric profile parameters in the three-dimensional gridded dataset and the angular data of the zenith angle and azimuth angle observed by the satellite; Calculating the atmospheric attenuation coefficient and scattering weight layer by layer based on the oblique path length and the atmospheric parameters of each layer in the three-dimensional gridded data set; The attenuation coefficients and scattering weights of all layers are integrated to generate correction parameters of the atmospheric radiation transmission path that vary with altitude and angle.

4. The method according to claim 1, characterized in that , by solving the radiation transfer equation to invert the vertical profile of atmospheric virtual temperature, and output the spatial distribution data of atmospheric virtual temperature field, including: According to the calibrated radiance observation value, the dry delay component data and the wet delay component data are solved by solving the radiative transfer equation using the atmospheric radiative transfer model; Based on the dry delay component and wet delay component data, combined with the atmospheric correction parameter table, calculate the contribution degree data of each atmospheric layer to the dry delay component and the wet delay component; Based on the contribution data, the virtual temperature value of each atmospheric layer is adjusted by iteratively solving the radiation transfer equation until the error between the calculated radiance value output by the atmospheric radiation transfer model and the calibrated radiance observation value is minimized, thereby generating an optimized vertical virtual temperature sequence; According to the optimized vertical virtual temperature sequence, atmospheric virtual temperature vertical profile data of a single geographical point is formed, and the atmospheric virtual temperature vertical profile data of all geographical points are integrated to obtain a continuous atmospheric virtual temperature field spatial distribution map.

5. The method according to claim 4, characterized in that According to the calibrated radiance observation value, the dry delay component data and the wet delay component data are solved by solving the radiative transfer equation using the atmospheric radiative transfer model, including: Establishing a radiation transfer relationship to express the calibrated radiance observation value as a linear superposition of a dry delay component and a wet delay component, wherein the dry delay component corresponds to the absorption effect of dry atmospheric components and the wet delay component corresponds to the absorption effect of water vapor components; Performing spectral decomposition on the calibrated radiance observation value, separating independent spectral response segments for different absorption characteristic bands, and obtaining radiance decomposition values ​​corresponding to each spectral response segment; Based on the radiation transfer relationship, a spectral response segment dominated by dry atmospheric components is selected, and the corresponding radiation brightness decomposition value is substituted into the independent expression of the dry delay component to generate dry delay component data; Based on the radiation transfer relationship, a spectral response segment dominated by water vapor components is selected, and the corresponding radiation brightness decomposition value is substituted into the independent expression of the wet delay component to generate wet delay component data.

6. The method according to claim 1, characterized in that Establishing an atmospheric radiation transfer model according to the correction parameters and inputting the calibrated radiance observation value into the atmospheric radiation transfer model includes: Constructing a vertically layered atmospheric radiation transfer model framework, dividing the atmosphere into multiple continuous layers according to altitude, with each layer corresponding to the attenuation coefficient and scattering weight in the correction parameters; In the model framework, the angle data is converted into the inclined path length of the electromagnetic wave propagation path in each layer, and the cumulative impact factor of the path is generated by combining the attenuation coefficient and scattering weight of each layer; The calibrated radiation brightness observation value is converted into an equivalent radiation energy input value, and the equivalent radiation energy input value is input into an atmospheric radiation transfer model containing the path cumulative impact factor.

7. The method according to claim 1, characterized in that The radiometric calibration coefficient is used to correct the atmospheric path radiation distortion of the radiometric brightness observation value, and the calibrated radiometric brightness observation value is generated by eliminating the atmospheric up-going radiation and down-going reflected radiation components, including: The radiance observation value is converted into an equivalent radiant energy value based on a radiation calibration coefficient, and the equivalent radiant energy value is separated into an atmospheric upward radiation component and a downward radiation component reflected by the surface, wherein the upward radiation component is calculated by a top-of-atmosphere radiation model, and the downward radiation component is determined by multiplying the surface reflectivity and the atmospheric downward radiation; The superposition value of the upgoing radiation component and the downgoing radiation component is subtracted from the equivalent radiation energy value to generate a calibrated radiation brightness observation value representing the target surface radiation.

8. An atmospheric model construction system based on virtual temperature correction, characterized in that: include: The acquisition module is used to obtain the radiance observation value of the satellite thermal infrared band, and synchronously record the angle data of the satellite observation zenith angle and azimuth angle, and obtain the surface cover type classification data, dew point temperature and corresponding height information; a calibration module for performing atmospheric path radiation distortion correction on the radiance observation value using a radiation calibration coefficient, and generating a calibrated radiance observation value by eliminating the atmospheric up-going radiation and down-going reflected radiation components; An inversion module is used to calculate the normalized vegetation cover index based on the reflectance data of the satellite red light band and the near-infrared band, and to construct a mapping relationship between the surface cover type and the emissivity in combination with the surface cover type classification data. The nonlinear relationship between the vegetation cover index and the emissivity is fitted through the mapping relationship to invert the spatial distribution data of the surface emissivity. A construction module is used to generate vertically layered atmospheric profile data based on the dew point temperature and the corresponding height information, fuse the spatial distribution data with the atmospheric profile data, and construct correction parameters of the atmospheric radiation transmission path in combination with the angle data; The output module is used to establish an atmospheric radiation transfer model based on the correction parameters, input the calibrated radiation brightness observation value into the atmospheric radiation transfer model, invert the atmospheric virtual temperature vertical profile by solving the radiation transfer equation, and output the atmospheric virtual temperature field spatial distribution data.

Images