Mathematical Model of a Surface Radiance Factor

Abstract

The article is devoted to the creation of a surface radiance factor mathematical model. The basis of the model is the solution of the boundary value problem of the radiative transfer equation (RTE). The surface is considered as a structure consisting of several turbid layers, each of which is characterized by its optical parameters. The top of the structure is randomly rough, uncorrelated, Fresnel. The lower boundary reflects perfectly diffusely. The complexity of solving the RTE boundary value problem for real layers is due to the fact that the suspended particles in each layer are always much longer than the wavelength. This leads to a strong anisotropy of the radiance angular distribution according to Mie theory. The solution comes down to a system of equations by the discrete ordinates method that consists of several hundred of differential equations. Subtraction of the anisotropic part from the solution based on an approximate analytical solution of the RTE allows avoiding this problem. The approximation is based on a slight decrease in the anisotropic part of the angular spectrum. The matrix-operator method determines the general solution for a complex multilayer structure. The calculation speed can be increased without compromising the accuracy of the solution with the help of the synthetic iterations method. The method consists of two stages: the first one repeats the described one with a small number of ordinates; on the second one the iteration of it is performed. The model is realised in the Matlab software.

Keywords

• radiative transfer equation
• matrix-operator method
• synthetic iterations method

1. Introduction

Reflective properties of surfaces play the key role in lighting calculations, and in most cases, calculations are impossible without knowing them. Some objects reflect both diffusely and specularly. Taking into account such surfaces is essential for correct calculations. Radiance factor is the classical value used to describe spatial reflective properties. In [14] a model of radiance factor of a uniform layer with diffuse bottom and Fresnel upper boundary was described. The present article is aimed to find the solution for a multi-layered structure and to optimize calculations (more fast calculations without lowering in accuracy).

2. Boundary value problem of the RTE for a uniform layer

A partially coherent wave reflection from a structure with a complex optical characteristics analysis shows [4] that the wave is dived into two components, after solution averaging: a coherent one determining the optical characteristics and a quasiuniform admitting the beam description and obeys the RTE. If the irradiating inhomogeneities of the medium are located in the Fraunhofer region from each other, which is most often realized in practice, then the optical characteristics of the medium are determined by summing them over the elementary volume of the medium.

Radiance factor is the relation between the radiance coming out from a surface in the specified direction and the radiance of an ideal diffuse surface in the same conditions. Therefore, the determination of the surface radiance factor is reduced to a boundary value problem for a flat layer. Let us analyze the algorithm for the numerical solution of the RTE boundary value problem for the case of irradiation of a turbid layer with optical thickness 0 by a plane unidirectional source (PU) in the direction ˆl0   1  02 ,0,0 , 0  (zˆ, ˆl0 )  cos 0 :  L(, ˆl)    L(, ˆl)    4  L(, ˆl) 

 L(, ˆl)x(ˆl, ˆl)dˆl, 0,0  (ˆl  ˆl0 ), L(, ˆl) 0,0  0; ( 1 ) where L(,,) is the radiance in the viewing direction ˆl   1  2 cos , 1  2 sin , ,   (zˆ, ˆl) at optical depth z    ()d  ; x(ˆl, ˆl) is indicatrix of scattering;  is single 0 scattering albedo. The problem ( 1 ) is defined in the coordinate system OXYZ, the axis OZ is directed perpendicularly down, zˆ is the unit vector along OZ. z = 0 on the upper boundary.

Numerical solution implies discretization of the RTE. The scattering integral should be replaced by a finite sum [2]. It is impossible, since there are singularities in the angular radiance distribution, which cannot be replaced by a finite series in any basis. The article [9] offers to single out the direct source radiation, which became an indispensable element of all RTE solving methods. However, a feature of all natural objects is the presence of suspended particles, which sizes significantly exceed the wavelength, which, in accordance with Mie Theory, leads to high angle scattering anisotropy, and the extraction of direct radiation is not enough efficient. In case of presence of strong anisotropy in the radiance angular distribution, replacing the integral with a finite sum can lead to significant uncertainties [15].

3. RTE discretization

The idea is of the accurate discretization is to analytically distinguish all the singularities and the anisotropic part La (, ˆl) , i.e. to represent the solution in the form [1]:

L(, ˆl)  La (, ˆl)  L(, ˆl) , here L(, ˆl) is the regular part of the solution (RPS), which is a smooth function. Such a smooth function can be represented by a finite basis of elements.

Considering ( 2 ) the boundary value problem ( 1 ) for L(, ˆl) is [2]:  L(, ˆl)    L(, ˆl )    4  L(, ˆl) 

 x(ˆl, ˆl)L(, ˆl) dˆl  (, ˆl);  0; L(, ˆl ) 0, 0 0 , 0  La (, ˆl), where the source function (, ˆl) is defined by the anisotropic part of the solution: (, ˆl)   dLa (, ˆl)  La (, ˆl) 

 d 4

As a numerical method implies finding an approximate value, it cannot exactly correspond to boundary conditions. That is why the second boundary condition is changed.  x(ˆl, ˆl)La (, ˆl)dˆl . ( 4 ) ( 2 ) ( 3 ) tion ( 8 ) with respect to flows exiting the layer. The solution is get in the form of scatterers [6]:  C (0)   F   R T   C (0)  , C (0 ) F  T R C (0 ) ( 9 ) where F   hJ, R T   h  u12 F  T R  e0 u22  e0 u11  ,  u21     u11 h    e0 u21

1 e0 u12  .

The expression ( 9 ) in the form of scatterers gives us the relation between the streams emerging from the layer and the incident ones and is a generalization of the radiance factor. The column F describes the inner radiation of the layer, the matrices R and T represent discrete values of the radiance factors for reflection and transmission.

5. Invariance of the solution

In the simplest case, a multilayer structure can be represented by only two layers: ( 10 ) ( 12 ) C1  F1   R1 T1  C1  C1   F1   T1+ R1+  C1  , CC22   FF22   RT22 RT22  CC22  , where the subscript defines whether the upper 1 or lower 2 belongs to. Vertical arrows indicate the direction of radiation incidence on the layer.

The radiation emerging from the first layer is the radiation entering the second layer and vice versa, that is:

C1  C2  C, C2  C1  C . ( 11 )

Solving the system with respect to the reflected radiation C1 and transmitted radiation C2 with respect to the incident radiations on the system, we obtain:

Since L(, ˆl) is a smooth function, it can be represented in a finite basis. For example, if the RTE is discretized by the discrete ordinate method (DOM), the representation is as follows:

M L(,i ,)   (2  0,m ) cos(m) Cm (,i ) , m0 ( 5 ) T where C ()  Cm (,i ), C()  C (), C () is a column vector of discrete values of the azimuthal expansion coefficients of radiance in a Forier series, j  0,5( j  1) , j are zeros of Gaussian quadrature of the order N/2 for zenithal discretization of the scattering integral. Index m is further absent due to the lack of need for it.

This allows to replace the scattering integral by a finite sum. The boundary value problem transforms to a boundary value problem for a matrix inhomogeneous linear differential equation of the first order with constant coefficients: d C()

d here wj are the weights of the Gaussian quadrature of the zen  BC()  M1 (), B  M1(1  AW) , ( 6 ) ithal order N/2, Pln () are associated polynomials Legendre, are polynomials

Legendre,  wi 0  , M  i 0  K

   , A  (2k  1) Pkm (i )xk Pkm (j ), 4  0 wi   0 i  k 0 Pl0 ()  Pl () W 

K x()  (2k 1)xk Pk (cos ) .

k 0

4. Propagator and scatterer

The solution of the matrix equation ( 6 ) can be presented [5] as the sum of the general solution of the homogeneous equation and the particular solution of the inhomogeneous: 0  C(0)  P(0, 0 ) C(0 )   P(0, ) M1 (,0 )d , ( 7 ) 0 where the function P(t, )  eB(t) is the solution of the homogeneous equation describing the relationship of the radiance distributions at two points t and  of the medium without internal sources. Such a function is called propagator.

There are problems with the solving of ( 7 ), because the propagator contains both positive and negative exponent functions. This physically corresponds to flows propagating top-tobottom and bottom-top. Thus, the matrix is poorly conditioned, and the calculations for fields with >1 become impossible. Scale conversion is proposed to avoid this effect [7]: 0 SU1 C(0)  HU1 C(0 )  J, J  S  etU1M1 (t)dt , ( 8 ) 0 where

U is the eigenvector matrix of the matrix B ;   diag( ,  ) is the eigenvalues matrix,    ; U1  u11 u12  , S  1 0  , H  e0 0 .

u21 u22  0 e0   0 1

The boundary value problem ( 1 ) and the corresponding system ( 6 ) are two-point problems, as the conditions on the boundaries are given, and the solution inside the layer needs to be found. Therefore, the solution ( 7 ) based on propagators is not complete [11].

Column vectors C (0), C (0 ) in ( 8 ) describe the flows falling on the layer and are defined by the boundary conditions. The vectors C (0), C (0 ) correspond to the radiation flows reflected from and transmitted the layer. Let us solve the equaC1    F1  T1  R 2F1  F2    C2  T2  F1  R1F2   F2     R1  T1  R 2 T1

T2  T1

T1  T2  C12  ,

R 2  T2  R1 T2  C  1 where   1  R2 R1  .

The expression ( 12 ) for two adjacent layers is completely equivalent to the expression for one layer ( 9 ) in the form of scatterers, but with effective parameters (coefficient) that can be obtained from the parameters of each layer separately. This shows the invariance property of the solution of scatterers, which allows us to logically go over to the invariant immersion of V.A. Ambartsymian [13].

On the other hand, invariance allows the calculation of a layer with an arbitrary vertical inhomogeneity, breaking it into an arbitrary number of homogeneous layers. In this case, two adjacent layers can be replaced by one layer described by expression ( 12 ). This approach in transport theory is called the matrix operator method [10] (MOM). The advantage of the obtained expression ( 12 ) is the allocation of the anisotropic part of the solution in arbitrary form ( 2 ).

6. Anisotropic part of the solution

If we talk about the complexity of calculations and, accordingly, the speed of calculations using some software, the key role for this has the sizes of the matrices included in the matrix solution ( 9 ) and ( 12 ), i.e. the values of the constants N, M, K, where N is the number of discrete ordinates, M is the number of azimuthal harmonics, K is the number of members of the decomposition of the indicatrix into polynomials Legendre.

In the general case of an arbitrary incidence angle, the values are approximately equal: M ≈ N ≈ K [2]. However, with a successful choice of the anisotropic part, it can be achieved that the angular dependence of the RPS will be close to isotropic. Thus, the calculation speed can be significantly increased because K >> N >> M.

The isotropic part was distinguished. However, it should be defined too. How can it be defined? The unequivocal answer in the spatial-angular representation is difficult, but the task is greatly simplified for the spectral representation of the angular distribution. The more anisotropic the angular distribution is, the more smooth is its spectrum. The radiance distribution can be represented as polynomial Legendre series:

 2k  1 La (,)   Zk () Pk () , ( 13 )

k 1 4 where the assumption is made that the anisotropy in the region of small angles is much stronger than the azimuthal asymmetry [2],   (ˆl, ˆl0 ) .

The spectrum Zk() of the anisotropic part of the solution slowly monotonically decreases from the index k. This allows us entering a continuous function Z(k,). Since the function slowly and monotonously decreases, the following expansion in Taylor series is valid:

Z (k  1, )

Zk ()  Z ( k, ) k

. kk

Substitution of ( 14 ) into an infinite system of differential equations for solving the RTE by the spherical harmonics method leads [3] to one equation of mathematical physics that allows an analytical solution

Zk ()  exp  (1  xk )  0  , which is called the small-angle modification of the spherical harmonics method (SHM).

In [8], a comparison was performed of the solution ( 9 ) with the separation of the anisotropic part based on the SHM (the MDOM program) with the main known programs: MDOM with the same accuracy in calculating the angular distribution of radiance exceeds other programs by 1-2 orders of magnitude in computational speed.

7. Reflection and refraction at the interface of two media

Let us consider a special case of a layer, in which the lower boundary reflects according to Lambert’s law and has a reflectance In this case, the boundary condition in ( 1 ) on the bottom is the following: ( 14 ) ( 15 ) (16) L(0 , ˆl) 0  

 L(0, ˆl)dˆl .

 (0)

The following matrix expressions can be obtained by substituting the azimuthal representation of radiance ( 5 ) and discrete ordinates, and the integral (16) by Gaussian quadrature: m  0 : C (0 )  2 R LC (0 ); m  0 : C (0 )  0 , (17) where the matrix of Lambert reflection R L consists of the same N/2 lines j wj .

In accordance with MOM ( 12 ), we obtain the following expression for the reflected component of zero harmonic: 1 C (0)  F  2T 1  2R LR  R LF , (18) where the index 1 related to the layer is omitted due to the absence of need. All other azimuthal harmonics m>0 are determined by the expression of a single-layer medium ( 9 ).

This approach does not work on the boundary with refraction, since the directions of the rays vary according to the Snell’s law:

n1 sin 1  n2 sin 2 , where n1, n2 are refractive indices of the medium, and the correspondence of ordinate directions is violated. A solution to this problem was proposed in [12].

Let us consider in more detail the refraction on the practically important case when n1 < n2. The first medium is called atmosphere for certainty (index a), na=1, and the second media is the ocean (index o) with no>1. The cosines of the rays with the axis OZ in both media according to (19) will be related to each other by the expression:

a  1  no2 (1  o2 ) .

An important feature of the problem and its further solution is the presence of a region of total internal reflection. It is seen from (20) that this region appears for o  t  1 1 no2 im the ocean medium: the rays do not exit the ocean, but are ideally reflected again into the ocean. Boundary conditions for the total internal reflection region are formulated without problems. Then the scattering integral can be represented as the sum of three integrals taking into account the total internal reflection region: 1 t  Qkm ()Cm (,)d   Qkm ()C m (,)d 1 1

t 1   Qkm ()C m (,)d   Qmk ()C m (,)d.

t t

The first and last integrals are related to the region of refraction, and the second relates to the total internal reflection region. For the second integral, we perform the transformation: t 1  Qkm ()Cm (,)d   Qmk ()Cm (,)d,   t , (22) 1 t which makes it possible to apply a double Gaussian quadrature with Nt nodes and subsequently move in this zone to two streams of ordinates Ct , Ct , which are connected by an ideal mirror reflection at the boundary.

For the first and last integrals in (21), we make the transformation of the integration variable to a by expression (20): o  1  (1  a2 ) no2 , do  ada no2  (1  a2 ) . (23) It is obvious that in the transition to discrete ordinates a complete correspondence is established between atmospheric ordinates Сa , Сa and ocean ordinates Сo , Сo . If we combine the vectors into one in accordance with the rules of the Matlab software Сocn  Сt ; Сo  , Сocn  Сo ; Сt  , then all the rela      tion in MOM will be valid for them. The introduced values also allow us to write down the condition at the atmosphere-ocean boundary: (19) (20) (21)  Ca    R Tao   Ca  , Cocn  Toa Roo  Cocn  (24)  0   0 1 where Tao  [T 0], Toa    , Roo    , R, T are Fresnel T R 0 reflective matrices.

8. Synthetic iterations method

Fig. 1 shows a comparison of the calculations in MDOM of radiance angular distribution reflected by a layer for different sets of parameters N and M. It is easy to see that the calculation of the almost complete distribution of the viewing angle is much faster than the calculation of individual small sharp peaks. The difference in calculation time t is more than 150 times. What is the reason?

The angular distribution of the RPS is actually close to isotropic, but with some small ripples. A fair question arises: how many discrete ordinates N are necessary to represent this small ripple? Since the RPS is a smooth function, its expansion into a series of polynomials Legendre has a finite number N: which corresponds to the Lagrange interpolation formula for the function L(,).

The latter relations are the analogue of the Nyquist-Shannon theorem on samples for the angular spectrum of the angular distribution over polynomials Legendre. Hence: 1. MDOM method provides average convergence; 2. All methods for isolating the anisotropic part are equivalent to each other in a uniform metric; 3. To achieve good convergence in a uniform metric, the sampling interval should correspond to the angular size of the smallest part, which should be reproduced on the radiance distribution.

When implementing a multilayer surface model taking into account diffuse reflection from the lower boundary and reflection from the upper boundary according to the Fresnel law in the Matlab software, the question arose of possible acceleration of calculations without loss of quality. The number of discrete ordinates N determines the size of the matrices with which the calculations are performed. That is, a decrease in N would speed up the calculations. Thus, the main question is: how not to lose in quality? The synthetic iterations method is the answer to this question.

The synthetic iterations (SI) method was proposed in nuclear physics [1]. In this case, the iteration splits into two stages. At the first stage, an approximate solution is sought that converges well in the average energy metric, and at the second, the usual iteration is performed, which significantly increases the convergence in a uniform angle metric. Since the developed method for solving MDOM has good convergence in the average metric, we should count on its significant increase in convergence after iteration.

A numerical comparison of the reflected radiance in the first iteration with MDOM is presented in Fig. 2, where t is the computation time.

a) full range of viewing angles b) range in the vicinity of gloria, the designations are the same as in Fig. 2a Fig. 2. Comparison of synthetic iteration (SI) with MDOM for the radiance of the reflected radiation of a layer

It is seen from the figure that the greatest difficulty for calculating in MDOM is the region near gloria (Fig. 2b). To calculate this region in the MDOM program, N = 801 and M = 256 are required, which corresponds to a sampling step of less than 0.5°. To achieve the same accuracy within the framework of synthetic iteration, only N = 11, M = 4 is necessary, which reduces the computation time by almost 60 times. Accordingly, the synthetic iteration from MDOM allows us to calculate the angular distribution of radiance with an accuracy in the uniform metric of no worse than 1% at a counting time of no more than 1 second.

This method allowed to significantly increase the speed of calculations carried out in the Matlab software using the created model.

9. Conclusion

The mathematical model of the luminance factor for a multilayer medium bounded above by the Fresnel and below Lambert surfaces is implemented. The calculation is optimized in terms of speed and accuracy of calculation. The obtained dependencies and characteristics qualitatively coincide with the expected ones. The model must be filled with parameters of real media, for which it is planned to experimentally test the model. In the future, it is also planned to take into account reflection from a randomly uneven border and polarization.

10. References

References

[1] Adams M.L. , Larsen E.W. Fast iterative methods for discrete-ordinates particle transport calculations // Progress in Nuclear Energy, 2002 . Vol. 40 , No.1. P.3- 159 .

[2] Budak V.P. , Klyuykov D.A. , Korkin S .V. Convergence acceleration of radiative transfer equation solution at strongly anisotropic scattering // In Light Scattering Reviews 5 . Single Light Scattering and Radiative Transfer / Ed. A.A. Kokhanovsky. Springer Praxis Books, 2010 . P. 147 - 204 .

[3] Budak V.P. , Korkin S .V. On the solution of a vectorial radiative transfer equation in an arbitrary threedimensional turbid medium with anisotropic scattering / / J. Quant. Spectrosc. Radial. Transfer., 2008 . Vol. 109 . P. 220 - 234 .

[4] Budak V.P. , Veklenko B.A. Boson peak, flickering noise, backscattering processes and radiative transfer in random media // J. Quant. Spectrosc. Radial. Transfer., 2011 . Vol. 112 . P. 864 - 875 .

[5] Flatau P.J. , Stephens G.L. On the Fundamental Solution of the Radiative Transfer Equation / / JGR, 1988 . V. 93 , No.D9. P. 11 , 037 - 11 , 050 .

[6] Fryer G.J. , Frazer L.N. Seismic waves in stratified anisotropic media // Geophys . J. R. Astr . Soc., 1984 . V.78, P. 691 - 698 .

[7] Karp A.H. , Greenstadt J. , Fillmore J.A. Radiative transfer through an arbitrarily thick , scattering atmosphere // J. Quant. Spectrosc. Radial. Transfer., 1980 . Vol. 24 . P. 391 - 406 .

[8] Kokhanovsky A.A. et al. Benchmark results in vector atmospheric radiative transfer / / J. Quant. Spectrosc. Radial. Transfer., 2010 . Vol. 111 . P. 1931 - 1946 .

[9] Milne E.A. The reflection effect of the eclipse binaries // Mon . Not. Roy. Astrophys. Soc. , 1926 . Vol. LXXXVII. P. 43 - 49 .

[10] Plass G.N. , Kattawar G.W. , Catchings F.E. Matrix Operator Theory of Radiative Transfer // Appl. Opt., 1973 . Vol. 12 . P. 314 - 326 .

[11] Redheffer R.M. Inequalities for a matrix Riccati Equation // Journal of Math. and Mechanics , 1959 . V.8. P. 349 - 367

[12] Tanaka M. , Nakajima T. Effects of oceanic turbidity and index refraction of hydrosols on the flux of solar radiation in the atmosphere-ocean system / / J. Quant. Spectrosc. Radial. Transfer., 1977 . Vol. 18 , No1. P. 93 - 111

[13] Ambartsumian V.A. K zadache o diffuznom otrazhenii sveta [To the problem of diffuse reflection of light] // JETF, 1943 . Vol. 132 , No. 9 - 10 , P. 323 - 334 (in Russian).

[14] Basov A.Yu. , Budak V.P. Model ' rasseivayushchego sloya s diffuznoj podlozhkoj i frenelevskoj granicej [Model of a scattering layer with a diffuse bottom and a Fresnel boundary] // GraphiCon 2018. Conference proceedings . Tomsk, Russia. P. 399 - 401 (in Russian).

[15] Krylov V.I. Priblizhennoe vychislenie integralov [Approximate calculation of integrals] // Nauka Publ., Moscow, 1967 (in Russian).