CN104880701A - Satellite-borne sensor imaging simulation method and device - Google Patents

Abstract

The embodiment of the invention discloses a satellite-borne sensor imaging simulation method and a satellite-borne sensor imaging simulation device. The satellite-borne sensor imaging simulation method comprises the steps of: acquiring apparent radiance parameters of the current atmosphere according to atmospheric path radiation, cross radiation effect weight of adjacent pixels and ground target reflectivity images; and acquiring entrance pupil radiance field images of a satellite-borne sensor according to an atmospheric point spread function and the apparent radiance parameters of the current atmosphere. The satellite-borne sensor imaging simulation method and the satellite-borne sensor imaging simulation device can improve the simulation precision of atmospheric radiation transmission.

Description

A kind of satellite borne sensor imaging simulation method and device

Technical field

The present invention relates to remote sensing satellite image technical field, be specifically related to a kind of satellite borne sensor imaging simulation method and device.

Background technology

Through country's giving more sustained attention and drop into emphasis in recent years; satellite remote sensing development rapidly; by the end of the year 2013; China has launched and remote sensing satellite in-orbit (is planned for the year two thousand twenty more than more than 40 and will reaches more than 100); China has begun to take shape the remote sensing land observation system of scale; achieve satellite remote sensing from the concept stage of popularization to the technology growth stage, then arrive the process of the industrialization stage of growth.But compared with the World Developed Countries such as the U.S., also there is very large gap from the aspect such as development scale, technical merit in China's remote sensing application and industrialization.This situation major reason is caused to be remote sensing full industrial chain integration scientific appraisal with Top-layer Design Method is not enough, remote sensing satellite develop and remote sensing application disconnection, remote sensing system comprehensive application ability is not enough, can not meet China's remote sensing industryization fast-developing to remote sensing application technical research and the new demand developing proposition.

In order to improve the image quality of China's High Resolution Remote Sensing Satellites, an important approach by the satellite imagery emulation before reflection, simultaneously to imaging simulation and the quality analysis of satellite in orbit, for the design of subsequent transmission satellite provides technical support, it is also a kind of practicable method.The research and development of the full link simulation platform of domestic remotely sensed image are started late, but the total body of the existing space flight Fifth Academy of Aerospace Industry, institute of the space flight Fifth Academy of Aerospace Industry 503 (comprising space flight sky to paint), institute of the space flight Fifth Academy of Aerospace Industry 508, University of Science and Technology, Central China, northern science and engineering, Harbin Institute of Technology, Bureau of Surveying and Mapping, remote sensing Suo Dengduojia unit of the Chinese Academy of Sciences independently have carried out the work of being correlated with in platform core texture and module composition.On the whole, by the complementation of constituent parts, the domestic research and development ability possessed each subsystem/module of platform, but there is no a unit and can develop complete remote-sensing flatform.

Atmospheric radiative transfer emulation is the important component part of the full link simulation platform of remotely sensed image, is also the link that influence factor is maximum, and simulation accuracy is difficult to ensure.

Summary of the invention

In view of this, the embodiment of the present invention provides a kind of satellite borne sensor imaging simulation method and device, to realize the emulation of high precision atmospheric radiative transfer.

The embodiment of the present invention is by the following technical solutions:

First aspect, embodiments provides a kind of satellite borne sensor imaging simulation method, comprising:

According to atmospheric path radiation, the cross radiance weighing factor of contiguous pixel, and terrain object albedo image obtains the apparent spoke luminance parameter of current atmospheric;

The entrance pupil spoke luminance field image of satellite borne sensor is obtained according to aerosol optical depth and the apparent spoke luminance parameter of described current atmospheric.

Second aspect, the embodiment of the present invention additionally provides a kind of satellite borne sensor imaging simulation device, comprising:

Current atmospheric apparent spoke luminance parameter acquiring unit, for according to atmospheric path radiation, be close to the cross radiance weighing factor of pixel, and terrain object albedo image obtains the apparent spoke luminance parameter of current atmospheric;

Entrance pupil spoke luminance field image acquisition unit, for obtaining the entrance pupil spoke luminance field image of satellite borne sensor according to aerosol optical depth and the apparent spoke luminance parameter of described current atmospheric.

The embodiment of the present invention is passed through, according to atmospheric path radiation, the cross radiance weighing factor of contiguous pixel, and terrain object albedo image obtains the apparent spoke luminance parameter of current atmospheric, the entrance pupil spoke luminance field image of satellite borne sensor is obtained according to aerosol optical depth and the apparent spoke luminance parameter of described current atmospheric, consider the cross radiance tripartite face sound of terrain object reflection, atmospheric path radiation, adjacent objects, ensure atmospheric radiative transfer simulation accuracy.

Accompanying drawing explanation

In order to be illustrated more clearly in the technical scheme in the embodiment of the present invention, below the accompanying drawing used required in describing the embodiment of the present invention is briefly described, apparently, accompanying drawing in the following describes is only some embodiments of the present invention, for those of ordinary skill in the art, under the prerequisite not paying creative work, other accompanying drawing can also be obtained according to the content of the embodiment of the present invention and these accompanying drawings.

Fig. 1 is the satellite borne sensor imaging simulation method flow diagram described in the specific embodiment of the invention one;

Fig. 2 is the satellite borne sensor imaging simulation method flow diagram described in the specific embodiment of the invention two;

Fig. 3 is the satellite borne sensor imaging simulation method flow diagram described in the specific embodiment of the invention three;

Fig. 4 is the satellite borne sensor imaging simulation method flow diagram described in the specific embodiment of the invention four;

Fig. 5 is the structured flowchart of the satellite borne sensor imaging simulation device described in the specific embodiment of the invention five;

Fig. 6 is satellite zenith angle geometric relationship schematic diagram;

Fig. 7 is satellite aximuth geometric relationship schematic diagram.

Embodiment

The technical matters solved for making the present invention, the technical scheme of employing and the technique effect that reaches are clearly, be described in further detail below in conjunction with the technical scheme of accompanying drawing to the embodiment of the present invention, obviously, described embodiment is only the present invention's part embodiment, instead of whole embodiments.Based on the embodiment in the present invention, those skilled in the art, not making the every other embodiment obtained under creative work prerequisite, belong to the scope of protection of the invention.

Technical scheme of the present invention is further illustrated by embodiment below in conjunction with accompanying drawing.

Embodiment one

Fig. 1 is the satellite borne sensor imaging simulation method flow diagram described in the present embodiment.As shown in Figure 1, the satellite borne sensor imaging simulation method described in the present embodiment comprises:

Step S101, according to atmospheric path radiation, the cross radiance weighing factor of contiguous pixel, and terrain object albedo image obtains the apparent spoke luminance parameter of current atmospheric.

The radiant quantity of ground scene arrives sensor entrance pupil place through atmospheric radiative transfer, the main contribution that there is following a few part: is to sun direct projection and Diffuse Light Reflectance effect from earth's surface target; Two be air to the scattering of sunshine and self thermal-radiating contribution, i.e. journey radiation; Three is impacts of bias light and diffused light, the cross radiance of ground close region.Therefore, when calculating above-mentioned three kinds of factors are all taken into account can obtain more comprehensively, objective result.

For the earth's surface of lambert's body reflection, satellite load camera observes that the brightness of entrance pupil spoke can represent with the publicity simplified below:

$L=G_t\frac{\mathrm{\ρ}}{1-{\mathrm{\ρ}}_{e}S}+G_b\frac{{\mathrm{\ρ}}_e}{1-{\mathrm{\ρ}}_{e}S}+L_0$

In formula, ρ is goal pels surface reflectivity; ρ _efor close region reflectivity; S is air spherical albedo; G _tfor the gain of goal pels; G _bfor the gain of backdrop pels; L0 is the journey radiation of air; L is the spoke brightness value that sensor entrance pupil place receives.

In above formula, ρ and ρ _ethe amount relevant to ground surface reflectance, and G _ag _bthese four parameters of S, L0 are determined by the sun, satellite observation direction and atmospheric conditions during imaging, therefore, when carrying out entrance pupil spoke brightness simulation calculation to a width albedo image, are determining G _a, G _b, look-up table between Reflectivity for Growing Season and the brightness of entrance pupil spoke can be set up after S, L0 fast.

In order to obtain G _a, G _b, S, L0 tetra-parameters value, can be that radiative transfer model calculating is carried out on 0%, 50% and 100% earth's surface to reflectivity respectively by MODTRAN software, and utilize the result of calculation on these three kinds of reflectivity earth's surfaces, solve G _a, G _b, S, L0 tetra-parameters, method for solving is shown below.

$S=\frac{{\mathrm{\ΔTOT}}_{100}-2{\mathrm{\ΔTOT}}_{50}}{\mathrm{\Δ}{\mathrm{TOT}}_{100}-{\mathrm{\ΔTOT}}_{50}}$

G _b＝ΔPATH ₁₀₀(1-S)

G _t＝ΔGRFL ₁₀₀(1-S)

L ₀＝PATH ₀

In above formula, Δ TOT ₁₀₀=TOT ₁₀₀-TOT ₀=TOT ₁₀₀-PATH ₀, Δ PATH ₁₀₀=PATH ₁₀₀-PATH ₀, Δ TOT ₅₀=TOT ₅₀-TOT ₀=TOT ₅₀-PATH ₀, Δ GRFL ₁₀₀=GRFL ₁₀₀-GRFL ₀=GRFL ₁₀₀.

Utilize PATH, GRFL and TOT tri-amount in MODTRAN software result of calculation, measure respectively for these three: PATH represents air path enters viewing field of camera journey radiation to sunshine scattering; GRFL enters the radiant quantity of viewing field of camera after representing earth surface reflection sunshine; TOT represents whole radiant quantity that camera entrance pupil place receives.

S102, obtain the entrance pupil spoke luminance field image of satellite borne sensor according to aerosol optical depth and the apparent spoke luminance parameter of described current atmospheric.

According to response wave band scope and the wavelength interval of satellite sensor, determine the number of wavelengths that will emulate, successively simulation calculation is carried out to the atmospheric radiative transfer process of each wavelength.First, read the aerosol optical depth on first wavelength, and determine the cross radiance weighing factor of contiguous pixel; Next, read terrain object albedo image, based on the computation model built in previous step, calculate the spoke brightness value of each pixel target after propagation in atmosphere in first wave length respectively; Repeat above-mentioned steps, the apparent spoke luminance picture successively in calculated response wave band on each wavelength.

The present embodiment passes through according to atmospheric path radiation, the cross radiance weighing factor of contiguous pixel, and terrain object albedo image obtains the apparent spoke luminance parameter of current atmospheric, the entrance pupil spoke luminance field image of satellite borne sensor is obtained according to aerosol optical depth and the apparent spoke luminance parameter of described current atmospheric, consider the cross radiance tripartite face sound of terrain object reflection, atmospheric path radiation, adjacent objects, ensure atmospheric radiative transfer simulation accuracy.

Embodiment two

Fig. 2 is the satellite borne sensor imaging simulation method flow diagram described in the present embodiment, and as shown in Figure 2, the satellite borne sensor imaging simulation method described in the present embodiment comprises:

Step S201, atmospheric path radiation according to history, the cross radiance weighing factor of contiguous pixel, and terrain object albedo image sets up air apparent spoke brightness look-up table.

For apparent spoke luminance parameter, mode adopts the computing formula shown in embodiment one to calculate in real time more intuitively, but the speed calculated in real time is due to hardware environment and data volume is different very big difference.In any case be optimized its algorithm, utilizing the mass data obtained to calculate apparent spoke luminance parameter in real time when emulating is all an operation more consuming time.In order to get apparent spoke luminance parameter at faster speed, the present embodiment adopts the mode of look-up table to replace real-time calculating.Because search operation only needs Data Matching simply, and without the need to a large amount of calculating, therefore fast doubly a lot of than calculating in real time.Adopt the mode of look-up table to obtain required apparent spoke luminance parameter in the present embodiment, specifically, the computing formula shown in embodiment one can be adopted, according to the atmospheric path radiation of history, the cross radiance weighing factor of contiguous pixel, and terrain object albedo image sets up air apparent spoke brightness look-up table.Wherein, described air apparent spoke brightness look-up table comprises the apparent spoke luminance parameter of air under history observation condition and this observation condition; Described observation condition comprises following at least one item: aerosol type, atmospherical model, month, solar zenith angle, view zenith angle and wavelength.These observation conditions reflect the pacing itemss such as atmospheric conditions, time, orientation substantially.

Step S202, acquisition Current observation condition.

Step S203, according to Current observation conditional search air apparent spoke brightness look-up table, obtain current atmospheric apparent spoke luminance parameter.

Obtain current observation condition, then in the apparent spoke brightness look-up table previously set up, the apparent spoke luminance parameter of the air corresponding with this observation condition can be gone out by Auto-matching according to this observation condition.Only need to provide Current observation condition, just can draw the spoke luminance parameter matched immediately, succinctly rapidly.

Step S204, obtain the entrance pupil spoke luminance field image of satellite borne sensor according to aerosol optical depth and the apparent spoke luminance parameter of described current atmospheric.

According to response wave band scope and the wavelength interval of satellite sensor, determine the number of wavelengths that will emulate, successively simulation calculation is carried out to the atmospheric radiative transfer process of each wavelength.First, read the aerosol optical depth on first wavelength, and determine the cross radiance weighing factor of contiguous pixel; Next, read terrain object albedo image, based on the computation model of a kind of structure of embodiment, calculate the spoke brightness value of each pixel target after propagation in atmosphere in first wave length respectively; Repeat above-mentioned steps, the apparent spoke luminance picture successively in calculated response wave band on each wavelength.

The present embodiment establishes sextuple air apparent spoke brightness look-up table, and the imaging simulation that can be applicable to global different time calculates.Because search operation is than calculating fast several order of magnitude in real time, therefore, the apparent spoke luminance parameter that the present embodiment is conformed to Current observation condition by Auto-matching from look-up table, significantly can improve simulation velocity.

Embodiment three

Fig. 3 is the satellite borne sensor imaging simulation method flow diagram described in the present embodiment, and as shown in Figure 3, the satellite borne sensor imaging simulation method described in the present embodiment comprises:

Step S301, the radiation of acquisition current atmospheric journey, obtain the cross radiance weighing factor of current contiguous pixel, and obtain Current terrestrial target reflectivity image.

According to the analysis in embodiment one, atmospheric path radiation, the contiguous cross radiance of pixel and the reflection of terrain object are the principal elements affecting the brightness of satellite borne sensor entrance pupil spoke, therefore, in order to carry out objective comprehensive emulation to the imaging of satellite borne sensor, the cross radiance weighing factor of the radiation of the present embodiment first-selection acquisition current atmospheric journey, contiguous pixel and terrain object albedo image.Wherein, the operation obtaining the radiation of current atmospheric journey specifically comprises: obtain aerosol optical depth according to current atmospheric optical property parameter, then obtains the radiation of current atmospheric journey according to aerosol optical depth.The operation obtaining the cross radiance weighing factor of current contiguous pixel specifically comprises: according to Determination of Aerosol Optical, Rayleigh scattering optics thickness, and absorbs other optical thickness, obtains the cross radiance weighing factor of contiguous pixel.

Step S302, according to the radiation of described current atmospheric journey, the cross radiance weighing factor of described current contiguous pixel, and the apparent spoke luminance parameter of described Current terrestrial target reflectivity image Real-time Obtaining current atmospheric.

The apparent spoke luminance parameter of air can adopt two kinds of diverse ways to obtain, and mode adopts the computing formula shown in embodiment one to calculate in real time more intuitively, but its computing velocity is slow, inefficiency; Another kind of mode adopts being obtained by Auto-matching apparent spoke brightness look-up table shown in embodiment two.This mode response speed is very fast, preferably adopts in this way; But, the prerequisite using look-up table to obtain apparent spoke luminance parameter previously must establish available look-up table according to historical data, if when without any historical data, still need to adopt the mode calculated in real time to obtain apparent spoke brightness value.In order to improve the speed of the apparent spoke brightness value of follow-up acquisition, can by related data as observation condition, apparent spoke brightness value result of calculation etc. be preserved in a lookup table while calculating in real time, thus to look-up table constantly accumulated history data upgrading, like this after after a while, the data in look-up table just can be utilized to implement Auto-matching, thus replace real-time account form.

Step S303, obtain the entrance pupil spoke luminance field image of satellite borne sensor according to aerosol optical depth and the apparent spoke luminance parameter of described current atmospheric.

According to response wave band scope and the wavelength interval of satellite borne sensor, determine the number of wavelengths that will emulate, successively simulation calculation is carried out to the atmospheric radiative transfer process of each wavelength.First, read the aerosol optical depth on first wavelength, and determine the cross radiance weighing factor of contiguous pixel; Next, read terrain object albedo image, based on the computation model of a kind of structure of embodiment, calculate the spoke brightness value of each pixel target after propagation in atmosphere in first wave length respectively; Repeat above-mentioned steps, the apparent spoke luminance picture successively in calculated response wave band on each wavelength.

Step S304, spectral response curve according to satellite borne sensor, be weighted the spectral radiance field picture on each wavelength in described response wave band, obtains all band entrance pupil spoke luminance field image of equivalence.

After the entrance pupil radiance data having calculated each wavelength, consider the spectral response curve of satellite borne sensor, according to the transmitance coefficient of each wavelength, calculate full wave equivalent entrance pupil spoke luminance field image.Following formula specifically can be adopted to calculate:

$L=\frac{\mathrm{\Σ}{L}_i*R_i}{\mathrm{\Σ}{R}_i}$

In formula, L represents the brightness of full wave equivalent entrance pupil spoke, L _irepresent the entrance pupil spoke brightness value of i-th wavelength, R _irepresent the transmitance value of satellite borne sensor on this wavelength i.By can obtain the equivalent entrance pupil spoke luminance field image in all band with upper type.

The present embodiment adopts wave spectrum segmentation mode to calculate atmospheric radiative transfer process, and consider the spectral response curve of satellite borne sensor, complete and multispectrally to calculate with the imaging simulation of panchromatic sensor, can radiation transmission impact in the whole response wave band of complete reflection, ensure the simulation accuracy of satellite borne sensor particularly full-color camera.

Embodiment four

Fig. 4 is the satellite borne sensor imaging simulation method flow diagram described in the present embodiment, and as shown in Figure 4, the satellite borne sensor imaging simulation method described in the present embodiment comprises:

The auxiliary data such as attitude, track, time of step S401, acquisition satellite.

The geometrical orientation of the sun, target, satellite in step S402, computer sim-ulation scene.

(1) based on simulating scenes longitude and latitude and simulation imaging time, calculate that sunshine enters zenith angle and the position angle of institute's ground scene.

Calculate solar zenith angle and azimuthal method as follows:

(1.1) according to revolution of earth and rotation rule, the sun can be calculated by solar declination, solar hour angle and observation station longitude and latitude relative to the zenith angle at any time of certain longitude and latitude place on the earth and position angle.

(1.2) computing method of solar declination are:

θ _δ＝0.3723-0.758cos(Γ)+23.2567sin(Γ)+0.3656cos(2Γ)+0.1149sin(2Γ)+0.0201cos(3Γ)-0.1712sin(3Γ)

Wherein, Γ=2 π t/365.2422, t=N-N ₀-1, N is day of year.

N ₀＝78.801+0.2422*(YYYY-1969)-INT((YYYY-1969)/4)

Wherein, the fraction of the year that YYYY representing, the integer being not more than this number is got in INT () expression.

(1.3) computing method of solar hour angle are as follows:

The local true solar time is:

θ _sh=15 × (t _{loc true sun}-12)

Wherein, t _{loc true sun}for the local true solar time can be tried to achieve by the correction of Greenwich mean solar time and geographic longitude:

T _{loc true sun}=t _{gMT is mean sun}+ θ _lon/ 15+E _t

Wherein, t _{gMT is mean sun}for the GMT mean solar time, θ _lonfor geographic longitude, being just eastwards, being westwards negative, E _tfor the time difference of true solar time and mean solar time is revised, can be calculated by many algorithms.

Lamm proposed with 4 years for calibration cycle calculates the time difference in 1981, and its algorithm is:

$E_t=\underset{k=0}{\overset{5}{\mathrm{\Σ}}}(A_k\cos \frac{2\mathrm{\πkN}}{365.25}+B_k\sin \frac{2\mathrm{\πkN}}{365.25})$

Table 1 Shicha algorithm coefficient and reference value, wherein: N is day of year, coefficient A _kand B _kvalue can look into following table obtain.

k	A k	B k
			0	2.0870e-4	0
1	9.2869e-3	-1.2229e-1
			2	-5.2258e-2	-1.5698e-1
3	-1.3077e-3	-5.1602e-3
			4	-2.1867e-3	-2.9823e-3
5	-1.5100e-4	-2.3463e-4

(1.4) solar zenith angle θ _stried to achieve by following formula

cosθ _s＝sinθ _Latsinδ _Lat+cosθ _Latcosδ _Latcosθ _st

Wherein θ _latfor observation station geographic latitude, δ _latfor solar declination, θ _stfor solar hour angle.

(1.5) solar azimuth computing method are:

Solar azimuth declination δ can be passed through _lat, solar hour angle θ _stwith sun altitude θ _hbe expressed as:

Sun altitude equals 90 ° and deducts solar zenith angle.

Solar azimuth also can be expressed as by sun altitude, observation station geographic latitude and declination:

(2) based on simulating scenes longitude and latitude and sea level elevation, satellite altitude, moonscope side-looking angle, zenith angle and the position angle of satellite sensor observation ground scene is calculated.

(2.1) according to the longitude and latitude of moonscope point, height above sea level, sub-satellite point longitude and latitude and moonscope side-looking angle, zenith angle and the position angle of moonscope can be calculated.

(2.2) the zenith angle computing method of moonscope are as follows:

As shown in Figure 6, wherein O is the earth's core to the geometrical orientation schematic diagram that moonscope zenith angle calculates, and R is the radius of a ball fifty-fifty, and L is observation station, and S is satellite position, and λ is moonscope elevation angle, and α is satellite observation direction angle.

Can obtain moonscope zenith angle by geometric relationship in figure is:

${\mathrm{\θ}}_v={\sin }^{-1}(\frac{R+h}{R}\×\sin \mathrm{\α})$

Wherein, R is average half axial length of the earth, and h is satellite orbital altitude, and α is satellite side-looking angle, and observed altitude angle is the complementary angle of zenith angle.

(2.3) azimuth calculation method of moonscope is as follows:

As shown in Figure 7, wherein O is the earth's core to the geometrical orientation schematic diagram that satellite aximuth calculates, and L is observation station, and S is satellite position, and B is observation station section and OS intersection point, and C, D are all intersection points of observation station section and x-axis and z.

By satellite S point, observation station L point by longitude and latitude 2000 ordinate transform under WGS84 coordinate system, S (x _s, y _s, z _s), L (x _l, y _l, z _l),

$\∠\mathrm{BOL}=\frac{\mathrm{\π}}{2}-{\mathrm{\θ}}_v-\mathrm{\α}$

$\mathrm{OB}=\frac{R}{\cos \∠\mathrm{BOL}}$

Solved by OB and OS conllinear: $B(x_b,y_b,z_b)=\frac{R}{(R+h)\sin (\mathrm{\λ}+\mathrm{\α})}(a,b,c),$ Vector BL,

It is vertical with vectorial CL by vectorial OL, can obtain:

$D(\mathrm{0,0},z_d)=D(\mathrm{0,0},\frac{R^2}{z_l}),C(x_c,y_c,0)=C(\frac{R^2{x}_l}{x_l^2+y_l^2},\frac{R^2{y}_l}{x_l^2+y_l^2},0)$ Solve to obtain vectorial CL, $\mathrm{Cl}=(-\frac{x_l{z}_l^2}{x_l^2+y_l^2},\frac{y_l{z}_l^2}{x_l^2+y_l^2},0).$

solve to obtain observed azimuth:

According to side-looking direction, satellite position judges that position angle makes a decision scope.

Step S403, carry out ground data parameter processing, obtain atmospheric optics characteristic parameter.

(1) the solar irradiation degrees of data arrived based on ground observation, simulating scenes longitude and latitude and simulation imaging time, the optical thickness that air is total is solved.

Air total optical thickness algorithm characterizes as follows:

According to Bouguer-Lamber law, the sun direct projection irradiance on earth's surface is

E(λ)＝E ₀(λ)R ^-2exp[-mτ(λ)]T _g

In formula, E ₀for the sun direct projection irradiance outside earth atmosphere when day, distance was 1, R is the solar distance measuring the moment, and m is air quality number, and τ (λ) is vertical whole atmosphere optical thickness, T _gabsorb gas permeation rate.

When utilizing CE318 to calculate atmosphere optical thickness, first select to process without the wave band of gas absorption, locate beyond 936nm wave band steam can not ignore the absorption of solar radiation, on its all band of CE318, the absorption of gas can be ignored, and thinks T _g=1, taken the logarithm in equation both sides above,

ln[E(λ)R ²]＝lnE ₀(λ)-mτ(λ)

Because air quality number m is the variable with time correlation, at short notice, can think that atmospheric conditions are stable, therefore, in during this period of time, to the repeatedly observation of sun direct projection illumination, the method for linear fit can be adopted, solve atmosphere optical thickness.

(2) based on the meteorological measuring in simulating scenes, as air pressure etc., Rayleigh diffuse optical thickness in imaging region air is calculated.

The opticalthicknessτ of Rayleigh scattering _r(λ) more stable, it is relevant with air pressure, sea level elevation and wavelength, can by formulae discovery,

${\mathrm{\τ}}_r(\mathrm{\λ},z)=0.00879{\mathrm{\λ}}^{-4.09}\frac{P\left(0\right)-P\left(z\right)}{P_n\left(0\right)}$

In above formula, P (0) is sea level pressure, P _n(0) standard atmospheric pressure on sea level when being 15 DEG C, the atmospheric pressure at P (z) to be atmospheric envelope heights of roofs be z place, earth air is mainly distributed in flow process, and atmospheric pressure P (z) of China's mid latitudes is about 100hPa.

(3) based on the ozone coefficient of imaging region air, the Ozone Absorption optical thickness in air is calculated.

When determining aerosol optical depth, the optical thickness relevant with ozone molecule calculates according to parameterized method usually.Ozone Absorption opticalthicknessτ O ₃with Ozone Absorption COEFFICIENT K O ₃with ozone content UO ₃pass between (unit DU) is

τO ₃＝KO ₃*UO ₃

Wherein the acquisition methods of ozone content can obtain from ozone product http://toms.gsfc.nasa.gov/ozone/o-zone_v8.html, the http://es-ee.tor.ec.gc.ca/cgi-bin/totalozone/ of OMI (global ground ozone measurement point); Ozone Absorption coefficient can with reference to existing coefficient table or the relation being calculated ozone transmitance and ozone content by Modtran software: KO ₃=-(lnTO ₃/ UO ₃).

(4) total according to simulating area aerosol optical depth, Rayleigh scattering optics thickness, calculate the aerosol optical depth without absorption bands 440nm and 870nm, and according to gasoloid relation equation, the aerosol optical depth of inverting 550nm wave band.

(4.1) Computing Principle of aerosol optical depth is:

Air can be divided into the Rayleigh scattering of gas molecule, gas molecules sorb and gasoloid to decay three parts on the impact of beam radia, and namely air total optical thickness can be expressed as

τ(λ)＝τ _r(λ)+τ _g(λ)+τ _a(λ)

In formula, τ _r(λ) optical thickness of gas molecule Rayleigh scattering is represented, τ _g(λ) optical thickness of gas molecules sorb is represented, τ _a(λ) aerocolloidal optical thickness is represented.

On the wavelength ignoring gas molecules sorb, obtain aerocolloidal optical thickness by following formula,

τ _a(λ)＝τ(λ)-τ _r(λ)

(4.2) usual, Determination of Aerosol Optical and meteorological visual level distance is characterized with the aerosol optical depth of 550nm wave band, according to the aerosol optical depth without absorption bands 440nm and 870nm wave band, utilize the relational expression of aerosol optical depth and wavelength, the aerosol optical depth of 550nm wave band can be calculated, relational expression is as follows

τ _a(λ)＝βλ ^-α

In formula, β is atmospheric turbidity parameter, and theory shows, when particle radii are 0.01-10.0 μm, the value of β has nothing to do without wavelength substantially; α is the power exponent of Jungle spectrum, is also Angstrom index.

Logarithm is asked to above formula both sides,

ln[τ _a(λ)]＝lnβ-αlnλ

According to above formula, utilize multiple aerosol optical depth calculated without absorption bands, linear fit goes out parameter beta and α, then calculates the aerosol optical depth on its all band, as the aerosol optical depth on 550nm wave band.

Step S404, according to atmospheric optics characteristic parameter, calculate aerosol optical depth.

(1) based on Determination of Aerosol Optical, Rayleigh scattering optics thickness and absorption other optical thickness, the radiation effect of different distance target in cross radiance coverage is calculated.

The approximate model of aerosol optical depth is that the geometric relationship formed according to sensor, goal pels and contiguous pixel describes the radiation flux reaching sensor from the diffuse reflection being close to pixel again through the primary event of air from principle, thus reaches the analytic approximation solution of aerosol optical depth.

According to atmospheric radiative transfer principle, the radiance that satellite borne sensor entrance pupil place receives comes from following several:

L _t＝L _r+L _a+L _ra+L _g+t _dL _s

On the right of above formula, first three items represents respectively and comes from Rayleigh scattering, gasoloid, the interactional contribution of Rayleigh scattering-gasoloid; L _gfor the contribution of mirror-reflection, moonscope avoids mirror-reflection usually, therefore L _gcan ignore; L _sfor the contribution of ground return or emitted energy, t _dfor air is to L _stransmitance, i.e. diffuse transmission.Diffuse transmission can calculate according to following formula,

t _d＝exp[-(0.5τ _r+τ _a(1-ω _aF)+τ _o)/cos(θ)]

In above formula, τ _rrepresent the optical thickness of Rayleigh scattering, τ _arepresent aerocolloidal optical thickness, τ _orepresent the optical thickness of Ozone Absorption, ω _arepresent gasoloid single albedo, F represents the probability of photon forward direction outgoing after aerosol particles scattering.

When calculating aerosol optical depth, only need consider gasoloid and the interactional contribution of Rayleigh scattering-gasoloid, gasoloid can be expressed as follows the contribution of aerosol optical depth,

t _aer＝t _d*exp(-0.5*τ _r)-t _b

T in formula _bfor light direct projection transmitance,

t _b＝exp[-(τ _r+τ _a+τ _o)/cos(θ)]

And Rayleigh scattering-gasoloid interaction can be expressed as follows the contribution of aerosol optical depth,

t _rra＝t _d-t _b-t _aer

The coverage of air point diffusion is by t _aerand t _rraapproximate treatment obtains, and aerosol optical depth approximate solution is:

$\mathrm{PSF}=[1/(2\mathrm{\πR}\×{10}^{-6})]*\frac{d}{\mathrm{dR}}\left[\mathrm{cum}\left(R\right)\right]$

(2) according to the radiation effect intensity of being close to pixel target and observed object, normalized, extrapolates the cross radiance weighing factor of coverage internal object.

After calculating the relation between aerosol optical depth and distance, determine that air point extends influence scope, calculate PSF distribution situation around object pixel, adopt normalized, acquisition closes on pixel weighing factor, and be saved in text, to use during air entrance pupil spoke brightness simulation calculation.

Step S405, judge whether to calculate atmospheric radiative transfer in real time, if then perform step S406, otherwise perform step S407.

When the apparent spoke luminance parameter of calculating air, there is two schemes: one is according to atmospheric radiation transmission, calculate apparent spoke luminance parameter in real time; Another kind calls air apparent spoke brightness look-up table storehouse, obtains the apparent spoke luminance parameter under respective conditions.Analyze above, the difference of the two is not repeat them here the difference in speed.

Step S406, according to atmospheric radiation mode, in real time calculate apparent spoke luminance parameter, perform step S408.

When the apparent spoke luminance parameter of selection real-time calculating air, call MODTRAN software to calculate propagation in atmosphere process, respectively in simulating scenes reflectivity be 0,0.5 and the atmospheric radiative transfer process of 1.0 3 kind of reflectivity calculate, the apparent spoke luminance parameter needed for acquisition.

The radiant quantity of ground scene arrives sensor entrance pupil place through atmospheric radiative transfer, the main contribution considering following a few part: is to sun direct projection and Diffuse Light Reflectance effect from earth's surface target; Two be air to the scattering of sunshine and self thermal-radiating contribution, i.e. journey radiation; Three is impacts of bias light and diffused light, the cross radiance of ground close region.

For the earth's surface of lambert's body reflection, satellite load camera observes that the brightness of entrance pupil spoke can represent with the publicity simplified below:

$L=G_t\frac{\mathrm{\ρ}}{1-{\mathrm{\ρ}}_{e}S}+G_b\frac{{\mathrm{\ρ}}_e}{1-{\mathrm{\ρ}}_{e}S}+L_0$

In formula, ρ is goal pels surface reflectivity; ρ _efor close region reflectivity; S is air spherical albedo; G _tfor the gain of goal pels; G _bfor the gain of backdrop pels; L ₀for the journey radiation of air; L is the spoke brightness value that sensor entrance pupil place receives.

In above formula, ρ and ρ _ethe amount relevant to ground surface reflectance, and G _a, G _b, S, L ₀these four parameters are determined by the sun, satellite observation direction and atmospheric conditions during imaging, therefore, when carrying out entrance pupil spoke brightness simulation calculation to a width albedo image, are determining G _a, G _b, S, L ₀after can set up look-up table between Reflectivity for Growing Season and the brightness of entrance pupil spoke fast.

In order to obtain G _a, G _b, S, L ₀the value of four parameters, can be that radiative transfer model calculating is carried out on 0%, 50% and 100% earth's surface to reflectivity respectively by MODTRAN software, and utilize the result of calculation on these three kinds of reflectivity earth's surfaces, solve G _a, G _b, S, L ₀four parameters, method for solving is shown below.

$S=\frac{{\mathrm{\ΔTOT}}_{100}-2{\mathrm{\ΔTOT}}_{50}}{\mathrm{\Δ}{\mathrm{TOT}}_{100}-{\mathrm{\ΔTOT}}_{50}}$

G _b＝ΔPATH ₁₀₀(1-S)

G _t＝ΔGRFL ₁₀₀(1-S)

L ₀＝PATH ₀

In above formula, Δ TOT ₁₀₀=TOT ₁₀₀-TOT ₀=TOT ₁₀₀-PATH ₀, Δ PATH ₁₀₀=PATH ₁₀₀-PATH ₀, Δ TOT ₅₀=TOT ₅₀-TOT ₀=TOT ₅₀-PATH ₀, Δ GRFL ₁₀₀=GRFL ₁₀₀-GRFL ₀=GRFL ₁₀₀.

Utilize PATH, GRFL and TOT tri-amount in MODTRAN software result of calculation, measure respectively for these three: PATH represents air path enters viewing field of camera journey radiation to sunshine scattering; GRFL enters the radiant quantity of viewing field of camera after representing earth surface reflection sunshine; TOT represents whole radiant quantity that camera entrance pupil place receives.

Step S407, call air apparent width brightness look-up table storehouse, obtain the apparent spoke luminance parameter under respective conditions, perform step S408.

When selecting to search the apparent spoke luminance parameter of air from look-up table storehouse, according to aerosol type, atmospherical model, in month, solar zenith angle, view zenith angle and wavelength six factors, Auto-matching is closest to the apparent spoke luminance parameter under condition.

The principle setting up air apparent spoke brightness look-up table is with to calculate air apparent spoke luminance parameter principle in real time the same, for difference emulation atmospheric conditions, set up the spoke luminance parameter of response, generate a sextuple look-up table storehouse, so that entrance pupil spoke brightness simulation algorithm model calls.

Step S408, according to aerosol optical depth, apparent spoke luminance parameter, calculate the spectral radiance field in entrance pupil place camera response wave band.

Inputted the albedo image of terrain object by step S411, aerosol optical depth, apparent spoke luminance parameter are also inputted in apparent spoke brightness calculation model, obtain the spectral radiance field in entrance pupil place camera response wave band.

(1) according to response wave band scope and the wavelength interval of satellite sensor, determine the number of wavelengths that will emulate, successively simulation calculation is carried out to the atmospheric radiative transfer process of each wavelength.

(2) read the aerosol optical depth on first wavelength, and determine the cross radiance weighing factor of contiguous pixel.

(3) read terrain object albedo image, based on the apparent spoke brightness calculation model built, calculate the spoke brightness value of each pixel target after propagation in atmosphere in first wave length respectively.

(4) (2) ~ (4) are repeated, the apparent spoke luminance picture successively in calculated response wave band on each wavelength.

Step S409, spectral response curve according to camera, calculate full wave equivalent entrance pupil spoke luminance field image.

After the entrance pupil radiance data having calculated each wavelength, consider the spectral response curve of camera, according to the transmitance coefficient of each wavelength, calculate full wave equivalent entrance pupil spoke luminance field image.Specific formula for calculation is as follows:

$L=\frac{\mathrm{\Σ}{L}_i*R_i}{\mathrm{\Σ}{R}_i}$

In formula, L _irepresent the entrance pupil spoke brightness value of certain wavelength i, R _irepresent the transmitance value of camera on this wavelength.

Step S410, acquisition entrance pupil spoke luminance field view data, terminate.

Can obtain star loaded camera entrance pupil spoke luminance field view data through above each step, simulation process terminates.

The imaging simulation process of the present embodiment to star loaded camera is described in detail, and has following technique effect compared to prior art:

(1) atmospheric radiation transmission that the present embodiment is set up considers the cross radiance tripartite face sound of terrain object reflection, atmospheric path radiation, adjacent objects, and call parameter required in MODTRAN software computation model, ensure atmospheric radiative transfer simulation accuracy.

(2) establish sextuple atmospheric radiative transfer Parameter lookup step, the imaging simulation that can be applicable to global different time calculates, and this 6 DOF comprises aerosol type, atmospherical model, month, solar zenith angle, view zenith angle and wavelength.

(3) establish aerosol optical depth computation model, according to atmospheric aerosol and atmospheric conditions, calculate the cross radiance ratio of contiguous pixel.

(4) wave spectrum segmentation mode is adopted to calculate atmospheric radiative transfer process, and adopt camera spectral response curve, complete and multispectrally to calculate with the imaging simulation of full-color camera, can radiation transmission impact in the whole response wave band of complete reflection, ensure the simulation accuracy of satellite sensor particularly full-color camera.

Embodiment five

Fig. 5 is the structured flowchart of the satellite borne sensor imaging simulation device described in the present embodiment, and as shown in Figure 5, the satellite borne sensor imaging simulation device described in the present embodiment comprises:

Air apparent spoke luminance parameter acquiring unit 501, for according to atmospheric path radiation, be close to the cross radiance weighing factor of pixel, and terrain object albedo image obtains the apparent spoke luminance parameter of current atmospheric;

Entrance pupil spoke luminance field image acquisition unit 502, for obtaining the entrance pupil spoke luminance field image of satellite borne sensor according to aerosol optical depth and the apparent spoke luminance parameter of described current atmospheric.

Further, described air apparent spoke luminance parameter acquiring unit 501 specifically for:

According to the atmospheric path radiation of history, the cross radiance weighing factor of contiguous pixel, and terrain object albedo image sets up air apparent spoke brightness look-up table, wherein said air apparent spoke brightness look-up table comprises the apparent spoke luminance parameter of air under history observation condition and this observation condition;

Obtain Current observation condition;

Air apparent spoke brightness look-up table according to Current observation conditional search, obtains the apparent spoke luminance parameter of current atmospheric.

Further, described observation condition comprises following at least one item: aerosol type, atmospherical model, month, solar zenith angle, view zenith angle and wavelength.

Or, described current atmospheric apparent spoke luminance parameter acquiring unit 501 specifically for:

Obtain the radiation of current atmospheric journey, obtain the cross radiance weighing factor of current contiguous pixel, and obtain Current terrestrial target reflectivity image;

According to the radiation of described current atmospheric journey, the cross radiance weighing factor of described current contiguous pixel, and described Current terrestrial target reflectivity image calculates the apparent spoke luminance parameter of current atmospheric in real time.Wherein, described real-time calculating is performed by MODTRAN software.

Further, obtain the radiation of current atmospheric journey specifically to comprise:

Obtain aerosol optical depth according to current atmospheric optical property parameter, obtain the radiation of current atmospheric journey according to aerosol optical depth.

Further, the cross radiance weighing factor obtaining current contiguous pixel specifically comprises:

According to Determination of Aerosol Optical, Rayleigh scattering optics thickness, and absorb other optical thickness, obtain the cross radiance weighing factor of contiguous pixel.

The satellite borne sensor imaging simulation device that the present embodiment provides can perform the satellite borne sensor imaging simulation method that the embodiment of the present invention one, embodiment two and embodiment three provide, and possesses the corresponding functional module of manner of execution and beneficial effect.

Further, multispectrally to calculate with the imaging simulation of panchromatic sensor to complete, radiation transmission impact in the whole response wave band of complete reflection, described device also comprises all band entrance pupil spoke luminance field image acquisition unit 503, for after obtaining the entrance pupil spoke luminance field image of satellite borne sensor according to the apparent spoke luminance parameter of described current atmospheric, according to the spectral response curve of satellite borne sensor, spectral radiance field picture on each wavelength in described response wave band is weighted, obtains all band entrance pupil spoke luminance field image of equivalence.Specific implementation is identical with embodiment one to four, does not repeat them here.The present embodiment adopts wave spectrum segmentation mode to calculate atmospheric radiative transfer process, and considers the spectral response curve of satellite borne sensor, has ensured the simulation accuracy of satellite borne sensor particularly full-color camera.

All or part of content in the technical scheme that above embodiment provides can be realized by software programming, and its software program is stored in the storage medium that can read, storage medium such as: the hard disk in computing machine, CD or floppy disk.

Note, above are only preferred embodiment of the present invention and institute's application technology principle.Skilled person in the art will appreciate that and the invention is not restricted to specific embodiment described here, various obvious change can be carried out for a person skilled in the art, readjust and substitute and can not protection scope of the present invention be departed from.Therefore, although be described in further detail invention has been by above embodiment, the present invention is not limited only to above embodiment, when not departing from the present invention's design, can also comprise other Equivalent embodiments more, and scope of the present invention is determined by appended right.

Claims (16)

1. a satellite borne sensor imaging simulation method, is characterized in that, comprising:

According to atmospheric path radiation, the cross radiance weighing factor of contiguous pixel, and terrain object albedo image obtains the apparent spoke luminance parameter of current atmospheric;

The entrance pupil spoke luminance field image of satellite borne sensor is obtained according to aerosol optical depth and the apparent spoke luminance parameter of described current atmospheric.

2. the method for claim 1, is characterized in that, according to atmospheric path radiation, and the cross radiance weighing factor of contiguous pixel, and the apparent spoke luminance parameter of terrain object albedo image acquisition current atmospheric specifically comprises:

According to the atmospheric path radiation of history, the cross radiance weighing factor of contiguous pixel, and terrain object albedo image sets up air apparent spoke brightness look-up table, wherein said air apparent spoke brightness look-up table comprises the apparent spoke luminance parameter of air under history observation condition and this observation condition;

Obtain Current observation condition;

Air apparent spoke brightness look-up table according to Current observation conditional search, obtains the apparent spoke luminance parameter of current atmospheric.

3. method as claimed in claim 1 or 2, it is characterized in that, described observation condition comprises following at least one item: aerosol type, atmospherical model, month, solar zenith angle, view zenith angle and wavelength.

4. the method for claim 1, is characterized in that, according to atmospheric path radiation, and the cross radiance weighing factor of contiguous pixel, and the apparent spoke luminance parameter of terrain object albedo image acquisition current atmospheric specifically comprises:

Obtain the radiation of current atmospheric journey, obtain the cross radiance weighing factor of current contiguous pixel, and obtain Current terrestrial target reflectivity image;

According to the radiation of described current atmospheric journey, the cross radiance weighing factor of described current contiguous pixel, and described Current terrestrial target reflectivity image calculates the apparent spoke luminance parameter of current atmospheric in real time.

5. method as claimed in claim 3, is characterized in that, described real-time calculating current atmospheric apparent spoke luminance parameter performs described real-time calculating for utilizing MODTRAN software.

6. the method as described in claim 3 or 4, is characterized in that, the operation obtaining the radiation of current atmospheric journey specifically comprises:

Obtain aerosol optical depth according to current atmospheric optical property parameter, obtain the radiation of current atmospheric journey according to aerosol optical depth.

7. the method as described in claim 3 or 4, is characterized in that, the operation obtaining the cross radiance weighing factor of current contiguous pixel specifically comprises:

According to Determination of Aerosol Optical, Rayleigh scattering optics thickness, and absorb other optical thickness, obtain the cross radiance weighing factor of contiguous pixel.

8. the method for claim 1, it is characterized in that, described method also comprises, after obtaining the entrance pupil spoke luminance field image of satellite borne sensor according to described aerosol optical depth and the apparent spoke luminance parameter of current atmospheric, according to the spectral response curve of satellite borne sensor, spectral radiance field picture on each wavelength in response wave band is weighted, obtains all band entrance pupil spoke luminance field image of equivalence.

9. a satellite borne sensor imaging simulation device, is characterized in that, comprising:

Air apparent spoke luminance parameter acquiring unit, for according to atmospheric path radiation, be close to the cross radiance weighing factor of pixel, and terrain object albedo image obtains the apparent spoke luminance parameter of current atmospheric;

Entrance pupil spoke luminance field image acquisition unit, for obtaining the entrance pupil spoke luminance field image of satellite borne sensor according to aerosol optical depth and the apparent spoke luminance parameter of described current atmospheric.

10. device as claimed in claim 9, is characterized in that, described current atmospheric apparent spoke luminance parameter acquiring unit specifically for:

According to the atmospheric path radiation of history, the cross radiance weighing factor of contiguous pixel, and terrain object albedo image sets up air apparent spoke brightness look-up table, wherein said air apparent spoke brightness look-up table comprises the apparent spoke luminance parameter of air under history observation condition and this observation condition;

Obtain Current observation condition;

Air apparent spoke brightness look-up table according to Current observation conditional search, obtains the apparent spoke luminance parameter of current atmospheric.

11. devices as described in claim 9 or 10, it is characterized in that, described observation condition comprises following at least one item: aerosol type, atmospherical model, month, solar zenith angle, view zenith angle and wavelength.

12. devices as claimed in claim 9, is characterized in that, described current atmospheric apparent spoke luminance parameter acquiring unit specifically for:

Obtain the radiation of current atmospheric journey, obtain the cross radiance weighing factor of current contiguous pixel, and obtain Current terrestrial target reflectivity image;

According to the radiation of described current atmospheric journey, the cross radiance weighing factor of described current contiguous pixel, and described Current terrestrial target reflectivity image calculates the apparent spoke luminance parameter of current atmospheric in real time.

13. devices as claimed in claim 12, is characterized in that, described real-time calculating current atmospheric apparent spoke luminance parameter performs described real-time calculating for utilizing MODTRAN software.

14. devices as described in claim 11 or 12, is characterized in that, obtain the radiation of current atmospheric journey and specifically comprise:

Obtain aerosol optical depth according to current atmospheric optical property parameter, obtain the radiation of current atmospheric journey according to aerosol optical depth.

15. devices as described in claim 11 or 12, it is characterized in that, the cross radiance weighing factor obtaining current contiguous pixel specifically comprises:

According to Determination of Aerosol Optical, Rayleigh scattering optics thickness, and absorb other optical thickness, obtain the cross radiance weighing factor of contiguous pixel.

16. devices as claimed in claim 9, it is characterized in that, described device also comprises all band entrance pupil spoke luminance field image acquisition unit, for after obtaining the entrance pupil spoke luminance field image of satellite borne sensor according to aerosol optical depth and the apparent spoke luminance parameter of described current atmospheric, according to the spectral response curve of satellite borne sensor, spectral radiance field picture on each wavelength in response wave band is weighted, obtains all band entrance pupil spoke luminance field image of equivalence.