CN103927455B - Land aerosol optical property retrieval method based on Gaofen-1 satellite - Google Patents

Abstract

The invention discloses a land aerosol optical property retrieval method based on the Gaofen-1 satellite. The land aerosol optical property retrieval method aims to solve the problem that an aerosol optical thickness retrieval method is limited by short wave infrared data, and technical supports are provided for air quality monitoring based on the domestic high-definition satellite. The land aerosol optical property retrieval method comprises the following steps that S1, an atmospheric optical parameter lookup table is generated; S2, a surface reflectance prior database is generated; S3, dark object pixels suitable for retrieval are selected; S4, the observation geometrical conditions of the dark object pixels are calculated; S5, the atmospheric optical parameters of the dark object pixels are obtained in an interpolation mode; S6, aerosol optical thickness retrieval is carried out based on radiation transfer. The method is efficient in operation, retrieval on the regional high-resolution atmospheric aerosol optical property is facilitated, and beneficial data sources are provided for particulate pollution monitoring.

Description

Land optical properties of aerosol inversion method based on a number satellite of high score

Technical field

The present invention relates to applications of atmospheric remote sensing techniques field, more particularly, to one kind cover ccd data based on number satellite width of high score Colloidal sol optical property inversion method.

Background technology

The aerocolloidal research of satellite remote sensing starts from last century middle nineteen seventies, in last more than 30 year of eighties of last century, Satellite remote sensing aerosol for studying and being formed businessization use is included in volcanic eruption aerosol monitoring, noaa (state of the U.S. Family ocean and Atmospheric Administration) that deep-sea overhead is derived from the colloidal sol floor optics that sandstorm and forest fire cause is thick for series of satellites Degree remote sensing, infrared remote sensing sand and dust aerosol remote sensing, total amount of ozone imaging spectrometer (toms) ultraviolet band are molten to absorbability gas Several aspect such as remote sensing of glue.The intermediate-resolution imaging spectral that U.S. earth observing system plan (eos) terra and aqua carries The distribution on global product of 10 km aerosol optical depths issued by instrument (modis) to the whole world；French promotion from 1996 Polarization multi-angle camera (polder) detects ocean and the aerocolloidal research in land.Up to the present, molten using satellite remote sensing gas Glue has defined certain global aerosol detection system.

Current more ripe and widely used land aerosol satellite remote sensing inversion algorithm is dark goal approach.The party Fado utilizes short infrared wave band (near 2.1 microns) to identify dark target and obtain Reflectivity for Growing Season, and then according to this wave band With the empirical relation of the Reflectivity for Growing Season of red, blue wave band, inverting aerosol optical depth.But the ccd number of the domestic satellite such as gf1 According to lack short infrared wave band it is difficult to be suitable for above-mentioned dark goal pels method.The present invention propose one kind be applied to high score one (under Abbreviation gf1) satellite width cover ccd data aerosol optical depth inversion method, the restriction of short-wave infrared data can be broken through, be Carrying out air quality monitoring based on domestic high-resolution satellite provides technical support.

Content of the invention

The technical problem to be solved in the present invention is: provide a kind of based on gf1 satellite width cover the operation of ccd data efficiently, The method of precision preferable region land aerosol remote-sensing inversion.

For solving the above problems, the invention provides a kind of land optical properties of aerosol inverting side based on gf1 satellite Method, the method comprising the steps of:

S1, generation atmospheric optical parameters look-up table: the multiple aerosol moulds being suitable for according to China's zones of different and under season Formula and Atmospheric models and gf1 satellite width cover that ccd data is red, the spectral response functions of blue wave band, are simulated with atmospheric radiative transfer Method calculates different aerosol optical depths, different observation geometrical condition and gf1 satellite width and covers red, the blue wave band of ccd data Corresponding atmospheric optical parameters, form multidimensional lookup table；

S2, generate Reflectivity for Growing Season prior data bank: based on other satellites obtain cover the whole nation, update in time big Region Reflectivity for Growing Season product, calculates red, blue wave band priori earth's surface luminance factor value, forms Reflectivity for Growing Season priori data Storehouse；

S3, the dark goal pels of the suitable inverting of selection: cover the radiation calibration coefficient meter of ccd data based on gf1 satellite width Red, blue wave band apparent reflectance；Based on cloud pixel criterion remove cloud pixel, based on water body pixel criterion remove water body pixel, Extract the dark goal pels of suitable aerosol optical depth inverting based on dark goal pels criterion；

S4, the dark goal pels of calculating observe geometrical condition: the auxiliary positioning data of ccd data is covered based on gf1 satellite width, Calculate the observation geometrical condition of each dark goal pels；

S5, interpolation obtain the atmospheric optical parameters of dark goal pels: the gf1 satellite width according to treating inverting covers ccd data The transit time of image and region, determine aerosol model and the Atmospheric models of coupling；Then it is directed to each dark goal pels, Show that different aerosol optical depths are corresponding in conjunction with its wave band, several how index slip part interpolation from look-up table described in s1 of observation Atmospheric optical parameters；

S6, be based on radiation transmission inverting aerosol optical depth: for each dark goal pels described in s3, by described in s3 Described in red, blue wave band apparent reflectance, s4, observation geometrical condition and atmospheric optical parameters described in s5 bring radiation transmission side into Journey, calculates red, blue wave band Reflectivity for Growing Season ratio under different aerosol optical depths；It is thick that this pixel aerosol optical is set up in matching The empirical relation of red, the blue wave band Reflectivity for Growing Season ratio of degree and Radiance transfer calculation；Extract described from data base described in s2 Described red, the blue wave band priori earth's surface luminance factor value of dark goal pels, obtains described dark target according to described empirical relation interpolation Pixel aerosol optical depth now.

In the above-mentioned methods it is preferable that described gf1 satellite width cover ccd data be through geometric approximate correction blue, green, Red, four wave bands of near-infrared, its assistance data includes: radiation calibration coefficient, transit time, space orientation and observation several are how believed Breath.

In the above-mentioned methods it is preferable that described multidimensional lookup table comprises aerosol model × Atmospheric models × observation wave band The various dimensions index condition of × observation geometry × optical thickness, and the corresponding atmospheric optical parameters of each group index condition, and deposit Storage is in the data file.

Preferably, the described big region Reflectivity for Growing Season product updating in time updates interlude every time less than 30 My god；Wave band involved by the Reflectivity for Growing Season product of described big region comprises red, the blue wave band of gf1 satellite, or with described gf1 satellite Red, blue wave band is close.Preferably, described Reflectivity for Growing Season prior data bank comprises the grid number with common projection coordinate According to described packet contains red, the blue wave band Reflectivity for Growing Season ratio of each grid.

For improving data processing quality, another program of this method comprises satellite data quality examination and pretreatment further Step, including, for gf1 satellite width cover ccd data, check satellite image disappearance or band impact situation, in image with Machine noise profile and degree, and whether have obvious geometric distortion；Based on the quality of data identify remove shortage of data, band and Other noise pel data.The method of the present invention is applied to the data characteristicses that gf1 satellite width covers ccd, runs efficient and inverting Effect preferably, is advantageously implemented the inverting of high-resolution region atmospheric aerosol optical property, and monitoring for Particulate Pollution provides Beneficial data source.

Brief description

Fig. 1 is the data processing stream of the land optical properties of aerosol inversion method covering ccd data based on gf1 satellite width Cheng Tu

Specific embodiment

The optical properties of aerosol inversion method covering ccd data based on gf1 satellite width proposed by the present invention, in conjunction with accompanying drawing And embodiment detailed description is as follows.

Step s1, generation atmospheric optical parameters look-up table: according to Various Seasonal, region, special based on aerosol physical chemistry Property etc. priori, determine and be applied to multiple aerosol models and the Atmospheric models of China；By atmospheric radiative transfer software mould Intend the different aerosol optical depths of calculating, different observation geometrical condition and gf-1 satellite width and cover red, the blue wave band pair of ccd data The atmospheric optical parameters answered, are stored in formation multidimensional lookup table in data file.Step s1 specifically includes:

S11, according to China zones of different aerosol source, atmospheric condition and its seasonal variations, in conjunction with aerosol physics, light Learn the prioris such as property, delimit China's zones of different and suitable aerosol model under season and Atmospheric models；Wherein, gas is molten Rubber moulding formula includes: continent type, ocean type, urban type, coal smoke type, sand and dust type, mixed type etc.；Atmospheric models include tropical, middle latitude Aestivate season, middle latitude winter, high latitude summer, high latitude winter etc.；Ccd data is covered according to gf1 satellite width red, blue wave band Spectral response functions, setting different observation geometrical conditions in the range of rational codomain (is view zenith angle, too in the present embodiment Positive zenith angle and relative bearing) and different aerosol optical depths, such as: 1) wave band that setting look-up table calculates is Gf1 satellite width covers red, the blue wave band of ccd；2) set different observation conditions: select 10 solar zenith angles (0 °～90 °, 10 ° of interval), 10 view zenith angle (0 °～90 °, be spaced 10 °), 19 relative bearings (0 °～180 °, be spaced 10 °)；3) Set different aerosol optical depths: with respect to the aerosol optical depth at 0.55 mum wavelength be set to 10 grades (0, 0.25th, 0.5,0.75,1,1.25,1.5,1.75,2.5 and 3)；

S12, from suitable radiative transmission mode (as second simulation of the satellite Signal in the solar spectrum, 6s), the index condition that parameters all kinds of in s11 and s12 are formed is simulated fortune Calculate, obtain the corresponding atmospheric optical parameters of each group index condition and (in the present embodiment, radiate ρ for air path₀, air flood double Journey transmitance t, three parameters of downward atmospheric hemispherical reflectance s), and store formation look-up table in the data file.

Step s2, generation Reflectivity for Growing Season prior data bank: produced based on the big region Reflectivity for Growing Season that other satellites obtain Product calculate and can cover the whole nation, the time updates higher red, the blue wave band Reflectivity for Growing Season ratio prior data bank of the frequency.First First, from suitable satellite Reflectivity for Growing Season product to meet: 1) completely can cover China, 2) update the frequency higher (if do not surpassed Spend 30 days), 3) comprise red, the blue wave band of gf-1 satellite or close with it (the Reflectivity for Growing Season product mod09 as from modis)； Then, the spatial orientation information according to selected satellite Reflectivity for Growing Season product, it is processed again and has common projection coordinate (such as Deng longitude and latitude projection) grid data (grid size be 0.01 degree)；Calculate red, the blue wave band priori earth surface reflection of each grid Rate ratio, forms data base.

Step s3, the dark goal pels of the suitable inverting of selection: ccd data is covered to gf-1 satellite width and carries out following quality inspection Look into and pretreatment；Extract the dark goal pels of suitable aerosol optical depth inverting based on empirical value.Step s3 is specifically wrapped Include:

S31, first, carries out quality examination to satellite image, if meeting: the 1) pixel of satellite data disappearance or band impact Number is less than the 20% of whole scape image, 2) no significant random noise in data, 3) data no significantly geometric distortion, then retain This scape data enters follow-up inverting link；Then, cover the assistance datas such as the quality identification of ccd data based on gf1 satellite width, go Except shortage of data or serious pixel affected by noise；Finally, each wave band is calculated based on the radiation calibration coefficient in assistance data Apparent reflectance: read radiation calibration coefficient g and l from assistance data₀, remote sensing image picture element brightness dn value is converted to apparent Spoke brightness l:

L=dn/g+l₀

Then described apparent spoke brightness l is converted into the apparent reflectance of this wave band:

$\mathrm{\ρ}=\frac{\mathrm{\π}l}{e_{\mathrm{\λ}}{\mathrm{cos\θ}}_s}$

Wherein, e_λFor the solar irradiance on atmosphere top, can be obtained by inquiring about solar constant；θ_sFor solar zenith angle, can Cover the assistance data of ccd from gf1 satellite width and obtain.

S32, according to specific empirical value, removed based on the multiband apparent reflectance data that gf1 satellite width covers ccd Cloud covers pixel, and cloud pixel criterion used is:

ρ_red< 0.25 and ρ_nir<0.2

Wherein, ρ_redAnd ρ_nirRepresent that gf1 satellite width covers that ccd is red, near-infrared two wave band apparent reflectance respectively；

Remove water body pixel, water body pixel criterion used is:

ρ_nir＜ ρ_red

The pixel retaining is considered as effective land pixel, will enter follow-up screening link.

S33, filter out the dark goal pels in land being applied to aerosol optical depth inverting based on dark goal pels criterion. As embodiment, using the normalized differential vegetation index (normalized sensitive to the dark goal pels such as dense vegetation Difference vegetation index, ndvi) carry out given threshold:

$ndvi=\frac{{\mathrm{\&rho;}}_{nir}-{\mathrm{\&rho;}}_{red}}{{\mathrm{\&rho;}}_{nir}+{\mathrm{\&rho;}}_{red}}$

Filter out the criterion of the dark goal pels in land being applied to aerosol optical depth inverting: ndvi and be more than 0.3.

Step s4, the dark goal pels of calculating observe geometrical condition: the auxiliary positioning data of ccd is covered based on gf1 satellite width, Calculate the observation geometrical condition of each dark goal pels；Step s4 specifically includes:

Obtain the longitude and latitude of each dark goal pels from assistance data, calculate solar zenith angle and the sun on this basis Azimuth.As solar zenith angle θ_sCan calculate according to following formula:

cosθ_s=sin (lat) sin δ+cos (lat) cos δ cost

Solar azimuthCan calculate according to following formula:

Wherein, lat is the latitude of this dark goal pels, and δ is the angle of sunlight and earth equatorial plane, and t is the sun Hour angle, local time defining latitude, 12 points of time is 0, is-pi/2 when 6 points, is pi/2 when 18 points；

Gf1 satellite width covers the observed azimuth in the assistance data of ccdAnd solar azimuth calculated above Relative bearing both calculating:

Step s5, interpolation obtain the atmospheric optical parameters of dark goal pels: the gf1 satellite width according to treating inverting covers ccd The transit time of image and region, determine aerosol model and the Atmospheric models of coupling；Then it is directed to each dark goal pels, Show that different aerosol optical depths are corresponding in conjunction with its wave band, several how index slip part interpolation from look-up table described in s3 of observation Atmospheric optical parameters.Step s5 specifically includes:

Ccd data cover region and transit time are covered according to place gf1 satellite width, in conjunction with the meteorological letter such as back trajectca-rles Breath, determines aerosol model and the Atmospheric models of suitable inverting；Then it is directed to each dark goal pels, according to wave band, the sun It zenith angle, view zenith angle and relative bearing, from look-up table described in s13, interpolation extracts each different aerosol The corresponding atmospheric optical parameters of optical thickness are (as the present embodiment includes air path radiation, air flood round trip transmitance and big The descending hemispherical reflectance of gas).

Step s6, be based on radiation transmission inverting aerosol optical depth: for each dark goal pels, by red described in s3, Blue wave band apparent reflectance and atmospheric optical parameters described in s5 bring radiation transfer equation into, calculate different aerosol optical depths Under red, blue wave band Reflectivity for Growing Season ratio；Extract this dark goal pels red, blue wave band priori earth's surface anti-from data base described in s2 Penetrate rate ratio, matching obtains this dark goal pels aerosol optical depth now.Step s6 specifically includes:

S61, it is directed to each dark goal pels, big according to obtain in blue, red wave band apparent reflectance in s32 and s5 Gas optical parametric, under plane-parallel atmosphere is assumed, calculates this pixel according to following formula blue, red under different aerosol optical depths The Reflectivity for Growing Season of wave band:

$\mathrm{\&rho;}={\mathrm{\&rho;}}_0({\mathrm{\&mu;}}_s,{\mathrm{\&mu;}}_v,\mathrm{\&delta;}\mathrm{\&phi;})+\frac{t\left({\mathrm{\&mu;}}_s\right)t\left({\mathrm{\&mu;}}_v\right)r}{\&lsqb;1-rs\&rsqb;}$

Wherein, μ_s=cos θ_s, μ_v=cos θ_v, θ_sWith θ_vIt is respectively solar zenith angle and view zenith angle, δ φ is phase Azimuthal；ρ is apparent reflectance, and r is this wave band Reflectivity for Growing Season, and s is the downward hemispherical reflectance of air, and t is that air is whole Layer round trip transmitance, ρ₀Path radiation term equivalent reflectivity for air；Calculate this pixel thick in different aerosol optical The Reflectivity for Growing Season of red under degree, blue wave band, and then try to achieve the two ratio:

k_i=r_red,i/r_blue,i

Wherein, r_red,i、r_blue,iIt is respectively i-th aerosol optical depth τ_iCorresponding red, blue wave band Reflectivity for Growing Season, k_iFor the two ratio.

S62, matching set up this pixel τ_iWith k_iEmpirical relation τ=f (k), selects this pixel institute from data base described in s22 Red, blue wave band priori earth's surface luminance factor value k in grid ', draw this pixel this moment in conjunction with this fitting empirical relation interpolation Opticalthicknessτ '.

Wherein, described gf1 satellite wide covering ccd data is blue, green, the red, near-infrared etc. four through geometric approximate correction Wave band, is stored as the image file of tiff form；Its assistance data includes: radiation calibration coefficient, transit time, space orientation and Observe several how information, be stored in xml document corresponding with each wave band image file.

Embodiment of above is merely to illustrate the present invention, and not limitation of the present invention, common about technical field Technical staff, without departing from the spirit and scope of the present invention, can also make a variety of changes and modification, therefore all Equivalent technical scheme falls within scope of the invention, and the scope of patent protection of the present invention should be defined by the claims.

Claims (10)

1. a kind of land optical properties of aerosol inversion method based on a number satellite of high score is it is characterised in that comprise following step Rapid:

Ccd number is covered according to China's zones of different and the multiple aerosol models being suitable under season and Atmospheric models, gf1 satellite width According to the spectral response functions of red, blue wave band, calculate different aerosol optical depths, different sight with atmospheric radiative transfer analogy method Survey geometrical condition and gf1 satellite width covers the corresponding atmospheric optical parameters of red, the blue wave band of ccd data, form multidimensional lookup table；

The big region Reflectivity for Growing Season product covering the whole nation, updating in time being obtained based on other satellites, calculates red, blue Wave band priori earth's surface luminance factor value, forms Reflectivity for Growing Season prior data bank；

Red, blue wave band apparent reflectance is calculated based on the radiation calibration coefficient that gf1 satellite width covers ccd data；Based on cloud pixel Criterion remove cloud pixel, remove water body pixel based on water body pixel criterion, to extract suitable gas based on dark goal pels criterion molten The dark goal pels of glue optical thickness inverting；

Cover the auxiliary positioning data of ccd data based on gf1 satellite width, calculate the observation geometrical condition of described dark goal pels；

Gf1 satellite width according to treating inverting covers transit time and the region of ccd image data, determines the aerosol model of coupling And Atmospheric models；Then it is directed to described dark goal pels, in conjunction with its wave band, observe several how index slip parts from described look-up table Interpolation draws the corresponding atmospheric optical parameters of different aerosol optical depths；

For described dark goal pels, by described red, blue wave band apparent reflectance, described observation geometrical condition and described big Gas optical parametric brings radiation transfer equation into, calculates red, blue wave band Reflectivity for Growing Season ratio under different aerosol optical depths；Intend Build the empirical relation of this pixel aerosol optical depth vertical and red, the blue wave band Reflectivity for Growing Season ratio of Radiance transfer calculation jointly； Described red, the blue wave band priori earth's surface luminance factor of described dark goal pels is extracted from described Reflectivity for Growing Season prior data bank Value, obtains the aerosol optical depth of described dark goal pels according to described empirical relation interpolation.

2. as claimed in claim 1 the land optical properties of aerosol inversion method based on a number satellite of high score it is characterised in that Described gf1 satellite wide covering ccd data is blue, green, red, four wave bands of near-infrared through geometric approximate correction, its assistance data Including: radiation calibration coefficient, transit time, space orientation and observation geological information.

3. as claimed in claim 1 the land optical properties of aerosol inversion method based on a number satellite of high score it is characterised in that Described multidimensional lookup table comprises the various dimensions of aerosol model × Atmospheric models × observation wave band × observation geometry × optical thickness Index condition, and the corresponding atmospheric optical parameters of each group index condition, and store in the data file.

4. as claimed in claim 1 the land optical properties of aerosol inversion method based on a number satellite of high score it is characterised in that The described big region Reflectivity for Growing Season product updating in time updates interlude every time and is less than 30 days；Described big region ground Wave band involved by table reflectance product comprises red, the blue wave band of gf1 satellite, or close with red, the blue wave band of described gf1 satellite.

5. as claimed in claim 1 the land optical properties of aerosol inversion method based on a number satellite of high score it is characterised in that

Described Reflectivity for Growing Season prior data bank comprises the grid data with common projection coordinate, and described packet contains each Red, the blue wave band Reflectivity for Growing Season ratio of grid；

Described common projection coordinate is to wait longitude and latitude projection, and described grid size is 0.01 degree.

6. as claimed in claim 1 the land optical properties of aerosol inversion method based on a number satellite of high score it is characterised in that The step also comprising satellite data quality examination and pretreatment: cover ccd data for gf1 satellite width, check that satellite image lacks Lose or band impact situation, random noise distribution and degree in image, and whether have obvious geometric distortion；Based on data matter Amount mark removes shortage of data, band and other noise pel data.

7. as claimed in claim 1 the land optical properties of aerosol inversion method based on a number satellite of high score it is characterised in that Described aerosol model includes: continent type, ocean type, urban type, coal smoke type, sand and dust type, mixed type；Described Atmospheric models include The torrid zone, middle latitude summer, middle latitude winter, high latitude summer, high latitude winter.

8. as claimed in claim 1 the land optical properties of aerosol inversion method based on a number satellite of high score it is characterised in that The condition of described observation geometry is: selected 10 solar zenith angles, 10 view zenith angle, 19 relative bearings, with respect to Aerosol optical depth at 0.55 mum wavelength is set to 10 grades.

9. as claimed in claim 1 the land optical properties of aerosol inversion method based on a number satellite of high score it is characterised in that Described big region Reflectivity for Growing Season product, is the Reflectivity for Growing Season product mod09 of modis.

10. the land optical properties of aerosol inversion method based on a number satellite of high score as claimed in claim 1, its feature exists In,

Described cloud pixel criterion is: ρ_red< 0.25 and ρ_nir< 0.2, wherein, ρ_redAnd ρ_nirRepresent that gf1 satellite width covers ccd respectively Red, near-infrared two wave band apparent reflectance；

Described water body pixel criterion is: ρ_nir＜ ρ_red, wherein, ρ_redAnd ρ_nirRepresent that gf1 satellite width covers ccd respectively red, closely red Outer two wave band apparent reflectances；

Described dark goal pels criterion is: ndvi is more than 0.3, and wherein, ndvi is normalized differential vegetation index.