WO2002082027A2

WO2002082027A2 - Method for obtaining a simulated digital representation of the radiative state of a heterogeneous three-dimensional scene

Info

Publication number: WO2002082027A2
Application number: PCT/FR2002/001181
Authority: WO
Inventors: Philippe Gastellu-Etchegorry
Original assignee: Universite Paul Sabatier Toulouse Iii
Priority date: 2001-04-06
Filing date: 2002-04-04
Publication date: 2002-10-17
Also published as: FR2823344A1; WO2002082027A3; FR2823344B1

Abstract

The invention concerns a method for obtaining a simulated digital representation of the radiative energy state of a heterogeneous three-dimensional scene cut out into a plurality of elementary cells (71, 72, 73). It consists in: defining a plurality of incidence sectors fewer in number than the incident directions; calculating cell by cell, an intercepted sector energy equal to the sum of intercepted energies of incident rays in the incidence sector; using a sector diffusion transfer function enabling to calculate a diffused sector energy at the original point of diffusion in the re-emission direction Φvj ; calculating the diffused volume energy as the sum of diffused sector energy.

Description

PROCESS FOR OBTAINING A REPRESENTATION

SIMULATED DIGITAL RADIATIVE STATE OF A SCENE

THREE-DIMENSIONAL HETEROGENEOUS

The invention relates to a method for obtaining a simulated digital representation of the radiative energy state, in the visible range and / or outside the visible range -in particular in the range understood from ultraviolet to thermal infrared- , of a heterogeneous three-dimensional scene, in particular a portion of space likely to extend in the field of vision of a satellite sensor above the ground of a planet such as Earth.

Such a simulated digital representation makes it possible either to determine the radiative energy state at any point on the scene, or to produce digital images of the scene each seen from a point outside the scene (for example a satellite or airborne sensor) located in one direction of observation (different simulated images seen in different directions of observation are generally produced simultaneously).

Obtaining such simulated digital representations capable of providing with high radiometric precision the radiative state of a three-dimensional scene (and not only the general visual aspect as in virtual image synthesis techniques) is important in certain applications, for example for the precise calculation of biophysical variables (biomass, study of vegetation ...) or biochemical from satellite measurements, or for the ground calibration of satellite sensors with high spatial resolution (can now reach l '' order of 1 m or less).

The publication "Modeling Radiative Transfer in Heterogeneous 3-D Vegetation Canopies" JP. Gastellu-Etchegorry et al Remote Sens. About. 58: 131-156 (1996) describes a method for obtaining such digital images. It is based on an iterative approach: at iteration n, there is diffusion of the energy W _out (Ω) according to the directions (Ω) of space, taking into account the spatial distribution of the intercepted energy W * _nt , and not absorbed, in the previous iteration n-1. The iteration n = l includes the illumination by an external source mono-directional (such as the sun), while iteration 2 includes the diffusion of the energy intercepted in iteration 1, and thermal emission, if applicable. This process suffers from two drawbacks:

1) As this publication indicates, the computing times remain extremely long and the memory capacities required are extremely large compared to the dimensions of the scene. Thus, any simulation applied to a 50mx50x20m forest landscape with 20cm cells, with 150 directions and with radiometric precision compatible with the usual precision of sensors (satellite and others) requires a calculation time of the order of a week, depending on the type of workstation used. This calculation time is essentially due to the method of spherical harmonics initially used to manage multiple diffusions. This method is indeed very long and is applied to any cell and to any radius, at each iteration.

2) Multiple diffusions are poorly taken into account, greatly affecting the accuracy of the model. The main cause of imprecision is the ignorance of the anisotropy of the incident radiation which is diffused. Indeed, any scattering depends on the angular distribution of the incident radiation. This is only well known in order 1, ie in iteration 1, because the solar direction is unique. For iterations n> l, the incident radiation can come from any direction and is a priori anisotropic. Memorizing the angular distribution of the incident radiation at each iteration, for each cell, requires very large memory volumes (of the order of a hundred gigabytes for a landscape of 50mx50mx20m). In theory, the decomposition of the scattered radiation in the form of spherical harmonics is supposed to solve this problem. In practice, in addition to the problems of computation time and storage memory volume mentioned above, this solution is imprecise, since it can only be stored a finite number of decomposition coefficients. In addition, the fact of multiplying a dimension of the scene by a factor m leads to multiplying the quantity of information to be processed by a factor m ⁶ . The drawbacks indicated above therefore prevent any practical industrial application. The invention therefore generally aims to overcome these drawbacks by proposing a method which makes it possible to obtain simulated digital representations of large scenes, having a large number of cells, with sufficient radiometric precision, and with calculation times. and storage capacities compatible with practical industrial operation.

The invention more particularly relates to such a method making it possible to simulate planetary remote sensing images - notably terrestrial - from space, for any experimental configuration (heterogeneous or non-heterogeneous landscape, any lighting and observation directions, etc. ).

The invention also more particularly relates to such a method making it possible to obtain radiometric images making it possible to calibrate modern satellite sensors on the ground before their launch.

The invention aims in particular to enable a better understanding of satellite measurements and consequently to improve their design, their exploitation and ultimately contribute to a better evaluation of the possibilities offered by satellite sensors.

The invention also aims to provide a simulated three-dimensional representation of the radiation balance of the studied landscape, essential information for many environmental applications such as the study of the functioning of vegetation (photosynthetic activity, gas flow, etc.).

To do this, the invention relates to a method for obtaining a simulated digital representation of the radiative energy state of a heterogeneous three-dimensional scene, in which:

the scene is divided into a plurality of elementary cells comprising cells, called turbid cells, having known radiative properties assumed to be homogeneous in the whole of each cell,

- transfer functions are memorized for each turbid cell making it possible to determine, in each direction of a discrete plurality of N _v directions of space, said directions of re-emission Ω _v> predetermined in an independent space reference cells, an energy re-emitted by the turbid cell as a function of an intercepted radiation due to an incident radiation which it receives formed of at least one incident ray defined by an incident direction Ω _sk belonging to a plurality N _s of directions of space, called incident directions Ω _S; and an incident energy Wi _n (Ω _Sk ), a fraction of which, called intercepted energy Wj _nt (Ω _sk ), is intercepted by the turbid cell,

- for at least one portion, called the portion of the ground, of the scene, a cell re-emitted energy W _0Ut (Ω _vj ) is calculated, cell after cell, by each turbid cell in each of the re-emission directions Ω _yj, this re-emitted energy W _out (Ω _vj ) being considered as coming from a diffused volume energy _d iff (Ω _vj ) in the direction of re-emission Ωy _j by the whole of the turbid cell, this volumetric energy diffused

being calculated by considering that it is emitted completely starting from a single point of the turbid cell, said point equivalent origin of diffusion M _{s> k} (Ω _vj ), characterized in that: - one defines in a reference mark of l space independent of the cells a plurality of M sectors, called incidence sectors ΔΩ ^* ι, i, in number M less than the incident directions Ω _s> and containing all the N _s incident directions Ω _s>

- to calculate, at least for the portion of soil, the volume energy scattered ^ _f (Ω _vj ) in each re-emission direction Ω _vj to each turbid cell:

. we calculate, for each incident sector ΔΩ, J, an energy, called intercepted energy, of sector Wi „ _t (ΔΩ ^* - _V | ^* _j j) _> equal to the sum of the intercepted energies

incident rays on the turbid cell in incident directions Ω _sk belonging to the incidence sector ΔΩ ^* _M , i,

. we use a function, called sector diffusion transfer function T (ΔΩM, i, Ω _j ) allowing to calculate for each incidence sector ΔΩ-y- ^{* *} and each re-emission direction Ω _vj , an energy, said scattered energy of sector W _d if ^ ΔΩ-M ^ Ω _yj ), scattered at the equivalent origin point broadcast

in the re-emission direction Ω _vj from intercepted sector energy

of the bearing sector

. for each re-emission direction Ω _{j (} we calculate the scattered volume energy W ^ Q _Vj ) by the turbid cell as the sum of the sector scattered energies

scattered by each incidence sector ΔΩ- -y in this re-emission direction Ω _Vj .

Advantageously and according to the invention, the sector diffusion transfer function T (ΔΩ _M) *, Ω _vj ) is a weighted average of diffusion transfer functions T (Ω _s , Ω _vj ) stored for each incident direction Ω _sk and each re-emission direction Ω _vj initially defined.

Advantageously and according to the invention, the sector diffusion transfer function T (ΔΩ _M! I, Ω _Vj ) is defined according to the following formula I:

Σ T (Ω _sk> Ω _vj ) .ΔΩ _sk

(I) T (ΔΩ _Mji , Ω _VJ ) = ^ Kd) k Σ = o ΔΩ _s ^s _k ^κ (i) where _] Σ ^ ΔΩ _Λ = ΔΩ _M , i,

and where K (i) is the number of incident directions Ω _Sk in the sector ΔΩ _Mji .

The functions T (ΔΩ _Mιi , Ω j) are pre-calculated in order to avoid repetitive calculation of these.

Advantageously and according to the invention, the number M of incidence sectors ΔΩM, i is greater than 2 and less than 20, in particular between 4 and 14, for example and preferably equal to 6, and much less than the number N _s of incident directions Ω _s which is greater than 100, in particular between 100 and 1000. In addition, advantageously and according to the invention, the possible incident directions Ω _s of the incident rays are determined and chosen from the same discrete plurality of directions of l space that the directions of re-emission Ω _v . This plurality of directions is determined absolutely in a reference frame of space, independently of the cells.

Advantageously and according to the invention, for each ray Ω _sk incident on the turbid cell, an equivalent point of origin is determined. M _{S) k} (Ω _yj ) of the scattering of the incident ray by the turbid cell in the re-emission direction Ω _VJ -. The subscript "k" indicates that the position of the point M _{s> k} (Ωy _j ) depends on the incident direction Ω _sk . The turbid cells being formed of polyhedra for each incident ray in a direction Ω _g and which enters a point Pin on one of its faces, called the entry face, the equivalent point of origin of the diffusion M _{s> k} (Ω _vj ) according to the re-emission direction Ω _vj is determined as the point located in the turbid cell, on the incident direction Ω _sk , and at a distance Δr _k (Ω _sk ) from the input face by the following formula II :

(H) Δr _k (Ω _sk ) =

G (Ω _y] ). | Μ _sk | .u _f

G (Ω _sk ). | Μ G (Ω _sk ) G (Ω _V1

^{aVeθ A [} [l-exp [-G (Ω _sk ) .u _f .Δl _k (Ω _3k )]] ^]

where θ _sk is the angle formed by the incident direction Ω _sk with a predetermined fixed direction of space, called the reference direction,

Θ _VJ is the angle formed by the re-emission direction Ω _vj with the reference direction,

and Δl (Ω _sk ) is the path of the incident ray through the cell, fg Ω _f ) G (Ω _sk ) =, 2. | Ω _s .Ω _f | .dΩ _f is a coefficient, called projection factor

perpendicular to the direction Ω _s . This projection applies to the material included in the turbid cell knowing that this material is oriented in the direction Ω _f , with Ω _f included in the sector 2π, and that the probability of having

the orientation Ω _f is ^st ir ^J1 . We have ^s i.

2π f fCΩf) G (Ω _Vj ) = 1 .jΩ _VJ .Ω _f | .dΩ _f is a projection factor perpendicular to the

2π direction Ωy _j .

In this way, the scattering point M _Sιk (Ω _Vj ) depends on the scattering directions Ω _vj and incident Ω _sk , which considerably improves the accuracy of the results.

Advantageously and according to the invention, one calculates and uses as an equivalent origin point of diffusion, for each incident ray, at least one average origin point of diffusion M _sk independent of the direction of re-emission Ω _vj) and one calculates the scattered energy of sector W _d i _f ï (ΔΩ _M7 ^* , Ω _j ) from each incident ray of incident direction Ω _sk as scattered at this point mean origin of scattering M _s - In this way, the calculation time and memory space required.

Advantageously and according to the invention, one calculates, and one uses as point of origin equivalent of diffusion, for each incident ray, two points of origin means of diffusion: a point M _sk (ΩÎ) associated with a radiation diffusion according to a first predetermined hemisphere of space, called upper hemisphere, and a point M ^ Ω-) associated with the scattering of radiation according to a second predetermined hemisphere of space, called lower hemisphere, the two upper and lower hemispheres being chosen so that the reference direction is orthogonal to the diametral plane separating these two hemispheres. We thus achieve a good compromise between on the one hand a good precision of results, and on the other hand computing times and a reasonable necessary memory space. Advantageously and according to the invention, the reference direction is the direction orthogonal to the diametral plane separating the upper and lower hemispheres.

In the case where the simulated digital representation comprises at least one simulated image seen from outside the scene in a direction of observation oriented towards an outside observer, fixed and predetermined with respect to the scene, advantageously and according to the invention, the direction of observation is oriented in the upper hemisphere. Advantageously and according to the invention, the direction of reference orthogonal to the diametral plane separating the two upper and lower hemispheres is the vertical of a ground zone of a planet, the scene extending from said ground zone of the planet. Advantageously and according to the invention, the method makes it possible to simultaneously obtain Nv images seen according to Nv directions of observations oriented towards the upper hemisphere.

Advantageously and according to the invention, one chooses

M _sk (Ωt) = M _{S) k} (Ω _vj (θ _vj = 20 °)) and M _sk (Ωi) = M _{s> k} (Ω _vj (θ _vj = 160 °)). These choices are generally very good compromises, especially for the simulation of remote sensing images.

In most practical applications, advantageously and according to the invention, two average origin points M _S (ΩÎ) and M _s (Ω-J _' ) are calculated respectively associated with the diffusion in the rising directions belonging to said upper hemisphere and respectively the descending directions belonging to said lower hemisphere, these two average origin points M _s (Ωf) and M _s (Ωl) being used as equivalent origin origin points for all the rays incident on the turbid cell. In this way, the calculation of the re-emitted energy W _0Ut (Ω _Vj ) remains simple even when the cells receive several incident rays, as is generally the case. For each turbid cell, the points M _S (ΩÎ) and M _s (Ω-l) are calculated respectively as energetic barycenters of the origin point of diffusion means M _sk (ΩÎ) and M _s (Ω-l- *), calculated on the different rays incident on the turbid cell, with the iterative relationships:

. _o , _W _int .t _s (Ωt) + W _int (Ω _sk ). (Ωt ^U Wint + intCΩd-

,. _{f QK} , VvWtJΩ-l) + W *, _t (Ω _sk ) .t, _k (ΩJ) and u "* JW _iπt + W _int (Ω _sk ) where t _s (ΩÎ) indicates the coordinates x, y and z point

M _S (ΩÎ), t _sk (ΩÎ) indicates the x, y and z coordinates of the point M _sk (Ωt), and W _int is the energy intercepted without counting the energy W _{] nt} (Ω _sk ) which has just been intercepted. After this last interception, the intercepted energy is therefore: Wi _nt + W _ιnt (Ω _sk ).

This calculation is performed sequentially during any interaction with a cell.

Advantageously and according to the invention, the method is further characterized in that before calculating the re-emitted energies W _0Ut (Ω _vj ) by each cell, each mean origin point M ₃ (ΩT) and M _S (Ω> 1 ') is replaced by an equivalent emission point which is the nearest point Pe chosen from a plurality of possible emission points Pe regularly distributed in each cell, for which all the possible trajectories of the rays emanating from all the points Pe according to all possible re-emission directions Ω _vj have been previously calculated and recorded.

Thus, advantageously and according to the invention, the prior calculation of all the possible trajectories of the rays in the scene is carried out for a grid of several points inside the origin cell, and not only for the center of the cell. So any ray path from the point

Pe of any cell with coordinates (X, Y, Z) in the scene, determined from an origin cell (0,0,0), is obtained by translation (X, Y, Z) of the path of a parallel radius from the corresponding point Pe in the origin cell (0,0,0).

During any diffusion, the two points Pe closest to the average origin points M _S (ΩÎ) and M _s (Ωi) are chosen.

Advantageously and according to the invention, the method is characterized in that one calculates by successive iterations a radiation re-emitted by each cell comprising at least one ray re-emitted in at least one direction of re-emission Ω _vj , by calculating and memorizing at each iteration:

- the radiation incident on the cell originating from rays re-emitted by the other adjacent cells as determined in the previous iteration,

- the radiation re-emitted by the cell originating in particular from the diffusion of the intercepted radiation due to the radiation incident on the cell. In addition, advantageously and according to the invention, the method is characterized in that one defines beforehand in the scene, independently of the cells, in particular by triangularization, geometric elements which represent different elements of the landscape, and the energy is calculated re-emitted W _out (Ω _vj ) taking into account these geometric elements as mainly reflecting and / or absorbing elements.

The invention thus makes it possible for the first time to produce simulated images taking into account simultaneously and with high radiometric precision, on the one hand, the translucent media (for example vegetation, grasses, leaves of trees ...) thanks to the division into turbid cells, but also the opaque environments (walls, roofs, floors, bodies of water, natural reliefs ...) thanks to the geometric elements. The diffusion, absorption and emission mechanisms are systematically taken into account using very precise physical modeling (transfer functions for translucent media and parametric models for opaque media).

Furthermore, in certain advantageous applications of the invention, the scene comprises a portion of soil from a planet, in particular the Earth, and a portion of space extending above this portion of soil. The invention then makes it possible to take into account with precision the three-dimensional nature of the atmosphere and of the planetary surfaces as regards the propagation of the radiation. This aspect is essential for obtaining images with good radiometric accuracy. Indeed, the "atmosphere-earth" interaction mechanisms are extremely complex, in particular due to the three-dimensional nature of the atmosphere and of terrestrial surfaces. Thus, advantageously and according to the invention, the scene comprising a portion of soil from a planet, an upper column of the scene, above the portion of soil, is defined as an atmospheric column divided into cells, called atmospheric cells having , like turbid cells, known radiative properties assumed to be homogeneous throughout each atmospheric cell. All atmospheric cells of the same altitude, that is to say located in the same plane perpendicular to the reference direction, have the same dimensions and physical properties (density, extinction coefficients, etc.). Advantageously and according to the invention, the dimensions of the atmospheric cells increase with altitude. Indeed, the gases and aerosols being less dense at altitude, the size of the cells can increase without compromising the radiometric precision. In addition, advantageously and according to the invention, the method is characterized in that the atmospheric column comprises a lower zone, called the lower atmosphere, in contact with the ground portion (defining a landscape zone up to the planetary ground), and an upper zone, called the upper and intermediate atmosphere. In this way, we note that we manage to take into account the radiative behavior (diffusion, absorption, thermal emission, transmittance) of the atmosphere with good precision, despite the complexity of its three-dimensional mechanisms. Thus, from a predetermined solar radiation defined by an incident ray on the upper face of each cell of the upper layer of the atmospheric column, the following successive executions are carried out:

1) determining the incident radiation on the upper face of the cells of the lower atmosphere from the solar radiation after propagation in the upper and intermediate atmosphere (in the field of long wavelengths, the incident radiation also includes atmospheric thermal radiation ).

2) the radiation (thermal radiation and reflected radiation due to the incident radiation determined during execution 1)) is determined, amount coming from the lower atmosphere,

3) the radiation backscattered by the upper and intermediate atmosphere towards the lower atmosphere is determined from the reflected radiation determined during execution 2),

4) the radiation reflected by the lower atmosphere to the upper and intermediate atmosphere is determined from the backscattered radiation determined at runtime 3), 5) the final radiation rising in the atmospheric column is determined - in particular at its face superior- from sum of the reflected radiations determined in executions 2) and 4) after propagation through the upper and intermediate atmosphere.

Preferably, the upper and intermediate atmosphere comprises an upper zone, called the upper atmosphere, which in embodiment 3) only diffuses radiation upwards, all the radiation arriving in this upper atmosphere being scattered upwards. In execution 3), only the intermediate atmosphere (zone of the upper and intermediate atmosphere located under the upper atmosphere) is therefore capable of diffusing radiation towards the lower atmosphere. The invention extends to a software product capable of being loaded into the random access memory of a computer system for the implementation of a method according to the invention so as to produce a computer device according to the invention.

This software product is a computer program suitable for loading into the internal RAM of a computer. It can be offered in the form of a recording medium or made available for download on an information transmission network (internet). The invention also extends to a recording medium adapted to be able to be read by a reader connected to a computer system, characterized in that it comprises a recorded program adapted to be able to be loaded into random access memory of the computer system and the program to implement a method according to the invention. A recording medium according to the invention therefore comprises a software product according to the invention and makes it possible, when read on a reader of a computer, the software product being loaded into random access memory, to obtain simulated digital representations - in particular simulated images - of the radiative energy state of a heterogeneous three-dimensional scene.

The invention extends to a method, a software product and a recording medium characterized in combination by all or some of the characteristics appearing above or below. Other objects, characteristics and advantages of the invention appear from the following description which refers to the appended figures in which: FIG. 1 is a general functional flow diagram of a method according to the invention,

- Figure 2 is a diagram illustrating an example of the three-dimensional scene of a method according to the invention, - Figure 3 is a diagram illustrating an example of a portion of soil and, partially, an atmospheric portion of a heterogeneous three-dimensional scene terrestrial cut into parallelepiped cells in a process according to the invention,

FIGS. 4a to 4c illustrate examples of incidence sectors defined in a method according to the invention,

FIG. 5 is a diagram illustrating the different stages of taking into account the atmosphere-soil coupling in a method according to the invention,

FIG. 6 is a block diagram illustrating the determination of an emission point,

FIGS. 7a and 7b are diagrams illustrating the determination of the means of origin points of diffusion of a turbid cell in a method according to the invention,

FIG. 8 is a general chronological flowchart of the successive steps of a method according to the invention,

FIG. 9 is a chronological flowchart of the calculation of emission points of the turbid cells during the first iteration of a method according to the invention,

FIG. 10 is a chronological flowchart of the calculation of an emission point of an opaque cell in a method according to the invention,

FIG. 11 is a chronological flowchart of the calculation of emission points of the turbid cells in subsequent iterations (other than the first) of the method according to the invention,

- Figure 12 is a chronological flowchart of the steps implemented by the atmosphere module of a method according to the invention.

FIG. 1 represents a diagram of the various functional modules of a method according to the invention. In this figure, we have represented in dotted the input modules 1, 2 and 3 which allow the user to enter the various configuration parameters according to the envisaged application.

The input module 1 makes it possible to define the geometric characteristics of the three-dimensional scene that one wishes to represent. In the example shown, the scene comprises a portion of soil 6 of a planet, in particular the Earth, and a portion of atmosphere 7 extending from the elements of the portion of soil to a predetermined altitude.

The input module 2 makes it possible to define the optical (leaf reflectance pg, leaf transmittance τ _fî , etc.) and thermal (T ^ of the three-dimensional scene that one wishes to represent) characteristics.

From the parameters entered in the input modules 1 and 2, a module, called the model module 4, creates a computer model 5 representative of the scene, represented in FIGS. 2 and 3. In this model 5, the scene is divided into a plurality elementary cells which in the example shown and preferably are cubic. Alternatively, these cells can be parallelepiped.

In the model 5, it is also defined by the input module 1, independently of the cells, of the geometric elements which represent different landscape elements mainly reflecting and / or absorbing -notably the relief of the ground (including the bare ground surfaces, aquatic surfaces and roads, but excluding plant surfaces), walls and roofs of buildings, and tree trunks. These geometric elements are defined by triangularization of the landscape, that is to say by a mesh representing the landscape as a juxtaposition of polygons (triangles and parallelograms) according to the coordinates of the extreme points of the geometric elements entered and according to the fineness of the spatial mesh chosen. in the input module 1.

The input module 3 makes it possible to define in a benchmark of the space independent of the cells the mode of discretization of different directions of the space according to which the radiation can propagate. At least a plurality of directions and their spatial distribution are defined. For each plurality of directions, the number N of directions Ω is defined, and possibly specific distribution densities of these directions in certain angular sectors. We can thus define N _s incident directions Ω _s , and

N _v re-emission directions Ω _v .

Preferably, the possible incident directions Ω _s of the incident rays are determined in the same discrete plurality of directions of space as the directions of re-emission Ω _v (N _S = N _V = N). In addition, an additional particular direction is defined corresponding to the direction, called the initial incident direction Ω °, of a main light source such as the sun. There are therefore Nt-1 defined directions. Consequently, any model 5 contains:

- turbid cells, each turbid cell having known radiative properties assumed to be homogeneous throughout the turbid cell; in particular, turbid cells, called leaf cells 71, represent the vegetation (grasses, crops, meadows, crowns of trees, etc.) and turbid cells, called atmospheric cells 72, represent the atmosphere; ,

- cells called opaque cells 73 containing at least part of a geometric element as described above, and possibly turbid material (vegetation or atmosphere). As can be seen in FIG. 2, the atmosphere portion advantageously comprises several discrete atmosphere levels:

a low atmosphere level 11, the cells 72 of which have the same dimension as the cells of the soil portion 6 and which extend above the leaf 71 or opaque cells 73 of the soil portion 6; an intermediate atmosphere level 12 extending above the low atmosphere 11, and the cells 72 of which have a larger dimension;

- A high atmosphere level 13 whose cells 73 have an even larger dimension - in particular each occupy an entire horizontal layer of the scene -.

The turbid cells 71, 72 are characterized by specific optical and geometric properties defined by means of the entry 1,2. Thus, for each type of vegetation, a leaf angular orientation is defined represented by a LAD coefficient and a leaf biomass represented by a LAI coefficient. For the atmosphere, it is defined for gases (exponent m) and aerosols (solid or liquid particles; exponent p) vertical profiles (varying according to altitude z) of albedo of simple diffusion ω ^m (z) and ω ^p (z), extinction coefficient α ™ (z) and α ^ z), absorption extinction coefficient α "(z) and α ^ (z), of diffusion transfer function T ^m ( Ω _sk , Ω _vj ) and T ^p (Ω _sk , Ω _vj ), as well as a temperature profile T (z).

The extinction coefficients α ™, α ™, a and α ^ are determined at any altitude z from knowledge of the total atmospheric optical thicknesses associated with gas absorption τ ™, with gas diffusion τ ^™ , with absorption of aerosols τ ^ and there diffusion of aerosols τ ^p . For this, the density profiles of gases and aerosols are assumed to decrease exponentially with altitude. Ozone is the only exception. Consequently: α (z) ^≈ 7T - .. eexxpp [[- gτT] | where the scale factor H specifies the type of exponential profile.

A scale factor H _{m is used} for gases uniformly mixed in the atmosphere, ie all gases other than ozone and water vapor. A scale factor H _{H∑O is used} for water vapor. An analytical law is used to represent the profile of ozone, since it is mainly centered around 20-30km of altitude. The uniformly mixed gases are essentially N ₂ , O ₂ , CO ₂ , CH and N ₂ O. On the other hand, a scale factor H _{p is used} to represent the density profile of aerosols.

The optical thicknesses τ ™ (λ) and τ ™ (λ) are calculated at any wavelength λ from a table of optical thicknesses τ ™ (λ) and τ ^ λ,) at wavelengths λi pre-defined. Values

are calculated with atmospheric models like 6-S (Nermote EF, Tanré D., Deuzé JL, Herman M. and Morcrette J.-J., "Second simulation of the satellite signal in the solar spectrum, 6S: An overview", IEEE Trans. Geosci. Remote Sensing, vol. 35, pp. 675-686, May 1997) for standard atmospheric conditions such as a tropical atmosphere. The user can give particular values, if necessary. The high spatial and temporal variability of atmospheric water vapor explains why the table stores the optical thickness τ ™ (λ *, d _k ) of H ₂ O for K values of water thickness that can be precipitated d _k . Consequently, the model module 4 calculates τ ™ (λ, d) by interpolation on the quantity of atmospheric water that can be precipitated d and on the wavelength λ.

The angular distribution of the gas diffusion is modeled with the Rayleigh phase function (Lenoble J., 1993, Atmospheric radiative transfer, ISBN 0-937194-21-2, Deepak Publishing, Hampton, Virginia, USA):

Pκ (Ψ _sv ) ≈ ^ ^1. (L + cos Ψ _sv ) +

≈0.7552 + 0J345.cos ² Ψ _sv

=> P _R (Ψ _sv ) ≈ |. (L + cos ² Ψ _sv )

where δ is the depolarization factor. It depends on the anisotropy of the molecules and is always small. For an average terrestrial atmosphere δ≈0.0279. the angle Ψgv is the angle between the directions Ω _s and Ω _v . It checks the relation: cosΨ _sv = [cosθ _v .cosθ _s + sinθ _v .sinθ _s .cos (Φ _v -Φ _s )]

The optical thickness of the aerosols is calculated at any λ wavelength λ by the empirical law of Angstrom: τjî (λ) = τ ^p _; (λ ₀ ). [r-] ^{'' | î} . The λ ₀ coefficient of Angstrom β depends on the type of aerosols present. It is generally between 0.5 and 2. On the other hand, we have: τ ^p (λ) = ω _p .τ _g (λ) and τ ^ (λ) = (l-ω _p ) .τ ^ (λ), where ω _p is the simple diffusion albedo. The calculation of the diffusion and absorption rate by atmospheric aerosols at any altitude z therefore requires the knowledge of 4 parameters: the optical thickness τ ^ (λ ₀ ) at a pre-defined wavelength λ ₀ , the coefficient β, the simple diffusion albedo ω _p and the scale factor H _p . The angular distribution of the aerosol diffusion is modeled with the sum of two P _P (Ψ) functions from Henyey-Greenstein (Lenoble, 1993):

where the asymmetry factor g represents the front diffusion surplus compared to the rear diffusion. It is therefore zero if P _p (Ψ) is symmetrical. It depends on the type of aerosol and the wavelength.

A geometric module 9 calculates and stores the incident directions Ω _s and the re-emission directions Ω _v as a function of the values entered in the input module 3. In addition, the geometric module 9 calculates and stores, in the reference frame of the space independent of the cells, a plurality of M volume sectors, called incidence sectors ΔΩ _ι j, in number M less than the N _s incident directions Ω _s but containing all the N _s incident directions Ω _S) complementary to each other in l solid solid angle and forming a partition of the solid solid angle.

Preferably, the number N _s of incident directions (the number N of directions in the case where N _S = N _V = N) is greater than 100, in particular between 100 and 1000.

Preferably, the number M of incidence sectors, entered using the entry module 3, is greater than 2 and less than 20. FIG. 4a represents an example with M = 6; FIG. 4b represents an example with M = 10; FIG. 4c represents an example with M = 14. The value M = 6 represents a good compromise for terrestrial scenes.

A module, called the optical module 8, calculates and stores for each type of turbid cell four transfer functions from the optical and geometric properties of the elements contained in the cell entered by the input modules 1, 2:

(1) a transmittance transfer function T (Ω), per unit of displacement and per unit of surface density of the material contained in the cell. This function calculates the energy transmitted W _tt (Ω) through the turbid cell, and therefore also the energy intercepted Wj „ _t (Ω) by the turbid cell, according to each of the N _s +1 pre-defined incident directions Ω _s , (cf. publication" Modeling Radiative Transfer in Heterogeneous 3-D Vegetation Canopies "JP. Gastellu-Etchegorry et al Remote Sens. Environ. 58: 131-156 (1996));

(2) an initial diffusion transfer function T (Ω _s , Ω _Vj ) memorized for the initial incident direction Ω ° and each re-emission direction Ω _VJ initially defined; this function indicates the energy that any voluminal element of the cell diffuses at the first iteration of the method according to each of the N _v directions of re-emission (Ω _j ), taking into account the energy Wi, _πt (Ω _Sk ) intercepted according to the initial incident direction (Ω °); this initial diffusion transfer function T (Ω °, Ω _vj ) is previously calculated by the optical module 8;

(3) a sector diffusion transfer function T (ΔΩ _M , i, Ω _Vj ) making it possible to calculate, for each incident sector ΔΩ _{M) i} and each re-emission direction Ω _vj , an energy, known as scattered energy sector ig (ΔΩ _M> i, Ω _j ), scattered at the equivalent origin point of scattering M _{S) k} (Ω _vj ) in the re-emission direction Ω _VJ from the intercepted sector energy Wi _n t ( ΔΩ _M , i) of the incident sector ΔΩ _M , i. This function T (ΔΩ, i, Ωv _j ) is calculated according to formula I above. The diffusion transfer functions T (Ω _sk , Ω _v *) are previously calculated by the optical module 8;

(4) a heat transfer function T (u _f , T _f , ε _f , Ω _vj ) ui makes it possible to calculate the thermal energy re-emitted by the cell in the re-emission direction, Ω _VJ , taking into account the quantity of material u _f present as well as the temperature T _f and the emissivity ε _f of this material (Thesis "Modeling of radiative and energy balances of plant cover", P. Guillevic, Paul Sabatier University, 13/12 / 1999. Paragraph m A3.2 p47 to 49).

It should be noted that in a method according to the invention, a transmittance transfer function T (Ω) is used which is the same for each type of cell and each direction Ω, and which does not depend in particular the path of the ray in the cell or the amount of material present in the cell.

For example, for leaf vegetation cells, the transmittance function defined by Beer-Lambert's law is written: T (Δli, Ω) = exp [-GG, Ω) .uf (i). Δli] where Δli is the distance of the ray path in the cell, uf (i) is the quantity of leaf area per unit of volume G (j, Ω) is a projection factor of the leaves perpendicular to the direction Ω considered, defining a type of vegetation. In the method of the invention, an approximation of this function is used, and only N values are stored for each type j of turbid cell, ie T (Ω) = exp [-G, Ω)]. For foliar vegetation cells with foliar volume density U (i), we therefore have: T (Δ1 ;, Ω) = [T (/ 2)] ^u ^ ^Λl

All these calculations being carried out, a module, called a radiative module 10, simulates the propagation of the radiation in the scene by successive iterations, and, at each iteration, by ray tracking through the cells of the scene, as explained below, and provides a simulated digital representation of the radiation state of the scene. This representation can include remote sensing images 11, for all rising Ω _VJ * re-emission directions and / or three-dimensional matrices 12 representative of the radiation balance of the scene.

FIG. 8 is a general flow diagram of the different steps implemented in a process for obtaining according to the invention, by the radiative module 10.

As illustrated diagrammatically in FIG. 6 in the initial step 14, the radiative module 10 predefines in a source cell with coordinates (0,0,0) (FIG. 3) a plurality of possible emission points Pe, regularly distributed within of the cell. In the example shown in FIG. 6, the cell is seen in elevation and comprises 5 ³ = 125 possible emission points Pe. In addition, the radiative module 10 calculates beforehand and records all the possible trajectories through the scene of the rays that can be emitted from each possible emission point Pe according to the different directions of re-emission Ω _v> by ray tracking. This initial calculation 14 therefore provides n ³ .N linked list which are stored, where n ³ is the number of points Pe and N is the total number of directions Ω other than the initial direction Ω °. A linked list therefore corresponds to a possible path in a given direction and from a given point Pe. Any linked list item is a structure that contains the coordinates of the entry point and the length of the path in the current cell as well as the index that indicates the next cell.

The solar illumination at the top of the atmosphere 7 is first of all determined. This solar illumination is defined by the initial incident direction Ω ° and by an initial incident energy Wu, (Ω °) as shown in FIG. 2.

A first step 16 of a module, called the atmosphere module, simulates the propagation of solar radiation ^" W _m (Ω °) through the atmosphere and determines 17 the direct solar illumination incident on the ground portion 6 (this is ie the energy resulting from the transmission of solar radiation through the atmosphere and incident on the ground portion 6). During this step 16 of the atmosphere module, it is also determined 18 the diffuse illumination of the landscape, that is to say the radiation resulting from the diffusion 19 after interception of the radiation incident on the various atmospheric cells, increased by the thermal emission 20 of the atmospheric cells calculated previously. The direct solar illumination of the portion of soil 6 is determined by simulation of the part of the initial solar radiation not intercepted by atmospheric cells.

During the subsequent step 22, the re-emitted energy W _0U t (Ω _vj ) is determined, cell after cell, by diffusion by each cell of the soil portion 6 in each of the re-emission directions Ω _vj. For each cell, first of all at least one emission point Pe is calculated. One calculates on the one hand two emission points Pe-1, Pet for the turbid cells, and on the other hand, and simultaneously, during step 24, an emission point Pe for each surface element of any opaque cell 73. FIG. 9 represents the different stages of the calculation of the points Pe-i, Pet. From direct solar radiation from step 17, 25 is calculated for the first turbid cell 71, 72 encountered by each incident ray in the solar direction Ω _s , an intercepted energy from sector Wi _nt (Ω _s ) equal to l energy intercepted by the turbid cell, and this by applying the transmittance transfer function T (Ω>: W _int (Ω °) - = W _in (Ω °) [1- | Y (Ω °)] ^U ^ - ^Δli ], where U _f (i) is for example the leaf density of a vegetation cell and where Δli is the length of the trajectory of the direction radius Ω ° in cell i. Then we calculate 26 two points of origin scattering means M ° (Ωt) and Mg (Ω-l ') respectively associated with the scattering of radiation in a first upper hemisphere of space, and, in a second lower hemisphere. To do this, apply the formula ( II) above and we choose for example and preferably M ° (Ωt) = M _S) (Ω _vj (θ _v = 20 ⁰ )) and M ° (Ω * i) = M _{S (k} ( Ω _vj (θ _v = 160 ⁰ )). During this first iteration Ω _s = Ω °.

Since the turbid cell 71, 72 can receive several incident rays of direction Ω _s , two average origin points M _s (Ωt) and M _s (Ω * i) are calculated as energy barycenters of the points M _sk (Ωt) and respectively M _s (Ω * l), and this during the subsequent step 27. This calculation is performed with the iterative relationships mentioned above.

FIG. 7a schematically represents the determination of the mean origin point of diffusion M _sk (Ωt) and FIG. 7b represents an example of calculation of a mean origin point M _s (Ωt) with three incident rays of incident energy W * _n ( Ω _sl ), W * _nt (Ω _s2 ), and W _int (Ω _s3 ) respectively in the general case of three distinct incident directions Ω _sl) Ω _s2> Ω _s3 . In the first iteration, the same calculation is made but all the incident directions are confused with Ω _s . More generally, and during subsequent iterations, for each incident ray, the three mean origin points of diffusion M _sl (Ωt), M _s2 (Ωt) and M _s3 (Ωt) are determined. We make the barycenter of these three points to find the mean origin point M _s (Ωt). It should be noted that FIGS. 7a, 7b, are plotted in the plan, the cell being seen in elevation, but must be considered in space, the cells being volumes.

These different steps 25 to 27 are repeated successively for all the turbid cells 71, 72 of the soil portion 6 as shown diagrammatically by the arrow 28 in FIG. 5.

We then calculate 29 the scattered energy W _d i _î (Ω °, Ω _vj ) scattered at the mean origin point M _s (Ωt) and M _S (Ω - *), depending on whether the direction of re-emission considered Ω _VJ is rising or descending, from intercepted sector energy

calculated in step 25 above and by applying the precalculated initial diffusion transfer function T (Ω °, Ω _Vj ), that is: W _dffl (Ω; A _j ) ≈ W _ιnt (Ω;). T (Ω, Ω _j ).

We then determine which are the possible emission points Pet and Pei closest to the average origin points M _s (Ωt) and M _s (Ω * i) as shown diagrammatically in FIG. 6. Indeed, we have previously calculated, for each cell, a plurality of possible emission points Pe, regularly distributed in the cell. Also, the radiative module 10 has previously calculated and recorded the possible trajectories of the different rays coming from the possible incident directions Ω _sk , and the possible re-emission directions Ω _VJ *, from the different possible emission points Pe of the source cell (0,0,0) (initial calculation 14). In any cell with coordinates (X, Y, Z), any ray path from a possible emission point Pe is determined from the linked lists stored for the origin cell (0,0,0) by translation ( X, Y, Z) of the path (linked list) of a parallel ray coming from the corresponding point Pe in the origin cell (0,0,0). Thus, in step 31, the energy scattered by the turbid cell is calculated in each of the re-emission directions Ω _j , that is to say the energy W _0Ut (Ωy *) assuming that this energy is re-emitted from the Pe point which is either the Pet point or the Pe * point! depending on whether the direction Ω * is rising or falling.

FIG. 10 represents the calculation 24 of the emission point Pe in an opaque cell from the direct solar illumination calculated in step 17, in the first iteration. We determine 33 if the incident ray is or not intercepted by a geometric element of the cell. If this is the case, the energy intercepted by this geometric element and the coordinates of the corresponding point of intersection during step 34 are determined from the known optical properties of the geometric element. Then the barycenter is calculated. energy of the different intersection points corresponding to the different rays incident on the surface element in the cell, as well as the corresponding total intercepted energy. These steps 33 to 35 are repeated for each surface element of each opaque cell and on all the opaque cells, as shown diagrammatically by the arrow 36. In the case where no geometrical element intercepts the incident ray, as determined during the step 33, we immediately pass to the subsequent geometric element by iteration as represented by the loop 37 in FIG. 10.

Once these iterations have been carried out, for each cell which contains at least one geometric element (surface), 38 for each geometric element whose interception energy is not zero is calculated the possible emission point Pe which is closest to the barycenter of the interception points and which is external to any volume circumscribed by geometric elements.

The input module 2 made it possible to define parametric laws making it possible to determine the diffusions originating from the geometric elements from the energies intercepted by these geometric elements in the different directions. It is therefore sufficient to apply these parametric laws to determine the energy scattered W _d * _fï (Ω _Vj ) at the point of emission Pe by each geometric element in each direction of re-emission Ω _VJ *. We then calculate 31 the re-emitted energy W _out (Ω _vj ) corresponding to this energy diffused by application of the transmittance transfer function T (Q, j), if the cell contains turbid material.

Furthermore, from the atmospheric thermal emission 20 making it possible to obtain the diffuse illumination of the landscape 18 and from the thermal emission of the portion of soil which can also be calculated 39, one obtains at the end of the first iteration II of the process, the radiation re-emitted by each of the cells of the model 5 from the initial solar radiation. In this first iteration II, therefore, multiple diffusion phenomena are not taken into account. In fact, unlike multiple scattering, order 1 scattering is only influenced by W _in radiation (Ω _s ) monodirectional in direction (Ω °). In addition, the importance of multiple broadcasts is much lower than that of order 1 broadcasts.

The process then continues by successive iterations by calculating and memorizing at each iteration In:

- during an initialization step 40 of the iteration, for each cell, the intercepted radiation W; _nt (Ω _sk ) resulting from the interception of the incident radiation Wi „(Ω _sk ), knowing that the latter is the re-emitted radiation W _0Ut (Ω _VJ ) by the other adjacent cells as determined in the previous iteration,

- the radiation re-emitted W _out (Ω _Vj ) by each cell originating in particular from the scattering of the radiation intercepted on the cell.

The energy intercepted W _int (ΔΩ _M , i) by each sector is initialized to zero at each step 40 of iteration initialization.

Here again, at least one point of emission of the re-emitted radiation is calculated for each cell, during each iteration In. As for the first iteration, there is a calculation 41 for the turbid cells and a calculation 42 for the opaque cells. The calculation 41 of the emission points Pet and Pe-i for the turbid cells is explained in more detail in FIG. 11. First of all, during step 43, a radius is chosen on the turbid cell according to the incident direction Ω _sk and with an incident energy Wi „(Ω _sk ). During the subsequent step 44, it is determined which sector ΔΩ ^* vι _, i belongs to the incident direction Ω _sk . We then calculate 45, the intercepted energy Wi _nt (Ω _s ) of the incident ray by the turbid cell by application of the transmittance transfer function T (Ω _& ) either: W _mt (Ω _sk ) = W _in (Ω _sk ) [l- [r (Ω _sk )] ^{Uf (i)} - ^Δl1 ], where U (i) is for example the leaf density of a vegetation cell and where Δl; is the length of the trajectory of the direction radius Ω _sk in cell i. We then increment 46 the intercepted energy W ^* jι _t (ΔΩ _M> *) for each sector ΔΩ ^* vι, i by adding to it the intercepted energy W ^* _nt (Ω _sk ) by each incident ray.

One then calculates 47 two average points of origin of diffusion M _s (Ωt) and M _sk (Ω-1) associated respectively with the diffusion of the radiation according to a first upper hemisphere of space, and, according to a second lower hemisphere. To do this, we apply the formula (II) above and we choose

When the turbid cell receives several incident rays, two average origin points M _s (Ωt) and M _S (Ω-I ') are calculated as energy barycenters of the points M _sk (Ωt) and respectively M _s (Ω-J' ), and this during the subsequent step 48. This calculation is performed with the iterative relationships mentioned above.

An iteration 49 is then made on the different incident rays M _s on the turbid cell.

Then, in step 50, the scattered energies of sector W _d ig (ΔΩ _M) i, Ωv _j ) scattered by each of the M sectors ΔΩ ^* _M , i are calculated in the re-emission directions Ω _Vj by the average origin points M _s (Ωt) and M _s (Ωl) determined in step 48, using the sector broadcast transfer function

calculated according to formula (I) above. By applying this sector diffusion transfer function to the intercepted energy of sector W _int (ΔΩM, i), one obtains, for each re-emission direction Ω _Vj , the energy scattered from sector W ifi (ΔΩ _M; i, Ω _vj ), that is:

W _diff (ΔΩ _{M; i} , Ω _V j) = W _int (ΔΩ _M , _i ) .T (ΔΩ _M) i, Ω _VJ ) The sum of these sector scattered energies is then calculated to determine the scattered volume energy W _d jfï (ΔΩv _j ) by the turbid cell in each re-emission direction Ω _vjt, i.e. the sum of the different sector energies diffused by the different impact sectors ΔΩM.L is performed from step 51, the two emission points are determined

Pet and Pe-i- closest to the average origin points M _s (Ωt) and M _S (Ω> 1--). In whereas the volume energy diffused

is emitted from the predetermined emission point Pet or Pei, and according to the energy intercepted by the turbid material between this emission point and the point Pout (Ω _Vj ) of exit of the corresponding radiation from the turbid cell, it is determined with the transmittance transfer function T (Ω ^) the re-emitted energy W _out (Ω _vj ). One reiterates the calculation 50 of the energies diffused and re-emitted for all the directions of re-emission Ωy _j .

In step 52, the energy re-emitted _out (Ω _Vj ) by the turbid cell in the direction of re-emission Ω _{vj is thus determined} from the two emission points Pet and? El depending on whether the direction of re-emission Ω _Vj is rising or falling. This calculation is repeated on all the turbid cells of the soil portion 6 of the model 5.

For an opaque cell 73, the calculation 42 of the emission point Pe is similar to the calculation 24 described in relation to FIG. 20 for the first iteration II. After this calculation, the re-emitted energy W _out (Ω _vj ) by the opaque cell 73 is also calculated during step 52 with a parametric expression as described previously in step 31 of the first iteration II. The difference is that this time, the parametric expression is adapted to the case of a more or less isotropic incident radiation and no longer monodirectional according to the direction Ωs. The whole calculation is reiterated for all the cells of the portion of soil 6 and this iteration is represented by the loop 53 in FIG. 4. To pass from one iteration to the next, it is considered in the initialization step 40 that the energy re-emitted during an iteration In corresponds to the incident energy of the subsequent iteration In + 1. After a predetermined number n of iterations, depending on the desired precision, the radiative state of the cells of the soil portion 6 is known.

It is then determined 54 whether an atmospheric backscatter should or should not be taken into account. If so, the atmosphere module performs a step 55 shown in more detail in FIG. 12. In fact, this atmospheric backscatter calculation 55 may not be carried out at all the executions of the calculation of the radiative module 10. As shown diagrammatically in FIG. 5, atmospheric diffusion is taken into account during a first execution El. During a second execution E2 of the radiative calculation, the radiative module 10 calculates the radiation reflected by the ground portion 6 towards the atmosphere portion 7 (corresponding to the upper and intermediate atmosphere 13) from the incident radiation on the upper face of the cells of the lower atmosphere 11 as determined during the first execution E1. During this second execution E2, the step 55 for calculating atmospheric backscatter is therefore not executed. On the other hand, during a third subsequent execution E3, the radiation backscattered by the upper 13 and intermediate atmosphere 12 (step 55) to the lower atmosphere 11 (and to the ground portion 6) is determined from the determined reflected radiation. during the second execution E2. During a fourth execution E4, the radiative module 10 determines the radiation reflected by the ground portion 6 and the lower atmosphere 11 towards the upper 13 and intermediate atmosphere 12, from the backscattered radiation as determined during the third execution E3. During the fourth execution E4, the calculation 55 of atmospheric backscatter is therefore not executed.

During a fifth subsequent execution E5, the final radiation rising in the atmospheric column is determined - in particular at the level of the upper face of the scene, that is to say of the atmosphere portion 7 - after propagation through the atmosphere. This sequence of executions of the method according to the invention is shown diagrammatically in FIG. 5.

During each backscatter calculation 55, as shown in FIG. 12, the reflection of the radiation by the ground portion 6 and the lower atmosphere 11 is first determined according to the rising directions of re-emission, that is to say say belonging to the upper hemisphere.

We then calculate 57 the atmospheric backscatter as a convolution of the rising energy from the lower atmosphere 11 with a transfer function noted FT _BA - _BA - The latter indicates the backscattered energy WB0A ('-j' - * Ω _s .) by the atmosphere on each pixel (i 'j') from the top of the lower atmosphere 11 and for any downward direction Ω _v |, knowing that we have a rising unit energy

in a rising direction Ω _v - _f , from of a pixel (ij). Consequently, the transfer function is noted FT _BA - _BA (Δi, Δj, Ω _v t, Ω _s j,), where Δi = i'-i and Δj = j-j '. The expression of the convolution is: _B O _A O'J'Λ-J,)

The transfer function is precomputed by the optical module 8 by simulating the energy W _B o _A (i'j ^l * .Ω _s .) Backscattered by the atmosphere, in any pixel (i'j ¹ ) and for any direction falling Ω _s i knowing that the only energy source is a pixel (ij) from the top of the lower atmosphere 11 which sends energy W _B o _A (i _» Ω _v ) in an upward direction Ω _v f. This preliminary calculation is carried out iteratively on all the rising directions Ω _v i. The use of the FT _BA - _{BA function} greatly reduces the calculation times, because it avoids simulating the atmospheric backscatter for each energy W _B o _A (iJ _» Ω _v * r) from any pixel at the top of the lower atmosphere 11 , according to each of the rising directions Ω _vf . The computation time is therefore appreciably reduced by a factor equal to the number of pixels of the lower atmosphere 11, a term generally greater than 500, or even much more.

The reflection by the ground portion 6 and the lower atmosphere 11 of the atmospheric backscatter previously calculated during step 57 is then calculated 58. During the reflection calculation 58, only the uplink re-emission directions are taken into account. Finally, in step 59, the radiation transmitted at all levels and in particular at the top of the atmosphere is calculated, according to the rising re-emission directions, that is to say towards the sensor 60 (observer), which can be, for example, an airborne or satellite sensor. The energy transmitted from the top of the lower atmosphere 11 to the sensor 60 is calculated using a transfer function denoted FT _BA - _sensor (Δi, Δj, Ω _v î, Ω _s î). This function indicates the rising energy W _BOA (Ï 'J', Ω _SÎ ) in any pixel (i'j ') at the altitude of the sensor 60, in any rising direction Ω _s τ, knowing that the only source of energy is due to a pixel (ij) of the lower atmosphere 11 which sends energy W _B o _A (ij _> Ω _v î) in the rising direction Ω _v τ. Consequently, the energy received at the level of the sensor 60 is W _sensor (i * j ', Ω _s t) = ∑i ∑j ∑Ωvî _B0A (ij, Ω _v t). FT _BA -sensor (ii 'j-j', Ω _v t, Ω _s t) where W _BOA J-Ω _VÎ ) is the total rising energy from any pixel (ij) in the lower atmosphere 11. This energy is therefore the sum of two terms which represent the energies that the portion of soil 6 diffuses upwards during steps 56 and 58. The invention can be the subject of numerous variants compared to the embodiment described above. It is possible to produce not only remote sensing images, but also a three-dimensional radiative energy balance of the scene. The remote sensing images can be of the satellite type or others. The invention also makes it possible to obtain simulated digital representations of scenes other than terrestrial scenes.

The various functions mentioned above can be easily implemented by computer programming, that is to say by computer programs. In particular, the method according to the invention has been implemented and programmed in C language on an IBM Risc 6000® workstation in a Linux environment. It was possible to simulate images representing natural and urban landscapes with a calculation time of around 4 hours for complex landscapes. These calculation times can be less than 15 min, and even much less, for relatively simple landscapes.

Claims

CLAIMS 1 / - Method for obtaining a simulated digital representation of the radiative energy state of a heterogeneous three-dimensional scene (5), in which: - the scene (5) is divided into a plurality of cells (71,

72, 73) elementary comprising cells, called turbid cells (71, 72), having known radiative properties assumed to be homogeneous in the whole of each cell,

- transfer functions are stored for each turbid cell (71, 72) making it possible to determine, in each direction of a discrete plurality of N _v directions of space, said directions of re-emission Ω _v> predetermined in a frame of the space independent of the cells, an energy re-emitted by the turbid cell as a function of an intercepted radiation due to an incident radiation which it receives formed of at least one incident ray defined by an incident direction Ω _s belonging to a discrete plurality of Ns directions of space, called incident directions Ω _s> and an incident energy W _n (Ω _s ] _c ), a fraction of which, called intercepted energy

is intercepted by the turbid cell,

- for at least one portion, called the portion of the ground (6), of the scene, one re- _emits energy, cell after cell, W _0llt (Ω _vj ) by each turbid cell (71, 72) in each of the directions of re-emission Ωy _j , this re-emitted energy W _or (Ωyj) being considered as coming from a diffused volume energy W _d i _ff Ωy _j ) in the direction of re-emission Ω _yj by the whole of the cell turbid, this scattered volumetric energy Wdi Ωy _j) being calculated by considering that it is emitted totally from a single point in the turbid cell, known as equivalent origin point of diffusion M _{s> k} (Ωy _j ), characterized in that than :

- a plurality of M sectors, called incidence sectors ΔΩ- y, are defined in a reference mark of the space independent of the cells (71, 72) number lower than the incident directions Ω _s and containing all the N _s incident directions Ω _s ,

- to calculate, at least for the portion of soil (6), the volume energy scattered W ^ f (Ω _v j) in each re-emission direction Ω _v j by each turbid cell:

. we calculate, for each incident sector ΔΩJ ^ J, an energy, called intercepted energy of sector W _int (ΔΩM, i), equal to the sum of the intercepted energies i _nt (Ω _s k _; i) of the incident rays on the turbid cell in incident directions Ω _s belonging to the incident sector ΔΩM, - ;,. a function is used, called the sector diffusion transfer function T (ΔΩ ^* vι, i, Ω _Vj ) making it possible to calculate for each incidence sector ΔΩ ^* v [, i and each re-emission direction Ω _v j, an energy, called energy diffused from sector W ₍ iiff (ΔΩ-- f, i- _> Ωyj), diffused at the equivalent origin point of diffusion M _Sj k (Ω _v j) in the re-emission direction Ωyj from l energy intercepted from sector Wint (ΔΩ, i) from the incident sector ΔΩM, i,

. for each re-emission direction Ω _v j _? we calculate the volume energy scattered Wdj (Ω _v j) by the turbid cell as the sum of the sector scattered energies

scattered by each incidence sector ΔΩM ^ in this re-emission direction Ωyj. 21 - Method according to claim 1, characterized in that the sector transfer transfer function T (ΔΩM, i, Ωyj) is a weighted average of diffusion transfer functions T (Ω _s k, Ω _v j) stored for each incident direction Ω _s , and each re-emission direction Ω _v j.

3 / - Method according to one of claims 1 or 2, characterized in that the number M of incidence sectors is greater than 2 and ώ lower than 20. A1 - Method according to one of claims 1 to 3, characterized in that the number N _s of incident directions is greater than 100, in particular between 100 and 1000.

5 / - Method according to one of claims 1 to 4, characterized in that the possible incident directions Ω _s of the incident rays are determined in the same discrete plurality of directions of space as the directions of re-emission Ωy.

61 - Method according to one of claims 1 to 5, characterized in that the turbid cells (71, 72) being formed of polyhedra, for each incident ray in a direction Ω _s and which enters at a point Pin of one from its faces, called the input face, the equivalent origin point M _Sιk (Ω _v j) of the scattering of the incident ray by the turbid cell (71, 72) is determined in the re-emission direction Ω _v j, as the point located in the turbid cell, on the incident direction Ω _sk , and at a distance Δr _k (Ω _sk ) from the input face determined by the following formula II:

G (Ω _5k ). | Μ _vi | G (Ω _sk ) G (Ω _V *)

_{. .} G (Ω _v . | _Μsk | ₊ G ₍ Ω _sk) . | _Μvi r ^[1 - ^{eXp [} - ^Uf - ^{ΔIk (Ωsk)} - '^ - ⁽ M ⁺ ] μ, | ^{aV6C L} [l-exp [-G (Ω _sk ) .u _f .Δl _k (Ω _sk )]] ^]

where θ _sk is the angle formed by the incident direction Ω _sk with a predetermined fixed direction of space, called the reference direction, θ _yj is the angle formed by the re-emission direction Ω _v j with the direction of reference, μ _sk = cosθ _sk , μ _vj = cosθ _yj , and Δl (Ω _sk ) is the path of the incident ray through the cell, G (Ω _sk ) is a projection factor of the material perpendicular to the direction Ω _sk predetermined for each type of cell,

G (Ω _Vj ) is a projection factor of the material perpendicular to the direction Ω _VJ * predetermined for each type of cell. 7 / - Method according to claim 6, characterized in that one calculates and uses as an equivalent origin point of diffusion, for each incident ray, at least one average origin point of diffusion M _sk independent of the direction of re- emission Ω _vj , and the scattered energy of sector W _d - ^ ΔΩ ^ Ω _yj ) is calculated from each incident ray of incident direction Ω _sk as scattered at this point of origin origin of scattering M _sk .

8 / - Method according to claim 7, characterized in that one calculates and uses as an equivalent point of origin of diffusion, for each incident ray, two points of origin of means of diffusion: a point M _sk (Ωt) associated with a scattering of radiation in a first predetermined hemisphere of space, called the upper hemisphere, and a point M _sk (Ω-l--) associated with the scattering of radiation in a second predetermined hemisphere of space, called lower hemisphere, the two upper and lower hemispheres being chosen so that the reference direction is orthogonal to the diametral plane separating these two hemispheres.

91 - Method according to claim 8, in which the simulated digital representation comprises at least one simulated image seen from outside the scene (5) in a direction of observation, fixed and predetermined with respect to the scene, characterized in that that the direction of observation is oriented towards the upper hemisphere.

10 / - Method according to one of claims 8 or 9, characterized in that the reference direction orthogonal to the diametral plane separating the two upper and lower hemispheres is the vertical of a ground area of a planet, the scene ( 5) extending from said soil area of the planet.

11 / - Method according to claim 6 and one of claims 8 to 10, characterized in that one chooses

12 / - Method according to one of claims 8 to 11, characterized in that two average origin points M _s (Ωt) and M _s (Ω-i) are respectively associated with the diffusion in the rising directions belonging in said upper hemisphere and the descending directions respectively belonging to said lower hemisphere, these two mean origin points M _s (Ωt) and M _S (Ω4) being common and used as equivalent origin points of diffusion for all the rays incident on the turbid cell (71, 72). 13 / - Method according to claim 12, characterized in that the mean origin points M _s (Ωt) and M _S (Ω * I) are calculated as energetic barycenters of the mean origin scattering points M _sk (Ωt) and respectively M _sk (Ω -), calculated on the different rays incident on the turbid cell. 14 / - Method according to claim 13, characterized in that before calculating the re-emitted energies W _0Ut (Ω _v j) by each cell (71, 72,

73), each mean origin point M _s (Ωt) and M _S (Ω * I) is replaced by an equivalent emission point which is the closest Pe point chosen from among a plurality of possible Pe emission points regularly distributed in each cell, for which all the possible trajectories of the rays coming from all the points

Pe according to all possible re-emission directions Ω _vj have been previously calculated and recorded.

15 / - Method according to claim 14, characterized in that to determine the possible emission points Pe, a cell (71, 72, 73) is divided into a plurality of similar zones and the media of these zones are chosen by as long as possible emission points Pe.

16 / - Method according to one of claims 1 to 15, characterized in that one calculates by successive iterations a radiation re-emitted by each cell (71, 72, 73) comprising at least one ray re-emitted in at least a re-emission direction Ωyj, by calculating and memorizing at each iteration:

- the radiation incident on the cell (71, 72, 73) from the rays re-emitted by the other adjacent cells (71, 72, 73), as determined in the previous iteration, - the radiation re-emitted by the cell resulting in particular from the scattering of intercepted radiation due to radiation incident on the cell. 17 / - Method according to one of claims 1 to 16, characterized in that one defines beforehand in the scene (5) independently of the cells (71, 72, 73), geometric elements which represent different landscape elements, and in that the re-emitted energy W _out (Ω _v j) is calculated by taking into account these geometrical elements as mainly reflecting and / or absorbing elements.

18 / - Method according to one of claims 1 to 17, characterized in that the scene (5) comprising a portion of soil (6) of a planet, and a portion of atmosphere (7) extending au- above this portion of the ground up to a predetermined altitude, an atmospheric column (7) divided into cells, called atmospheric cells (72), is defined above the scene, having known radiative properties assumed to be homogeneous in the whole of each atmospheric cell (72), all atmospheric cells (72) of the same altitude having the same dimensions and radiative properties.

19 / - Method according to claim 18, characterized in that the dimensions of the atmospheric cells (72) increase with altitude.

20 / - Method according to one of claims 18 or 19, characterized in that the atmospheric column (7) comprises a lower zone, called the lower atmosphere, in contact with the scene and an upper zone, called the upper and intermediate atmosphere, and in that from a predetermined solar radiation defined by an incident ray on the upper face of each cell (72) of the upper layer of the atmospheric column (7), the following successive executions are carried out: 1) one determines the incident radiation on the upper face of the cells of the lower atmosphere from solar radiation after propagation in the upper and intermediate atmosphere,

2) the radiation reflected by the landscape area and the lower atmosphere in the upper and intermediate atmosphere is determined from the incident radiation determined during execution 1), 3) the radiation backscattered by the upper and intermediate atmosphere to the lower atmosphere is determined from the reflected radiation determined in execution 2),

4) the radiation reflected by the landscape zone and the lower atmosphere towards the upper and intermediate atmosphere is determined from the backscattered radiation determined during execution 3),

5) the final radiation rising in the atmospheric column is determined - notably at the level of its upper face - resulting from the sum of the reflected radiation determined in executions 2) and 4) after propagation through the upper and intermediate atmosphere