Abstract

The article is devoted to modelling the reflection of radiation from real surfaces. Most of the existing models consider only the near-surface processes, although in reality the processes occurring in the volume of the layer also have a fairly significant contribution to the final radiance. The authors propose a mathematical model that takes into account volume scattering. The model is based on the radiative transfer equation, the numerical solution of which makes it possible to find the reflection and transmission radiance factors. The paper describes the features of the implementation of the model in the case of a multilayer medium, a method is proposed that makes it possible to consider the effect of a randomly uneven surface. To validate the model, calculations of the radiance of the reflected radiation from the asphalt concrete pavement were carried out: for this, the corresponding parameters were selected that describe the optical properties of the medium. A comparison of the simulation results with the results of measuring the radiative characteristics of an asphalt concrete pavement sample was carried out, which showed that the model gives sufficiently reliable results not only qualitatively, but also quantitatively. The created model allows high-speed calculations of the radiative characteristics of various surfaces for different angles of incidence and observation, which can be used both in lighting calculations and in the formation of realistic images in computer graphics. road lighting metrics Radiance factors, radiative transfer equation, matrix-operator method, Monte-Carlo methods, GraphiCon 2021: 31st International Conference on Computer Graphics and Vision, September 27-30, 2021, Nizhny Novgorod, Russia ORCID: 0000-0002-9670-5413 (A. Basov); 0000-0001-9725-4266 (G. Boos); 0000-0003-4750-0160 (V. Budak); 0000-0002-1253-7687 (A. Grimailo)

1. Introduction

Spatial-angular characteristics of reflection are usually characterized by a bidirectional reflectance distribution function. This function allows calculating the radiance of a surface in a certain direction, created when radiation falls on this surface at a certain angle. Such a task is especially relevant in calculations, for example, in the street and architectural lighting, where the standardized value is exactly the luminance and not illuminance. The existing reflection models, which underlie most programs for performing lighting calculations, have several disadvantages and do not allow obtaining reliable results. The authors have developed and implemented a mathematical model of reflection, taking into account the processes occurring in the volume. To validate the model, it is necessary to select the input data, and compare the modelling results with the measurement results. 2. Boundary value problem of radiative transfer equation for a homogeneous layer

The model is based on the radiative transfer equation (RTE), the numerical solution of which allows finding the reflection and transmission radiance factors. The radiance factor, by definition, is the ratio of the radiance of the medium in a given direction to the radiance of an ideal diffuse surface under the

2021 Copyright for this paper by its authors.   μ  L(τ, ˆl)   L(τ, ˆl) τ τ0,μ0 Λ 4π  L(τ, ˆl) 

 L(τ, ˆl)x(ˆl, ˆl)dˆl,  δ(ˆl  ˆl0 ), L(τ, ˆl) ττ0 ,μ0  0; where L(,,) is radiance of the light field in the viewing direction ˆl   1  μ2 cos φ, 1  μ2 sin φ,μ

0 Boundary value problem ( 1 ) is defined in the Cartesian coordinate system OXYZ, in which the OZ axis is directed downward perpendicular to the layer boundary, zˆ is a unit vector along OZ. The upper boundary of the layer is z = 0. Unit vectors will be marked with a “^” sign, column vectors will be marked with a right arrow, row vectors with a left arrow, and matrices will be marked with a double arrow above the corresponding symbols.

For the numerical solution of the RTE, it is necessary to accurately discretize that means analytically isolate from the solution all the features and the strongly anisotropic part of the solution La (τ, ˆl) , i.e. present the solution in the form [1]:

L(τ, ˆl)  La (τ, ˆl)  L(τ, ˆl) , where L(τ, ˆl) is the regular part of the solution, for which the boundary value problem ( 1 ), taking into account ( 2 ), is transformed into the system [2]: the direction ˆl0   1  μ02 , 0,μ0, μ0  ( zˆ, ˆl0 )  cos0 : same radiation conditions. Let us consider the solution of the RTE boundary value problem for the case of irradiation of a layer of a turbid medium of optical thickness 0 by a flat monodirectional source in ( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 ) ( 6 ) ( 7 ) where the source function Δ(τ, ˆl) is defined:

dτ the solution of which [3] has the form: d C(τ)   BC(τ)  M1 Δ(τ),

B  M1(1  AW) , τ0  C(0)  eBτ0 C(τ0 )   eBτM1 Δ(τ,μ0 )dτ .

A smooth function L(τ, ˆl) can be represented in a finite basis using the method of discrete ordinates (MDO):

M L(τ,μi ,φ)   (2  δ0,m ) cos(mφ) Cm (τ,μi ) ,

m0 where C (τ)  Cm (τ,μi ), C(τ)  C (τ), C (τ)

T radiation expansion coefficients in a Fourier series in azimuth, μj  0,5(ζ j 1) , j are zeros of Gaussian quadrature of order N/2 for discretizing the scattering integral over the zenith angle. Because of the obviousness, we omit the superscript m everywhere.

Replacing the scattering integral with a finite sum, we transform the boundary value problem ( 3 ) into a boundary value problem for a matrix inhomogeneous linear differential equation of the first order with constant coefficients: is a column vector of discrete values of the  L(τ, ˆl) μ  τ L(τ, ˆl )  Δ(τ, ˆl)  μ τ0, μ0 dLa (τ, ˆl) dτ Λ 4π  0; L(τ, ˆl )  L(τ, ˆl)   x(ˆl, ˆl)L(τ, ˆl) dˆl  Δ(τ, ˆl)  La (τ, ˆl)

The conditionality of the matrix of the system decreases rapidly with the thickness of the medium. To eliminate this, we apply a scale transformation:

 u21   eΓτ0 u21 u22 

Expression ( 9 ) gives the connection between the flows leaving the layer and the incident ones and is a generalization of the radiance factors. In this case, the column F describes the intrinsic radiation of the layer, and the matrices R and T represent discrete values of the reflection and transmission radiance factors.

1 eΓτ0 u12  .

When implementing the model in the MATLAB environment, special attention should be paid to the transformation of the source function ( 4 ), considering the multilayer.

The problem in calculating ( 4 ) is the presence of the -singularity in the angle in the anisotropic part of the solution, which is the forward radiation attenuated according to Bouguer. To eliminate this, we select the direct radiation from the anisotropic part of the solution:

La (τ1  τ,μ, φ)  La (τ1  τ,μ, φ)  e(τ1τ) μ0 δ(ˆl  ˆl0 ) , where τ1 is the optical thickness of the layers above the layer τ.

Taking into account the discrepancy of the direct radiation containing the -function, the source function takes the form: ( 8 ) ( 9 ) ( 11 ) (12)

3. Source function for a multilayer medium

T2 α  F1  R1F2   F2  

As shown in [5], in the case of a multilayer medium, the solution is completely equivalent to the expression for one layer ( 9 ) in the form of scatterers, but with effective parameters, which reflects the invariance property of the solution: CC12    F1  T1 α  R2F1  F2    R1 TT21ααTR12 T1 R2 TT12ααTR21 T2  CC12  , ( 10 ) 1 where α  1  R2 R1  . where U is a matrix of eigenvectors of the matrix B ; Γ  diag(Γ, Γ ) is a matrix of eigenvalues  Γ  Γ and besides U1   u11 u12  , S  1 0  , H  eΓτ0 0 .

u21 u22  0 eΓτ0   0 1

To obtain a complete solution, let us solve equation ( 8 ) with respect to the fluxes leaving the layer, and obtain a solution in the form of scatterers [4]: Δ(τ1  τ,μ, φ)  μ dLa (τ1  τ,μ, φ) dτ

 La (τ1  τ,μ, φ)  Λ  x(ˆl, ˆl)La (τ1  τ,μ, φ)dˆl 

Λ e(τ1τ) μ0 x(ˆl0, ˆl) 4π 4π The anisotropic part of the solution from (12) can be expanded by the spherical harmonics method: La (τ1  τ,μ,φ)   k 2k  1 edk (τ1τ) μ0  e(τ1τ) μ0 Qkm (μ0 ) Qkm (μ) eimφ , (13) k0 mk 4π

The rest of the transformations (before integration in ( 8 )) are similar to the transformations for the source function for one layer, taking into account the replacement   1+

For ( 8 ), considering the new source function:

τ F(τ0 )  S 0 eΓtU1M1 Δ(τ1  t,μ0 )dt 

0 K 2k  1  km 4π

τ dk qk SeΓτ1 0 eΓt Zk (τ1  t)dt U1 1 / μ0  M1 Qkm 

0 0 Further transformations (14) lead to the final expression:

τ K 2k  1 Λ SeΓτ1 0 eΓtet μ0 dt U1 M1 1 / μ0  

xk qk Qkm km

4π

K 2k  1 F(τ0 )   km 4π

qk (edk τ1 μ0 Jk  eΓτ1J0 ) U1 1  μ0 M1 Qkm 

K 2k  1 Λ  xk qk edk τ1 μ0 J k U1 1  μ0 M1 Qkm

km 4π

To show the correctness of the solution and the operation of the program, a calculation was carried out for one layer with an optical depth equal to and for three layers with the optical thickness /4, /2 and /4 (Figure 1). (14) (15) 4. Formation of t he randomly rough boundary and account for slope correlation

Modelling of the randomly rough surface is a part of great importance of the whole model of the reeflctive surface. It inuflences both incident ray s and rays re-entering the outer medium. uThs, the randomly rough boundary has an immediate and intense impact on all processes in the scattering layer and hence the radiance factor.

Having been treated with a tool or aefcted by natural forces, facets of a real reeflctive surface are always correlated. In this regard, it appears to be highly substantial to account for slope correlation on the boundary “outer medium – scattering layer”. In [6] the authors showed that slope correlation taken into account can lead to the statistical lens emergence (the state of a wavy surface that acts like a real optical lens).

Currently, one can distinguish two approaches to the randomly rough surface formation. eTh rfist of them is the method of mathematical expectations [7]. iThs is a promising way which essence consists in modelling only points of rays’ intersection with the surface instead of modelling the surface itself. ehT method seems quite attractive since it provides a considerable decrease in computational time. Nevertheless, current development level of this method does not appear to sufice its implementation. ehT authors did not manage to nfid robust research on this way of a randomly rough surface modelling and particularly slope correlation account. iThs fact f orced the authors to pass to the spectral representation method.

Detailed description of this way is shown in [6]. eTh method is based on Fourier -transformation of a random efild and a correlation function. eThn realizations of the random surface can be go t by dint of Monte-Carlo Methods.

In case of an isotropic surface and correlation function

B(r)  exp(k0r) , (16) where k0 is constant and r2 = x2 + y2 one can obtain the point height Ψ(x, y) above the XOY-plane assuming DΨ(x, y) ≠ 1 as: Ψ( x, y)  σ M 2 ln αm cos k0 β2m ( x cos 2πγm  y sin 2πγm )  2παm  ,

M m1  1  β2m  where σ is standard deviation; α , α β, γ are independent and uniformly distributed in ( 0, 1 ). (17)

5. Description of the measuring system

To validate the created model, it was decided to compare the modelling results with the results of measurements.

Road safety depends on the luminance of the road surfaces created by street lighting. Asphalt concrete coatings, which are most often used as materials for highways, have both a diffuse and a specular component in the reflected radiation. Therefore, their reflection characteristics, coupled with the quality of the coating itself [8], are of great importance for road traffic.

Measurements of the luminance characteristics of asphalt concrete pavements are carried out in Russian Lighting Research Institute named after Sergey Vavilov (VNISI). The results are obtained on a measuring device created by the staff of the Institute. The setup is based on a goniophotometer, the kinematic measurement scheme of which is based on the B-β photometry system (Figure 2). To measure the spatial distribution of radiation reflected from the road surface sample (luminance), an illuminator with low beam divergence, fixed on a lyre, is introduced into the measuring circuit. To determine the luminance of the samples in a certain direction, created when light is incident on the sample at a certain angle, a luminance meter with a small measurement angle (1°) is used.

The test sample of asphalt concrete pavement in the form of a disk is installed so that its center is aligned with the photometric center of the installation. The sample is rotated relative to the vertical and horizontal axes of the installation. To switch from the angles of rotation of the installation Uh and Uv to the angles α and β, which determine the direction of observation on the road (Figure 3), the following formulas are used:

cos θ  cosUv cosUh , sin β 

sinUh 1  cosUv cosUh  2 where θ is the angle between the direction of observation and the normal to the sample, i.e. complementary to the angle α: α = 90o - θ. At Uh = 0 and Uv = 0 the angle β is taken equal to 0.

The change in the angle of incidence of light on the sample ε is performed by rotating the illuminator about the horizontal axis. In this case, the rotation of the illuminator does not depend on the rotation of the sample about the same axis. The radiation receiver is located at a certain distance from the photometric center of the setup along its photometric axis, which is perpendicular to the rotation axes.

6. Comparison of modelling results with measurement results

First, it was decided to compare the simulation results with the measurement results at the installation described in section 5. The measurements were carried out with a sample of asphalt pavement removed from the pedestrian zone. The sample lacked large crushed stone, which is typical for highways.

The measurement results are shown in Figure 4 and Figure 5. At normal incidence of light on the sample, the light is reflected symmetrically relative to the normal viewing direction (Figure 4). At small angles of incidence, light is reflected primarily in the direction coinciding with the direction of incidence (Figure 5). As the angle of incidence increases, a peak appears in a direction that is specular with respect to the direction of incidence of the light.

To calculate the light transfer in the sample under study and find its luminance, the model included the following parameters characterizing its optical properties: τ0 = ∞, Λ = 0.9, the coefficients for representing the scattering indicatrix in the form of the Henyey-Greenstein phase function g1 = 0.8, g2 = 0.3, a = 0.7 — the model of large absorbing balls is assumed — forward diffraction and significant backward reflection (scattering), which seems reasonable for asphalt and soil.

As can be seen in Figure 6, the modelling results coincide with the measurement results (taking into account the measurement error ~ 20 %). Thus, with correctly selected parameters describing the optical properties of the layer, the model created in the MATLAB environment makes it possible to calculate reflection and transmission from real surfaces at a sufficiently high speed (~ 0.2 s).

It was also decided to simulate the reflection from the snow and compare it with the measurement results [9]. In [9] luminance factors are given, but the developed model considers the luminance. Therefore, for comparison with the results of modelling, the values must be multiplied by the cosine of the angle of incidence and divided by π. Information about the optical properties of snow was taken from [10]. The results of modelling and measurements were close (within 10 %), taking into account the spread of measurements and variations in the albedo of a single scattering (Figure 7).

In the future, it is planned to expand the model: to solve the problem by the Monte Carlo method, taking into account polarization, which will allow for comparison with other measurement results [11].

7. Conclusion

The proposed model has shown its suitability for modelling reflection from real objects. In the future, it is planned to carry out a comparison with a large number of measurements of the luminance coefficients on different samples and for different angles of incidence and observation. Perhaps in order to obtain correct results of modelling the luminance coefficients of road surfaces for a small observation angle of 1°, it will be necessary to make some changes in the parameters: if for such an almost opaque formation the role of reflection from the substrate is small, then the reflection from a randomly uneven upper boundary can strongly affect the results.

8. Acknowledgements

The authors are grateful to V. M. Pyatigorskiy and A. A. Korobko from Russian Lighting Research Institute named after Sergey Vavilov (VNISI) for providing measuring system and samples of asphalt concrete pavements for measurements. 9. References

References

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

[2] V.P. Budak , D.A. Klyuykov , S.V. Korkin , 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 , pp. 147 - 204 .

[3] P.J. Flatau , G.L. Stephens , On the Fundamental Solution of the Radiative Transfer Equation , JGR 93 ( 9 ) ( 1988 ) 11037 - 11050 .

[4] G.J. Fryer , L.N. Frazer , Seismic waves in stratified anisotropic media, Geophys . J. R. Astr . Soc. 78 ( 1984 ) 691 - 698 .

[5] A.Y. Basov , V.P. Budak , Mathematical model of a surface radiance factor , CEUR Workshop Proceedings , 2019 , Vol. 2485 , pp. 43 - 47 . doi: 10 .30987/graphicon-2019-2- 43 -47.

[6] V.P. Budak , A.V. Grimailo , Light reflection from real surfaces: Probabilistic model of the layer radiance factor , CEUR Workshop Proceedings , Vol. 2744 , 2020 . doi: 10 .51130/graphicon-2020- 2-4-37.

[7] B.A. Kargin , K.B. Rakimgulov , A weighting Monte Carlo method for modelling the optical radiation fi eld in the ocean-atmosphere system , Russ. J. Numer. Anal. Math. Model 7 ( 3 ) ( 1992 ) 221 - 240 .

[8] H. Bowden , M.J. Almond , W. Hayes , C. Browne , S. McRobbie , The use of diffuse reflectance infrared spectroscopy to monitor the oxidation of UV irradiated and naturally aged bitumen and asphalt , Road Materials and Pavement Design 22 ( 6 ) ( 2021 ) 1254 - 1267 . doi: 10 .1080/14680629. 2019 .1684982

[9] Yu. Lv , Z. Sun , The reflectance and negative polarization of light scattered from snow surfaces with different grain size in backward direction , Journal of Quantitative Spectroscopy and Radiative Transfer 133 ( 2014 ) 472 - 481 .

[10] Yu.M. Timofeev , A.V. Vasiliev , Osnovy teoreticheskoj atmosfernoj optiki [Fundamentals of theoretical atmospheric optics] , St Petersburg University publ., 2007 , 152 p.

[11] Z. Sun , J. Zhang , Z. Tong , Yu. Zhao, Particle size effects on the reflectance and negative polarization of light backscattered from natural surface particulate medium: Soil and sand , Journal of Quantitative Spectroscopy and Radiative Transfer 133 ( 2014 ) 1 - 12 .