Prof Zbyněk Malenovský

Dr Nicolas Lauret

Eric Chavanon

Prof Jordan Guilleux

Prof Abdelaziz

Omar Ronan Paugam

Yingjie Regaieg

Xuebo Wang

Zhu Yang

Paul Boitard Tao

Madam Peiying

Tianqi Cang

Frédéric Baret

J Gastellu-Etchegorry

Y Wang

O Regaieg

T Yin

Z Malenovsky

Z Zhen

X Yang

Z Tao

As the creator of the Discrete Anisotropic Radiative Transfer (DART) model, he gives me the greatest guidance on the usage and development of DART. Despite his very busy days, he still patiently teaches me the meaning of each parameter in the DART model and the relevant mathematical and physical formulas. He carefully writes each formula down on the paper in case I cannot follow him. He holds weekly seminars to identify and solve our problems as quickly as possible. I am deeply touched by his insistence and dedication to research. He is always very enthusiastic about helping others. He keeps running the DART model Google forum (https://groups.google.com/forum/#!forum/dart-cesbio), where people from every corner of the world have the opportunity to discuss with others and share their questions and answers efficiently.

In addition to the scientific research, he always shows me the best food and culture (i.e., language, music and so on) of France with pride and provides the greatest help for my life in France. All of these are the most precious treasures of my life that I will never forget. I hope the DART model achieves more success under his leadership in the future. I also thank my co-supervisor, Professor Shengbo Chen. Under his guidance, I have actively participated in research projects and accumulated solid basic knowledge. In life, I am

of leaves in urban areas is achieved by introducing an innovative calibration of DART. Finally, we analyze the vegetation isolines behaviours and propose the intersection point right shift phenomenon based on the DART simulation data to mitigate the soil-adjusted vegetation index (SAVI) saturation effect in high FC areas.

In order to achieve a 3D radiative transfer simulation capability in the continuous-time phase, a static 3D maize modelling model constructed from the extended L-system (ELSYS) is coupled with a dynamic maize growth equation using degree days as the growth factor and development rules describing canopy structure throughout the plant growing season from seedling emergence to the male flowering stage. Maize canopy reflectance for the corresponding scenes is simulated using various 3D radiative transfer models for cross-validation, and a good agreement among them is achieved. The simulated canopy reflectances from the 3D radiative transfer models are compared with the 1D radiative transfer model under the same leaf area index (LAI) conditions as true values.

Results show that the homogeneity assumptions of the 1D radiative transfer model under the same LAI conditions resulted in an approximately 1.5 times overestimation of FC, leading to a significant overestimation of the reflectance in the Near-InfraRed (NIR) in the nadir direction. This overestimation may be because the homogeneity assumption fails to consider that the vegetation seen at the nadir is the smallest part of the scene and that multiple scattering is mainly from vegetation. Considering that the nadir direction is a vital observation direction for the sensor, and the NIR is an essential band for vegetation monitoring, we believe that the overestimation due to the homogeneity assumption is not negligible and result in an underestimation of the amount of vegetation (LAI, etc.) in the remote sensing inversions based on the homogeneity assumption.

In order to invert the optical properties of leaves in urban areas, the DART calibration is introduced. DART calibration first separates the single scattered reflectance of the canopy using a linear spectral mixture model and iteratively corrects the optical properties of the input leaves to simulate the canopy reflectance close to the mixed pixel separation. The inversion results are evaluated for accuracy assessment and sensitivity analysis. The mean relative errors of the ideal noise-free simulation experiment are 0.013, 0.005, 0.027, 0.297, and 0.250 for ground, roof, water, tree, and shrub in all bands. Under noise interference experimental conditions (pixel shifting, geometric accuracy of the 3D scene and modulation transfer function with some deviation from the true value), significant errors in the inversion are observed: for ground, roof, water, trees and shrubs, the mean relative errors are 0.233, 0.507, 3.088, 0.834, and 1.256 respectively. The parameters that have the most significant influence on the inversion accuracy of the optical properties of urban matter are the SZA, the spatial resolution of the satellite images, the pixel shifting, the inaccuracy of the 3D urban scene modelling and the modulation transfer function, in descending order of importance.

In order to alleviate the saturation effect of SAVI in areas with high FC, the vegetation isolines are analyzed, and a right-shift phenomenon is proposed for the vegetation isolines intersection point.

The right-shift phenomenon shows that as the FC increases in the homogeneous canopy (defined in this paper as a canopy with a clumping index equal to 1), the vegetation isolines and soil line's intersection points gradually move towards the positive direction of the red band axis. When the intercept of the vegetation isolines is smaller than that of the soil line, the final intersection point can reach the positive region of the red band axis. The right-shift phenomenon successfully solves the two primary debates in the current academic community and achieves a dialectical unification from a new perspective. Considering that the optimal soil adjustment factor is the negative value of the abscissa of the intersection of the vegetation isolines and the soil line, based on the right-shift phenomenon, the hypothesis that the optimal soil adjusted factor should be negative in high FC areas is put forward. Results show that the optimal soil adjustment factor is approximately equal to -0.148 when the average LAI equal 5.35; the optimal soil adjustment factor is approximately equal to -0.183 when the average LAI equal 6.72. The hypothesis can mitigate the saturation effect of SAVI and improve the accuracy of LAI estimation in the high FC area.

In this thesis, vegetation reflectance simulation and bio-optical property inversion based on a three-dimensional radiative transfer model is the research object in terms of both simulation and inversion. The first half of this thesis focuses on the modelling and reflectance simulation of a 3D vegetation scene with a coupled growth model. The second half focuses on the inversion of vegetation bio-optical parameters using the 3D radiative transfer model. Potential applications include providing high-quality analytical validation data for sensor design and adequate data support for quantitative remote sensing inversion modelling, spatial and temporal scale conversion, and multi-source data assimilation to achieve the objective of "quantitative simulation, quantitative evaluation".

Résumé

En tant que composante de l'écosystème terrestre, la végétation peut maintenir l'équilibre écologique et améliorer l'environnement. La végétation est indissociable de l'homme, jouant ainsi un rôle essentiel dans sa survie et son développement. Les données de télédétection sont de plus en plus utilisées, en particulier via les méthodes d'inversion, pour la surveillance de la végétation en raison de l'amélioration de la précision des mesures et de la résolution spatiale/spectrale/temporelle, ainsi que des progrès des méthodes d'interprétation des données de télédétection. Les modèles de transfert radiatif unidimensionnels simulent souvent avec une mauvaise précision la réflectance de la végétation, si bien que leur utilisation dans des procédures d'inversion des mesures de télédétection en termes de paramètres bio-optiques conduit à des paramètres estimés qui peuvent être très imprécis. Les modèles de transfert radiatif tridimensionnels sont beaucoup plus précis, car ils peuvent prendre en compte de manière réaliste l'architecture des couverts végétaux (FC), si bien qu'ils peuvent modéliser avec précision les mécanismes radiatifs complexes qui surviennent au sein des FC. Cependant, l'emploi des modèles 3D, et en particulier leur inversion, se heurte à trois problèmes majeurs. (1) Les modèles de simulation actuels ne permettent pas de simuler de manière continue l'évolution de la végétation du fait de la méconnaissance de la variation spatiale et temporelle des paramètres clés du sol. Cette difficulté est une contrainte importante pour la mise en oeuvre de l'inversion de séries temporelles d'images de télédétection. [START_REF] Li | Spatio-temporal Characteristics of Area Coverage and Observation Geometry of the MISR Land-surface BRF Product: A Case Study of the Central Part of Northeast Asia[END_REF] Les méthodes d'inversion spectrale foliaire existantes sont principalement applicables aux couverts homogènes tels que chaque unité spatiale d'inversion comprend plusieurs éléments de paysages (e.g., arbres). Cependant, l'observation de scènes qui comportent de nombreuses composantes dont la dimension est du même ordre de grandeur ou inférieure à la résolution spatiale du capteur, conduit à l'apparition de nombreux pixels mixtes. C'est en particulier le cas des villes, où la précision de l'inversion des propriétés optiques foliaires des arbres tend à être imprécise, voire très imprécise. (3) Les indices de végétation existants sont facilement saturés dans les zones à forte couverture végétale, ce qui limite fortement leur capacité d'inversion. Pour résoudre les problèmes mentionnés ci-dessus concernant l'inversion des images de télédétection, nous avons couplé le modèle DART (Discrete Anisotropic Radiative Transfer) de transfert radiatif 3D avec un modèle de croissance. Cette approche permet de simuler avec une meilleure précision l'évolution de la réflectance des couverts en fonction de leur développement. L'inversion de la végétation en milieu urbain est particulièrement complexe du fait de sa distribution éparse et de l'impact de l'architecture urbaine.

Nous avons résolu ce problème via le développement de la méthode appelée "Etalonnage de DART".

Nous avons aussi analysé les comportements des isolignes de végétation et proposé un décalage vers la droite du point d'intersection des isolignes calculé à partir de simulations DART pour atténuer l'effet de saturation de l'indice de végétation ajusté au sol (SAVI) dans les zones à fort FC. Afin de simuler l'évolution temporelle du transfert radiatif 3D dans les couverts végétaux, un modèle de modélisation 3D statique du maïs construit à partir du système L étendu (ELSYS) a été couplé à une équation de croissance dynamique du maïs utilisant les degrés-jours comme facteur de croissance et des règles de développement décrivant la structure de la canopée tout au long de la saison de croissance de la plante, de l'émergence des semis au stade de la floraison. La réflectance de couvert de maïs a ainsi été simulée avec plusieurs modèles de transfert radiatif 3D pour une validation croisée, et un bon accord entre eux est obtenu. Les réflectances de la canopée simulées avec plusieurs modèles de transfert radiatif 3D sont comparées à un modèle de transfert radiatif 1D dans les mêmes conditions d'indice de surface foliaire (LAI) que les valeurs réelles. Les résultats montrent que les hypothèses d'homogénéité du modèle de transfert radiatif 1D, utilisé avec le même indice de surface foliaire, entraînent une surestimation d'environ 1,5 fois de la FC, ce qui conduit à nettement surestimer la réflectance dans le proche infrarouge (NIR), pour la direction d'observation du nadir. Cette surestimation peut être expliquée par l'hypothèse d'homogénéité, car celle-ci ne considère pas le fait que la section efficace de la végétation est en général minimale si la direction d'observation est le nadir, alors que les diffusions multiples sont surtout dues à la végétation.

Considérant que la direction du nadir est la direction d'observation la plus commune et que le NIR est une bande essentielle pour la surveillance de la végétation, nous pensons que la surestimation due à l'hypothèse d'homogénéité n'est pas négligeable et entraîne une sous-estimation de la quantité de végétation (LAI, etc.) dans les inversions de mesures de télédétection.

Afin d'inverser les propriétés optiques des feuilles dans les zones urbaines, la calibration DART est introduite. Cette calibration décompose tout d'abord la réflectance d'ordre 1 de la canopée à l'aide d'un modèle de mélange spectral linéaire et corrige de manière itérative les propriétés optiques foliaires de manière à ce que la réflectance de la canopée simulée soit proche de la somme de ses composantes. Les résultats de l'inversion sont analysés via une évaluation de leur précision et une analyse de sensibilité. Les erreurs relatives moyennes de l'expérience de simulation idéale sans bruit sont de 0,013, 0,005, 0,027, 0,297 et 0,250 pour le sol, le toit, l'eau, les arbres et les arbustes dans toutes les bandes. Pour des conditions expérimentales réalistes, c'est-à-dire en présence de différents bruits (décalage des pixels, précision géométrique de la scène 3D et fonction de transfert de modulation avec une certaine déviation de la valeur réelle), des erreurs significatives dans l'inversion sont observées : pour le sol, le toit, l'eau, les arbres et les arbustes, les erreurs relatives moyennes sont respectivement de 0,233, 0,507, 3,088, 0,834 et 1,256. Les paramètres qui ont la plus grande influence sur la précision de l'inversion des propriétés optiques des matériaux urbains sont la SZA, la résolution spatiale des images satellites, le décalage des pixels, l'imprécision de la modélisation de la scène urbaine en 3D et la fonction de transfert de modulation, par ordre décroissant d'importance. Afin d'atténuer l'effet de saturation de SAVI dans les zones à FC élevé, les isolignes de végétation sont analysées, et un phénomène de décalage vers la droite est proposé pour le point d'intersection des isolignes de végétation. Le phénomène de décalage vers la droite montre que lorsque le FC augmente dans la canopée homogène (définie dans cet article comme une canopée avec un indice d'agglutination égal à 1), les points d'intersection des isolignes de végétation et de la ligne de sol se déplacent progressivement vers la direction positive de l'axe de la bande rouge. Quand l'interception des isolignes de végétation est plus petite que celle de la ligne de sol, le point d'intersection final peut atteindre la région positive de l'axe de la bande rouge. Ce décalage vers la droite résout avec succès deux controverses majeures dans ce domaine. Il réalise une unification Résumé IX dialectique à partir d'une nouvelle perspective. Considérant que le facteur optimal d'ajustement du sol est la valeur négative de l'abscisse de l'intersection des isolignes de végétation et de la ligne de sol, sur la base du phénomène de décalage vers la droite, l'hypothèse selon laquelle le facteur optimal d'ajustement du sol devrait être négatif dans les zones à FC élevé est avancée. Les résultats montrent que le facteur optimal d'ajustement du sol est approximativement égal à -0.148 si le LAI moyen est égal à 5.35 ; le facteur optimal d'ajustement du sol est approximativement égal à -0.183 si le LAI moyen est égal à 6.72. L'hypothèse peut atténuer l'effet de saturation de SAVI et améliorer la précision de l'estimation du LAI dans les zones à FC élevé.

Dans cette thèse, la simulation de la réflectance de la végétation et l'inversion des propriétés bio-optiques basées sur un modèle de transfert radiatif tridimensionnel constituent l'objet de recherche tant en termes de simulation que d'inversion. La première moitié de cette thèse se concentre sur la modélisation et la simulation de la réflectance d'une scène de végétation en 3D avec un modèle de croissance couplé. La seconde moitié se concentre sur l'inversion des paramètres biooptiques de la végétation en utilisant le modèle de transfert radiatif 3D DART. Les applications potentielles comprennent la fourniture de données de validation analytique de haute qualité pour la conception de capteurs et un support de données adéquat pour la modélisation d'inversion de télédétection quantitative, la conversion d'échelle spatiale et temporelle, et l'assimilation de données multi-sources pour atteindre l'objectif de "simulation quantitative, évaluation quantitative". Therefore, many studies have been carried out to monitor vegetation change. However, human activities (e.g., deforestation) put enormous pressure on the environment, causing significant land cover changes, affecting the regional and global climate system. These changes and the resulting impacts need to be harnessed on a regional or global scale. Remote sensing is an essential tool to monitor global vegetation changes as they allow for rapid, extensive scale observations in outer space, offering the possibility to study human activities and climate changes on a large scale.

List of Figures

XV

List of Figures

List of Tables

List of Acronyms

Motivation

Remote sensing data are increasingly used because of improvements in the accuracy of sensor radiometry and remote sensing imagery in terms of spatial/spectral/temporal resolution and advances in remote sensing data interpretation methods [START_REF] Gastellu-Etchegorry | Why To Model Remote Sensing Measurements In 3d? Recent Advances In Dart: Atmosphere, Topography, Large Landscape, Chlorophyll Fluorescence And Satellite Image Inversion[END_REF][START_REF] Li | Spatio-temporal Characteristics of Area Coverage and Observation Geometry of the MISR Land-surface BRF Product: A Case Study of the Central Part of Northeast Asia[END_REF]. These methods typically use remote sensing models to simulate the surface's bi-directional reflectance factor (BRF). In addition, realistic 3D structural models (terrain, trees, buildings, etc.) are increasingly used due to the growing demand for satellite inversion information and radiative budget accuracy. Remote sensing data are increasingly used to obtain vegetation optical and biophysical properties due to advantages such as high-frequency re-entry cycles and wide observational ranges [START_REF] Zhen | Potentials and Limits of Vegetation Indices With BRDF Signatures for Soil-Noise Resistance and Estimation of Leaf Area Index[END_REF][START_REF] Zhen | Using the Negative Soil Adjustment Factor of Soil Adjusted Vegetation Index (SAVI) to Resist Saturation Effects and Estimate Leaf Area Index (LAI) in Dense Vegetation Areas[END_REF][START_REF] Zhu | A regional maize yield hierarchical linear model combining landsat 8 vegetative indices and meteorological data: Case study in jilin province[END_REF]. Inferring the state of the land surface through the interpretation of remote sensing electromagnetic signals (reflections or emissions) is the "inversion problem" in remote sensing. In general, inversion aims to determine the model's input parameters to match the available measurements. In the noise-free ideal case, the input and output are mathematically related. However, the scene structure is often difficult to describe in terms of simple statistical parameters due to its complexity. It explains why scenes are often simplified in remote sensing models. For example, radiative transfer (RT) models based on the homogeneity assumption model the vegetation canopy as a horizontal layer with randomly distributed leaves. Based on this scene simplification, many one-dimensional models relying on atmospheric RT theory have been developed [6]. For example, canopy RT models are usually developed by formulating relationships between the BRF and vegetation parameters such as leaf area index (LAI), leaf angle distribution (LAD), and leaf optical properties (OPs). So far, 1D RT models, for example, Scattering by Arbitrary Inclined Leaves (SAIL) model and Soil Canopy Observation, Photochemistry and Energy fluxes (SCOPE) model, have been widely used for parametric inversion of remote sensing. However, due to their high abstraction of the earth's surface, the simulation accuracy of 1D models is usually too low. For example, they usually fail to consider the significant gaps between tree canopies and trees of different heights. Therefore, threedimensional (3D) RT models are needed to accurately invert remote sensing data because they can present complex heterogeneous remote sensing scenes.

3D RT models use schematic geometric objects (e.g., ellipsoids), triangular meshes and/or voxels with turbid media to describe 3D remote sensing scenes [7] and often employ ray tracing or radiosity methods in order to solve the RT equations [START_REF] Disney | Monte Carlo ray tracing in optical canopy reflectance modelling[END_REF][START_REF] Gastellu-Etchegorry | DART: a 3D model for simulating satellite images and studying surface radiation budget[END_REF]. They can simulate remote sensing data under arbitrary conditions, which is essential for relating remote sensing signals to realistic scene structures. Three-dimensional radiative transfer models allow the maximum use of information from multi-source remote sensing data and alleviate the 'intractable' problems remote sensing inversions faced for many years [START_REF] Liang | Recent progress in land surface quantitative remote sensing[END_REF]. Today, with the increasing availability of 3D spatial data because of advances in computer technology, the potentials of 3D RT models are growing, and their developments have two main implications [START_REF] Yang | Review of Optical Remote Sensing Imaging Simulation of Farmland Radiation Transfer Process[END_REF].

( The canopy reflectance simulation model describes the propagation and interaction (scattering, transmission, absorption, etc.) of energy between the various components of the canopy (leaves, soil, branches, etc.). Canopy reflectance simulation models allow the scaling up of leaf OPs to canopy level using canopy structural parameters. Besides, they are the base for the leaf OPs inversion. There are three main types of canopy reflectance simulation models commonly used today.

(1) Empirical and semi-empirical models Empirical models do not attempt to explain the physical mechanisms between biophysical parameters and the BRF; they describe mathematically the patterns observed in the BRF data set.

Semi-empirical models rely on the physical principles of simplified geometric optical (GO) models and radiative transfer theory. For example, the kernel-driven model [START_REF] Roujean | A bidirectional reflectance model of the Earth's surface for the correction of remote sensing data[END_REF][START_REF] Wanner | On the derivation of kernels for kernel-driven models of bidirectional reflectance[END_REF][START_REF] Liu | The angular and spectral kernel model for BRDF and albedo retrieval[END_REF] describes the BRF as the sum of isotropic and anisotropic scattering functions (i.e., kernel functions) characterizing voxel and surface scattering. Widely used semi-empirical models benefit from the use of numerical solutions with few parameters. For example, the Rahman, Pinty Verstraete (RPV) model [START_REF] Rahman | Coupled surface-atmosphere reflectance (CSAR) model: 2. Semiempirical surface model usable with NOAA advanced very high resolution radiometer data[END_REF] and several subsequent versions of accelerated inversion, such as the Modified RPV model (MRPV) [START_REF] Martonchik | Determination of aerosol optical depth and land surface directional reflectances using multiangle imagery[END_REF] and the Exponential Modified Rahman, Pinty Verstraete (EMRPV) model [START_REF] Engelsen | Parametric bidirectional reflectance factor models: Evaluation, improvements and applications[END_REF] Exhibiting Reflectance and Transmittance Yields (LIBERTY) [START_REF] Dawson | LIBERTY-Modeling the effects of leaf biochemical concentration on reflectance spectra[END_REF], which simulates the OPs of leaves.

( Hapke model [START_REF] Hapke | Bidirectional reflectance spectroscopy: 1. Theory[END_REF]), leaf models (e.g., the PROSPECT model [START_REF] Feret | PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments[END_REF]), and canopy models (e.g., the SAIL model [START_REF] Verhoef | Light scattering by leaf layers with application to canopy reflectance modeling: The SAIL model[END_REF]). In addition, some models such as DART use the output of soil and leaf radiative transfer models to simulate the radiative transfer of the entire "Earth-atmosphere system".

One of the main difficulties in radiative transfer modelling is the accurate calculation of multiple scattering under the precondition of energy conservation. Different methods of calculation are available, and four main ones are described.

(i) Discrete ordinate method

The discrete ordinate method assumes that the radiation propagates along N discrete directions and is therefore based on N radiative transfer equations. For example, the SAIL model [START_REF] Verhoef | Light scattering by leaf layers with application to canopy reflectance modeling: The SAIL model[END_REF] relies on a system of 4 differential equations with 4 fluxes / 4 directions applied to horizontally homogeneous scenes: one solar flux, two upward and downward isotropic fluxes, and one flux along the direction of the sensor observation. For better consideration of anisotropic radiation, the number of fluxes can be much larger (e.g., 100 and more) [START_REF] Yin | A new approach of direction discretization and oversampling for 3D anisotropic radiative transfer modeling[END_REF], as in the DART model [START_REF] Gastellu-Etchegorry | Modeling radiative transfer in heterogeneous 3-D vegetation canopies[END_REF][START_REF] Gastellu-Etchegorry | Recent Improvements in the Dart Model for Atmosphere, Topography, Large Landscape, Chlorophyll Fluorescence, Satellite Image Inversion[END_REF]. A classical approach [START_REF] Kimes | Radiative transfer model for heterogeneous 3-D scenes[END_REF][START_REF] Myneni | A three-dimensional radiative transfer method for optical remote sensing of vegetated land surfaces[END_REF] is to discretize the spatial variables into a set of spatial nodes to consider the threedimensional heterogeneity of the scene (i.e., voxelization).

(ii) Successive orders of scattering (SOS)

Successive orders of scattering is the oldest and conceptually most straightforward methods for solving multiple scattering problems. It uses an iterative approach to calculate the continuous order of scattering: the total irradiance vector is the sum of the contributions after several times photon scattering. For example, the Successive Orders of Scattering Vector Radiative Transfer (SOSVRT) model [START_REF] Duan | A polarized radiative transfer model based on successive order of scattering[END_REF] simulates polarised radiative transfer in a vertically inhomogeneous parallel planar medium.

(iii) Monte Carlo

Monte Carlo simulates a series of scattering events of photons between a light source and a receiver. It has the advantage that only a single scattering characteristic needs to be accurately modelled [START_REF] Disney | Monte Carlo ray tracing in optical canopy reflectance modelling[END_REF]. However, calculation time is its major limitation. 3-D Forest LIGHT model (FLIGHT) [START_REF] North | Three-dimensional forest light interaction model using a Monte Carlo method[END_REF], Drat [START_REF] Lewis | Three-dimensional plant modelling for remote sensing simulation studies using the Botanical Plant Modelling System[END_REF] and Raytran [START_REF] Govaerts | Raytran: A Monte Carlo ray-tracing model to compute light scattering in three-dimensional heterogeneous media[END_REF]

(iv) Radiosity

Radiosity is based on the radiative balance equation for a finite number of N discrete scattering facets [START_REF] Borel | The radiosity method in optical remote sensing of structured 3-D surfaces[END_REF]. Unlike the radiative transfer, radiosity is based on the radiative balance of the voxels in three dimensions. It requires solving the inverse matrix of an N × N matrix. It needs to calculate the view factor between all N facets to do this. This step can be very time-consuming and complex if N is large, especially for scene elements such as trees. Therefore, one objective of our study is to accurate the computation of the view factor of the radiosity based model.

In addition to radiative transfer tracking, classification can also be based on simplifying remote sensing scenes: homogeneous and realistic landscapes. Stacking homogeneous and horizontal layers of turbid media creates homogenous landscapes (i.e., random distribution of infinitely small facets).

With this assumption, "basic" models may replicate trends as crop BRF changes in response to changes in LAI. Landscape architecture is frequently just given a sliver of a picture. The realistic landscape is widely used in RS and radiative budget (RB) models. There are two main approaches:

discretizing the spatial variable into voxels filled with turbid media and/or simulating each individual scene element as a juxtaposition of facets. It should be able to prove a more accurate simulation result on the condition that an accurate landscape is provided [START_REF] Tao | Effects of plant and scene modeling on canopy NDVI simulation: a case study on Phragmites australis and Spartina alterniflora[END_REF]. Therefore, one objective of our study is to simulate the 3D maize scene accurately.

Inversion of leaf optical properties study

The OPs of leaves reflect their combination of biochemical, morphological and physiological properties and play an essential role in many ecological and earth system processes. There are three principal methods for obtaining OPs of leaves: instrumental measurements, simulations of leaf biochemical parameters, and canopy spectral inversion.

(

1) Optical instrument measurement

Remote sensing simulations or inversions of leaf OPs require spectral measurements to verify accuracy. There are currently three main methods for measuring the OPs of leaves [START_REF] Jacquemoud | Measurement of Leaf Optical Properties[END_REF]: (1) The instrument's optical fibre is connected to a port on an external integrating sphere, such as the LI-1800-12 or ASD RTS-3ZC, with an internal halogen light source for measuring the diffuse reflectance or transmittance. (2) The instrument's optical fibre is connected to unique accessories (plant probes and leaf clips) with an internal halogen light source. Using a double-sided rotating head containing a recessed background plate, reflectance (black panel, reflectance < 5%) and transmittance (Spectralon white panel, reflectance > 99% in the visible to the near-infrared band) measurements can be carried out simultaneously. This measurement method has the advantage of reducing external disturbances such as angle of incidence or unstable light sources. (3) The instrument's fibre optic is connected to a gun (Pistol Grip), and the leaf is observed vertically, for example, in the laboratory with light as the light source or the field with the sun at any angle of incidence as the light source. Then the bi-directional reflectance factor is measured, assuming that sensors and sunlight are assumed to be mono-directional.

Calibrating these data can be a daunting task, and most of the time, the calibration is not performed correctly. Few studies compare the reflectance of the same leaf recorded using different systems or spectrometers, each with a unique spectral resolution and sampling intervals. Castro-Esau et al. [START_REF] Castro-Esau | Comparison of spectral indices obtained using multiple spectroradiometers[END_REF] perform an instrument comparison and standardization of sampling procedures.

Depending on the instrument, they observe slight to severe differences in shape and amplitude between spectra of the same leaf. For some applications, such differences may not be significant at a practical level, as the calculated spectral indices are similar. For other applications, such as the use of physical models, such differences may be critical. Therefore, inverting PROSPECT model parameters using data measured by ASD probes or leaf clamps may lead to biased estimates of leaf biochemical parameters.

(2) Simulation of leaf radiative transfer

Variations of leaf reflectance and transmittance are modelled by analyzing the interaction processes of electromagnetic waves with the biochemical components within the leaf and the leaf structure. Light propagation in plant leaves is mainly governed by absorption and scattering interactions. In parallel with advances in experimental measurements of the OPs of leaves, deterministic methods using different representations of the interaction of light with plant leaves have been developed. These models differ depending on the underlying physics and the complexity of the leaf. The simplest models treat the leaf as a single scattering and absorbing layer, while in the most complex models, all cells' shape, size, location, and biochemical content are described in detail.

Regardless of the approach used, these models have improved understanding of the interaction of light with plant leaves. Baranoski and Rokne [START_REF] Baranoski | Light interaction with plants[END_REF], Ustin et al. [START_REF] Ustin | Remote sensing of environment: State of the science and new directions[END_REF], and Jacquemoud et al. [START_REF] Jacquemoud | Compréhension et modélisation de la couleur des feuilles[END_REF] have extensively reviewed computer-based leaf models that have improved the understanding of the interaction between light and plant leaves from the late 1960s to the present. These models are divided into different categories and arranged in increasing complexity [START_REF] Hecht | The interpretation of diffuse reflectance spectra[END_REF].

(i) Flat model

The flat plate model, first proposed by Allen et al. [START_REF] Allen | Interaction of isotropic light with a compact plant leaf[END_REF] in 1969, treats a dense leaf as a homogeneous flat plate with a rough surface. Light incident on such a leaf surface is reflected and transmitted several times. It is partially reflected and partially transmitted at the first interface, and the transmitted part is then reflected back and forth between the two interfaces. The total reflectance of the plate can be obtained by summing the amplitudes of successive reflections and refractions.

However, this model does not apply to non-compact leaves. After Allen et al. [START_REF] Allen | Mean effective optical constants of cotton leaves[END_REF], Breece and Holmes [START_REF] Breece | Bidirectional scattering characteristics of healthy green soybean and corn leaves in vivo[END_REF] extended the flat plate model to non-compact leaves by introducing the generalized flat plate model, subdividing the leaf into N homogeneous compact plates separated by N-1 air spaces. In the 20 th century, N values were extended from integers to the real number domain. The PROSPECT model, now widely used in the remote sensing community, was developed from the flat model [START_REF] Jacquemoud | PROSPECT: A model of leaf optical properties spectra[END_REF]. It is one of the first radiative transfer codes to accurately model the hemispheric reflection and transmission of various plant leaves and conditions (monocotyledons, dicotyledons of healthy or senescent leaves) across the solar spectrum from 400 nm to 2500 nm. The input parameters for the early PROSPECT models are: structural parameter N, green pigment and water content to simulate fresh leaves. If dry leaf spectra need to be simulated, protein and cellulose + lignin content must be added. However, the plate model cannot be applied to needle leaves because they cannot be treated as discrete parallel plates.

(ii) Compact spherical particle model

Dawson [START_REF] Dawson | Modelling the spectral response of coniferous leaf structures for the estimation of biochemical concentrations[END_REF][START_REF] Dawson | A new technique for interpolating the reflectance red edge position[END_REF] modified the equation for the interaction of light with spherical particles and designed a degree-of-freedom model precisely to calculate the OPs of dried and fresh slash pine (Pinus elliottii) needles to model the OPs of needles. The model treats needle cell structures as spherical cells with diameters and air gaps according to the laws of geometrical optics and then combines leaf thickness, absorption coefficients of water, chlorophyll, cellulose, lignin and nitrogenous compounds in the leaf to correct for the OPs of the needles.

(iii) N flux model

The N-flux model is derived from the Kubelka-Munk (KM) theory, which considers leaves plates filled with absorption and scattering coefficients [START_REF] Kubelka | Ein Beitrag zur Optik der Farbanstriche (Contribution to the optic of paint)[END_REF]. In order to eliminate edge effects, lateral extension and boundary reflections below the plate are assumed to be absent, and the OPs are expressed as a function of the absorption coefficient, the scattering coefficient and the thickness of the leaf. The absorption coefficient and scattering coefficient are expressed in later versions as a function of the leaf thickness. Radiative transfer within the leaf is modelled using paired fluxes, such as two-and four-flux models. N-flux model has the advantage that the pigment content can be inverted without damage. The complete leaf biochemistry is described by Conel et al. [START_REF] Conel | Application of a Two-Stream Radiative Transfer Model For Leaf Lignin and Cellulose Concentrations from Spectra Reflectance Measurements (Part 1)[END_REF], who used a two-flux model to study the effects of water, protein, cellulose, lignin and starch on midinfrared radiation in leaves. However, they did not validate their model. Finally, note that the parameter N describing the internal structure of the leaf in the generalized flat plate model plays a role similar to that of the scattering coefficient in the KM model. One drawback of this approach is that it cannot consider the multiple scattering of leaf flesh objects of a size comparable to the wavelength of the incident radiation (cells, organelles, bubbles and others) and cannot characterize them.

(iv) Radiative transfer equation

The primary electromagnetic theory can describe the propagation of light through plant leaves.

The leaf can be considered a random medium with a spatially varying dielectric constant, and the electric field variation can be described by Maxwell's theory. Thus, the problem is simplified as the one-dimensional energy flow through the medium. In contrast to canopy radiative transfer, only a few models directly use the radiative transfer equations at the leaf scale. The mathematical complexity of the algorithm and the lack of information on the leaf's internal structure and biochemical distribution lead to substantial simplifications, making this approach less efficient than more powerful formulations [START_REF] Jacquemoud | Leaf optical properties[END_REF].

(v) Stochastic models

Tucker and Garratt [START_REF] Tucker | Leaf optical system modeled as a stochastic process[END_REF] present a primitive stochastic model in which the radiative transfer is modelled by a Markov chain, which is a stochastic process. They divide the leaf into two different tissues (palisade parenchyma and spongy mesophyll) and define states of photons (incident solar radiation, specularly reflected radiation, diffuse reflected radiation, diffuse transmitted radiation, absorbed radiation and scattered radiation, in each tissue). The random photon states can only take on discrete values, such as "absorbed in the palisade parenchyma" or "scattered in the spongy mesophyll". The initial vector sets the incident radiation, and the OPs of the leaf is obtained when iterative state shifts are made until they smooth out.

(vi) Light tracing model

The advantage of ray-tracing models is that they can describe the complexity of the internal structure of a leaf (e.g., individual cells and their unique arrangement within the tissue). By defining the OPs of the leaf material (cell walls, cytoplasm, pigments, stomata, etc.), the propagation of individual photons incident on the leaf surface can be simulated using the laws of reflection, refraction and absorption. Once a sufficient number of rays have been simulated, a statistically valid estimate of the radiative transfer in the leaf can be deduced. The technique has been applied to many variants. The earliest studies were carried out at the cellular level. Senn [START_REF] Senn | Die Gestalts-und Lageveränderung der Pflanzen-Chromatophoren: mit einer Beilage: Die Lichtbrechung der lebenden Pflanzenzelle[END_REF], Haberlandt [START_REF] Haberlandt | Optical sense-organs[END_REF] and more recently, Gabrys-Mizera [START_REF] Gabrys-Mizera | Model considerations of the light conditions in noncylindrical plant cells[END_REF] and Bone et al. [START_REF] Bone | Epidermal cells functioning as lenses in leaves of tropical rain-forest shade plants[END_REF] construct geometric models of the passage of light through cross-sections of plant cells of different shapes, particularly epidermal cells, whose shape may influence the path of the incident light. Allen et al. [START_REF] Allen | Willstätter-Stoll theory of leaf reflectance evaluated by ray tracing[END_REF] treat leaves as cell walls and cell gaps and model the spectrum of maple leaves using ray tracing. However, this method requires a significant computational load and has been commonly applied to validate simple models and understand the light transmission processes inside the leaf.

(

3) Inversion of canopy radiative transfer

Taking the inversion of the OPs of leaves as an example, as the leaf and canopy radiative transfer models continue to be studied in-depth, researchers have found that the radiation signals observed by sensors can be used to interpret the biophysical state of the vegetation. The canopy radiation signal captured by the sensor can be simulated by a priori knowledge of the vegetation (canopy structure, leaf OPs, etc.) using a suitable canopy radiative transfer model, a process known as forwarding simulation. In contrast, the bio-optics of the leaf or the canopy's geometry can be obtained by using a suitable canopy radiative transfer model based on the a priori knowledge of the radiation signal captured by the sensor. For example, Otterman [START_REF] Otterman | Inferring spectral reflectances of plant elements by simple inversion of bidirectional reflectance measurements[END_REF] assumed a canopy structure and inverted a protrusion model to determine the reflectance of the leaf in the nadir observation direction and several bands. Otterman [START_REF] Otterman | Inferring parameters for canopies nonuniform in azimuth by model inversion[END_REF] later attempted to invert the leaf orientation but found it difficult to isolate the leaf reflectance from the leaf orientation or LAI. Kuusk [START_REF] Kuusk | Determination of vegetation canopy parameters from optical measurements[END_REF] and Privette et al. [START_REF] Privette | Invertibility of a 1-D discrete ordinates canopy reflectance model[END_REF] invert the OPs of leaves and the spatial distribution of scatters in the canopy by inversion of a bidirectional reflectance analysis model.

A major problem limiting the accuracy of the inversion is multi-solvability. An accurate canopy radiative transfer model is difficult to invert the input parameters from the output parameters, and it is difficult or even impossible to find the inversion function corresponding to the radiative transfer equation. Inverting the model is to find the parameter that minimizes the difference between the measured and simulated data. The critical factor for the success of leaf OPs inversion is choosing a suitable canopy radiative transfer model and inversion method. Commonly used canopy radiative transfer models and their applicability are listed in section 1.2.1, and standard inversion methods are described in detail below, including iterative optimization algorithms [START_REF] Jacquemoud | PROSPECT: A model of leaf optical properties spectra[END_REF][START_REF] Jacquemoud | Comparison of four radiative transfer models to simulate plant canopies reflectance: Direct and inverse mode[END_REF], Look-Up Tables algorithms (LUTs) [START_REF] Atzberger | Comparative analysis of different retrieval methods for mapping grassland leaf area index using airborne imaging spectroscopy[END_REF][START_REF] Schlerf | Inversion of a forest reflectance model to estimate structural canopy variables from hyperspectral remote sensing data[END_REF], and machine learning algorithms such as neural networks tool (NNT)

and support vector machines [START_REF] Bacour | Neural network estimation of LAI, fAPAR, fCover and LAI× Cab, from top of canopy MERIS reflectance data: Principles and validation[END_REF][START_REF] Fang | Retrieving leaf area index with a neural network method: Simulation and validation[END_REF][START_REF] Durbha | Support vector machines regression for retrieval of leaf area index from multiangle imaging spectroradiometer[END_REF].

(i) Iterative optimization algorithm

An Many authors used iterative optimization to obtain biochemical parameters for various vegetation types [START_REF] Jacquemoud | PROSPECT: A model of leaf optical properties spectra[END_REF][START_REF] Jacquemoud | Comparison of four radiative transfer models to simulate plant canopies reflectance: Direct and inverse mode[END_REF][START_REF] Zhen | Quantitative Analysis of DART Calibration Accuracy for Retrieving Spectral Signatures Over Urban Area[END_REF].

(ii) Look-Up Tables method

The use of iterative optimization algorithms for bio-optical property inversion with remote sensing from large vegetation areas requires many cost functions and a high operational load, which is a problem that can be well avoided with LUTs. LUTs are an efficient and straightforward inversion method. LUTs can improve the problem of machine learning algorithms falling into local optima when the training sample is small. In addition, LUTs can reduce the outliers between the simulated spectra and the observed spectra of the sensors. Moreover, LUTs do not rely on ground sampling sites for modelling, so they can invert vegetation parameters based on radiative transfer models when the number of ground sampling sites is insufficient, effectively eliminating the instability of statistical methods. However, the LUT size can be proliferated if the step size of the input parameters is too small, leading to non-representative results. Therefore, a parameter sensitivity analysis should first be conducted to set reasonable parameters' range and step size before building a LUT to improve the efficiency and accuracy of the inversion. Darvishzadeh and Matkan [START_REF] Darvishzadeh | Inversion of a radiative transfer model for estimation of rice canopy chlorophyll content using a lookup-table approach[END_REF] obtain canopy chlorophyll content using a LUT approach by inverting the PROSPECT + SAIL (PROSAIL) vegetation radiative transfer model. Dorigo [START_REF] Dorigo | Improving the robustness of cotton status characterisation by radiative transfer model inversion of multi-angular CHRIS/PROBA data[END_REF] uses Compact High Resolution Imaging Spectrometer (CHRIS) / Project For On-Board Autonomy (PROBA) multi-angle remote sensing images to obtain the chlorophyll content of the canopy.

(iii) Machine learning methods

The goal of machine learning is to draw lessons from training samples. The main machine learning methods for the inversion of vegetation bio-OPs are currently artificial neural networks and support vector machines. The inversion based on artificial neural networks consists of two main steps: [START_REF] Gastellu-Etchegorry | Why To Model Remote Sensing Measurements In 3d? Recent Advances In Dart: Atmosphere, Topography, Large Landscape, Chlorophyll Fluorescence And Satellite Image Inversion[END_REF] training the neural network using many training samples; and (2) using the training network to invert the canopy bio-optical property. Neural networks have a solid non-linear mapping capability and are suitable for complex internal mechanisms problems. However, neural networks require a large number of training samples in the training process and are prone to overfitting. The support vector machine approach is mainly applied in remote sensing for classification and less frequently in vegetation bio-optical parameter inversion. Durbha [START_REF] Durbha | Support vector machines regression for retrieval of leaf area index from multiangle imaging spectroradiometer[END_REF] uses support vector machines and the PROSAIL radiative transfer model to invert the leaf area index from Multi-angle Imaging SpectroRadiometer (MISR) images with reasonable accuracy. Therefore, the applications of machine learning methods to invert vegetation biochemical parameters needs further research.

A short comparison between the look-up tables method and machine learning methods is that machine learning methods do not need to know the underlying data distribution explicitly. Hence, machine learning methods are developed without assuming a particular probability density distribution, which is why they work well with all kinds of data types [START_REF] Verrelst | Quantifying vegetation biophysical variables from imaging spectroscopy data: a review on retrieval methods[END_REF]. However, for the lookup table methods, the data distribution must be considered when creating the lookup table. Machine learning methods also offer the possibility of incorporating prior knowledge and the flexibility to include different data types in the analysis. The machine learning methods use all available training data. As a result, all data take effect. However, for the LUT method, not all data take effect. The main problem for machine learning is the collinearity complicating regression. However, noise removal is still a standard and much-needed pre-processing step for the lookup table method.

Although various inversion methods are used, at present, the research on the inversion of leaf spectral characteristics is still based on the assumption that the scene is a homogeneous canopy. This assumption may be reasonable for dense vegetation cover areas such as forests. However, this assumption may not be feasible for urban plants, especially in the current commonly used satellite images (i.e., 10 m for Sentinel-2). Therefore, one objective of our study is to retrieve leaf optical property (OP) in the urban areas based on DART radiative transfer model and iterative optimization algorithm.

Inversion of leaf area index study

LAI is a critical biophysical parameter in studies of ecophysiology, atmosphere-ecosystem interactions and global climate change [START_REF] Xiao | Long-Time-Series Global Land Surface Satellite Leaf Area Index Product Derived From MODIS and AVHRR Surface Reflectance[END_REF][START_REF] Yan | Scale Effect in Indirect Measurement of Leaf Area Index[END_REF][START_REF] Zheng | Assessing the Contribution of Woody Materials to Forest Angular Gap Fraction and Effective Leaf Area Index Using Terrestrial Laser Scanning Data[END_REF][START_REF] Kandasamy | Assessment of Three Methods for Near Real-Time Estimation of Leaf Area Index From AVHRR Data[END_REF][START_REF] Ma | Simultaneous Estimation of Leaf Area Index, Fraction of Absorbed Photosynthetically Active Radiation, and Surface Albedo From Multiple-Satellite Data[END_REF][START_REF] Hu | Using Airborne Laser Scanner and Path Length Distribution Model to Quantify Clumping Effect and Estimate Leaf Area Index[END_REF][START_REF] Yan | Generating Global Products of LAI and FPAR From SNPP-VIIRS Data: Theoretical Background and Implementation[END_REF][START_REF] He | Reduced Prediction Saturation and View Effects for Estimating the Leaf Area Index of Winter Wheat[END_REF][START_REF] Qu | Retrieval of a Temporal High-Resolution Leaf Area Index (LAI) by Combining MODIS LAI and ASTER Reflectance Data[END_REF][START_REF] Qu | Remote Sensing of Leaf Area Index from LiDAR Height Percentile Metrics and Comparison with MODIS Product in a Selectively Logged Tropical Forest Area in Eastern Amazonia[END_REF]. Since measuring LAI in the field is time consuming and laborious, remote sensing methods are now widely used to invert LAI over extensive areas [START_REF] Kang | How Universal Is the Relationship between Remotely Sensed Vegetation Indices and Crop Leaf Area Index? A Global Assessment[END_REF][START_REF] Zhang | Development of a New BRDF-Resistant Vegetation Index for Improving the Estimation of Leaf Area Index[END_REF][START_REF] Wang | Analysis of Differences in Phenology Extracted from the Enhanced Vegetation Index and the Leaf Area Index[END_REF][START_REF] Sun | Crop Leaf Area Index Retrieval Based on Inverted Difference Vegetation Index and NDVI[END_REF]. Inversion of LAI using remote sensing observations can be broadly classified into physical and empirical methods [START_REF] Myneni | Leaf area index and fraction of absorbed PAR products from terra and aqua MODIS sensors: Analysis, validation, and refinement[END_REF]. Physical radiative transfer models have the advantage of being based on radiative transfer model mechanisms and vegetation ecology theory. Therefore, they are not susceptible to the influence of vegetation type and non-vegetation factors [START_REF] Gascon | Retrieval of forest biophysical variables by inverting a 3-D radiative transfer model and using high and very high resolution imagery[END_REF]. Commonly used physical radiative transfer models are described in section 1.2.1, while common inversion methods are described in section 1.2.2. The procedure for physical model-based LAI inversion is similar to that for leaf spectral inversion and is not repeated here. However, the a priori parameters required for inversion using physical models are complicated. Moreover, in many cases, it is not possible to obtain all required inversion parameters precisely. For example, Goel and Thompson [START_REF] Goel | Inversion of vegetation canopy reflectance models for estimating agronomic variables. IV. Total inversion of the SAIL model[END_REF][START_REF] Goel | Inversion of vegetation canopy reflectance models for estimating agronomic variables. V. Estimation of leaf area index and average leaf angle using measured canopy reflectances[END_REF] and

Deering [START_REF] Deering | Shinnery oak bidirectional reflectance properties and canopy model inversion[END_REF] uses directional reflectance data to invert the Leaf Area Index (LAI) and lower accuracy average leaf inclination angle (ALA), respectively, from SAIL [START_REF] Verhoef | Light scattering by leaf layers with application to canopy reflectance modeling: The SAIL model[END_REF] and Three-dimensional Radiation Interaction Model (TRIM) [START_REF] Goel | Inversion of canopy reflectance models for estimation of biophysical parameters from reflectance data[END_REF] models. However, they must know the reflectance and transmittance of the leaves, the reflectance of the soil, and the skylight ratio in advance. In addition, physical models may fail due to uncertainties in input reflectance data such as clouds and aerosols when empirical models are often used as backup algorithms [START_REF] Myneni | Leaf area index and fraction of absorbed PAR products from terra and aqua MODIS sensors: Analysis, validation, and refinement[END_REF]. The most common method used in empirical models is vegetation indices (VIs) [START_REF] Kang | How Universal Is the Relationship between Remotely Sensed Vegetation Indices and Crop Leaf Area Index? A Global Assessment[END_REF][START_REF] Goel | Influences of canopy architecture on relationships between various vegetation indices and LAI and FPAR[END_REF].

In order to invert LAI over large areas, a variety of VIs have been devised over the last few decades [START_REF] Xue | Significant remote sensing vegetation indices: A review of developments and applications[END_REF]. The applications of these indices can be broadly classified into two categories [START_REF] Rondeaux | Optimization of soil-adjusted vegetation indices[END_REF].

(

1) Single angular index

According to their application, single-angle VIs can be broadly classified into three categories:

intrinsic VIs, atmosphere-adjusted VIs, and soil-adjusted VIs [START_REF] Rondeaux | Optimization of soil-adjusted vegetation indices[END_REF]. The intrinsic VIs are the most commonly used and contain only spectral reflectance with no other parameters. They are widely used because of their simplicity. However, they are susceptible to the atmosphere and soil OPs [START_REF] Rondeaux | Optimization of soil-adjusted vegetation indices[END_REF][START_REF] Huete | Spectral response of a plant canopy with different soil backgrounds[END_REF][START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF][START_REF] Qi | A modified soil adjusted vegetation index[END_REF][START_REF] Gilabert | A generalized soil-adjusted vegetation index[END_REF]. This shortcoming has led to the emergence of atmosphere-adjusted VIs to minimize the atmospheric impact. However, with the development of atmospheric radiative transfer models, atmosphere-corrected images are increasingly available as remote sensing images [START_REF] Berk | MODTRAN: A moderate resolution model for LOWTRAN[END_REF][START_REF] Myneni | Atmospheric effects and spectral vegetation indices[END_REF][START_REF] Qin | A fast, accurate algorithm to account for non-Lambertian surface effects on TOA radiance[END_REF][START_REF] Yin | Simulation of satellite, airborne and terrestrial LiDAR with DART (II): ALS and TLS multi-pulse acquisitions, photon counting, and solar noise[END_REF][START_REF] Gastellu-Etchegorry | DART: recent advances in remote sensing data modeling with atmosphere, polarization, and chlorophyll fluorescence[END_REF][START_REF] Yin | Atmospheric correction of ground-based thermal infrared camera through dart model[END_REF][START_REF] Jiang | Development of a two band enhanced vegetation index without a blue band[END_REF]. On the other hand, the soil effect in the vegetation mixed-image meta-spectrum can hardly be eliminated [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF][START_REF] Qi | A modified soil adjusted vegetation index[END_REF][START_REF] Baret | Potentials and limits of vegetation indices for LAI and APAR assessment[END_REF], which explains the development of soil-adjusted VIs.

(i) Inherent vegetation index

One of the first VIs was proposed by Jordan and Carl [START_REF] Jordan | Derivation of leaf-area index from quality of light on the forest floor[END_REF] in 1969 and named the Ratio Vegetation Index (RVI), which is based on the OP of vegetation. The leaves absorb more red band light than the Near-InfraRed (NIR) band, so plants have a lower reflectance in the red band.

nir red RVI R R = ••••••••••••••••••••••••••••••••••••• (1.1)
where Rnir and Rred are the atmospheres calibrated surface reflectance in the red and NIR bands, respectively.

The reflectance of vegetation shows a high correlation with LAI, Leaf Dry Biomass Matter (LDBM) and chlorophyll content of leaves [START_REF] Quan | Eco-environment variable estimation from remote sensed data and eco-environment assessment: models and system[END_REF]. RVI is widely used to estimate and monitor green biomass, especially in high-density Foliage Coverage (FC), as this index is susceptible to vegetation and correlates well with plant biomass. However, RVI is susceptible to atmosphere influences and becomes sensitive if FC is sparse (less than 50% coverage), and its correlation with biomass weakens. Although very sensitive to soil spectral variation, the Differential Vegetation Index (DVI) can also monitor vegetation ecology.

The most widely used Vegetation Index (VI), Normalised Difference Vegetation Index (NDVI) is proposed by Rouse et al. [START_REF] Rouse | Monitoring vegetation systems in the Great Plains with ERTS[END_REF] as

nir red nir red NDVI R R R R - = + ••••••••••••••••••••••••••••••••••• (1.2)
As the NDVI values are calculated by normalization, the index ranges from 0 to 1. The NDVI is sensitive to green vegetation even in areas of low FC. The index is often used in regional and global vegetation assessment studies and is related to canopy structures such as LAI and canopy photosynthesis [START_REF] Gamon | Relationships between NDVI, canopy structure, and photosynthesis in three Californian vegetation types[END_REF][START_REF] Grace | Can we measure terrestrial photosynthesis from space directly, using spectral reflectance and fluorescence?[END_REF]. However, NDVI is sensitive to soil brightness, soil colour, atmosphere, clouds and canopy shadows and requires remote sensing calibration. It has been shown that when soil background brightness increases, NDVI also increases systematically.

(ii) Soil adjusted vegetation index

In order to mitigate the interference of soil noise, many scholars have worked on soil-adjusted VIs. Richardson and Wiegand [START_REF] Richardson | Distinguishing vegetation from soil background information[END_REF] propose a method to distinguish vegetation from soil background spectra by analyzing soil lines: soil lines describe the linear relationship between soil reflectance in the red-NIR two-dimensional plane. The spectral response of the soil (soil line) is oblique in the red-NIR two-dimensional plane, as the soil shows a higher spectral response in the red-NIR plane. Therefore, it is considered an integrated description of a large amount of soil spectral information in a number of environments [START_REF] Baret | The soil line concept in remote sensing[END_REF]. Soil-adjusted VIs, such as the perpendicular vegetation index (PVI), are mainly designed based on the soil line. For example, PVI is defined as the distance between the reflectance of the vegetation pixels and the soil line. In this way, it effectively filters the influence of the soil background. At the same time, PVI is also less sensitive to the atmosphere, and it is mainly used for the inversion of vegetation parameters (LAI, grass production, chlorophyll content), vegetation identification and classification [START_REF] Kaufman | Atmospherically resistant vegetation index (ARVI) for EOS-MODIS[END_REF].

Due to some shortcomings of NDVI and PVI in describing the OP of vegetation and soil background, Huete [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF] develops a Soil Adjusted Vegetation Index (SAVI) to reduce the sensitivity of NDVI to soil background noise by adding a soil adjustment factor, XSAVI, and assuming that vegetation isolines intersect on the third quadrant angle bisector.

( )

nir red SAVI nir red SAVI (1 ) SAVI R R X R R X - + = + + •••••••••••••••••••••••••••• (1.3)
where XSAVI is the soil adjustment factor of SAVI, which is suggested to be equal to 0.5 in its original paper [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF]. In contrast to NDVI, SAVI uses a soil adjustment factor XSAVI to account for first-order soil background variability. When XSAVI = 0, SAVI is equivalent to NDVI. The advantage of SAVI is that the appropriate value of XSAVI can be chosen autonomously with prior knowledge of the vegetation density. XSAVI is the soil adjustment factor; its recommended value for general FC conditions is (XSAVI = 0.5) [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF], and the optimum soil adjustment factor should be close or equal to 0 under dense FC and close or equal to 1 under sparse FC.

Compared to RVI, SAVI is much less sensitive to background noise due to soil colour or soil surface moisture content. Three new versions of SAVI (SAVI2, SAVI3 and SAVI4) have been developed based on the spectral theory of wet and dry soils [START_REF] Major | A ratio vegetation index adjusted for soil brightness[END_REF]. SAVI2, SAVI3 and SAVI4 reduce background soil noise by considering the effects of changes in the angle of solar incidence and changes in the soil's physical structure.

Based on the changing rule of XSAVI, the Modified Soil Adjusted Vegetation Index (MSAVI) [START_REF] Qi | A modified soil adjusted vegetation index[END_REF] replaces the constant XSAVI in the SAVI of Eq.(1.3) with an inductive soil adjustment function.

In this way, MSAVI reduces the disturbance of SAVI by soil noise.

According to Eq.(1.3), we assume any seed value X 0 SAVI(-∞,+∞), which can minimize the soil noise ( )

SAVI SAVI 0 nir red 0 nir red (1 ) SAVI R R X R R X - + = + + ••••••••••••••••••••••••••• (1.4)
Based on the rule of optimal XSAVI, we can get an inductive soil adjustment function to minimize the soil noise further

1 1 MSAVI u u L - = - •••••••••••••••••••••••••••••••••• (1.5)
where u is the iteration number.

These iterations continue U times until soil noise cannot reduce further, namely MSAVIU= MSAVIU-1, then we have ( )

nir red nir red MSAVI 2 MSAVI 1 MSAVI U U U R R R R - = - + + - •••••••••••••••• (1.6)
One of the two solutions for Eq.(1.6) within the range of 0 and 1 is

( ) ( ) 2 nir ni ni re r r d 2 1 2 1 8 MSAVI 2 R R R R + - + - - = •••••••••••••••••• (1.7)
Based on MSAVI, Richardson and Wiegand [START_REF] Richardson | Distinguishing vegetation from soil background information[END_REF] propose a modified secondary soil-adjusted vegetation index (MSAVI2). MSAVI2 does not rely on the soil line principle, and its algorithm is simple; it is mainly used for plant growth analysis, desertification studies, grassland yield inversion, LAI assessment, soil organic matter analysis, drought monitoring, soil erosion analysis, etc. Some VIs are designed based on the concept of the soil line nir,soil red,soil

R a R b = ⋅ + •••••••••••••••••••••••••••••••• (1.8)
where Rred,soil and Rnir,soil represent the reflectance of soil pixels in the red and NIR band, and a and b represent the slope and intercept of the soil line.

The Perpendicular Vegetation Index (PVI) is the distance between vegetation reflectance point and soil line in the red-NIR plane.

nir red 2 PVI 1 R a R b a -⋅ - = + ••••••••••••••••••••••••••••••• (1.9)
The Transformed Soil Adjusted Vegetation Index (TSAVI) [START_REF] Baret | Potentials and limits of vegetation indices for LAI and APAR assessment[END_REF] considers both soil adjustment factors and soil line parameters and is free from SAVI restriction that vegetation isolines must intersect on the angle parallel.

( ) ( )

nir red 2 nir red TSAVI TSAVI 1 a R a R b a R R a b X a ⋅ -⋅ - = ⋅ + -⋅ + ⋅ + ••••••••••••••••••• (1.10)
where a and b are the slope and intercept of the soil line, set to 1.2 and 0.04, respectively; they are considered the global soil line parameters [START_REF] Baret | Potentials and limits of vegetation indices for LAI and APAR assessment[END_REF]. XTSAVI is the soil adjustment factor for TSAVI and is suggested to take a value of 0.08 in its original paper [START_REF] Baret | Potentials and limits of vegetation indices for LAI and APAR assessment[END_REF].

Simulated data from the SAIL model suggest that TSAVI appears to be the most reliable VI when the leaf inclination angle is known. However, it is not widely used due to the many parameters required (soil line parameters and soil adjustment factors). Rondeaux et al. [START_REF] Rondeaux | Optimization of soil-adjusted vegetation indices[END_REF] study the sensitivity of five VIs, including NDVI, SAVI, MSAVI, TSAVI and Global Environmental Monitoring Index (GEMI), to the soil background. They simulate the performance of the VIs under different soil textures, moisture and roughness by using the SAIL model. They determine an optimal soil adjustment factor to reduce the effect of soil background and then propose the Optimised Soil Adjusted Vegetation Index (OSAVI).

nir red

nir red OSAVI OSAVI R R R R X - = + + ••••••••••••••••••••••••• (1.11)
where XOSAVI is the soil adjustment factor and equals 0.16, it is easy to observe that OSAVI is a special case of TSAVI with soil line parameters a equal 1 and b equal 0.

OSAVI is independent of soil line and can effectively remove the effect of soil background.

However, OSAVI is not widely used, and it is mainly used to calculate above-ground biomass, leaf nitrogen content and chlorophyll content [START_REF] Wenlong | Vegetation index controlling the influence of soil reflection[END_REF].

(iii) Atmosphere-adjusted vegetation index

Since the atmosphere effect has a much more significant effect on the red band than the NIR band, Kaufman and Tanre [START_REF] Kaufman | Atmospherically resistant vegetation index (ARVI) for EOS-MODIS[END_REF] modified the red reflectance to differentiate between the blue and red bands. Thus, Soil and Atmosphere Resistant Vegetation Index (SARVI) can effectively reduce the effects of the atmosphere effect, especially aerosols. However, this VI still requires the remote sensing images undergoing atmosphere calibration by the 5S [START_REF] Tanred | Description of a computer code to simulate the satellite signal in the solar spectrum: The 5s code[END_REF] atmospheric transport model. In addition, to achieve the calibration, the corresponding atmospheric parameters are required, which are difficult to obtain. Without atmospheric calibration of remote sensing images, this index is not expected to perform better than the NDVI.

Huete et al. [START_REF] Huete | A comparison of vegetation indices over a global set of TM images for EOS-MODIS[END_REF] further proposed an Enhanced Vegetation Index (EVI) for the elimination of both atmosphere and soil noises nir red

nir 1 red 2 blue EVI EVI R R G R C R C R X - = + - + ••••••••••••••••••••••• (1.12)
where Rblue is the atmosphere corrected surface reflectance in the blue band. G is the gain factor and is equal to 2.5. C1 and C2 are aerosol adjustment factors (C1 = 6, C2 = 7.5) that correct the aerosol effect for the red band reflectance using the blue band reflectance. XEVI is the soil adjustment factor initially equal to 1 [START_REF] Huete | A comparison of vegetation indices over a global set of TM images for EOS-MODIS[END_REF].

Jiang et al. [START_REF] Jiang | Development of a two band enhanced vegetation index without a blue band[END_REF] established a two-band EVI without the blue band (EVI2) as an alternative to EVI with three bands, defined as

nir red nir 1 2 red EVI EVI2 ( / ) R R G R C C c R X - = + - + ••••••••••••••••••••• (1.13)
where C1, C2 and XEVI have the same values as in EVI. c is the band correlation coefficient designed to achieve the best similarity between EVI and EVI2, c = Rred / Rblue. Jiang et al. [START_REF] Jiang | Development of a two band enhanced vegetation index without a blue band[END_REF] derive c equal to 2.08 by fitting EVI2 from MODIS data to EVI. (2) Directionality Index

The studies show that multi-angle remote sensing is more valuable than fixed-angle remote sensing in obtaining three-dimensional structural information of vegetation. For example, Baret et al. [START_REF] Baret | GAI estimates of row crops from downward looking digital photos taken perpendicular to rows at 57.5 zenith angle: Theoretical considerations based on 3D architecture models and application to wheat crops[END_REF] found that the highest accuracy of LAI is obtained when the zenith angle is 53.5°. Stavros et al. [START_REF] Stagakis | Monitoring canopy biophysical and biochemical parameters in ecosystem scale using satellite hyperspectral imagery: An application on a Phlomis fruticosa Mediterranean ecosystem using multiangular CHRIS/PROBA observations[END_REF] use CHRIS / PROBA data and found that the VIs calculated using reflectance obtained from off-nadir angles are better than those obtained from the nadir direction. Breunig et al. [START_REF] Breunig | Directional effects on NDVI and LAI retrievals from MODIS: A case study in Brazil with soybean[END_REF] used MODIS data to invert the LAI and found that the inversion accuracy of the pixels in the backward observation area is higher than that of the pixels in the MOD15 alternate algorithm (LAI-NDVI empirical model) for the inversion of the LAI. The common LAI-NDVI empirical model used by MODIS products. Therefore, many studies have been carried out to calculate VIs using the reflectance obtained from off-nadir observation directions, which led to the development of directionality indices.

(i) Angular indices

The reflectances of the canopy have a clear directional character. For example, in the backward direction with the same view zenith angle as sun zenith angle, very high reflectance is observed and is named as "hot spot"; in the forward direction with the same view zenith angle as sun zenith angle, a very low reflectance is observed and is named as "dark spot" (Figure 1.1). The hotspot is observed because when the observation direction is the same as the sun direction, the sensor observed the least proportion of the shadows; The darkspot is observed because when the observation direction is opposite to the sun direction, the sensor observed the largest proportion of the shadows. HDS [START_REF] Lacaze | Retrieval of vegetation clumping index using hot spot signatures measured by POLDER instrument[END_REF] is defined as

HS DS DS HDS R R R - = •••••••••••••••••••••••••••••••••• (1.14)
where RHS is the reflectance at the hotspot and RDS is the reflectance at the darkspot. As hotspot and darkspot are challenging to sample by sensors routinely, the model can also be used to fit a bidirectional reflectance coefficient to existing observed data and calculate the hotspot and darkspot values [START_REF] Lacaze | Retrieval of vegetation clumping index using hot spot signatures measured by POLDER instrument[END_REF][START_REF] Leblanc | Mapping vegetation clumping index from directional satellite measurements[END_REF][START_REF] Chen | Multi-angular optical remote sensing for assessing vegetation structure and carbon absorption[END_REF].

NDHD [START_REF] Chen | Multi-angular optical remote sensing for assessing vegetation structure and carbon absorption[END_REF][START_REF] Chen | Predicting leaf area index in wheat using an improved empirical model[END_REF] is defined as

HS DS HS DS NDHD R R R R - = + •••••••••••••••••••••••••••••••••• (1.15)
Reflectance tends to be minimal at darkspot, where shadows from vegetation can be observed to the greatest extent possible [START_REF] Lacaze | Retrieval of vegetation clumping index using hot spot signatures measured by POLDER instrument[END_REF]. The reflectance in the red band depends on the physiological conditions of the green leaves, especially at darkspot, where reflectance is low. It explains the potentially significant errors in calculating NDHD [START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF]. Therefore, we use the reflectance in the NIR band to calculate NDHD. The main difference between HDS and NDHD is that the linear correlation with NDHD and the clumping index is better than that of HDS, suggesting that NDHD may be a better index than HDS [START_REF] Chen | Multi-angular optical remote sensing for assessing vegetation structure and carbon absorption[END_REF].

(ii) Hotspot signature vegetation index

As a type of directionality VIs, the hotspot signature VI takes the form of a multiplication of the nadir observed VI and the angular index. The best known of these is the normalised hotspot signature vegetation index (NHVI) [START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF], which is defined as

NHVI = NDVI HDS × ••••••••••••••••••••••••••••••• (1.16)
The reason for considering the clumping index in VIs is that clumping indirectly affects the relationship between LAI and VIs through the fraction of absorbed photosynthetically active radiation (FAPAR). The angular indices (HDS, NDHD) are negatively linear to the clumping index, whereas FAPAR shows an excellent linear relationship with various VIs (NDVI, SAVI, ...) [START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF][START_REF] Chen | Estimation of winter wheat LAI using hotspot-signature vegetation indices[END_REF]. In addition, the angular indices show good sensitivity to LAI even in high LAI regions [START_REF] Chen | Predicting leaf area index in wheat using an improved empirical model[END_REF].

Thus, based on measured data in boreal forests, researchers have found that NHVI, generated by the product of NDVI and HDS, better improves linear correlation with LAI [START_REF] Chen | Predicting leaf area index in wheat using an improved empirical model[END_REF][START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF][START_REF] Chen | Estimation of winter wheat LAI using hotspot-signature vegetation indices[END_REF].

Although the formations of vegetation indices are enormous, they are suffering from the saturation effect to varying degrees. Therefore, studies should be carried out to alleviate the saturation effect of vegetation indices. It is one of our study objects. Besides, even though vegetation indices are designed to enhance the optical signal of vegetation, they are inevitably affected by environmental factors. Therefore, these influencing factors should be quantitatively analyzed before using VIs.

Problems

The three-dimensional radiative transfer model is the theoretical basis and a powerful tool for performing a component spectral inversion. It also points out three issues facing the simulation and inversion studies.

(1) How can the 3D vegetation remote sensing scenes be accurately modelled? Threedimensional remote sensing scenes are the necessary basis for modelling reflectance using threedimensional radiative transfer models. Current research on remote sensing modelling lacks knowledge of the spatial and temporal variation of critical parameters of the ground and cannot achieve continuous-time phases simulation. For example, in China's extensive crop growing areas, crop growth leads to significant changes in plant structure and biochemical parameters and significant changes in the spatial distribution of crops. Therefore, maize growth processes and models must be introduced to achieve continuous temporal phase data simulation.

(2) How can we accurately invert the OPs of leaves in areas with high heterogeneity and complex component distribution? Much of the current research assumes that the canopy is homogeneous and that the pixel contains only two components, the leaf and the soil. This assumption applies to areas of high FC. However, the spectra that make up the pixel's reflectance may originate from the vegetation, soil and the surrounding artificial structures in the urban area. In this case, the current inversion method would no longer be applicable, and spectral unmixing of the pixels would need to be considered.

Outline of the thesis

Based on the above problems, three study contents are proposed below: (

2) Leaf optical property inversion based on radiative transfer model at the sub-pixel scale

We proposed to use DART calibration to accurately invert the OPs of leaves in areas with high heterogeneity and complex component distribution. Leaf OP inversion from the mixed pixel in the urban area is achieved based on DART calibration. Based on the DART calibration, inversion of leaf OPs at sub-pixel scale in urban areas is achieved. DART simulated images are used for validation with two cases, including ideal and non-ideal cases. Only the OPs are assumed to be unknown in the ideal scene, and all other radiative transfer simulation parameters are known. In this case, an inversion accuracy analysis is performed. In the non-ideal scene, a certain amount of error is assumed to exist in the model input parameters, and sensitivity analysis of the error and an inversion accuracy assessment are performed. This method is also applied to the Basel area to obtain spectral data of leaves, soils and buildings.

( The relationships between the above chapters are visualized in Figure 1.2.

Thesis structure

The thesis is divided into 7 chapters: For the inversion part, it can then be divided into leaf OP inversion and LAI inversion part.

Chapter 2 Basics of vegetation canopy radiative transfer models

The syntax rules of the ELSYS are introduced, and the programming ideas for modelling real structural trees using the ELSYS based on fractal typology are described.

(i) Leaf area index

When radiation is transmitted in the medium, the influence is greatly related to the density distribution of scatters and absorbers. For vegetation, it is the leaf.

The leaf area index is defined as the one-sided area of a leaf per unit background area, i.e. the integral of the density distribution function ul(z) of the leaf area in the vertical direction (H)

1 0 LAI ( )d H u z z = ∫ ••••••••••••••••••••••••••••••••••• (2.1)
where the upper limit of integral H is the depth of vegetation canopy and the orientation of z is downward (z = 0 is the upper bound of vegetation and z = H is the lower bound of vegetation). Leaf area density refers to the sum of leaf (single-sided) area per unit volume. Its spatial distribution is called leaf area density distribution, and the unit is m -1 .

There are many models to describe ul(z). For example, under the assumption of plane parallel distribution of vegetation, the leaf area density distribution of many canopies tends to be higher and closer to the top of the canopy, and the simplest assumption is that the density function is constant along with the canopy height [such as ul(z) = L / H].

(ii) Leaf projection function

The radiation transfer process of vegetation has a great relationship with the scattering and absorption medium-leaf orientation, unavailable in other fields. Therefore, geometric functions G(Ω) is defined to describe the average projection of the normal direction of leaves in the Ω direction.

( )

2π 1 1 1 1 1 0 0 1 ( ) d 2π G g Ω = Ω Ω ⋅ Ω Ω ∫ ∫ ••••••••••••••••••••••••• (2.2)
where gl(Ωl) is the probability density of leaf normal distribution, representing the probability of a blade in a unit solid angle near Ωl in its normal direction (taking the unilateral normal direction of its upper hemisphere space). gl(Ωl) has normalization condition

l l l 2π 1 ( )d 1 2π g + Ω Ω = ∫ •••••••••••••••••••••••••••••••• (2.
3)

The integral region 2π + is the upper hemisphere space because the blade can only calculate one side.

Based on Eq.( 2.3), we get

π 2π l l l 2 l l l l l l 2π 0 0 ( , ) 1 ( )d sin d 1 2π 2π g g d θ ϕ ϕ θ θ + Ω Ω = = ∫ ∫ ∫ •••••••••••••••• (2.4)
From Eq. (2.4), we get

π 2 l l l l l 0 ( , ) sin d 1 g θ ϕ θ θ = ∫ ••••••••••••••••••••••••••••••• (2.5)
It is generally assumed that the zenith angle and azimuth angle of the normal distribution of the leaf are not correlated and uniformly distributed along the azimuth direction, so

l l l l l ( , ) ( ) g g θ ϕ θ = ••••••••••••••••••••••••••••••••••• (2.6)
Based on Eq.(2.5) and Eq.(2.6), we get

π 2 l l l l 0 ( ) sin d 1 g θ θ θ = ∫ ••••••••••••••••••••••••••••••••• (2.7)
Define g l * (θ l )=g l (θ l )sinθ l as LAD. The choice of LAD depends on the vegetation type. Some common forms of LAD are listed in Table 2.1.

(iii) Scattering phase function

Like other radiative transfer theories, scattering phase function is also defined in vegetation, denoted as function. The function is also related to the leaf orientation at the scattering point and is not normalized.

Firstly, the scattering phase function γl(Ωl, Ω'→Ω ) is introduced, and the ratio of scattering to Ω direction is expressed when the radiation in Ω' direction is incident to the blade with normal orientation Ωl.

If the scattering characteristics of the blade can be seen as two different radius reflection and transmission hemispheres, namely:

Incoming flux E + can be written as:

0 =π cos ' E L α + •••••••••••••••••••••••••••••••••••• (2.8)
Projection of reflection radiance Lr(Ωl) is ( )

l 0 l cos ' c s = | o | r L L r α α Ω •••••••••••••••••••••••••••••• (2.9) Projection of transmission radiance Lt(Ωl) is ( ) l 0 l co | s | s ' co t L L t α α = Ω ••••••••••••••••••••••••••••• (2.10)
Based on Eq. (2.8), Eq. (2.9) and Eq. (2.10)

l l l l l ( ) 1 cos , cos cos 0 π ( ) 1 cos , cos cos 0 π r t L r E L t E α α α γ α α α + + Ω  ′ = <   =  Ω  ′ = >   •••••••••••••••••• (2.11)
where α'=cos -1 (Ω'•Ωl) is incidence angle, α=cos -1 (Ω•Ωl) is exit angle, rl is the leaf reflectance, tl is the leaf transmittance.

In order to characterize the scattering characteristics of leaf groups, a function Г must be introduced. At the scattering position z, the differential probability of the blade with normal Ωl is gl(z, Ωl) • dΩl. When incidents with Ω', the incident intensity must also multiply by the factor |cosα'|.

Therefore, we introduce:

l l l l l 2π 1 1 ( , ' ) ( , ) cos ' ( , ' )d π 2π z g z α γ + Γ Ω → Ω = Ω Ω Ω → Ω Ω ∫ ••••••••••• (2.12)
If the leaf has double-hemisphere scattering characteristics, the group scattering phase function

is l l l l l l l l 1 1 ( , ' ) ( , ) cos cos d ( , ) cos cos d 2π 2π z g z t g z r α α α α + - Ω Ω ′ ′ Γ Ω→ Ω= Ω Ω- Ω Ω ∫ ∫ (2.13)
where the integral region Ω ± satisfies ±cosα cosα' > 0, and Ω + + Ω -= 2π + .

(iv) Leaves size

In addition to LAI and LAD, another feature to be considered in canopy radiative transfer modelling is the size of the leaves, which is often used as a driver of the hotspot correction function.

It is because the hotspot effect is mainly explained based on shadow shading theory, whereas statistically speaking, the number of gapes is the same as the number of leaves, and the size and number of leaves uniquely determines the size of the gapes.

( The system of SAIL equations takes the following form:

S S dE KE dz = ••••••••••••••••••••••••••••••••••• (2.14) S ' dE aE bE c E dz - - - = - - ••••••••••••••••••••••••••••• (2.15) S dE bE aE cE dz + - + = - + •••••••••••••••••••••••••••••• (2.16) 0 S 0

dE uE vE wE kE dz

+ - = + + - ••••••••••••••••••••••••••• (2.17

Discrete Anisotropic Radiative Transfer Model

The Scattering in a turbid medium causes a change in the initial direction of the light striking the medium.

Whether using the turbid medium or triangular facets to fill the canopy depends primarily on the size of the trees relative to the resolution of the simulated scene (i.e., the voxels size that make up the scene). The triangular facets can be independent of the unit voxels that make up the scene, whereas the turbid medium depends on the unit voxels. For example, if the size of the simulated vegetation is much smaller than the resolution of the simulated image (i.e., the size of the constituent scene voxels), then it is preferable to use triangulated surfaces for filling. If the size of the simulated vegetation is much larger than the resolution of the simulated image, then both are possible. When the scene resolution is too large, and the canopy size is too small, care needs to be taken when filling with turbid media because no corresponding vegetation voxel can be created when the foliage makes up less than 50% of the unit voxel.

DART supports three approaches to modelling the vegetation in a scene. The first is an entirely homogeneous ideal voxel consisting of a turbid medium or triangular surface, referred to as a homogeneous assumption canopy. The second is a tree with a simple geometric shape (which includes a canopy modelled by a voxel filled with a turbid medium or a discrete triangular surface), and the vegetation scene constructed in this way is said to be a simplified heterogeneous canopy.

The third is a three-dimensional real structural canopy composed of facets, and the vegetation scenes constructed in this way are called real structural canopies. All three vegetation modelling approaches are addressed in this study.

The simplified heterogeneous canopy scene creates a tree consisting of two parts: the trunk and the crown. The concept of "species" is used to manage the different tree types. Trees of the same species share some of the same characteristics, such as the way the canopy is modelled, the OP of the tree, the vertical hierarchical stratification of the defined canopy, the density of the leaf area volume, the size and shape of the canopy. The trunks are vertically stacked in an octagonal shape, both within or below the canopy. The crown can be defined in different shapes: ellipsoid, compound ellipsoid, truncated cone, trapezoid, and compound cone. For canopies with hierarchical layers in the vertical direction, the relative height of the layer, the relative leaf area voxel density of the layer, the relative trunk diameter, the OP of the leaves and trunk, the proportion of clumping in the horizontal direction, and so on can be defined on each hierarchical layer. There are two types of infill methods within the canopy. The first method is the turbid medium canopy layer, which is composed of a turbid medium and void voxel. In the DART model, a voxel is created as a turbid medium voxel only when the percentage of leaves within the voxel reaches 50% or more of the total voxel. The user specifies the parameters that make up the turbid medium voxels, such as the leaf area index or leaf area voxel density, the branch area index, the leaf inclination index, etc. In addition, the proportion of the whole leaf voxels occupied by the leaf cannot exceed 100%, and any portion exceeding 100% will be preserved in the surrounding neighbourhood of other voxels. The second infill method is a discrete isosceles triangular facet, where the triangle's height is twice the height of the bottom edge, and the user can specify the total area or the total number of triangles within the canopy. It has been shown that whether or not short branches are simulated within the canopy has little effect on the total BRF [START_REF] Banskota | An LUT-based inversion of DART model to estimate forest LAI from hyperspectral data[END_REF]. Similar filling rules can also be applied to the creation of vegetation plots.

DART's graphical user interface (GUI) provides three modes for users to set the spatial distribution, size and species of trees [START_REF] Gastellu-Etchegorry | DART: recent advances in remote sensing data modeling with atmosphere, polarization, and chlorophyll fluorescence[END_REF][START_REF] Gastellu-Etchegorry | Simulation of satellite, airborne and terrestrial LiDAR with DART (I): Waveform simulation with quasi-Monte Carlo ray tracing[END_REF][START_REF] Gastelluetchegorry | Modeling specular reflectance and polarization in DART model for simulating remote sensing images of natural and urban landscapes[END_REF][START_REF] Gastelluetchegorry | Time series of urban radiative budget maps derived from EO satellites using a physical remote sensing model[END_REF][START_REF] Gastelluetchegorry | Dart: A Tool For Studying Earth Surfaces -Time Series of Urban Radiative Budget From Eo Satellites[END_REF][START_REF] Gastelluetchegorry | Simulation of Chlorophyll Fluorescence for Sun-and Shade-Adapted Leaves of 3D Canopies with the Dart Model[END_REF]. The first requires the user to enter the exact location and the exact size of each tree. Input parameters (i.e., location and size) need to be written in a text file, and this text file can be specified in the GUI of the DART. An example file called "trees.txt" is available in the DART's database folder and can be used as a basis for user modifications. In addition, the user can use some of the output parameters of the forest growth model to obtain the required geometry parameters in "trees.txt" and use them as input parameters to the DART model for the forest scene simulation. The second method requires the input of exact locations and random dimensions, where the location fields are entered as in the first method using a file and only the "Species_ID, POS_X, POS_Y" fields are retained, and the mean and standard deviation of the dimension fields are set in the DART GUI. The third method uses random positions and random dimensions. In this mode, no files are used, and the user can specify a small area (ellipse or rectangle) or the whole scene as a tessellation and define the intervals Δx and Δy between neighbouring nodes. The simplified heterogeneous forest canopy simulated in section 2.2.2 was created from the DART model using this method.

In addition to the tree models created by the GUI, DART also supports 3D realistic structural plants constructed by other plants' modelling software using facets, and the 3D plants composed of facets are also resolution independent from the simulated scene.

( In addition to these four modules, DART includes a simulation sequencer and many tools, such as the ability to simulate broadband satellite imagery, generate schematic Digital Terrain Model (DTM), integrated bidirectional surface reflectance models from Hapke [START_REF] Hapke | Bidirectional reflectance spectroscopy: 1. Theory[END_REF] and RPV [START_REF] Rahman | Coupled surface-atmosphere reflectance (CSAR) model: 2. Semiempirical surface model usable with NOAA advanced very high resolution radiometer data[END_REF], Fluspect and PROSPECT [START_REF] Jacquemoud | PROSPECT: A model of leaf optical properties spectra[END_REF] optical leaf models, etc. The format of the input data depends on the nature of the data: Structured Query Language (SQL) databases of OPs of materials, solar irradiance and physical properties of the atmosphere, obj files of 3D scenes (e.g., city geometry databases, trees, etc.), simple text files of solar irradiance measured in the field, etc. On the other hand, the DART model is controlled via a GUI or directly via scripts written in Python.

The N-flux tracking method used by DART solves the radiative transfer equations along all previously defined discrete directions. The general equation describing the propagation of a singleband stationary electromagnetic wave of wavelength λ (natural light is the sum of monochromatic waves), at position r along direction Ω, with luminosity L(r, Ω) in a medium of temperature T is described as (for simplicity of notation, no spectral dependence is specified)

( ) ( ) e a B d 4π d d d ( ,Ω) d d d , Ω Ω ( ,Ω) ( ,Ω) ( ,Ω) ( ) ( ,Ω) , Ω dΩ 4π L r x y z P r r L r r L T r L r ξ η µ α α α   + + =     → - + ′ ′ + ′ ∫ •••• (2.18)
where ξ, η and µ are the directional cosines of the propagation direction along the x, y and z axes. αa(r,Ω) is the absorption coefficient of the medium. αd(r,Ω) is the diffusion coefficient of the medium. αe(r,Ω) is the extinction coefficient of the medium (diffusion + absorption). P(r,Ω'→Ω) is the normalised phase function. LB(T) is the Planck illuminance. In the case of discrete directions, this equation can be written as [START_REF] Peddle | Structural change detection in a disturbed conifer forest using a geometric optical reflectance model in multipleforward mode[END_REF] where n and m are both directions, and Ndir is the total number of directions, L(r,Ωn) is the radiance of the point r in the direction of Ωn. To ensure a suitable discretisation of the equation, the value of the solid angle ΔΩn must be sufficiently low, depending on the degree of anisotropy of the phase function and the inhomogeneity of the medium under consideration. The flux-tracking mode consists of light rays that carry a stream of spectral energy in a limited number of directions. In heterogeneous 3D scenes [START_REF] Gastellu-Etchegorry | Modeling radiative transfer in heterogeneous 3-D vegetation canopies[END_REF] and atmospheres [START_REF] Grau | Radiative transfer modeling in the Earth-Atmosphere system with DART model[END_REF], it is based on an exact kernel and discrete longitudinal coordinate approach using iterative and convergent methods. All or part of the radiation intercepted by the scene elements in the i th iteration is diffused in the next i+1 st iteration. The iterative process stops when the relative difference in scene exitance between two successive iterations is less than a previously specified threshold. In addition, if the angular energy of the light is less than the average angular energy of the scene diffused in the first iteration multiplied by a factor chosen by the user, the light is stopped.

( ) ( ) ( ) ( ) ( ) ( ) ( ) dir e a B d 0 d d d ,Ω ,Ω ,Ω ,Ω ( ) d d d ,Ω Ω ,Ω ,Ω ΔΩ 4π n n n n N m n n m m m L r r L r r L T x y z P r r L r ξ η µ α α α =   + + = - +     → + ∑ (2.
The general scheme for modelling remote sensing sensor measurements includes five main steps: (1) Calculation of scene illuminance: the transmission and scattering of solar radiation and thermal emission from the atmosphere give a direct and diffuse illuminance of the scene under consideration (at BOA, Top of the atmosphere (TOA) and other specified sensor levels in the atmosphere). ( 2) Calculation of the radiative transfer in the model: the model and its elements iteratively scatter the radiation they intercept and the thermal radiation they emit. The radiation balance and images can be stored at the end of each iteration. (3) Calculation of the coupling irradiance: the exitance of the landscape at the BOA level results in backscattered energy from the atmosphere. (4) Calculation of the Earth-Atmosphere Coupling: the energy backscattered by the atmosphere is scattered back by the Earth scene. This scattering is calculated with a single iteration, but an extrapolation based on step 2 allows multiple scatterings to be considered. The Earth-Atmosphere coupling has been successfully cross-comparisons with simulations of the MODTRAN atmospheric radiative transfer model. ( 5) Calculation of radiation at TOA and sensor levels: transfer of upwelling BOA radiation at TOA and sensor levels into the atmosphere. A more detailed theoretical presentation of the DART model can be found in the literature [START_REF] Gastellu-Etchegorry | Modeling radiative transfer in heterogeneous 3-D vegetation canopies[END_REF].

The DART model was validated by comparison with measurements [START_REF] Gastellu-Etchegorry | Modeling BRF and radiation regime of boreal and tropical forests: I. BRF[END_REF] and other 3D reflection models such as FLIGHT [START_REF] North | Three-dimensional forest light interaction model using a Monte Carlo method[END_REF], Spreading Of Photons For Radiation Interception (SPRINT) [START_REF] Thompson | Two models for rapidly calculating bidirectional reflectance of complex vegetation scenes: Photon spread (PS) model and statistical photon spread (SPS) model[END_REF], Raytran [START_REF] Govaerts | Raytran: A Monte Carlo ray-tracing model to compute light scattering in three-dimensional heterogeneous media[END_REF] in Radiative Transfer Model Inter Comparison (RAMI) studies. It is being and has been used successfully in many scientific applications, including the development of satellite or inversion methods of aircraft reflectance images [START_REF] Gascon | Retrieval of forest biophysical variables by inverting a 3-D radiative transfer model and using high and very high resolution imagery[END_REF][START_REF] Banskota | Investigating the utility of wavelet transforms for inverting a 3-D radiative transfer model using hyperspectral data to retrieve forest LAI[END_REF], satellite sensor design such as NASA DES Dynl, CNES Pleiades, CNES LIDAR missions [START_REF] Durrieu | Preliminary studies for a vegetation ladar/lidar space mission in france[END_REF], estimating the effect of canopy structure on satellite image texture [START_REF] Gastellu-Etchegorry | Modeling BRF and radiation regime of boreal and tropical forests: I. BRF[END_REF][START_REF] Bruniquel-Pinel | Sensitivity of texture of high resolution images of forest to biophysical and acquisition parameters[END_REF], modelling the 3D distribution of photosynthesis in plant canopies [START_REF] Malenovský | Retrieval of spruce leaf chlorophyll content from airborne image data using continuum removal and radiative transfer[END_REF][START_REF] Malenovský | Influence of woody elements of a Norway spruce canopy on nadir reflectance simulated by the DART model at very high spatial resolution[END_REF][START_REF] Malenovský | Discrete anisotropic radiative transfer modelling of solarinduced chlorophyll fluorescence: Structural impacts in geometrically explicit vegetation canopies[END_REF], designing of new chlorophyll estimation indices for coniferous forests [START_REF] Malenovský | Retrieval of spruce leaf chlorophyll content from airborne image data using continuum removal and radiative transfer[END_REF],

in tropical forest research [START_REF] Barbier | The variation of apparent crown size and canopy heterogeneity across lowland Amazonian forests[END_REF][START_REF] Barbier | Linking canopy images to forest structural parameters: potential of a modeling framework[END_REF][START_REF] Proisy | Biomass prediction in tropical forests: the canopy grain approach[END_REF], etc.

Radiosity-graphics based model (1) Radiosity-graphics based model scene building

RGM uses the ELSYS to build remote sensing scenes. The ELSYS originated at the State University of New York (SUNY), is a language interpretation system written in C. It is similar to the compilers widely used in programming. ELSYS accepts source files with the suffix ".lsy" for compilation. The outputs are files with the suffix ". in", containing the vertex 3D coordinates of the facets that make up the remote sensing scene. ELSYS is an executable program in a Disk Operating System (DOS) environment and therefore has no code editing interface. Any text editing software (Textbook, Notepad++, etc.) can be used to edit the input file (*.lsy) and output files (*.in). Besides, it supports compilation errors. The compilation process is similar to that of other compilers and consists of: word scanning → assignment statement identification → syntax pre-analysis → syntax analysis → variable separation → instruction sequence generation → instruction execution.

ELSYS uses facets to build realistic structures of trees and therefore needs to model the realistic 3D structure of trees using its specific syntax rules through fractal principles. The grammar rules of ELSYS are presented in Appendix.A. Here the basic programming idea is explained using the example code provided by ELSYS for a broadleaf tree: the structure of a tree is fractal in nature, with the tree consisting of roughly two parts: the crown and the trunk. The crown is made up of different levels of branches, and the branches are made up of poles, twigs and leaves. Therefore, the programming is based on the idea of moving from the general to the local, abstracting the whole scene into major parts to build the scene, and then replacing the general part with details to subdivide the scene. As shown in Figure 2.3, the scene is divided into the ground part and the tree part. The ground part is made up of many small squares. A predecessor "soil" statement is used to determine the ground space position, and then the "soil" statement is replaced with a subsequent "quarter" statement to determine precisely how many small squares consist of the ground. Finally, the "quarter" statement is replaced with a subsequent "tile" statement to draw each small square. The tree part is determined in the "tree" statement by determining the tree's exact location in space and then dividing the tree into two major parts: the "crown" statement and the "trunk" statement. The trunk is first determined in the predecessor "trunk" statement by determining the number of octagons is stacked vertically, then the "trunk" statement is replaced by the successor "layer" statement to specify each of the trapezoidal sides of the octagons, all of which are enclosed together to form the bark. In contrast, the bottom and top sides of the octagons are not drawn. Finally, the "layer" statement is replaced by a subsequent "facet" statement to draw each trapezoid. The "crown" statement divides the crown into three layers based on the vertical height of the trunk: top, middle and bottom. The branches inclination angle varies at each level: the branches at the top and bottom of the crown extend along the zenith direction, while the branches in the middle of the crown extend horizontally. The "crown" statement is replaced by a subsequent "branch" statement to specify the length of each first-level branch in each of the three tiers crown. The middle of the crown is longer than the top and bottom of the crown. The succeeding "branch1" statement replaces the "branch" statement to implement the combination of primary and secondary branches (a secondary branch "twig" statement grows outward between every two primary branch poles), the "pole" statement and the "twig" statement are used to draw the square to represent the branches. In addition, the "twig" statement is required to determine the position of the leaves in space. The subsequent "leaf" statement replaces the "twig" statement to draw the leaves on the secondary branch. The leaf consists of two parts: represented by a rectangle (the "stem" statement) and consisting of triangles and trapezoids (the "leaf2" statement).

In addition, the primary branches grow on the main trunk, and secondary branches grow on the primary trunk, with all leaves growing on the secondary trunk only. The leaves grow on the twig, which grows on the secondary trunk. All the branches are consisted of rectangles, with a secondary branch growing between every two primary branches. In summary, the scene can be simplified to the form shown in Figure 2.4. The single tree remote sensing scene in section 2.2.1 was modelled using this method based on ELSYS. The radiance Bi for a single differential facet dSi can be defined as

2 d d 1 d d d 1 ,
, 2, , 2

j i N i i i i i j S S j S j j B S E S B F S i N χ → = = + ∑ = … ∫ •••••••••••••••• (2.20)
where Bi is the radiance (W•m -2 ) per unit facet dSi. N is the total number of facets that make up the scene, and since each facet has both positive and negative faces, the total number of faces to be counted is 2N. χi represents the OPs per unit facet i, which can be expressed as reflectance ρi or transmittance τi, depending on the relative positions of facet i and facet j ( ) ( )

,if 0 ,if 0 i i j i i i j ρ χ τ  ⋅ <  =  ⋅ >   n n n n ••••••••••••••••••••••••••••••• (2.21)
where ni is the normal vector of facet i, and nj is the normal vector of facet j. Ei is the emitted flux density that the sunlight irradiates to facet i such that the facet has, which is defined as

1 ( , )d d , 1, , 2 i i i i S E E x y x y i N S = = … ∫ ••••••••••••••••••••••••• (2.22)
where Si is the total area of the facet i, (x,y) is a point in the local coordinate system defined in terms of facet i, and Ei(x,y) is the emission flux density at that point, which is defined as 0 ( , ) HID( , ) , 1, ,2

i i i E x y x y E i N χ = ⋅ = … n s ••••••••••••••••••• (2.23)
where HID(x,y) is the occlusion coefficient, which is 0 if there is another facet between position (x,y) of facet i and the sun, and 1 otherwise (no occlusion). E0 is the incident solar flux density. s is the direction of solar incidence.

d d j i S S
F → is the view factor, which represents the proportion of the radiation flux density (unitless) leaving the infinitesimal surface dSj and reaching the infinitesimal facet dSi, and is calculated as

d d 2 cos cos d d 1 π j j i i j i j i S j S S S S S F S r θ θ → = ∫∫ ••••••••••••••••••••••••• (2.24)
where θi and θj are the angles formed by the line between the normal and facet i and facet j, and r is the distance between the two facets. A more detailed description of the principle of the radiosity model and its derivation can be found in the literature [START_REF] Borel | The radiosity method in optical remote sensing of structured 3-D surfaces[END_REF].

Simulating Canopy BRF with Radiosity-Graphics based Model (RGM) at Pixel Scale based on OpenACC

As one type of computer simulation model, the facet-based radiosity fully preserves the spatial correlations of the canopy, the gap fraction and clumping effect and so on within the canopy in simulating the radiative transfer between facets, which is an excellent solution to the shortcomings of traditional 1D radiative transfer models [START_REF] Gerstl | Principles of the radiosity method versus radiative transfer for canopy reflectance modeling[END_REF]. However, the view factor calculation has been the main bottleneck limiting the radiosity models. Since each facet needs to calculate the view factor with all other facets in the scene, if the number of facets making up the scene is too large (50,000 facets), its calculation will be very slow or even crash [START_REF] Huang | RAPID: A Radiosity Applicable to Porous IndiviDual Objects for directional reflectance over complex vegetated scenes[END_REF]. It directly limits the application of radiosity models in the field of remote sensing. For this reason, Huang et al. [START_REF] Huang | RAPID: A Radiosity Applicable to Porous IndiviDual Objects for directional reflectance over complex vegetated scenes[END_REF] proposed the Radiosity Applicable to Porous Individual Objects (RAPID) model by simplifying the facets that make up the scene and using porous individual thin objects as a basic unit instead of facets. This model has been successfully applied to simulate canopy reflectance at the pixel scale [START_REF] Zhen | Simulation and Sensibility Analysis of Earth Surface Bidirectional Reflectance Factor Based on Radiosity[END_REF].

However, RAPID cannot be applied to real 3D structural remote sensing scenes. Based on the GPU acceleration technology, the calculation part of the view factors of the RGM model [START_REF] Xie | Influences of leaf-specular reflection on canopy BRF characteristics: A case study of real maize canopies with a 3-D scene BRDF model[END_REF] was rewritten in parallel to increase the simulation speed of the RGM model while increasing the maximum number of computable facets. Based on both real and simplified canopy scenes, the simulated reflectance was verified in both the principal plane (PP) and the cross-principal plane (CP) observation directions.

Reflectance simulation of realistic structural single tree canopy

The purpose of the real structural single tree validation was to evaluate the computational speed of the improved RGM model. Due to the limited computational power of the original RGM model, we selected a previously used single tree constructed with a certain number of facets as the research target. The previously used single young elm tree in Chaoyang Park, Changchun, Jilin Province [START_REF] Meng | Simulation and Analysis of Tree Canopy BRF Based on the Real Scene[END_REF], was selected as the study subject (Figure 2.6.a), and the detailed data collection and processing flows can be found in the original literature.

The basic geometric parameters of the elm were obtained using multiple measurements for averaging, including basic parameters such as total tree height, trunk height, trunk radius, crown length, crown width, number of levels, branch length and branch inclination, number of twig leaves, leaf size and shape (Table 2.2). The length units were measured using a leather ruler and straightedge, and the angle units were measured using a protractor and compass. LAI and ALA were measured using the LAI-2000 plant canopy analyzer. Three-dimensional modelling of a single elm was achieved using the collected geometric structure parameters as input variables to the ELSYS (Figure

Reflectance simulation of simplified heterogeneous canopy scene

The improved RGM model has substantially increased its computational power, allowing for simulations at the pixel scale. Therefore, the DART model, also computationally powerful, was chosen to validate the simulation results. The pixel-scale validation was based on the official DART work package (WP6B). The work package can be downloaded via the website (https://dart.omp.eu/#/doc). The OPs of the components were derived from the corresponding components in the DART database (Table 2.5). The input parameters were slightly modified to fit the study objectives (Table 2.6). The 3D scene at the pixel scale was a mixed scene consisting of trees and soil, where the trees consist of trunks and canopies. The trunk was a square octagon, and the crown was filled with small triangles distributed homogeneously. All trees had identical geometry and were randomly distributed in the horizontal scene (Figure 2.9). The entire 3D scene was generated from the DART model and inputted into the RGM model.

Since the total number of facets that make up the scene was 580817, which was approximately eight times the amount of operations that the original RGM could handle [START_REF] Huang | RAPID: A Radiosity Applicable to Porous IndiviDual Objects for directional reflectance over complex vegetated scenes[END_REF]; therefore, the DART model was used to validate the pixel scale scene. From Figure 2.10, there was a substantial similarity between the RGM simulation results and the DART simulation results, with the overall shape of the curves being nearly the same and the trend of the BRF with the observed azimuth being consistent with each other overall. The best agreement was found in the NIR band, and the worst was found in the blue band. In the blue band, the BRF simulated by the RGM model was slightly lower than that of the DART model.

Summary

As a basis for the subsequent chapters, three radiative transfer models (SAIL, RGM and DART) are introduced briefly. Two methods to build the 3D remote sensing scene, including ELSYS and DART GUI, are introduced briefly. One improvement method of RGM using OpenACC to lift calculated load and running speed is introduced.

(1) The RGM model uses the ELSYS platform to build the 3D scene, an L-system-based 3D modelling platform similar to a compiler but without a GUI. The RGM considers the flux density leaving the facets consisting of the reflected flux density, the transmitted flux density, and the emitted flux density after sun irradiation and calculates the radiative transfer.

(2) The DART model supports using a voxel or a facet as the basic unit to build a scene.

Several simple 3D scenes can be generated in the GUI of the DART model. In addition, the DART model supports the import of 3D scenes composed of facets from external sources. The DART model uses the N-flux tracking method to calculate the radiative transfer.

( this real structure scene and the homogeneous assumption scene over the same growing period (same LAI) were compared, the corresponding canopy reflectance was simulated using SAIL, RGM and DART and analyzed for comparison.

3D maize scene modelling

The maize growth equation [START_REF] López-Lozano | Sensitivity of gap fraction to maize architectural characteristics based on 4D model simulations[END_REF] uses a simple geometric description of plant structure and a small number of canopy variables. The equation uses maximum leaf area per plant, the maximum number of leaves, plant density, row spacing, maximum plant height and thermal time-driven plant growth stage as input parameters to simulate the maize canopy. The model can be used in various situations without having extensive a priori knowledge of the geometric structural parameters of the canopy. The time variables are based on the maize canopy leaf sequence N proposed in the literature [START_REF] Grant | Simulation of maize phenology[END_REF]. N is calculated based on two variables: cumulative degree day (Td) since planting and thermal spacing (Td interval between the appearance of two consecutive leaves). Leaf stages can then be converted to calendar dates. Only the necessary information is provided here. For a detailed description, see the literature [START_REF] López-Lozano | Sensitivity of gap fraction to maize architectural characteristics based on 4D model simulations[END_REF][START_REF] Espana | Modeling maize canopy 3D architecture: application to reflectance simulation[END_REF].

For a leaf with order N, the area of the leaf SN is

d max 1 1 sin 2 a 2 N N S N π π           = + ⋅ -             •••••••••••••••••••••••••• (3.1)
where a and d are fixed parameters 0.59526 and 1.4158. Nmax is a user-specific input parameter that means the maximum number of leaves a single maize plant has. After that, the leaf area needs to be normalized to SN' ( )

max sum N N N S S S S ′ ⋅ = •••••••••••••••••••••••••••••••••••• (3.2)
where Smax is a user specifies input parameter that means the maximum total leaf area produced by a plant. The basal width of the leaf WN is width,length

N N W S R ′ = ⋅ ••••••••••••••••••••••••••••••••• (3.3)
where Rwidth,length is the leaf base to leaf length ratio and is a user-defined input parameter (here, we set it as 0.35). So, the leaf length LN is width,length

N N S L R ′ = •••••••••••••••••••••••••••••••••••• (3.4)
Some maize leaves may break. For leaf sequence numbers less than 12, there is a 15% probability of breaking; for leaf sequence numbers greater than or equal to 12, there is a 25% probability of breaking. If a leaf breaks, the breakage is randomly located between 0.25 and 0.95 times the LN, and the probability distribution follows a uniform distribution [START_REF] Espana | Modeling maize canopy 3D architecture: application to reflectance simulation[END_REF].

The inclination angle of the leaf θN is ( )

max rand randn N N S θ θ θ θ ′′ = ⋅ ∆ + + •••••••••••••••••••••••••• (3.5)
where Δθ is the difference between the first leaf's inclination and the leaf's inclination with the largest leaf area (fixed at 20°). θmax is the maximum leaf inclination (fixed at 45°), and the randn(θrand) term introduces some randomness (-θrand / 2 < rand(θrand) < θrand / 2, θrand fixed at 5°). S '' N is the normalized difference between the largest and the smallest leaf, calculated as

ma max x N N N N S S S S ′ ′ ′′ ′ - = ••••••••••••••••••••••••••••••••••• (3.6)
The leaves grow alternately on both sides of the plant, and the distribution of azimuth angle has a Gaussian standard deviation of φrand (fixed at 5°) at random.

The vertical height HN of a leaf with leaf order N is [START_REF] Espana | Modeling maize canopy 3D architecture: application to reflectance simulation[END_REF] ( ) ( )

max 1 max 1 1 1 c N c H H H N H N - = ⋅ -+ - ••••••••••••••••••••••••••• (3.7)
where H1 is the height of the first leaf of the maize, which takes a fixed value of 0.015 m, and Hmax is the user-defined maximum plant height. The coefficient c is calculated as follows max max 3.7 0.1 0.36

c N S = -⋅ - ⋅ ••••••••••••••••••••••••••••• (3.8)
The width of the base stem (RN) at order N is calculated as a linear function:

max 1 1 max ( 1) 1 N N R R R R N N   - = + ⋅ -   -   •••••••••••••••••••••••••••• (3.9)
where R1 is the width of the stem in the first stage of the leaf, which is equal to 0.002, and RNmax is the maximum width of the stem proportional to Hmax, and is calculated as

max max 0.008 N R H = ⋅ ••••••••••••••••••••••••••••••••• (3.10)
The number of leaf sequences N present in the plant was determined based on the degree days required to grow the first leaf, the degree days required to grow the second leaf and the degree days the leaf could survive. According to Eq. (3.1) to Eq. (3.10), the relevant geometric parameters were calculated, and these geometric parameters were applied to a static 3D maize modelling model developed based on ELSYS.

ELSYS [START_REF] Goel | A High-Level Language for L-systems and Its Applications[END_REF] is described in Appendix. A. For maize models developed using ELSYS, the input parameters were geometric parameters, including the maximum height of the plant, the radius Angle (VAA), were input into the radiative transfer (RT) model to create remote sensing scenes to simulate canopy bi-direction reflectance factor (BRF). See Table 3.1 for the meaning of abbreviations.

ELSYS output files were stored in a specific format and can be used by RGM or Thermal Radiosity-Graphics based Model (TRGM) models. Here, the format was converted programmatically to a standard 3D structure file with the suffix ".obj" and applied to other radiative transfer models like DART.

The 3D maize structure was modelled using the ELSYS platform. FC witnessed a considerable difference during the same growth period even though homogeneous and real structure maize canopies shared the same true LAI. For example, at the growth stage (Td = 100 to Td = 1100), the absolute difference between homogeneous and real structure maize canopy ranged from 0 to 0.25. Then, during the stabilization stage (Td = 1150 to Td = 1600), the absolute difference ranged from 0.26 to 0.33. Finally, during the decay stage (Td = 1600

to Td = 1950), the absolute difference ranged from 0.27 to 0.33. Moreover, the FC of the homogeneous canopy was approximately 1.5 times greater than that of the real structure maize canopy for the same true LAI. There was a significant difference in FC between homogeneous and real structure maize canopies, especially for fully developed plants with high LAI (Figure 3.2.g). It implied that the homogeneity assumption leads to an overestimation of FC.

Multi-temporal simulation analysis of maize canopy reflectance

Three vegetation radiative transfer models were used: The Scattering by Arbitrarily Inclined Leaves (SAIL) model [START_REF] Verhoef | Light scattering by leaf layers with application to canopy reflectance modeling: The SAIL model[END_REF] simulated BRFs based on a homogeneous canopy, and the accelerated RGM model (for simplicity, we call it the RGM model after this text) simulated BRFs based on a 3D real structural canopy. BRFs based on homogeneous and real structural canopies were simulated using the DART model and compared with the SAIL and RGM simulations. The geometrical parameters and soil OPs of the scenes were taken from the literature [START_REF] Duthoit | Assessing the effects of the clumping phenomenon on BRDF of a maize crop based on 3D numerical scenes using DART model[END_REF]. Leaf OPs were obtained using reflectance and transmittance data of maize leaves from the Leaf OPs Experiment 93 (LOPEX93) database. The size of the voxels along the x, y and z axes was set to 0.1 m to ensure that the true LAI ("sum of the true LAI" / "number of voxel layers") of each voxel in the homogeneous scene was always less than 0.5 [START_REF] Zhen | Potentials and Limits of Vegetation Indices With BRDF Signatures for Soil-Noise Resistance and Estimation of Leaf Area Index[END_REF]. There were many definitions of LAI, such as "effective LAI" or "true LAI". Here, LAI refers to "true LAI", defined as "the area of green leaves per unit area on one side". The LAD was spherical in the homogeneous canopy simulations, a typical value for describing maize canopy [START_REF] Duthoit | Assessing the effects of the clumping phenomenon on BRDF of a maize crop based on 3D numerical scenes using DART model[END_REF]. At the same T d, the true LAI of a homogeneous canopy was the same as that of a real structure canopy. real structure) maize scene simulations in the red and NIR bands along with the three observation angles (hotspot, darkspot and nadir direction) that were considered to provide richer information on the canopy structure. There was a good agreement between SAIL (red line) and DART homogeneous canopy (blue line) for the three observation angles in the red band. However, the RGM (light green)

and DART heterogeneous canopy (dark green) were less consistent than the SAIL and DART homogeneous canopy in the red band, especially in the nadir direction. This inconsistency may be due to the different radiative transfer modules of the different models. The consistency in the NIR band is not as good as the red band, which may be due to the strong multiple scattering effects in the NIR band.

There was a significant similarity between homogeneous and heterogeneous canopies in the red band's hotspot observation direction (solid line). However, there were minor differences between homogeneous and heterogeneous canopies in the darkspot observation direction (dashed line) and the nadir observation direction (dotted line). The differences between homogeneous and heterogeneous canopies in the NIR band were minor in the hotspot observation direction (solid line) and slightly different in the darkspot observation direction (dashed line). However, there were significant differences between homogeneous and heterogeneous canopies in the NIR band's nadir observation direction (dotted line).

In the NIR band, the darkspot reflectance was greater than the nadir reflectance (Figure 3.3).

Thus, our results confirmed that the lowest reflectance does not always occur in the darkspot direction, especially in the NIR band, where multiple scattering was strong [START_REF] Zhen | Potentials and Limits of Vegetation Indices With BRDF Signatures for Soil-Noise Resistance and Estimation of Leaf Area Index[END_REF][START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF]. Moreover, in the NIR band, the nadir direction reflectance of homogeneous canopies was much greater than the nadir reflectance of real maize canopies, implying that the homogeneity theory assumption would lead to overestimating reflectance. Other studies of maize also showed that the homogeneity theory assumption could lead to an overestimation of NIR band reflectance in the nadir direction [START_REF] Espana | Modeling maize canopy 3D architecture: application to reflectance simulation[END_REF].

They attributed this phenomenon to the fact that multiple scattering was mainly from vegetation.

The vegetation fraction seen at the nadir was minimal; however, the homogeneity assumption was not well accounted it. Most VIs used the NIR band, and the nadir observation direction was essential in remote sensing studies. Therefore, the error caused by the homogeneous canopy assumption was not negligible for maize canopies.

Summary

The maize growth equation [START_REF] López-Lozano | Sensitivity of gap fraction to maize architectural characteristics based on 4D model simulations[END_REF] is introduced into a 3D static maize model developed using the ELSYS platform to achieve multi-temporal maize scene simulations and radiative transfer simulations. The reflectance of three vegetation canopy scenes with homogeneous assumptions, simplified heterogeneity and real structures are simulated with SAIL, DART and RGM radiative transfer models for comparisons and cross-validations.

(1) The reflectance of the darkspot is not necessarily the lowest in all observation directions.

In the red band, the reflectance of the darkspot observations is the lowest. However, the nadir observations have the lowest reflectance in the NIR band rather than the darkspot.

(2) For the simulation of BRFs with row crops (i.e., maize), the canopy homogeneity assumption leads to an overestimated nadir reflectance in the NIR band. It is because multiple scattering is mainly from vegetation, and the proportion of vegetation seen at the nadir is minimal, which is not well considered by the homogeneity assumption. demonstrating the potential of the method to invert optical properties in complex environments with the presence of mixed pixels.

The chapter is presented in the paper:

(1) Z. Zhen, J-P Gastellu-Etchegorry, S. Chen, T.

Inversion of leaf optical properties using simulated images

A virtual 3D city structure scene was constructed using the DART model, and a corresponding satellite remote sensing image was simulated. This image replaced the real satellite image to assess the accuracy of the inversion-derived leaf OPs.

The study area was a synthetic 128 m × 128 m 3D urban scene with five components: ground, water, buildings, shrubs and trees (Figure 4.1). Three 15 m water surfaces and nine buildings with three types of roofs were evenly distributed in the scene. The roofs were triangular prism on the left side of the scene, flat in the middle and frustum of a square pyramid on the right side, all with a slope equal to 20°. The length, width and height of the buildings were all 15 m. The length and width of the ground between the buildings were also 15 m. These values were typical of European cities [START_REF] Meister | Large-scale bidirectional reflectance model for urban areas[END_REF]. The scene contained two shrubs with a length and width of 15 m and a height of 1 m. The shrubs were modelled as the homogeneous canopy filled with leaves with a spherical LAD and a leaf area density equal to 0.5 m 2 /m 3 . The scene also contained four trees. These trees were created directly by DART to simplify the analysis of the inversion procedure. In addition, the leaves in the canopy and the shrubs were simulated with fractional triangles rather than turbid medium. The trunks of the trees created with DART were replaced with regular octahedrons, and the crowns had specific vertical and horizontal LAD profiles and leaf area voxel density. Here, all trees had the same geometric parameters: an elliptical crown with a height of 10 m and a diameter of 10 m at the mid-height, uniformly filled with leaves with a leaf density of 0.5 m 2 /m 3 . The urban components were spatially distributed to conform to conventional light conditions, and the shading and mixed pixels effect could be considered.

DART calibration is an iterative inversion method applied to satellite images to obtain OP maps with urban anthropogenic heat flux for each urban component. First, the DART calibration continuously corrects the input to invert the OP of the components by iteratively comparing simulated and satellite images. This two-dimensional component OP map allows one to compute anthropogenic radiation flux maps. A detailed description of the DART calibration can be found in the literature [START_REF] Landier | Modélisation 3D du bilan radiatif des milieux urbains par inversion d'images satellites en cartes de réflectance et de température des matériaux urbains[END_REF][START_REF] Landier | 3D modeling of radiative transfer and energy balance in urban canopies combined to remote sensing acquisitions[END_REF][START_REF] Landier | Calibration of urban canopies albedo and 3D shortwave radiative budget using remote-sensing data and the DART model[END_REF]. The inversion procedure consists of two major steps.

The first step is to provide a first-order OP map using a linear spectral mixture analysis (LSMA) method, where DART estimates the "amount" of urban components in a pixel by simulating reflectance images and reflectance images for each type of component in the city (called "component reflectance images"). The term "amount" indicates the spectral contribution of each component rather than its relative area. The inversion starts with a user-defined OP value, which can be selected from the DART database. The first-order OP map is calculated by assuming that all urban components of the same type share the same OP within a certain distance. The urban components with a small "amount" are ignored because their OPs are challenging to be estimated.

Here, a threshold of 1% is set. The estimation of the first-order reflectance relies on a linear mixed pixel decomposition

( ) ( ) ( ) ( ) DART , 0 1 1 DART 0 1 , , , , n N n n n x y x y x y x y ρ ρ ρ ρ = = ⋅ ∑ ••••••••••••••••••••••• (4.1)
where ρ 0 n(x,y) represents the OPs of component n at coordinate (x,y) for the initial iteration 0. ρ 1 n(x,y) is the OP solved for in the first iteration and is unknown here. ρ 0 DART,n(x,y) is the reflectance value of the component image simulated by the DART model using ρ 0 n(x,y). ρ 1 DART(x,y) is the reflectance value simulated by DART in the first iteration. Because we want the DART simulated ρ 1 DART(x,y) to be infinitely close to the satellite reflectance value at the first iteration, which leads to the following equation ( ) ( )

1 sat DART , , x y x y ρ ρ = ••••••••••••••••••••••••••••••• (4.2)
From the Eq., (4.1) only one equation can be written For a single pixel. However, if the number of components within a pixel is larger than one, this equation system is unsolvable. Therefore, it is necessary to assume that the same components n within adjacent pixels have the same ρ 0 n(x,y), increasing the number of equations to make the equation system solvable.

The second step is a series of iterations in which the final OPs maps of each type of urban component converge at the spatial resolution of the satellite images. It combines the dichotomous and Newtonian methods. The OPs are the independent variables, and the component reflectance images are the dependent variables. If the pixel reflectance of the DART simulated images is less than the reflectance of the satellite image, the OP values of the components in this pixel increase, and vice versa. The iterative process is carried out on a pixel-wise basis. If the relative difference between the reflectance of the DART and the satellite image is less than a user-defined threshold, the OP values of the components in the pixel are not updated in the following iterations; here, this threshold is set to 1%. This step implicitly considers the multiple scattering between components.

At the end of the procedure, the DART simulated image with the OP map obtained by inversion is almost equal to the satellite image (i.e., the mean relative error is less than 10 -3 ).

For DART simulated remote sensing images with one-pixel reflectance ρn k-1 (x,y) and ρn k (x,y), if the simulated pixel reflectance is lower than the remote sensing image reflectance in two consecutive iterations k-1 and k, the OPs of component n at the next iteration k+1 can be approximated using Newtonian iteration

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 DART, DART, 1 1 DART, DART, , , , , , , , ,
k k k k n n n n k k n n k k n n x y x y x y x y x y x y x y x y ρ ρ ρ ρ ρ ρ ρ ρ + - + -     - ⋅ -     = ± - • (4.3)
We want the DART simulated image's pixel reflectance to be very close to the satellite image in the next iteration. Therefore, we have

( ) DART 1 sat ( , ) , k x y x y ρ ρ + = •••••••••••••••••••••••••••••••• (4.4)
Combining Eq.(4.3) and Eq.(4.4), we have

( ) ( ) ( ) ( ) ( ) 1 1 sat DART, 1 DART, DART, DART , , , 1 ,
,

k k n n k k k n n n k k k n n x y x y x y x y x y ρ ρ ρ ρ ρ ρ ρ ρ ρ - + - -   =± ⋅ ⋅ -   -   ••••••• (4.5)
For DART simulated remote sensing images with one-pixel reflectance ρn k-1 (x,y) and ρn k (x,y), if in two consecutive iterations (k-1 and k), the simulated image reflectance is one lower and one higher than the satellite image, the OPs of the component n at the next iteration (k+1) can be approximated by a weighted dichotomous method ( ) ( ) ( )

1 1 , , , k k k k n n n n x y x y x y ρ ρ α ρ ρ + - = ± - ••••••••••••••••••••• (4.6)
where α is the weight and is calculated as

1 DART, DART, 1 DART, DART, k k n n k k n n ρ ρ α ρ ρ + - - = - •••••••••••••••••••••••••••••••• (4.7)
The accuracy of the DART calibration method was assessed using the DART-simulated pseudo-satellite BOA reflectance image. DART, therefore, calculated two products: the pseudosatellite image and the so-called initial image, which converged to the pseudo-satellite image as the calibration method iteratively changed the component OPs image. The relative difference between the true and inversion-derived OP values was used to indicate the accuracy of the DART calibration.

Six spectral bands were considered (three visible bands, NIR and two shortwave infrared bands). It was important to note that the DART calibration handled all bands independently; it could even handle single-band images. In addition, the DART atmosphere radiative transfer module [START_REF] Gascon | Radiative transfer model for simulating high-resolution satellite images[END_REF][START_REF] Wang | DART: Improvement of thermal infrared radiative transfer modelling for simulating top of atmosphere radiance[END_REF] was not used because the satellite images were considered to be atmosphere corrected.

Noiseless ideal case for leaf optical properties inversion

As mentioned previously, two base simulation cases were investigated: the noise-free ideal case and the noise-interference case. In the noise-free ideal case, the OP values were the only unknowns (Table 4.1). ρ and τ referred to reflectance and transmittance, respectively. The OP of the tree was equal to reflectance + transmittance, and SWIR referred to the shortwave infrared band. All other input parameters were the same in the base and satellite simulations. SZA equal 30°, a perfect sensor Modulation Transfer Function (MTF), a 4 m spatial resolution of the satellite images, a perfect geometric calibration between the initial image and the satellite images, etc. (Table 4.2). The SZA greatly influenced the bidirectional reflectance effect [START_REF] Gastellu-Etchegorry | Modeling BRF and radiation regime of boreal and tropical forests: I. BRF[END_REF][START_REF] Guillevic | Modeling BRF and radiation regime of boreal and tropical forest: II. PAR regime[END_REF][START_REF] Li | Analysis of Accuracy of MODIS BRDF Product (MCD43 C6) Based on MISR Land Surface BRF Product-A Case Study of The Central Part of Northeast Asia[END_REF] and the distribution of shadows in the images. Furthermore, together with the roof slope, they determined the angle of incidence of light striking the roof. MTF defined how the sensor retained much contrast for the original target. It defined the genuine transition of the spatial frequency content of the object in the image. It could cause neighbouring pixels to interfere with each other's spectra, thus affecting the accuracy of OP inversion [START_REF] Xu | Simulation of the Impact of a Sensor's PSF on Mixed Pixel Decomposition: 1. Nonuniformity Effect[END_REF]. The spatial resolution referred to the size of a pixel on the ground. |ε| ��� mixed tree =0.3806. In addition, the accuracy of the OP was very similar in sunlight and shade, and the initial confusion of the OP spectra had a negligible effect on the results ε̅ component . Ground was excellent, even with shaded ground, except for the mixed pixels with vegetation. The results for roofs and water surfaces were excellent, even in the shaded areas, except for the mixed water surface pixels containing vegetation. The shape of the roof had a negligible effect on the results. However, the accuracy of vegetation, which included shrubs and trees, was relatively unsatisfactory compared to other components.

Artificially noise case for leaf optical properties inversion

The additional noise experiment considered the complexity of the satellite images. Several and the 3D geometric accuracy of the urban scene.

In the noise interference case, the SZA and spatial resolution of satellite images were the same in pseudo-satellite and initial images, as they were usually known. In contrast, the MTF, the The first-order and total-order sensitivity indices were calculated using the Fourier decomposition term of the model output [START_REF] Tarantola | SimLab 2.2 Reference Manual[END_REF]. The first-order sensitivity index is a direct measure of sensitivity based on variance and contributes to the output variance of the main effects of the input parameters. As such, it measures the effect of variation in only one input parameter but averages over the variation in the other input parameters. The total order sensitivity index measures the contribution of each input parameter to the output variance, including all the variance resulting from its interaction with other input variables in any order.

A sensitivity study was carried out in the presence of noisy interference. MTF: MTF had a minimal effect on calibration accuracy. The highest mean first-order sensitivity (0.095) and total order sensitivity (0.658) of all bands were for shrubs. The lowest mean first-order sensitivity (0.018) and total order sensitivity (0.093) of all bands were for the ground.

The mean first-order and total-order sensitivity of the MTF over all bands were 0.056 and 0.346.

SZA: The highest mean first-order sensitivity (0.721) and mean total order sensitivity (0.861)

were ground in all bands. Water had the lowest mean first-order sensitivity (0.047) and mean total order sensitivity (0.348) over all bands. Overall, the average first-order and total-order sensitivities for all bands of the SZA were 0.253 and 0.660, respectively.

Satellite resolution: water had the highest mean first-order sensitivity (0.327), and shrubs had the highest mean total-order sensitivity (0.891) over all bands. The ground had the lowest mean first-order sensitivity (0.036) and total-order sensitivity (0.137) over all bands. The mean first-order and total-order sensitivities for spatial resolution over all bands were 0.222 and 0.628.

X-axis shifting: the highest mean first-order sensitivity (0.145) and total order sensitivity (0.806) over all bands was for shrubs. The lowest mean first-order sensitivity (0.033) and total-order sensitivity (0.129) over all bands was ground. The mean first-order and total-order sensitivity of the

x-shift over all bands for all components were 0.096 and 0.441, respectively.

Y-axis shifting: the roof had the highest mean first-order sensitivity (0.105) over all bands, and the shrub had the highest mean total-order sensitivity (0.727) over all bands. The ground had the lowest mean first-order sensitivity (0.050) and total-order sensitivity (0.172) over all bands. The mean y-shift first-order and total-order sensitivity for all components were 0.085 and 0.437 over all bands.

Geometric accuracy of the 3D city database: the highest mean first-order sensitivity (0.094) and total-order sensitivity (0.758) over all bands was for shrubs. Roofs had the lowest mean firstorder sensitivity (0.015) over all bands, while the ground had the lowest mean total-order sensitivity (0.088) over all bands. The mean first-order and total-order sensitivities of the 3D scene model accuracy over all bands over all components were 0.054 and 0.407.

In general, the mean sensitivity values for all components of the six bands were close to each other: 0.317 for the blue band, 0.305 for the green band, 0.307 for the red band, 0.264 for the NIR, 0.301 for the shortwave infrared 1 (SWIR1) band, and 0.342 for the shortwave infrared 2 (SWIR2) band. In addition, the mean total order sensitivity (0.485) was four times greater than the mean firstorder sensitivity (0.127) over all bands and all components.

In decreasing order, the DART calibration accuracy was affected by SZA, satellite spatial resolution, pixel shifting, 3D city database inaccuracy, and MTF. SZA significantly affected the apparent ground reflectance through shadowing and bi-directional reflectance effects. In addition, the SZA and the roof angle determined the angle of incidence for roofs. The low sensitivity of the SZA to water maybe because the water was treated as a surface without topography here.

The spatial resolution of satellite images was essential in urban studies because it determined the distribution of pure pixels in satellite images [START_REF] Franke | Hierarchical multiple endmember spectral mixture analysis (MESMA) of hyperspectral imagery for urban environments[END_REF]. Welch et al. [START_REF] Welch | Spatial resolution requirements for urban studies[END_REF] used averaged city maps to demonstrate that high-frequency detail of urban scene components was necessary for remote sensing data with a spatial resolution of 0.5-10 m. A spatial resolution of at least 5 m was required to capture applications of urban structure adequately [START_REF] Small | High spatial resolution spectral mixture analysis of urban reflectance[END_REF]. The spatial resolution was also a constraint on the DART calibration, as it assumed that the OPs of the components were constant within that pixel. Subclassification of urban components could be a potential solution. For example, roofs in the same pixel could be classified as sub-components if additional information (e.g., type of material) was available. Furthermore, this assumption could also be mitigated by using higherresolution remote sensing imagery.

The pixel shifting and the accuracy of the 3D model had little impact on the inversion accuracy of continuous and homogeneous components such as the ground, as they had a similar geometric impact on the problem of estimating the "amount" for each component inside a mixed pixel. MTF had less impact on large and homogeneous scene components such as the ground than on small components such as roofs, trees, water and shrubs. In summary, small and highly heterogeneous components were much more sensitive than large and highly homogeneous ones. The DART calibration was validated only by comparing the time series of the DART-simulated shortwave Radiative budget with the time series measured in the flux tower in Basel, Switzerland [START_REF] Landier | Calibration of urban canopies albedo and 3D shortwave radiative budget using remote-sensing data and the DART model[END_REF]. Over one year, the mean relative difference was 2.7%, which was a very encouraging result.

However, this comparison was only a single-pixel validation. Here, the validation was extended to each pixel by using pseudo-satellite images from the DART simulation.

The inversion results (Figure 4. Other work has highlighted the poor performance of vegetation in unmixing models [START_REF] Franke | Hierarchical multiple endmember spectral mixture analysis (MESMA) of hyperspectral imagery for urban environments[END_REF],

particularly for mixed pixels containing vegetation. The complexity of vegetation structure, including multiple scattering, leads to a complex non-linear relationship between leaf albedo and canopy reflectance [START_REF] Zhen | Potentials and Limits of Vegetation Indices With BRDF Signatures for Soil-Noise Resistance and Estimation of Leaf Area Index[END_REF][START_REF] Malenovský | Variability and uncertainty challenges in scaling imaging spectroscopy retrievals and validations from leaves up to vegetation canopies[END_REF]. Consistent with these conclusions, DART calibration provides a less accurate map of vegetation's OPs than urban surface components, such as ground and roof. Similarly, the low reflectance of the water complicates the inversion of its OPs [START_REF] Franke | Hierarchical multiple endmember spectral mixture analysis (MESMA) of hyperspectral imagery for urban environments[END_REF]. Here, however, the results for water are excellent, even in the shaded areas. It is because the water is treated as an opaque surface in this work (no transmission is considered).

Shadows are often considered a significant factor in the inaccuracy of pixel unmixing [START_REF] Meister | Large scale multispectral BRDF of an urban area[END_REF].

For example, shadows cast by buildings in dense urban areas can significantly obscure local information in an image, leading to results that may be corrupted or misinterpreted [START_REF] Dare | Shadow analysis in high-resolution satellite imagery of urban areas[END_REF]. Methods like Multiple Endmember Spectral Mixture Analysis (MESMA) solve this problem by treating shadows as one component [START_REF] Franke | Hierarchical multiple endmember spectral mixture analysis (MESMA) of hyperspectral imagery for urban environments[END_REF]. Here, the small relative error of shadows on the inversion of OPs suggests that the DART calibration algorithm is well suited for dealing with urban areas with shadows.

The presence of noisy interference situations emphasises that inaccurate DART input parameters can seriously affect calibration accuracy. For example, inaccurate geometric calibrations of satellite images lead to inaccurate estimates of the "amount" and, finally, inaccurate OPs.

Compared to some traditional unmixing methods, the DART calibration does not estimate the abundance calculated from the OPs of the endmembers, which is a significant advantage as it allows the OPs of the endmembers to vary between pixels. However, it requires a very precise coregistration of the satellite image and the 3D city database. In addition, some components with low OPs (e.g., water and vegetation) produce high relative errors and variance. It is because the error is assessed using a relative error. If the pixel shifting problem results in a high reflectance pixel to a low optical property component, the denominator in the relative error is too small, resulting in high relative error and variance.

Inversion of leaf optical properties using satellite image 4.2.1 Study area overview and PlanetScope data pre-processing

The DART calibration was applied to real satellite data to validate its accuracy. The study area was chosen as Basel, one of the research targets of the Urban Anthropogenic Heat FLUX from Earth Observation Satellites (URBANFLUXES) project provided by the Earth Observation Satellite [START_REF] Chrysoulakis | Urban energy exchanges monitoring from space[END_REF].

Basel is on the Rhine River in north-western Switzerland, where the Swiss, French and German borders meet [START_REF] Landier | Modélisation 3D du bilan radiatif des milieux urbains par inversion d'images satellites en cartes de réflectance et de température des matériaux urbains[END_REF][START_REF] Landier | Calibration of urban canopies albedo and 3D shortwave radiative budget using remote-sensing data and the DART model[END_REF]. The study area encompassed almost the entire city of Basel, with latitudes ranging from 47.53° N to 47.57° N and longitudes from 7.55° E to 7.62° E. Its geometric database was very detailed (Figure 4.7). Basel is a representative city in Europe, so the results can be extended to other European cities. Basel's 3D realistic remote sensing scenes were constructed and used by Lucas et al. [START_REF] Landier | Modélisation 3D du bilan radiatif des milieux urbains par inversion d'images satellites en cartes de réflectance et de température des matériaux urbains[END_REF][START_REF] Landier | 3D modeling of radiative transfer and energy balance in urban canopies combined to remote sensing acquisitions[END_REF][START_REF] Landier | Calibration of urban canopies albedo and 3D shortwave radiative budget using remote-sensing data and the DART model[END_REF]. According to the URBANFLUXES project [START_REF] Feigenwinter | Spatial distribution of sensible and latent heat flux in the city of Basel (Switzerland)[END_REF][START_REF] Feigenwinter | Spatial distribution of sensible and latent heat flux in the URBANFLUXES case study city Basel (Switzerland)[END_REF], land cover categories in Basel included buildings, impervious cover, water, low vegetation, bare soil, deciduous and evergreen forest. In our study, some categories were regrouped. The tree category contained deciduous and evergreen; the ground category contained impervious layer, low vegetation and bare soil; and the building and water categories were the same as that in the URBANFLUXES project. Sun-Synchronous Orbit (SSO) solar synchrotron orbit is 3.7 m. More details can be found in the document [START_REF] Marta | Planet Imagery Product Specifications[END_REF]. A 4-band multispectral surface reflectance image with 3 m resolution was downloaded from the Planet website (https://www.planet.com/). This image was orthorectified and atmosphere corrected, and it was selected because it had an off-nadir angle of exactly 0.0°. Table 4.3 shows the detailed parameters of this image. The data processing consisted of a co-registration between the satellite image and the DART simulated image. The satellite images were first clipped according to geographic latitude and longitude, retaining only the remote sensing images within the study area. The 3D urban database was then shifted along the x and y axes to warp to the satellite images. The geometric accuracy of the DART simulated image was much higher than that of the satellite image due to observation angle and image distortion factors. Therefore, the DART simulated image was used as the base image, and the satellite image was warped to the DART simulated image for co-registration using GeFolki [START_REF] Plyer | A new coregistration algorithm for recent applications on urban SAR images[END_REF][START_REF] Brigot | Adaptation and evaluation of an optical flow method applied to coregistration of forest remote sensing images[END_REF], developed by the French Aerospace Laboratory. Gefolki is a module that performs coregistration of two remote sensing images with a geometry error of less than 0.1 pixels. Finally, the warped satellite image was used for inversion with the DART simulated image. The Subarctic Summer (SAS) was chosen as the atmospheric model in the DART simulations.

This atmospheric model was selected based on the total water vapour column and atmospheric temperature [START_REF] Cooley | FLAASH, a MODTRAN4-based atmospheric correction algorithm, its application and validation[END_REF]. These data were derived from water vapour and surface temperature data from the National Centers for Environmental Prediction and National Center for Atmospheric Research (NCEP/NCAR) reanalysis data [START_REF] Kalnay | The NCEP/NCAR 40-year reanalysis project[END_REF]. The aerosol model was chosen directly from the study area as the urban model [START_REF] Cooley | FLAASH, a MODTRAN4-based atmospheric correction algorithm, its application and validation[END_REF]. Atmosphere parameters had a minimal influence on the results (only the skylight ratio distribution was affected) as atmosphere corrected surface reflectance data were used, assuming atmosphere effects were removed. The input parameters of the spectral band, spatial resolution and solar-target-satellite geometry remained consistent with the satellite images (Table 4.3) in the DART simulations. the OP of some bridges was assigned to the water component. In addition, some ground pixels showed typical vegetation spectral features (Figure 4.9). It is because there was no distinction between ground and grass in the urban 3D database. As a result, the grass pixels were also considered an opaque surface during the DART radiative transfer simulations. Overall, there were apparent spectral differences between the components, indicating good spectral unmixing. It should be noted that a satellite image with an off-nadir angle of 0.0° was used in this study.

Inversion of optical properties of leaves from multispectral images

The major difficulty in using off-nadir images was the inconsistent coordinate system between the satellite images and the DART simulated images. Most surface reflectance products (e.g., Sentinel2

L1C products dated before March 2021, https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2msi/definitions) used DEM to do orthorectification. This process was known as "orthorectification".

As a result, some buildings were skewed in the image (the walls are visible in the image). However, the DART model used the Digital Surface Model (DSM) for orthorectification, known as "true orthorectification ". As a result, all buildings were vertical (there are no wall elements in the simulated image). The coordinates of the corresponding image points in the two images were different. However, geometric accuracy was crucial for DART calibration, especially for scenes with strong heterogeneity. The DART model is currently being investigated to provide an "orthorectification" image in addition to the "true orthorectification" image.

Using negative soil adjustment factor of SAVI to mitigate vegetation index saturation effect

According to the original SAVI literature [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF], the optimum soil adjustment factor should be close or equal to zero under dense vegetation area and close or equal to one under sparse vegetation area. However, Ren et al. [START_REF] Ren | Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands[END_REF] find a good performance of negative soil adjustment factor using field-measured data in arid grasslands where LAI is very low. In order to study the different results of optimal adjustment factors, we deduce the geometric interpretation of vegetation indices. Then, vegetation isoline behaviours are studied based on the DART simulated data. Finally, the assumption that using negative soil adjustment factor to mitigate vegetation index saturation effect is proposed and validated.

Study of index isoline and vegetation isoline

The variations of the VI value of constant canopy optical properties induced by canopy background may be simply expressed in the red-near infrared (NIR) reflectance space as the difference between vegetation isolines and index isolines [START_REF] Yoshioka | Analysis of vegetation isolines in red-NIR reflectance space[END_REF].

(1) Index isoline Vegetation indices are usually designed based on the isoline concept. Every VI has its isoline, and usually, their isolines cannot fully overlap on the vegetation isolines. The ideal VI should be designed with the maximum overlap between isoline and vegetation isoline. Here we take RVI, NDVI, SAVI, PVI and TSAVI as an example to illustrate their isolines, and their corresponding geometry meanings are shown in Figure 5.1.

(i) RVI

RVI is the slope of isoline RVI tan( )

α = •••••••••••••••••••••••••••••••••••• (5.1) (ii)
1 π NDVI tan( ) 4 1 R R R R R R R R R R R R R R α - - - = + = = = - + + ••••••••••••• (5.2)
As a result, NDVI is a tan value of α-π/4, and it is a positive correlation with α when -π/2 < απ/4 < π/2. The index isolines of NDVI is depicted in Figure 5.2.a.

(iii) SAVI

Based on Eq.(1.3), we can get ( )

nir red SAVI SAVI SAVI nir red (1 ) SAVI= 2 2 R R X X X R R - ⋅ + + + + ••••••••••••••••••••••••• (5.
- - - - - -( 1 ) ( a R b a R a R b a R R a b X a X X X X X X a a R a R a R b a R b a a R R X b X a R a a ⋅ + ⋅ ⋅ ⋅ + + ⋅ -⋅ - = ⋅ + -⋅ + ⋅ + = ⋅ + ⋅ + ⋅ - -- + ⋅ - ⋅ ⋅ -- = + •••••••••••••• (5.4) 1 ( ) ( ) ( ) ) -
Because XTSAVI is the negative abscissa of soil adjustment factor, by assuming the slope of vegetation isoline as kv and based on Eq.(5.4), we get

v v TSAVI= 1 tan( + ) tan( ) k a a k a a a β γ γ β - ⋅ + ⋅ = ⋅ - = ⋅ ••••••••••••••••••••••••••••••• (5.5)
It should be noted that when β = 90°, TSAVI equal ∞; When 90° < β < 180°, TSAVI is negative. (2) Vegetation isoline Vegetation isolines are a series of canopy reflectance points with the same optical and structural properties (e.g. constant leaf reflectance, leaf transmittance, LAI and green cover) but with different background brightness [START_REF] Yoshioka | Analysis of vegetation isolines in red-NIR reflectance space[END_REF].

Most of the two-band VIs were designed based on the geometry of the convergence points of all vegetation isolines in the red-NIR plane or the intersection of vegetation isolines with soil lines.

Understanding vegetation isolines behaviour is essential because the desired VI should yield a constant value for the reflectance points on the same vegetation isolines where vegetation conditions are considered constant. Therefore, it is essential to study the behaviours of vegetation isolines.

BRFs of homogeneous and heterogeneous canopies simulated by the DART model was used to investigate the vegetation isolines behaviours. The DART model simulated the reflectance of homogeneous and heterogeneous vegetation canopy scenes and compared the differences between different levels of simplification of the remote sensing scene. This model can simulate the energy radiation balance and remote sensing images of an Earth scene with an atmosphere and simulate both homogeneous and non-homogeneous canopy scenes [START_REF] Gastellu-Etchegorry | Modeling radiative transfer in heterogeneous 3-D vegetation canopies[END_REF][START_REF] Gascon | Retrieval of forest biophysical variables by inverting a 3-D radiative transfer model and using high and very high resolution imagery[END_REF][START_REF] Gastellu-Etchegorry | Dart: Radiative Transfer modeling for simulating terrain, airborne and satellite spectroradiometer and LIDAR acquisitions and 3D radiative budget of natural and urban landscapes[END_REF][START_REF] Yin | DAta simulation and fusion of imaging spectrometer and LiDAR multi-sensor system through dart model[END_REF][START_REF] Gascon | Using multi-directional high-resolution imagery from POLDER sensor to retrieve leaf area index[END_REF][START_REF] Gastellu-Etchegorry | Discrete anisotropic radiative transfer (DART 5) for modeling airborne and satellite spectroradiometer and LIDAR acquisitions of natural and urban landscapes[END_REF][START_REF] Grau | Estimation of 3D vegetation density with Terrestrial Laser Scanning data using voxels. A sensitivity analysis of influencing parameters[END_REF][START_REF] Fawcett | Advancing retrievals of surface reflectance and vegetation indices over forest ecosystems by combining imaging spectroscopy, digital object models, and 3D canopy modelling[END_REF]. A variety of LAI, ALA, and canopy bi-directional reflectance factors for soil OPs was simulated using DART.

Most of the input parameters (Table 5.1) were taken from the literature [START_REF] Baret | Potentials and limits of vegetation indices for LAI and APAR assessment[END_REF], where a well-known soil adjusted family VI was proposed, called TSAVI. It was designed to improve SAVI. Several parameters were modified to fit the study objectives better: two SZAs of 30° and 60° were added to examine the effects of SZA. Hotspot directions and darkspot directions were also considered in calculating hotspot signature VIs. As the maximum LAI field measurement was 6.3, the maximum LAI value (12.8) was removed to keep the LAI range between the simulated and measured data not too large for the LAI inversion. The OPs of leaves and soils were referenced from the literature [START_REF] Baret | Potentials and limits of vegetation indices for LAI and APAR assessment[END_REF]. and horizontal layers (i.e., a random distribution of elements in infinitely small planes). Therefore, the components of the scene included soil and homogeneous vegetation as the substratum. The voxel sizes were set to ensure that the true LAI ("sum of the true LAI" / "number of voxel layers") of each voxel in the scene was always less than 0.5, even at the maximum LAI (6.4), to ensure an accurate simulation of the multiple scattering within the voxel [START_REF] Zhen | Potentials and Limits of Vegetation Indices With BRDF Signatures for Soil-Noise Resistance and Estimation of Leaf Area Index[END_REF]. In the heterogeneous scene (Figure 5.3.b) simulated by DART (Table 5.2), trees consisted of a trunk and a crown filled by an isosceles triangle (twice the height of the bottom edge). All trees with random dimensions were randomly positioned in the scene using the DART "random tree positioning" model [START_REF] Gastelluetchegorry | Modeling specular reflectance and polarization in DART model for simulating remote sensing images of natural and urban landscapes[END_REF][START_REF] Gastelluetchegorry | Time series of urban radiative budget maps derived from EO satellites using a physical remote sensing model[END_REF][START_REF] Gastelluetchegorry | Dart: A Tool For Studying Earth Surfaces -Time Series of Urban Radiative Budget From Eo Satellites[END_REF][START_REF] Gastelluetchegorry | Simulation of Chlorophyll Fluorescence for Sun-and Shade-Adapted Leaves of 3D Canopies with the Dart Model[END_REF][START_REF] Grau | Estimation of 3D vegetation density with Terrestrial Laser Scanning data using voxels. A sensitivity analysis of influencing parameters[END_REF]: trees were randomly set in a disk (radius r) centred on a Δx / Δy grid node on the scene, with a probability of existence.

In this model, the geometric properties of the tree were defined as the mean and standard deviation.

The standard deviation was set to 20% of its value for each parameter. Each trunk consisted of 8 trapezoids (horizontal segments = octagonal). Each crown was ellipsoidal, and the interior was randomly filled with isosceles triangles. The scene included only vegetation and soil. In the DART model, tree LAI was defined as the ratio of "total leaf area of all trees of the selected species" / "scene area". In this experiment, only one tree species was used.

Table 5.2 Additional input parameters for the heterogeneous scenes.

Parameters Value

Scenes

Number of trees 72

Probability of tree presence at nodes 0.90 Due to the significant variation in LAI, the range of variation in the reflectance of the simulated canopy in red and NIR was very large (Figure 5.4). However, when the LAI increased above 3, the reflectance no longer varied dramatically, indicating saturation of the VIs. Reflectance also varied with the observation direction due to anisotropic scattering by vegetation. The dynamic range of the NIR reflectance was significantly greater than that of the red reflectance, as described in the literature [START_REF] Huete | A comparison of vegetation indices over a global set of TM images for EOS-MODIS[END_REF]. Besides, hotspot reflectance was highest in the PP, while darkspot reflectance was not always the lowest, as described in the literature [START_REF] Chen | Predicting leaf area index in wheat using an improved empirical model[END_REF][START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF]. In addition, the highest values of NIR band reflectance were smaller in the heterogeneous canopy than in the homogeneous canopy because soils with low NIR occupied a large part of the scene. Also, as the SZA increased, the dynamic range of red band reflectance was much smaller than that of the NIR band, which may increase VIs values, as described in the literature [START_REF] Walter-Shea | Relations between directional spectral vegetation indices and leaf area and absorbed radiation in alfalfa[END_REF].

X
Here, vegetation condition was summarised in terms of "FC" (i.e., the proportion of the scene covered by foliage) as a variable to characterise vegetation isolines [START_REF] Rondeaux | Optimization of soil-adjusted vegetation indices[END_REF]. As the FC of turbid media cannot be calculated using its definition, LAI and ALA [START_REF] Campbell | Extinction coefficients for radiation in plant canopies calculated using an ellipsoidal inclination angle distribution[END_REF] were used to calculate it, as used in the literature [START_REF] Rondeaux | Optimization of soil-adjusted vegetation indices[END_REF]. Figure 5.5 shows the vegetation isolines figures for the direction of nadir observation with SZA equal to 45° and ALA equal to 70°. why VIs such as PVI, which assume parallel vegetation isolines, perform better in areas with low FC. In contrast, VIs such as SAVI, which assume a convergence of vegetation isolines to a common point, perform better in areas with high FC. Furthermore, there was a clear tendency for the intersections of the vegetation isolines with the soil lines to shift to the right, and these intersections can even move from the second or third quadrant towards the first quadrant. Thus, the results further confirmed the conclusion from Ren et al. [START_REF] Ren | Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands[END_REF] that the intersections of vegetation isolines with the soil line can be located in the first quadrant. Previous conclusions suggested that the intersections of vegetation isolines moved close to the origin as FC increased [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF][START_REF] Qi | A modified soil adjusted vegetation index[END_REF][START_REF] Gilabert | A generalized soil-adjusted vegetation index[END_REF]. In contrast, Ren et al. [START_REF] Ren | Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands[END_REF] noted that the intersection of vegetation isolines moved away from the origin as FC increased, thus contradicting the previous findings [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF][START_REF] Qi | A modified soil adjusted vegetation index[END_REF][START_REF] Gilabert | A generalized soil-adjusted vegetation index[END_REF]. This contradiction was that vegetation isolines in previous studies intersected the soil line in the second or third quadrant, whereas in Ren et al. [START_REF] Ren | Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands[END_REF] case, they intersected in the first quadrant.

Therefore, here we proposed to describe the intersections movement using the abscissa rather than the absolute distance from the origin and to propose a phenomenon of the right shift of the intersections, which was not only consistent with Ren et al. [START_REF] Ren | Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands[END_REF] but also with the previous findings [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF][START_REF] Qi | A modified soil adjusted vegetation index[END_REF][START_REF] Gilabert | A generalized soil-adjusted vegetation index[END_REF].

It is not hard to prove that a decrease in XSAVI implies an increase in SAVI. Thus, it must be shown that y= )( )

1 2 1 2 1 2 ( 1) 1 0 1 a b x x y y a b x x + - - - = > + + + ••••••••••••••••••••••••••••• (5.6)
It implied: y1> y2. Thus, a decrease in XSAVI leads to an increase in the SAVI value.

A decrease in XSAVI led to a right shift in the intersection of the SAVI isoline with the soil line [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF] and an increase in the slope of the SAVI isoline, which was why the optimal XSAVI should decrease as LAI increases. It was worth noting that this conclusion applied both when XSAVI was positive and negative. Therefore, there was no significant difference between positive and negative values of XSAVI. The sign of the intersection abscissa depended strongly on the soil line's intercept and the vegetation isolines' intercept. It was negative when the intercept of the vegetation isoline was greater than that of the soil line. When the intercept of the vegetation isoline was smaller than that of the soil line, the abscissa of the intersection was positive. If the intersection of the vegetation isoline lines with the lowest LAI had a positive transverse coordinate, then most of the intersection points had a positive transverse coordinate.

Compared to the homogeneous canopy, the dynamic slope range in the heterogeneous canopy was much less with increasing LAI (Figure 5.5.b). In addition, there was no significant right shift in the position of the intersection (Figure 5.6.b). This conclusion contradicted the results of homogeneous canopy studies [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF][START_REF] Qi | A modified soil adjusted vegetation index[END_REF][START_REF] Gilabert | A generalized soil-adjusted vegetation index[END_REF][START_REF] Ren | Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands[END_REF] and previous studies, suggesting that the conclusions drawn from homogeneous canopies did not necessarily apply to heterogeneous canopies.

Relationship between vegetation indices and LAI (1) The empirical relationship between vegetation index and leaf area index

For convenience, a linear relationship is normally used to descript the relationship between VIs and LAI. That is to say, the maximum value of VIs corresponds to the maximum value of LAI, and LA1 can be inferred from standardized VI values min max max min

VI VI LAI LAI

VI VI

- = + •••••••••••••••••••••••••••• (5.7)
However, the relationship between VI and LAI is not always linear. For example, the VI always tends to be asymptotically saturated within the range of LAI from 2 to 6. It also depends on vegetation cover type and environmental conditions. When the linear relationship is no longer satisfied, the saturation effect occurs. Therefore, the linear relationship is a major focus in our antisaturation study.

The semi-empirical relationship for saturation case between VIs and LAI can be expressed in terms of a modified Beer's law to descript the non-linear relationship [START_REF] Baret | Potentials and limits of vegetation indices for LAI and APAR assessment[END_REF] ( ) ( )

max min max VI VI=VI + VI VI exp LAI K - ⋅ - ⋅ ⋅ Ω •••••••••••••••••••• (5.8)
where KVI is the coefficient controlling the slope (equivalent to the extinction coefficient), VImax is the maximum value of the VI (equivalent to the asymptotic value of VI as LAI tends to infinity), and VImin is the minimum value of VI (equivalent to the VI of bare soil).

(

2) Time series decomposition analysis of remote sensing LAI data

Time series analysis is a statistical method for predicting, controlling and understanding characteristics of time series data [START_REF] Hamilton | Time series analysis[END_REF]. The ability of the VIs to invert LAI in areas of high FC was evaluated. Seasonal decomposition of the time-series leaf area index (LAI) was first performed to test the data stability. The VIs were then calculated for the period corresponding to the LAI to estimate trends in LAI variation.

The study area (Figure 5.7) comprised the Apiacás region (latitude: -7.35° ~ -9.82°N, longitude:

57.04° ~ 58.57°W) located in the northernmost intact part of Mato Grosso, Brazil. There were five land cover types in the study area: water body (blue), evergreen broadleaf vegetation (dark green), deciduous broadleaf vegetation (yellow), annual grass vegetation (yellow-green), and urban and built-up land (dark grey). The area was chosen for its dense vegetation. It contains part of the 19,582 square kilometres of the Juruena National Park, one of the largest conservation units in Brazil. Its elevation is around 200 m [START_REF] Conservação | Unidade de Conservação: Parque Nacional do Juruena[END_REF].

The multispectral satellite data and MODIS LAI data were provided and processed by GEE [START_REF] Mutanga | Google Earth Engine Applications[END_REF][START_REF] Kumar | Google Earth Engine applications since inception: Usage, trends, and potential[END_REF][START_REF] Gorelick | Google Earth Engine: Planetary-scale geospatial analysis for everyone[END_REF] (https://earthengine.google.com/). They were used to validate the LAI of inversions of LAI data were obtained from the MODIS 4-day global 500 m LAI product (MCD15A3HV6 Level 4). The global MODIS 4-day 500 m LAI product [START_REF] Myneni | MODIS collection 6 (C6) LAI/FPAR product user's guide[END_REF] was used. The MCD15A3H V6 Level 4, Fraction of Photosynthetically Active Radiation (FPAR) and LAI product is a composite dataset with a temporal resolution of 4 days and a spatial resolution of 500 m. Over four days, the algorithm selects the "optimal" available pixels for inversion from all observations from the two MODIS sensors located on NASA's Terra and Aqua satellites. LUT is generated using the three-dimensional radiative transfer equation [START_REF] Knyazikhin | Synergistic algorithm for estimating vegetation canopy leaf area index and fraction of absorbed photosynthetically active radiation from MODIS and MISR data[END_REF]. Inputs to the algorithm include (i) the type of vegetation structure;

(ii) the geometry of the sun and sensor; (iii) the BRF in the red (648 nm) and NIR (858 nm) spectral bands and (vi) its uncertainty. Each pixel compares the observed and modelled spectral BRFs for a set of vegetation structures and soil patterns that represent the expected range of typical conditions of a given biome type. All canopy/soil patterns and corresponding FPAR values are acceptable solutions if the modelled and observed BRFs differed within a given level of uncertainty. Inversion products included LAI, mean FPAR, dispersion, Standard Deviation (STD) LAI and STD FPAR, and uncertainties [START_REF] Knyazikhin | Synergistic algorithm for estimating vegetation canopy leaf area index and fraction of absorbed photosynthetically active radiation from MODIS and MISR data[END_REF]. In areas of high FC, the optical reflectance saturates, and therefore the sensitivity of the spectral to canopy properties becomes weak. The reliability of inversion parameters under saturated conditions is low because of the large dispersion of the solution distribution. Such inversion results will be flagged in the quality detection layer. A backup method is utilised when the LUT method does not yield a solution. The type of inversion algorithm used is archived in the Quality Assurance (QA) layer. Analysis of the performance of the algorithms shows that the best quality and most accurate inversion results are obtained from the main algorithm [START_REF] Yang | Analysis of prototype Collection 5 products of leaf area index from Terra and Aqua MODIS sensors[END_REF][START_REF] Yang | MODIS leaf area index products: From validation to algorithm improvement[END_REF]. Therefore, the algorithm used for the inversion is an essential quality indicator. The algorithm has interfaces with the MODIS surface reflectance product (MOD09GA) and the MODIS land cover product (MCD12Q1). Technical details of the algorithm can be found in the Algorithmic Theoretical Basis Document (ATBD) [START_REF] Cover | MODIS land cover product algorithm theoretical basis document (ATBD) version 5.0[END_REF].

The multispectral images were first processed. All available S2 and L8 Tier 1 surface reflectance datasets in the GEE by 10 November 2020 were selected, and the datasets were first filtered by selecting images that covered the study area and images with less than 20% average cloud coverage. Each remaining image was then clipped to ensure that only the pixels in the study area were retained. Afterwards, two boundary masks were used in a pre-processing stage. The first was a cloud mask, and the second was a vegetation mask. GEE provided cloud masks for S2 and L8.

Pixels in the cloud mask were removed in the subsequent data processing. The vegetation mask for S2 was obtained by classifying NDVI and GRVI [START_REF] Henrich | The Index Database[END_REF] images using the thresholds proposed in the technical report for S2 [START_REF] Richter | Sentinel-2 MSI-level 2A products algorithm theoretical basis document[END_REF]: classifying (NDVI > 0.40) or (0.36 ≤ NDVI ≤ 0.40 and GRVI > 2.50) pixels were labelled as vegetation pixels. The same rules were used to create the L8 vegetation mask.

The VIs and reflectance values were calculated pixel-wise, and the final VIs and reflectance for each satellite image were obtained by aggregating the VIs and reflectance calculated for each image. In the next step, LAI extraction was performed on each satellite image. Depending on the type and location of the sensor, the study area could correspond to multiple remote sensing images of the day. In this case, the VIs, reflectance and LAI were averaged using the images acquired on that day. The reason for treating the whole study area as one object rather than pixels was to attenuate inconsistencies due to differences in the resolution of remote sensing data. A time-series analysis of the LAI of the study area was carried out to ensure the stability of the vegetation condition in the study area. The independent variable was the date in the time series analysis, and the dependent variable was LAI. The time series of LAI can be divided into trends, seasons and residuals. The trend was the increasing (or decreasing) value in the data series.

Seasonality was the recurring short-term cycle in the data series. After the trend and seasonal components were removed, residual residuals were the time series data. The decomposition was achieved by sliding average using the classical seasonal decomposition [START_REF] Hirsch | A nonparametric trend test for seasonal data with serial dependence[END_REF]. The seasonal multiplicative component was selected for a base-level static analysis using the rolling mean and Dickey-Fuller test. A 99% significance level was applied to the null hypothesis (i.e., no trend) as a threshold for analysis. This test explored trends from one quarter of the year to the following year.

The results of the LAI time series decomposition analysis can be seen in Figure 5.8. All LAI data (marked as observations in the graph) were averaged using quarters (winter, spring, summer and autumn). The LAI for the whole year was first divided into four quarters. The average of the quarters was then calculated and used for the time series analysis. However, there were very few quarters where the corresponding LAI was not available due to cloud cover, in which case an average was assigned between the previous and the next value of this null value. All quarterly LAI data were then decomposed into the trend, seasonal and residual components at a quarterly frequency. The smoothness of the LAI residuals was verified using the rolling mean (Figure 5.9) and the Dickey-Fuller test (Table 5.3). According to the results of the Dickey-Fuller test (Table 5.3), the statistic was much smaller than the critical value at 1%, implying that the data could be considered stable at the 99% confidence level.

The time-series variation of the VIs and LAI was shown in Figure 5.10, where the trend in VI with leaf area index can be seen. VIs were calculated using the default soil adjustment factors (XSAVI = 0.5, XTSAVI = 0.08, XEVI = 1). LAI values ranged from 2 to 6.5. NDVI had the highest values throughout the period (i.e., 0.7 to 0.9), followed by EVI (i.e., 0.4 to 0.7). The range of SAVI and TSAVI was almost the same, from 0.2 to 0.6. These high values indicated a high level of FC in the area. In addition, Figure 5.10 shows the time series of VIs (left vertical axis) and LAI (right vertical axis). It suggests that VIs may or may not be an indicator of trends in LAI over time. Solid circles highlight a period when the LAI is trending against all VIs; the dot-dash circles highlight a period when the LAI is trending in line with the NDVI and against the other VIs; the dotted circles highlight a period when the LAI is trending in line with all VIs. It can be seen that the saturation effect severely limits the ability of VIs to invert the LAI.

(3) Using vegetation indices to estimate LAI Since the optimum soil adjustment factor is the negative value of the abscissa of the intersection of the vegetation isolines and the soil line, the optimum soil adjustment factor of SAVI in high FC areas can be taken as a negative value according to the right shift phenomenon. The ability of the negative soil adjustment factor to resist saturation was verified using linear regression with remote sensing data and global LAI measurements. We found that some of these data were from measurements of different biomes on further examination. Therefore, these data were removed as it was impossible to establish a one-to-one correspondence between LAI and VI. ( 5) Matching of LAI with satellite reflectance imagery: As the temporal accuracy of all measured LAI data is one month, all surface reflectance data from one month of Landsat-5 Tier 1 data were searched and averaged to obtain the final surface reflectance data. Due to cloud cover, some of the reflectance data could not be inverted, so the corresponding LAI data were excluded.

Table 5.4 shows the LAI data and corresponding station information that remained after the above steps. These data were used to evaluate the accuracy of the VIs of LAI inversion. The LAI values ranged from 4.06 to 10.59 and were mainly distributed in Canada, Japan and the USA.

Linear regression was used to establish the relationship between VIs and LAI. Linear regression is a method of establishing a linear relationship between a dependent variable and one or more independent variables [START_REF] Montgomery | Introduction to linear regression analysis[END_REF]. Using a univariate linear regression with VI as the independent variable and LAI as the dependent variable, the linear regression equation is defined as

LAI= VIs c d ⋅ + •••••••••••••••••••••••••••••••••••• (5.9)
where c and d are the intercept and slope of the fitting line, linear regression is widely used to invert LAI from VIs because of its simplicity [START_REF] Zhen | Potentials and Limits of Vegetation Indices With BRDF Signatures for Soil-Noise Resistance and Estimation of Leaf Area Index[END_REF][START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF]. When LAI increases to a certain threshold, saturation effects are observed and no longer maintain a linear relationship with VIs. For discrete problems for which no valid solution method is known, it may be necessary to perform a sequential test for each possibility to determine whether the result is a solution. This exhaustive test of all possibilities is called an exhaustive search [START_REF] Weisstein | Exhaustive Search[END_REF]. The exhaustive search was used to evaluate the performance of a SAVI calculated with the exhaustively optimal soil adjustment factor XSAVI. Here the independent variable was XSAVI, and the dependent variables were R 2 , slope and P of the linear regression equation. Work by Ren et al. [START_REF] Ren | Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands[END_REF] demonstrated that a negative XSAVI value (e.g., -0.2) is acceptable, which allows the specific interval of XSAVI to be set to [-0.3,1] and the bound to be divided into 2000 by 0.001 intervals sections. The ability of SAVI to invert the LAI was investigated throughout the study period using XSAVI sequences. A linear regression model between LAI and SAVI calculated using the XSAVI sequence was used for inversion. MODIS and field measurement LAI data were used here. Finally, the linear regression performance of the various SAVI values calculated from positive or negative XSAVI value was evaluated using R 2 and P.

As can be seen from the two LAI data sets, the best results occurred mainly at XSAVI between -0.18 and 0, indicating that negative XSAVI values were highly desirable and performed well compared to positive values. In the MODIS LAI data set (Figure 5.12.a), XSAVI appeared with a negative R 2 and slope in the range of 0.225 to 1, indicating that the correlation between SAVI and LAI was extremely weak. Furthermore, the negative slope implied a negative correlation between LAI and SAVI. A maximum R 2 (0.2472) and a minimum P (0.0003) occur at XSAVI = -0.148. In the field measurement data, XSAVI shows a negative R 2 and slope between -0.3 and -0. According to an earlier study [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF], the value of the optimal soil adjustment factor was related to the general condition of vegetation: XSAVI should be close but equal or lesser than 1 for sparse vegetation and close but equal or greater than 0 for dense vegetation. However, Ren et al. [START_REF] Ren | Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands[END_REF] and Zhen et al. [START_REF] Zhen | Potentials and Limits of Vegetation Indices With BRDF Signatures for Soil-Noise Resistance and Estimation of Leaf Area Index[END_REF][START_REF] Zhen | A Modified Transformed Soil Adjusted Vegetation Index for Cropland in Jilin Province, China[END_REF] questioned the zero point of the lower boundary for dense vegetation based on field-measured data from arid grasslands and DART simulated data, respectively. Similar to Ren et al. [START_REF] Ren | Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands[END_REF], the negative soil adjustment factor performed optimally in this study.

However, the vegetation condition differed from Ren et al. [START_REF] Ren | Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands[END_REF], where NDVI in the dry grassland was around 0.2, indicating very sparse vegetation. In our study area, NDVI values were mainly distributed around 0.8, indicating very dense vegetation. It is assumed that the reason for the better performance of the negative soil adjustment factor in both cases was that the intersections of the soil line with the vegetation isolines lied in the first quadrant, despite the huge differences in vegetation density.

Field measurements and remote sensing data (Figure 5.12) further confirmed that the lower boundary of the soil adjusted factor could be negative [START_REF] Zhen | Potentials and Limits of Vegetation Indices With BRDF Signatures for Soil-Noise Resistance and Estimation of Leaf Area Index[END_REF]. The locations of the intersections between the soil line and the vegetation isolines are determined by the intercept of the vegetation isoline and that of the soil isoline. If the vegetation isoline intercept is greater than the soil line intercept, the intersection is in the second or third quadrant of the red-NIR plane. The intersection is located in the first quadrant if the vegetation isoline intercept is smaller than the soil line intercept. Furthermore, the intercept depends heavily on the OPs of the leaves and soil, in addition to canopy structural parameters such as LAI. This study further confirms this conclusion using field measurements and remote sensing data. Furthermore, it is noted that NDVI has better performance than SAVI, even though SAVI is considered an improved version of NDVI. Similar suboptimal performance of SAVI has been reported in the literature [START_REF] Ren | Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands[END_REF][START_REF] Ren | Estimation of green aboveground biomass of desert steppe in Inner Mongolia based on red-edge reflectance curve area method[END_REF]. The reason for this poor performance may be that NDVI assumes that the vegetation isolines converge at the origin, and SAVI assumes that the vegetation isolines converge at a point located in the third quadrant (-0.5, -0.5). In this study (a very densely vegetated area), the vegetation isolines converge to a point in the first quadrant. Therefore, the origin-based assumption is better than the point (-0.5, -0.5). This phenomenon can also be seen from the slightly better performance of TSAVI over SAVI since TSAVI assumes that the intersection of the vegetation isolines lies between the origin and the point (-0.5, -0.5). Another possibility is the use of NDVI in the inversion backup algorithm for MODIS LAI products if the radiative transfer inversion method fails, which may also affect the results.

Using hotspot signature vegetation indices to estimate LAI

Hotspot signature vegetation indices are superior to single angular indices for LAI estimation [START_REF] Chen | Predicting leaf area index in wheat using an improved empirical model[END_REF][START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF][START_REF] Chen | Estimation of winter wheat LAI using hotspot-signature vegetation indices[END_REF]. To evaluate the performance of hotspot signature vegetation indices for LAI estimation, NHVI [START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF], together with other kinds of hotspot signature VIs proposed below, are validated for LAI estimation using both simulated and field measured data. Similar to NHVI, see Eq. (1.16), here we defined Hotspot-signature Soil-adjusted Vegetation Index (HSVI) as the product of the multispectral index SAVI and the angular index NDHD

HSVI = SAVI NDHD × ••••••••••••••••••••••••••••••• (5.10)
where SAVI [START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF] is the soil-adjusted vegetation index, expressed by Eq.(1.3).

We replace HDS in the Eq. (1.14) with NDHD and propose NHVI2, which is defined as (5.11) where NDVI [START_REF] Rouse | Monitoring vegetation systems in the Great Plains with ERTS[END_REF] is the normalised difference vegetation index, expressed by Eq.(1.2).

NHVI2 = NDVI NDHD × •••••••••••••••••••••••••••••
In order to evaluate the ability of hotspot signature VIs for soil noise resistance and LAI estimation, the Hotspot-signature 2-band Enhanced Vegetation Index (HEVI2) is also defined as (5.12) where EVI2 [START_REF] Jiang | Development of a two band enhanced vegetation index without a blue band[END_REF] is the two-band atmospheric enhanced vegetation index, expressed by Eq. (1.13).

HEVI2 = EVI2 NDHD × ••••••••••••••••••••••••••••••
Reflectances tend to be minimal at the darkspot because the field of view is maximally filled with shadows from vegetation in the darkspot observation direction [START_REF] Lacaze | Retrieval of vegetation clumping index using hot spot signatures measured by POLDER instrument[END_REF]. The reflectance in the red band depends on the physiological conditions of the green leaves and is particularly low at darkspot, which explains the potentially significant errors in calculating NDHD [START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF], so the reflectance in the NIR band is used to calculate NDHD and HDS. A fixed viewing zenith angle (57.5°) is also added because it is a long-established view angle to give LAI estimation close to actual value [START_REF] Baret | GAI estimates of row crops from downward looking digital photos taken perpendicular to rows at 57.5 zenith angle: Theoretical considerations based on 3D architecture models and application to wheat crops[END_REF][START_REF] Weiss | Methods for in situ leaf area index measurement, Part II: From gap fraction to leaf area index: Retrieval methods and sampling strategies[END_REF][START_REF] Weiss | Review of methods for in situ leaf area index (LAI) determination: Part II. Estimation of LAI, errors and sampling[END_REF][START_REF] Frederic | LAI, fAPAR and fCover CYCLOPES global products derived from VEGETATION: Part 1: Principles of the algorithm[END_REF][START_REF] Weiss | LAI and fAPAR CYCLOPES global products derived from VEGETATION. Part 2: validation and comparison with MODIS collection 4 products[END_REF][START_REF] Comar | A semi-automatic system for high throughput phenotyping wheat cultivars in-field conditions: description and first results[END_REF]. In this paper, the VIs calculated using this view angel are called VIs57.5.

Two data sets were used, including DART simulations (section 6.1) and field measurements.

Field measurements were from the Boreal Ecosystem-Atmosphere Study (BOREAS) as a heterogeneous canopy [START_REF] Chen | Leaf area index of boreal forests: Theory, techniques, and measurements[END_REF][START_REF] Gamon | Remote sensing in BOREAS: Lessons learned[END_REF]. Data were downloaded from the website (https://daac.ornl.gov/cgi-bin/dataset_lister.pl?p=2). For accuracy, LAI data were used rather than LAIe data. In situ data used in this study were from BOREAS [START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF][START_REF] Chen | Leaf area index of boreal forests: Theory, techniques, and measurements[END_REF][START_REF] Gamon | Remote sensing in BOREAS: Lessons learned[END_REF]. This area was chosen because boreal forests have a more complex canopy structure [START_REF] Chen | Predicting leaf area index in wheat using an improved empirical model[END_REF]. Field measurements were taken from two specific study areas of the BOREAS project (Figure 5.13), the Southern Study Area (SSA) near Prince Albert and Candle Lake, Saskatchewan and Thompson, and the Northern Study Area (NSA) near Nelson House, Manitoba [START_REF] Brandt | The extent of the North American boreal zone[END_REF]. Bi-directional reflectance factors were acquired by the POLDER [START_REF] Gamon | Remote sensing in BOREAS: Lessons learned[END_REF] instrument on the helicopter and tower. The instrument had five bands: 443 nm (blue), 550 nm (green), 670 nm (red), 864 nm (NIR1) and 910 nm (NIR2). At the NSA sampling point, multi-angle reflectance data sets were obtained by varying the angle of observation in 1° increment. Data were randomly varied at the SSA sampling site in the observed zenith and azimuth angles. The BOREAS project [START_REF] Chen | Leaf area index of boreal forests: Theory, techniques, and measurements[END_REF] measured LAI intensively along three transects at 10 m step intervals. Hemispheric photographs were taken on some of the transects. A detailed description of the field measurement site can be found in the literature [START_REF] Chen | Leaf area index of boreal forests: Theory, techniques, and measurements[END_REF]. All data were carefully examined, and six sampling sites were found to satisfy the conditions for calculating the hotspot signature VI while having LAI data. The Ross-Li model [START_REF] Wanner | On the derivation of kernels for kernel-driven models of bidirectional reflectance[END_REF] was used to fit the bidirectional reflectance distribution function (BRDF) curves to the available data points and to find reflectance values for the specific observed directions to calculate the VIs [START_REF] Lacaze | Retrieval of vegetation clumping index using hot spot signatures measured by POLDER instrument[END_REF][START_REF] Leblanc | Mapping vegetation clumping index from directional satellite measurements[END_REF][START_REF] Chen | Multi-angular optical remote sensing for assessing vegetation structure and carbon absorption[END_REF]. The exponential fit performed the best of all fits for both types of vegetation scenes due to the "saturation phenomenon", as shown in Eq. (5.8). Multispectral VIs performed well in homogeneous simulated data (Figure 5.15), with the highest R 2 and lowest RMSE, followed by VIs57.5. The hotspot signature VIs did not perform better than the other two VIs in the homogeneous scene, which Hasegawa et al. [START_REF] Hasegawa | Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures[END_REF] ignored. However, it can be noted that HEVI2, HSVI and NHVI2 outperform NHVI. Good performance of the hotspot signature VIs can be seen in both the heterogeneous simulated data (Figure 5.16), and the field measured data (Figure 5.17). Furthermore, the performances of VIs57.5 were unstable and outperformed the multispectral VIs in the heterogeneous simulated data when the SZA was small (SZA = 30°). However, when SZA was large, only NDVI57.5 outperformed its multispectral counterpart VIs (SZA = 45° and SZA = 60°). In the field measurements, VIs57.5 outperformed the multispectral VIs overall.

Evaluation of inversion of leaf area index without noise interference

Evaluation of inversion of leaf area index with random noise interference

As the VIs saturate with LAI over 3, 14 different sets of random errors were added to the reflectance to compare the robustness of the VIs. The maximum values of these random errors ranged from 0 to 0.007 in steps of 0.0005, with the expectation always being zero. The absolute value of the maximum error of the random noise was set according to the minimum value of all reflectances to ensure a constant positive reflectance. The R 2 and RMSE of the exponential fit model were used to compare the performance of the VIs considered.

For the homogeneous simulated data (Figure 5.18), all three multispectral VIs (NDVI, SAVI and EVI2) performed well on the LAI inversions, and HEVI2 and HSVI performed about the same as SAVI57.5 and EVI257.5, while NHVI and NHVI2 did not perform as well. Furthermore, HEVI2

and HSVI performed similarly to SAVI57.5 and EVI257.5, while NHVI and NHVI2 did not perform well. For both the heterogeneous simulated data (Figure 5.19) and the field measurements (Figure 5.20), the hotspot signature VIs performed best in the LAI inversion, having the highest R 2 and lowest RMSE, resulting in the most reliable LAI estimation. VIs57.5 could give overall reasonable LAI estimates, especially NDVI57.5. On the other hand, SAVI57.5 and EVI257.5 performed unsteadily.

= - ∫ ••••••••••••••••••••••• (6.1)
where VI(LAI) ����������� is the mean value of the VIs at the corresponding LAI. The "signal" compares the mean value of the VI corresponding to the maximum and minimum LAI canopy. The "noise" is measured by the area between the highest and the lowest curves (i.e., the product of the change in the index due to the change in the OPs of the soil and the LAI range value). Since a desirable index is susceptible to vegetation spectral while unaffected by soil spectral, a larger SNR indicates a better VI.

The second criterion for evaluating VI is based on the LAI dependent parameter (T) [START_REF] Gilabert | Designing a generalized soiladjusted vegetation index (GESAVI)[END_REF],

which is defined as follows

( ) LAI LAI T σ σ = •••••••••••••••••••••••••••••••••••••• (6.2)
where σLAI is the standard deviation of VIs corresponding to a specific LAI value and 𝜎𝜎 � is the standard deviation of VIs corresponding to the whole LAI variation range. In contrast to SNR, this parameter decreases as the VI efficiency increases.

The hotspot signature VIs performed best under both types of canopies compared to the multispectral VIs and VIs57.5. Of all the hotspot signature VIs, HSVI had the highest SNR ratio (Figure 6.1), which averages the SNR of the nine sets of ALA and is presented here. HSVI had the lowest T value (Figure 6.2), which averaged the values of the T corresponding to the nine sets of ALA and was presented here. VIs57.5 was superior to the multispectral VIs, especially when SZA was high (SZA=45° and SZA=60°). Of all VIs57.5, SAVI57.

Effects of atmosphere on vegetation indices

Absorption and scattering of the atmosphere can lead to changes in wavelength radiance [START_REF] Gastellu-Etchegorry | DART: recent advances in remote sensing data modeling with atmosphere, polarization, and chlorophyll fluorescence[END_REF][START_REF] Wang | DART: Improvement of thermal infrared radiative transfer modelling for simulating top of atmosphere radiance[END_REF], which can cause changes in band reflectance. Four VIs (NDVI, SAVI, EVI, EVI2) were used as examples to explore the atmosphere effect. VIs calculated with TOA and SR data from L8 and S2 sensors were compared to assess the atmosphere effect on VIs.

Here we used a similar pre-processing step as in section 6.3. TOA images were needed to evaluate the atmosphere effect on VIs. However, it was found that TOA and SR images did not always exist in pairs in the GEE. Therefore, an additional filter was needed for the SR images: all SR images without a corresponding TOA image were also removed, and this filter was only applied for the study of atmosphere effect (see section 6.2) and not for the EVI2 cross-sensory calibration (see section 6.3).

The VIs calculated using the SR data were considered true values, while the VIs calculated using the TOA reflectance images were considered the measured values. The relative error E atmosphere,VI(x,y) of atmosphere effect for each pixel (x,y) was then defined as [START_REF] Zhang | Approach for a vegetation index resistant to atmospheric effect (in Chinese with English abstract)[END_REF] ( ) ( ) ( ) ( )

- = × ••••••••••••••• (6.3)
where VISR(x,y) was the VIs calculated with SR images (i.e., NDVI, SAVI, EVI, EVI2), and VITOA(x,y) was the VIs calculated with TOA reflectance image. The Eatmosphere,VI(x,y) of all pixels in an image was then summed and averaged to the Eatmosphere,VI,i for each image i. The VIs calculated using TOA and SR data were compared separately to evaluate the atmosphere effect. Figure 6.3 shows the Eatmosphere,VI distribution for the four VIs, including the EVI, EVI2, NDVI and SAVI. Table 6.1 depicts the statistical distribution of Eatmosphere,VI. It should be noted that EVI2 and EVI produced opposite results when influenced by the atmosphere.

Eatmosphere,NDVI, Eatmosphere,SAVI and Eatmosphere,EVI2 were mainly distributed in the negative region, indicating that NDVI, SAVI and EVI2 were underestimated when influenced by the atmosphere.

However, Eatmosphere,EVI was mainly located in the positive region, suggesting that atmospherically influenced EVI was overestimated. In addition, S2 had the smallest median value of Eatmosphere,EVI, indicating good atmospheric resistance in most cases. However, the maximum and minimum values of Eatmosphere,EVI were also non-negligible (Table 6.1), indicating that EVI may cause poor results in some extreme cases. The atmosphere effect was approximately the same for EVI2, NDVI and SAVI.

The literature [START_REF] Myneni | Atmospheric effects and spectral vegetation indices[END_REF][START_REF] Zhang | Approach for a vegetation index resistant to atmospheric effect (in Chinese with English abstract)[END_REF] results that the atmosphere effect can underestimate common dual-band VIs were confirmed by analysing all available global L8 Tier 1 and S2 SR data up to 2021 (Figure 6.3), with a larger temporal and spatial resolution scales data. Aerosols can reduce the difference between NIR and red reflectance captured by orbital sensors because the atmosphere effect has a net increase in reflectance in the red band [START_REF] Myneni | Atmospheric effects and spectral vegetation indices[END_REF]. The results suggest that for the threeband VIs like EVI, the atmosphere effect may have different impacts than for two-band VIs. EVI2, NDVI and SAVI are underestimated by the atmosphere effect (and potentially lead to an underestimation of inversion of LAI), while EVI is overestimated (and potentially leads to an overestimation of LAI inversion). It is suspected that this overestimation may be due to the positive net atmosphere effect in the blue band. As a result, the denominator of EVI becomes smaller, and the fractional value becomes large. Consequently, an opposite reaction to the atmosphere effect would inevitably increase the difference between EVI and EVI2. In addition, EVI2, NDVI and SAVI show minor differences in terms of atmosphere effect. Therefore, the study emphasises the importance of considering the different reactions of VIs to the atmosphere effect.

Despite its low median, E atmosphere,EVI had much high maximum and minimum values than the other three double-band VIs (e.g., EVI2, NDVI, SAVI). All images with a very large Eatmosphere,EVI were carefully examined and found that they were concentrated mainly in the Greenland region.

The S2 classified images were then examined, and it was found that the number of pixels classified as the vegetation was minimal in the image and that the larger Eatmosphere,EVI was mainly from these pixels. Due to experimental limitations, it is challenging to determine whether the vegetation pixels corresponding to the S2 classification bands in this region correspond to actual vegetation or whether they are misclassified pixels. Further research is needed on the accuracy of S2 classification products.

Effect of spectral response functions on vegetation indices

The Enhanced Vegetation Index (EVI) [START_REF] Huete | A comparison of vegetation indices over a global set of TM images for EOS-MODIS[END_REF] represents an optimized estimation of landscape vegetative conditions in space, spectrum, and radiometry evaluated by the two regular Terra and Aqua Modis products. Compared with NDVI, EVI can weaken the effects of atmosphere and soil background on vegetation index [START_REF] Huete | A comparison of vegetation indices over a global set of TM images for EOS-MODIS[END_REF]. Besides, EVI has more excellent biophysical linearity with indicators for vegetation and a more comprehensive range of values that estimates the leaf area index (LAI) [START_REF] Kang | How Universal Is the Relationship between Remotely Sensed Vegetation Indices and Crop Leaf Area Index? A Global Assessment[END_REF]. NDVI starts to reach saturation thresholds in cropland during the growing season, while EVI is more robust during development [START_REF] Wardlow | Analysis of time-series MODIS 250 m vegetation index data for crop classification in the US Central Great Plains[END_REF]. EVI is less prone to saturation in forests [START_REF] Xiao | Satellite-based modeling of gross primary production in an evergreen needleleaf forest[END_REF]. The EVI has also been shown to effectively monitor, detect, and assess seasonal changes in evergreen forests [START_REF] Deng | Seasonal variation of MODIS vegetation indexes and their statistical relationship with climate over the subtropic evergreen forest in Zhejiang, China[END_REF].

However, apart from the red and NIR bands, EVI requires the sensor systems to be equipped with a blue band, making EVI less commonly used than NDVI. For example, the generation of EVI products is not possible for sensors without a blue band, such as the Advanced Very High Resolution Radiometer (AVHRR). In addition, the spatial resolution of MODIS' blue band is 500 m, lowering the resolution of MODIS EVI products from 250 m to 500 m. Therefore, to solve these issues, the 2-band enhanced vegetation index (EVI2) [START_REF] Jiang | Development of a two band enhanced vegetation index without a blue band[END_REF] is developed without requiring the blue band. EVI2

shows good agreements with EVI in MODIS data [START_REF] Jiang | Development of a two band enhanced vegetation index without a blue band[END_REF][START_REF] Motohka | Evaluation of Sub-Pixel Cloud Noises on MODIS Daily Spectral Indices Based on in situ Measurements[END_REF]. However, for some non-MODIS sensors (i.e., PlanetScope), it is observed that the performance of EVI and EVI2 differ significantly [START_REF] Myers | Assessing the Impact of Satellite Revisit Rate on Estimation of Corn Phenological Transition Timing through Shape Model Fitting[END_REF]. Although many studies [START_REF] Fan | A global study of NDVI difference among moderate-resolution satellite sensors[END_REF][START_REF] Steven | Intercalibration of vegetation indices from different sensor systems[END_REF] for each satellite are divided into three categories (primary, secondary and real-time).

(i) Tier 1 (T1): data meeting geometric and radiometric quality requirements.

(ii) Tier 2 (T2): data not meeting the requirements of Tier 1.

(iii) Real-time: data that has not yet been evaluated (takes as long as one month).

At the same time, Landsat offers images of three products:

(i) Raw images: Digital Number (DN) values, representing the scaled and calibrated sensor radiance images.

(ii) Top of the atmosphere (TOA): A calibrated reflectance image of the TOA.

Band correlation coefficient c in EVI2 was obtained by fitting MODIS data to assess the influence of the spectral response function on the VIs. SR data from L8 and S2 were used to perform a linear regression between the blue and red bands. As SR products were used, the atmosphere effect was assumed to be almost completely removed. Therefore, the atmosphere effect was not considered.

Pre-calibration error assessment

All processes of remote sensing data were done on Google Earth Engine (GEE), a planetarystar platform for earth science data and analysis [START_REF] Gorelick | Google Earth Engine: Planetary-scale geospatial analysis for everyone[END_REF]. Earth Engine's public data archive includes over forty years of historical imagery and scientific datasets updated and expanded daily. Preprocessing of the data consisted of two main steps: filtering and masking. First, all available global L8 Tier 1 SR data and S2 SR data until 1 January 2021 were selected. These data were then filtered based on cloudiness to ensure that each view image had less than 20% cloudiness. After filtering, the cloud-covered areas were masked. The cloud mask provided by GEE was first implemented on the images to remove cloud pixels. The vegetation mask was then used to extract the vegetation pixels. S2's vegetation mask was created by using thresholds for NDVI and the simple ratio NIR / green ratio vegetation index (GRVI) [START_REF] Henrich | The Index Database[END_REF], a land cover classification rule and corresponding thresholds provided by S2's technical report [START_REF] Richter | Sentinel-2 MSI-level 2A products algorithm theoretical basis document[END_REF]: with (NDVI > 0.40) or (0.36≤ NDVI ≤ 0.40 and GRVI > 2.50) were classified as vegetation pixels. These rules were also used to create the L8 vegetation mask [START_REF] Zhen | Using the Negative Soil Adjustment Factor of Soil Adjusted Vegetation Index (SAVI) to Resist Saturation Effects and Estimate Leaf Area Index (LAI) in Dense Vegetation Areas[END_REF].

As EVI2 is designed to replace EVI, the EVI2 value should be sufficiently close to EVI.

Therefore, the relative error (EEVI2) between EVI2 and EVI was used to assess accuracy by treating EVI as "real data" and EVI2 as "measured data".

( ) ( ) ( ) ( ) EVI2 EVI2 , EVI , , 100% EVI , x y x y E x y x y - = × ••••••••••••••••••• (6.4)
EEVI2(x,y) was calculated pixel-wise. Obviously, the smaller the EEVI2,i was, the better performance of EVI2, and vice versa.

Figure 6.5 shows the EEVI2 before calibration on L8 and S2, where EVI is considered the "real" data and EVI2 is considered the "measured" data. Table 6.2 depicts the statistical distribution of EEVI2. The EVI2 calculated from L8 and S2 mainly was underestimated. Furthermore, it could be observed that the EVI2 calculated from L8 performed better than that of S2.

Band correlation coefficient calibration

The cross-sensory calibration aimed to make the EVI and EVI2 calculated from the same satellite sensors as equal as possible. The spectral correlation coefficients (a and b) between the blue and red band reflectances were used to calibrate the spectral correlation coefficients c for EVI2 calculated for S2 and L8, respectively. A linear equation was fitted between the blue band reflectance Rblue and the red band reflectance Rred for the S2 and L8 sensors

red blue R a R b = × + •••••••••••••••••••••••••••••••••• (6.5)
where a and b are the linear regression coefficients for slope and intercept, respectively, and the linear regression is performed pixel-wise, with one image corresponding to a pair of a and b values. The coefficient with the highest frequency was finally determined as the final value and used to calibrate c in EVI2. The calibrated band correlation coefficients cL8 and cS2 were applied to the L8 and S2 sensors. In addition to using remote sensing images, the vegetation spectral library of Chapter 5 from the USGS Spectral Library, Version 7 [START_REF] Kokaly | USGS Spectral Library Version 7: US Geological Survey Data Services 1035[END_REF], was used to validate the improvements of the calibrated EVI2 on L8 and S2 sensors. Field measurements were less affected by soil, topography, and bidirectional reflectance effects than remote sensing data. Therefore, field-measured data were used for validation. The spectral library contained spectral measurements using laboratory, field Two vegetation spectral libraries from the USGS Spectral Library, Version 7, Chapter 5

Vegetation, were used. The first was s07LSAT8. The spectral resolution of this spectral library has been resampled to fit the SRF of the OLI L8. The seven-band SRF of this sensor covered visible and short-wave infrared wavelengths, and the ENVI 5.3 software provided the SRF. The second was the s07SNTL2 spectral library, whose spectral resolution was resampled into the SRF of the S2 MSI.

Careful examination of all the data revealed negative reflectance values in the blue and red bands in several spectral curves. Therefore, these spectral curves were simply removed. The confidence interval is estimated using a bootstrap. sensing data using hexagonal binning maps to calculate pixel-wise coefficients. In the S2 data, the most frequent point was located at (1.60, 0.00), which was closer to the value derived from the fieldmeasured data (1.262, 0.021) than the original EVI2 value (2.08, 0.00). In the L8 data, most of the fitted parameters were located at the point (1.80, 0.00). This point was closer to the point (1.262, 0.021) derived from the field measured data than the original point (2.08, 0.00). Therefore, it was assumed that cL8 = 1.80 and cS2 = 1.60 may be a better proxy than c = 2.08 with S2 and L8 when calculating EVI2 using S2 or L8 surface reflectance data. To distinguish the calibrated EVI2 from the original EVI2, EVI2L8 and EVI2S2 represented the EVI2 calculated with cL8 and cS2. Good agreement between EVI and EVI2 has been reported in MODIS data [START_REF] Motohka | Evaluation of Sub-Pixel Cloud Noises on MODIS Daily Spectral Indices Based on in situ Measurements[END_REF]. However, EVI2 calculated using PlanetScope data was in suboptimal agreement with EVI [START_REF] Myers | Assessing the Impact of Satellite Revisit Rate on Estimation of Corn Phenological Transition Timing through Shape Model Fitting[END_REF]. Moreover, Myers et al. [START_REF] Myers | Assessing the Impact of Satellite Revisit Rate on Estimation of Corn Phenological Transition Timing through Shape Model Fitting[END_REF] attributed this inconsistency to undetected atmospheric water vapour or aerosols affecting the reflected spectra of the inversion. However, it is suspected that the inconsistency between EVI and EVI2 varies from sensor to sensor. Kang et al. [START_REF] Kang | How Universal Is the Relationship between Remotely Sensed Vegetation Indices and Crop Leaf Area Index? A Global Assessment[END_REF] found that EVI2 performed better in the L8 sensor than S2 (Figure 6.5). Furthermore, the linear regression results (Figure 6.7) also indicated that cL8 was closer to c than cS2. Therefore, it is speculated that the better performance of L8 on the cross-sensory calibration may be due to c being closer to cL8 than cS2 (Figure 6.7). Thus, the results suggest that the difference between SRFs may also be one of the main reasons for the variation in EVI and EVI2 performance. Based on the spectral correlation theory, the red and blue band correlation coefficients of EVI2 were calibrated on the L8 and S2 sensors, respectively, and the results (Figure 6.8) show that the calibrated EVI2L8 and EVI2S2 were closer to the EVI than the values given by EVI2 [START_REF] Jiang | Development of a two band enhanced vegetation index without a blue band[END_REF]. Differences in the regression parameters between laboratory measurements and remote sensing data may be due to external noise, such as vegetation, soil and atmosphere. Current cross-sensor studies of VIs [START_REF] Steven | Intercalibration of vegetation indices from different sensor systems[END_REF] have mainly been devoted to NDVI and EVI to make VIs obtained from different sensors comparable. However, because EVI2 was designed to replace 3-band EVI, we aimed to make EVI2 closer to EVI from the S2 and L8 sensors. It was where the research approach differs from another cross-sensory study [START_REF] Steven | Intercalibration of vegetation indices from different sensor systems[END_REF]. Furthermore, unlike Jiang et al.

Post-calibration error assessment

[106], the difference between EVI and EVI2 was compared using relative errors rather than absolute errors. The reason for this was that the overall values of the VIs in the low FC areas were small, so absolute errors may create the illusion of some minor errors in the low FC areas, whereas the use of relative errors could avoid this problem.

Summary

Spectral response function, atmosphere effect and soil noise may impair the vegetation signal expressed using vegetation indices. We compare the soil noise resistance among hotspot-signature vegetation indices and directionality indices, study the atmosphere effect on NDVI, SAVI, EVI and EVI2, calibrate the spectral correlation coefficient of EVI2 for S2 and L8 sensors.

(1) Hotspot-signature vegetation indices can well resist the influence of soil noise.

(2) The atmosphere effect causes the common two-band VIs (i.e., NDVI, SAVI, EVI2) to be underestimated. However, the atmospheric effect causes the EVI to be overestimated, which stems from the fact that the atmosphere effect increases the reflectance of the blue band. Thus, the atmosphere effect can increase the difference between EVI and EVI2.

(3) The spectral correlation coefficient c for EVI2 is calibrated (cL8 = 1.80, cS2 = 1.60) to produce an EVI2 closer to the EVI with L8 and S2 satellites. The new EVI2 calculated using cL8 and cS2 are named EVI2L8 and EVI2S2. After calibration, the R 2 for S2 increases by 6%, and the RMSE decreases by 25%. Meanwhile, the R 2 for L8 increases by 3%, and the RMSE decreases by 17% based on the USGS vegetation Spectral Library.

(4) For the S2 and L8 sensors, the EVI2 values are slightly smaller than EVI, with the median value of EEVI2 equal to -0.025 (L8) and -0.063 (S2). Furthermore, the EVI2 calculated from L8 is closer to the EVI than S2 because cL8 is closer to c than cS2.

for an average LAI equal to 5.35) and more generally in areas with larger LAI values (e.g., the optimal soil adjustment factor takes a value of -0.183 for an average LAI equal to 6.72).

(4) The influences of soil, atmosphere and SRF on VIs are analysed. The spectral correlation coefficients between S2 and L8 red and blue bands are calibrated and applied to the EVI2 calculation.

The calibration coefficients are validated using the USGS vegetation spectral library. The results of the atmosphere effect show that the atmosphere effect exacerbates the differences between EVI and EVI2. As a result, atmospherically influenced EVI is overestimated, and atmospherically influenced EVI2 is underestimated. The results of the quantitative analysis show that EVI2 is slightly smaller than EVI in the S2 and L8 sensors. Assuming that EVI is the true value and EVI2 is the measured value, the median relative error equals -0.025 for L8 and -0.063 for S2. The minor relative error in EVI2 for L8 than for S2 is that c L8 is larger than cS2 and is closer to the original c value of 2.08. The proposed two new spectral band correlation coefficients, cL8 = 1.80 and cS2 = 1.60, are estimated from SR data for S2 and L8. The new EVI2 calculated with cL8 and cS2 is named EVI2L8 and EVI2S2.

Validations based on the USGS vegetation spectral library show a 6% increase in R 2 between EVI and EVI2S2 and a 25% decrease in RMSE for S2 for the relative error. There is a 3% increase in R 2 between EVI and EVI2L8 and a 17% decrease in RMSE for L8 for the relative error.

Innovations

The main points of innovations are: (2) Based on the DART calibration, the inversion of the leaf optical properties can be achieved at the sub-pixel scale.

Based on the DART calibration, canopy spectral unmixing and leaf OP inversion is performed based on the DART model to invert leaf OPs from a mixed spectrum containing three components (vegetation, soil, urban) or more. Compared to current leaf inversion methods based on the assumption of a homogeneous canopy, this method can consider the strong heterogeneity and nonplanarity of urban areas and the spectral-spatial variability of the endmembers. In addition, the method requires a low number of spectral bands (at least one band) to achieve the inversion and therefore has few requirements on sensors.

(3) The phenomenon of the right shift of the intersections is proposed to mitigate the saturation effect of SAVI in high FC areas.

The "right shift phenomenon" is proposed to describe the movement of the intersection points between the soil line and the vegetation isolines in the red-NIR plane with the FC, and the formula is derived to prove it. The right-shift phenomenon shows that, as the FC increases, the intersections between the vegetation isolines and the soil line gradually move towards the positive red band coordinate direction. When the intercept of the vegetation isoline is smaller than the intercept of the soil line, the intersection point can eventually reach the positive region. The right-shift phenomenon successfully solves the two contradictory debates on describing the movement of intersection points between vegetation isolines and soil lines in the current academic community and achieves a dialectical unification in a new perspective. Because the optimal soil adjustment factor is the negative value of the abscissa of the intersection of the vegetation isolines and the soil line, based on the right-shift phenomenon, the hypothesis that the optimal soil adjustment factor should be negative in the high FC area is proposed and verified. This hypothesis can mitigate the saturation effect of SAVI and improve the accuracy of LAI inversion.

Shortcomings and perspectives

Reflectance simulations based on the 3D radiative transfer models are an essential theoretical basis for understanding 3D radiative transfer simulations and an effective tool for the inversion of vegetation bio-optical parameters. Researches have been carried out in both modelling and applications, but there are still a number of shortcomings that need further improvements.

(1) Reflectance based on 3D maize structure simulations is generated under minimal experimental conditions. A more comprehensive range of experimental conditions should be investigated. Also, the simulated BRF has not been validated with field measurement data. It is the next stage of work. The currently considered growth control factors are relatively few. Several environmental factors that significantly impact plant growth, such as nitrogen, phosphorus, and potassium, will be added to the future growth model. The potential of this growth model is enormous because it can bridge the knowledge of plant morphology and remote sensing rather than the previously used simplified canopy structure parameters.

(2) A number of uncontrollable external factors may significantly affect the accuracy of leaf spectral inversions, such as geometric co-registration errors or user-defined ratios between leaf reflectance and transmittance. The noise resistance for the inversion methods of leaf OPs should be adequately improved. Besides, even though this thesis only involves the spectral from visible to short wave infrared, the potentials of DART calibration are far more from this. The a priori abundance information can well solve the mixed pixels problem in the thermal infrared band. In addition, combining the temperature and emissivity separation algorithm and DART calibration allows for the temperature emissivity decomposition of non-isothermal pixels. The consideration of spatial heterogeneity and non-planarity in the DART calibration makes it suitable for use in complex study areas such as cities. However, one limitation of the DART inversion is the use of the 3D urban geometry database, which is still not widely open access for many cities. However, with the development of surveying and mapping technology, we have reason to believe that 3D urban data structure data will become more and more common. The potential application of the inverted optical properties of urban materials is enormous, such as estimating energy budget (urban climate/heat island, etc.). It is vital to the urban planning (e.g., to reduce or prevent turbulent sensible heat flux and anthropogenic heat flux hot spots), health studies (e.g., to estimate the impact on thermal comfort) and future-proofing (e.g., to plan and implement interventions to reduce heat emissions).

(3) The accuracy of the proposed inversion of LAI using the negative soil adjustment factor in SAVI is much higher than using a positive soil adjustment factor in the high FC area, but the accuracy still falls short compared to methods such as RT-based inversion or active remote sensing.

Therefore, the proposed method is better suited to using common optical satellites with a single (4) Determination of L-system rules: This section specifies the iteration rules of the L-system, i.e., the predecessor and successors. The rule layer (n, h, delta) creates a series of n polygons ( by repeatedly calling the rule "facet". )

Their bases are engraved on a circle. This rule is applied to construct a "polygonal approximation" of a cone or cylinder.

The rule trunk(n, m) creates a polygonal approximation of a cylinder with a height of m layers and has a circular base that approximates a symmetric n-sided polygon

The rule object(n, m, h, delta) constructs a flowering plant with stem composed of rule trunk (n, m). The flower (constructed using regular "layers") has length h and grow from the trunk with angle delta. By changing the parameters of h and delta, the different stages of flowering process can be obtained. As the flower becomes larger, it expands from the stem with an increasing angle.

The rule animation(kx,ky, n,m, h1,h2, delta1,delta2, nsteps) accomplish all procedure by continually constructing a field of kx times ky flowers (kx flowers in the x direction and ky in the y direction) and saving various phases of development into files an.OOl, an.002, an.003, and so on.

The built-in operator "setname" is used to set the structure file name. In addition to kx and ky, the rule also describes the parameters n and m of the stem (such as in the trunk rule), the parameters h1 and delta1 of the initial growth stage of the flower, the parameters h2 and delta2 of the final growth stage of the flower, and the integer n steps that represent the number of steps required for animation from the initial stage to the final stage. The polygon that describes the flower field is generated in a loop by rule animation, where the parameters h and delta of the rule "object" are linearly interpolated between the values of the initial and final stages of growth. When a growth phase is completed, the resulting polygon is recorded in the file, the file counter increases, and the number of polygons is reset to 0. The built-in operator "flush" executes these three purposes.

Figure 1 . 1 15 Figure 1 . 2 19 Figure 2 . 1 22 Figure 2 . 2 25 Figure 2 . 3

 111512192122222523 Figure 1.1 Example of changes in the reflectance of red and near-infrared bands according to the view zenith angle changes. Data were simulated using the DART model with homogeneous canopy and LAI equal 0.8............................................................................................................................ 15 Figure 1.2 Thesis structure. The whole thesis can be divided into simulation and inversion parts. For the inversion part, it can then be divided into leaf OP inversion and LAI inversion part. ........ 19 Figure 2.1 Schematic diagram of vegetation canopy in one-dimensional mixed media. The leaves are evenly and randomly distributed in the whole space. ................................................................ 22 Figure 2.2 Flux tracking theory of SAIL model. It is based on the four radiative fluxes method. . 25 Figure 2.3 Iterative diagram of the tree structure. The whole scene comprises trees and soil, and trees comprise trunk and crown. ..................................................................................................... 30

Figure 2 . 4 31 Figure 2 . 5 31 Figure 2 . 6 33 Figure 2 . 7

 24312531263327 Figure 2.4 Tree structure diagram created using ELSYS. .............................................................. 31 Figure 2.5 Schematic diagram of the radiosity model principle. The green arrow represents the incident energy, and the yellow arrow represents the outgoing energy. .......................................... 31 Figure 2.6 Modelling of a real structural mono-tree scene using ELSYS. (a) Real structure single tree; (b) ELSYS modelling of the single tree. ................................................................................. 33 Figure 2.7 Comparison of the simulated BRF of real structural single tree scene using RGM before and after acceleration. (a) Blue band in the principal plane direction (PP); (b) Green band in the PP direction; (c) Red band in the PP direction; (d) NIR band in the PP direction; (e) Blue band in the cross-principal plane (CP) direction; (f) Green band in the CP direction; (g) Red band in the CP direction; (h) NIR band in the CP direction .................................................................................... 35

Figure 2 . 8 36 Figure 2 . 9 37 Figure 2 . 37 Figure 2 . 11 38 Figure 3 . 1 43 Figure 3 . 2 44 Figure 3 . 3

 28362937237211383143324433 Figure 2.8 Accuracy validation of post-accelerate RGM model for real structural single tree scene. (a) Blue band; (b) Green band; (c) Red band; (d) NIR band ........................................................... 36 Figure 2.9 Simplified heterogeneous canopy scene. Trees are randomly distributed in the scene, and small triangles fill the crown with uniform random distribution. ................................................... 37 Figure 2.10 Simulated BRF of RGM and DART model with simplified heterogeneous canopy. (a)

Figure 4 . 1

 41 Figure 4.1 Schematic urban scene: ground, water, shrubs, trees, and buildings with the flat, triangular prism and frustum of a square pyramid roof. (a) 3D view. (b) Top view. ....................... 50

Figure 4 . 2 52 Figure 4 . 3

 425243 Figure 4.2 (a) DART pseudo satellite image of the 3-D scene. It has some spectral confusions (e.g., the different reflectance for nine roofs and the same reflectance of three roofs as ground reflectance). (b) DART image simulated with constant SS per urban component. (c) DART image after the calibration procedure. All images are RGB colour composites. ..................................................... 52 Figure 4.3 Evolutions of the mean absolute value of the relative error of the component SS with the iterations in the green band. The mean absolute value of the relative error at iteration 0 is due to the selected initial SS value. Relative errors are calculated pixel-wise for each component, and a mean statistical value of the absolute value of the relative error is calculated. Similar trends are found in the other bands. ............................................................................................................................... 56

Figure 4 . 4

 44 Figure 4.4 Final calibration results of the error analysis for ground (a and b), roof (c and d), water (e and f), shrub (g and h), and trees (i and j) in the green band in the spatial and frequency domain. Relative errors are calculated pixel-wise. Blank in the image means there is no corresponding component in that pixel. Violin range is limited to within the range of the relative error. A similar phenomenon was also observed for other bands. ............................................................................ 57

Figure 4 . 5

 45 Figure 4.5 Sensitivity analysis results. The eFAST method is used for sensitivity analysis based on 420 cases. Sensitivity indices are calculated per component per band and stacked over all bands. All subplots share the same x-axis. In descending order of magnitude, the effect magnitudes are SZA, satellite spatial resolution, pixel-shift, landscape model inaccuracy, and MTF. ............................. 61

Figure 4 . 6

 46 Figure 4.6 Mean absolute value of the relative errors on the retrieved SS of scene components for the blue, green, red, NIR, SWIR1, and SWIR2 bands. The relative errors were calculated for each component in terms of pixels. The absolute value of the relative error was the average of all 420 cases. The violin range was constrained to the boundary of the mean relative error. ..................... 62

Figure 4 . 7

 47 Figure 4.7 Basel 3-D scene. (a) Top view. (b) Side view of the red square in (a). Buildings (roofs: dark, walls: light grey), trees (green), river (blue), and ground (yellow). ....................................... 63

Figure 4 . 8

 48 Figure 4.8 (a) PlanetScope image in the Basel area. (b) DART image simulated with constant SS per urban component. (c) DART image after the calibration procedure. All images are RGB colour composites. ...................................................................................................................................... 64

Figure 4 . 9 65 Figure 4 . 10 66 Figure 5 . 1 71 Figure 5 . 2 72 Figure 5 . 3 73 Figure 5 . 4 74 Figure 5 . 5 75 Figure 5 . 6 76 Figure 5 . 7 79 Figure 5 . 8 80 Figure 5 . 9 81 Figure 5 . 10 82 Figure 5 . 11 84 Figure 5 . 12 85 Figure 5 . 13 87 Figure 5 . 14 88 Figure 5 . 15 89 Figure 5 . 16 90 Figure 5 . 17 90 Figure 5 . 18 97 Figure 5 . 19 98 Figure 5 . 20 98 Figure 6 . 1 102 Figure 6 . 2 103 Figure 6 . 3 104 Figure 6 . 4 108 Figure 6 . 5 109 Figure 6 . 6 111 Figure 6 . 7 111 Figure 6 . 8

 49654106651715272537354745575567657795880598151082511845128551387514885158951690517905189751998520986110262103631046410865109661116711168 Figure 4.9 Retrieved SS of vegetation, river, roof, and ground by DART calibration. All images are RGB colour composites. The black pixels indicate no corresponding landcover in that pixel. ...... 65 Figure 4.10 Frequency distribution of SS obtained from DART calibration in (a) blue, (b) green, (c) red, and (d) NIR bands. The cross symbol indicates the SS in iteration 0 of DART calibration in Fig. 8(b). ................................................................................................................................................. 66 Figure 5.1 The geometry explanation of vegetation indices. A vegetation pixel is depicted as (Rnir, Rred), and the isoline of TSAVI (upper) and RVI (lower) go through it. A soil pixel is depicted as (Rnir,soil, Rred,soil), and the soil line goes through it. ........................................................................... 71 Figure 5.2 Concepts of vegetation isolines for various vegetation indices. (a) NDVI; (b) SAVI; (c) TSAVI; (d) PVI. The convergence point of NDVI is located in the original point; the convergence point of SAVI is located in the angle dividing line of the first and third quadrants; the convergence point of TSAVI is located in the soil line; the isolines of PVI is parallel to the soil line. ............... 72 Figure 5.3 Canopy scene created by DART for reflectance simulation. (a) Homogeneous scene; (b) heterogeneous scene. ....................................................................................................................... 73 Figure 5.4 DART simulated BRF in red and NIR spectral bands for homogeneous and heterogeneous canopy as a function of LAI with ALA equal to 40 • and soil reflectance in red band equal to 0.25 when SZA equal to (a) 30 • , (b) 45 • , and (c) 60 • . Similar trends are also found for other ALA and soil optical property values. The dynamic ranges of the homogeneous canopy are slightly larger than the heterogeneous canopy because the low-NIR-reflectance soil occupies a large fraction of the heterogeneous canopy. .......................................................................................................... 74 Figure 5.5 Vegetation isolines with SZA equal to 45 • and ALA equal to 70 • for (a) homogeneous and (b) heterogeneous canopy from nadir observation. In the homogeneous canopy, the highest LAI vegetation isoline intersects with the soil line in the first quarter, whereas in the heterogeneous

)

  Radiative transfer modelRadiative transfer models, also known as physical models, simulate radiation propagation based on physical mechanisms. Radiative transfer models are based on the basic radiative transfer equation, which relates changes in the brightness of radiation along the light path to local absorption and scattering. Radiative transfer models are potentially robust and accurate because they can be applied to three-dimensional realistic structural scenes. In general, simulations of the bottom of atmosphere (BOA) BRFs involve three main radiative transfer components: soil models (e.g., the

Figure 1 . 1

 11 Figure 1.1 Example of changes in the reflectance of red and near-infrared bands according to the view zenith angle changes. Data were simulated using the DART model with homogeneous canopy and LAI equal 0.8.

( 3 )

 3 How to improve the robustness of vegetation indices? Vegetation indices are widely used for vegetation monitor and LAI estimation. The vegetation indices are designed to enhance the spectral information of vegetation. However, vegetation indices are fragile to environmental noise, such as spectral response functions, atmosphere effect, soil noise, and saturation effect. Therefore, studies should be carried out to study the impact of these factors on vegetation indices and propose a corresponding improvement approach to enhance the ability of vegetation indices.

1. 3 . 1 ( 1 )

 311 Study content Simulation of multi-temporal three-dimensional radiative transfer in maize canopies by coupling growth modelsWe propose using the ELSYS model, combined with the maize dynamic growth model, to model the 3D vegetation remote sensing scenes accurately. Maize time-controlled growth parameters are input into the maize growth equation to obtain maize geometric parameters, and the geometric parameters are input into the ELSYS model to obtain a time-series simulation of a threedimensional maize scene with cumulative degree days as a variable is achieved by coupling the maize growth equation. The time-series maize canopy reflectance is simulated using a Radiosity-Graphics based Model (RGM) after parallel acceleration and DART, respectively, and a good agreement is achieved. The simulation results are compared with 1D radiative transfer simulations to analyse the shortcomings of the homogeneous canopy-based assumption.

  ) Vegetation indices study for noise immunity and LAI estimation We propose to study the behaviours of vegetation isoline for LAI saturation studies. The behaviours of vegetation isolines are analysed, and the intersection right-shift phenomenon is proposed. The use of negative soil adjustment factors to mitigate SAVI saturation is deduced from the intersection right-shift phenomenon and validated. We propose to study the impact of environmental factors (i.e. SRF, soil noise, atmosphere) on vegetation indices to improve the robustness of vegetation indices. The effect of spectral response functions and atmosphere on vegetation indices are studied. The spectral correlation coefficient of EVI2 is calibrated for Sentinel-2 and Landsat-8 sensors. The ability of soil noise resistance is compared between single angular index and hotspot signature vegetation indices, together with LAI estimation.

Chapter 1

 1 introduces the background of the study, gives the purpose and significance of the research and summarises the research background. Chapter 2 introduces the physical basis of quantitative remote sensing. Three common vegetation radiative transfer models, SAIL, DART and RGM models, are introduced for subsequent chapters. The simplified assumptions for the canopy of them are demonstrated and the radiative transfer theory. A parallel accelerated RGM version based on OpenACC is proposed, and the resultsshow that acceleration and power increase is possible without loss of quality of the simulated signal.Chapter 3 focuses on the canopy reflectance simulation by coupling the maize growth equation with the cumulative degree day to achieve a three-dimensional scene of maize, and the canopy reflectances of this scene are simulated using RGM and DART for cross-validation. The corresponding homogeneous canopy reflectance is simulated using SAIL to compare the differences between the one-dimensional radiative transfer model and the three-dimensional radiative transfer model.Chapter 4 focuses on the inversion of OPs of leaves by performing mixed pixel decomposition and physical inversion of leaf OPs from the mixed pixels in urban areas. The inversion results are validated using simulated and satellite images; Accuracy assessment and sensitivity analysis are carried out.Chapter 5 evaluate the performance of VIs to estimate the leaf area index. We propose several new hotspot signature vegetation indices and compare their performance with the single angular index for leaf area index estimation. We propose a hypothesis that a negative soil adjustment factor can alleviate the saturation effect of SAVI based on the shift rule of vegetation isolines and soil line. This hypothesis is well observed in both remote sensing and field measured data sets. Chapter 6 analyses the robustness of the VIs as influenced by the spectral response function, atmosphere effect and soil noise. The spectral correlation coefficient of EVI2 is calibrated for Sentinel-2 and Landsat-8 sensors. The effect of the atmosphere on various vegetation indices are studied, and the soil noise resistance of several vegetation indices are compared. Chapter 7 summarises the whole text, points out the innovations and provides an outlook for future research work.

Figure 1 . 2

 12 Figure 1.2 Thesis structure. The whole thesis can be divided into simulation and inversion parts.

2. 1

 1 Radiative transfer models 2.1.1 Scattering by Arbitrary Inclined Leaves (SAIL) model (1) Canopy structureCanopy radiation characteristics depend primarily on canopy structure. Different modelling techniques reflect the structure of the canopy at different scales. In the radiative transfer model, the vegetation canopy can be assumed to be a one-dimensional mixed medium (e.g. atmosphere and water), and the canopy units are randomly distributed, as shown in Figure2.1. The canopy is infinitely uniform horizontally for the one-dimensional case with a limited vertical variation. The substrate of the canopy may be soil or another canopy (e.g. grass), and for simplicity, the canopy may be considered foliage only because the main photosynthetic tissue in the vegetation is leaves.When radiation is transmitted in the vegetation, it is more likely to interact with the leaves to change the radiation characteristics. The average effect of proportional distribution can summarise the true canopy consisting of leaves and branches. The physical properties of the leaf include leaf size, leaf orientation, leaf surface roughness, and leaf optical properties such as reflectivity, transmittance, and absorptivity. Since we pay more attention to the overall nature of the leaves, it is necessary to define some vegetation group characteristic parameters. They are a refined description of the vegetation canopy structure and optical properties and are the results of the statistical average of the distribution of all leaves. The main structural parameters of the canopy include the leaf area index (LAI) and the leaf inclination angle distribution function (LAD).

Figure 2 . 1

 21 Figure 2.1 Schematic diagram of vegetation canopy in one-dimensional mixed media. The leaves are evenly and randomly distributed in the whole space.

  ) SAIL model radiative transfer simulation SAIL (light scattering by arbitrarily inclined leaves) model is extended from Suits model, which assumes that the canopy consists only of vertical and horizontal leaves and parameterizes canopy structure, sun, and observation geometry. The main difference between the SAIL model and the Suit model is that leaves replace the horizontal and vertical projections of the Suit model at arbitrary angles close to reality. Results show that the odd inflexion point disappears compared with the Suits model. The input parameters of the SAIL model include optical properties of leaves (i.e., reflectance and transmittance) and soil (i.e., reflectance), view geometry (sun zenith angle, view zenith angle, and relative azimuth angle between sun azimuth angle and view azimuth angle), Leaf area index (LAI), the average leaf angle (ALA). The output parameter is the canopy reflectance. The computer simulation code is open-source and provided by several language versions on the website (http://teledetection.ipgp.jussieu.fr/prosail/).The continuous vegetation radiative transfer equation is a physical condition describing the value of the radiation radiance in any direction at any point within the vegetation layer that should be satisfied, together with the boundary conditions. This differential-integral equation is in principle solvable, requiring that the solution of such an equation must satisfy: (i) the determination of the phase function; and (ii) the determination of the boundary conditions. So far, it has not been possible to find a strictly analytical solution of the RT equation, only a variety of approximate solutions, each of which has its own adaptation conditions. SAIL (Scattering by Arbitrary Inclined Leaves) is one of the earliest canopy reflectance models to simulate the bidirectional reflectance factor of turbid medium plant canopies by solving the scattering and absorption of four upward / downward radiative fluxes. SAIL actually provides all four-stream optical properties (nine in total) of the canopy layer at the output Figure2.2.

  ) where z is the radiation transfer direction, K and k represent the extinction coefficient of light propagating along with the incident and observed radiation directions, respectively. ES is the topdown transmitted direct radiation flux density, E0 is the radiation flux density in the direction of observation, E + and E -represent upstream and downstream radiation flux densities, respectively; a is the extinction coefficient, b is the backscattering coefficient, c' is the scattering coefficient for forwarding scattering direct radiation and c is the scattering coefficient for backward scattering direct radiation. u, v and w represent the scattering coefficients of the radiance transmitted from E + , E -and ES to the observed direction, respectively.

Figure 2 . 2

 22 Figure 2.2 Flux tracking theory of SAIL model. It is based on the four radiative fluxes method.

  ) DART radiative transfer simulation DART consists of four main parametric modules: (1) the direction module, where the direction of radiation propagation in spherical space is pre-calculated; (2) the phase module, which precalculates the OPs of facet and voxel of the existing earth surface and atmosphere (e.g., simulating the phase function of vegetation under turbid conditions) in order to minimise the computational time required to simulate radiative transfer; (3) the maket module, which simulates the optical and thermal properties of 3D buildings and scenes; (4) the DART model, which simulates the radiative transfer.

Figure 2 . 3

 23 Figure 2.3 Iterative diagram of the tree structure. The whole scene comprises trees and soil, and trees comprise trunk and crown.

Figure 2 . 4

 24 Figure 2.4 Tree structure diagram created using ELSYS.

( 2 )

 2 Radiosity-graphics based model radiative transfer simulation Radiosity has been widely used in reflectance simulations of the 3D real scenes, and its theoretical core is the radiosity equation. The radiosity equation describes the radiative energy balance within a scene containing N discrete facets, assuming that all facets comprising the scene are Lambertian and have specific OPs (reflection and transmission). The radiant flux density (radiosity) leaving a facet consists of three components: the reflected flux density, the transmitted flux density and the emitted flux density obtained after being irradiated by the sun (Figure 2.5).

Figure 2 . 5

 25 Figure 2.5 Schematic diagram of the radiosity model principle. The green arrow represents the incident energy, and the yellow arrow represents the outgoing energy.

Figure 2 . 6

 26 Figure 2.6 Modelling of a real structural mono-tree scene using ELSYS. (a) Real structure single tree; (b) ELSYS modelling of the single tree.

Figure 2 . 7

 27 Figure 2.7 Comparison of the simulated BRF of real structural single tree scene using RGM before and after acceleration. (a) Blue band in the principal plane direction (PP); (b) Green band in the PP direction; (c) Red band in the PP direction; (d) NIR band in the PP direction; (e) Blue band in the cross-principal plane (CP) direction; (f) Green band in the CP direction; (g) Red band in the CP direction; (h) NIR band in the CP direction

Figure 2 .

 2 Figure 2.8 further analyses the simulated reflectance difference before and after the parallel computing improvements. Scatter plots were drawn, and coefficient of determination (R 2 ) and RootMean Square Error (RMSE) were calculated using the reflectance of the pre-improved simulation as the horizontal axis and the post-improved simulation as the vertical axis. The results indicate that the differences were minor.

Figure 2 . 8

 28 Figure 2.8 Accuracy validation of post-accelerate RGM model for real structural single tree scene. (a) Blue band; (b) Green band; (c) Red band; (d) NIR band

Figure 2 . 9

 29 Figure 2.9 Simplified heterogeneous canopy scene. Trees are randomly distributed in the scene, and small triangles fill the crown with uniform random distribution.

Figure 2 .

 2 Figure 2.11 further analyses the differences in simulated reflectance between the RGM and DART. The scatter plot was plotted, and R 2 and RMSE were calculated using the DART simulated reflectance as the horizontal axis and the improved RGM simulated reflectance as the vertical axis. The correlation coefficients of canopy BRF values in the visible and NIR bands were 0.38, 0.45, 0.83 and 0.96, respectively, with RMSEs of 0.0033, 0.0050, 0.0041 and 0.0141, respectively. The results showed that the simulated reflectances between the RGM and DART models were in good agreement.

Figure 2 .

 2 Figure 2.10 Simulated BRF of RGM and DART model with simplified heterogeneous canopy. (a) Principal plane (PP) direction in the blue band; (b) PP direction in the green band; (c) PP direction in the red band; (d) PP direction in NIR band; (e) Cross-principal Plane (CP) direction in the blue band; (f) CP direction in the green band; (g) CP direction in the red band; (h) CP direction in NIR band.

Figure 2 . 11

 211 Figure 2.11 Accuracy validation of the RGM model for simplified heterogeneous canopy. (a) Blue band; (b) Green band; (c) Red band; (d) NIR band.

  ) Take RGM as an example; it shows that the efficiency increase is possible without losing the quality of the simulated signal for radiative transfer simulation. Compared to the pre-acceleration RGM model, the parallel accelerated RGM has an R 2 of up to 0.97 in the solar PP and CP in a realistic structural single tree scene. Also, the parallel accelerated RGM model has an R 2 of up to 0.96 in the NIR band compared to the DART model in a simplified heterogeneous forest scene. The RGM model takes 91.6% less time to simulate and can tolerate a minimum of 840% increase in the maximum number of facets.The maize growth equation is coupled with a 3D realistic structural maize scene simulation based on the ELSYS platform. A static 3D maize modelling model developed based on the ELSYS platform was driven by coupling the maize growth equations presented in the literature[START_REF] López-Lozano | Sensitivity of gap fraction to maize architectural characteristics based on 4D model simulations[END_REF][START_REF] Espana | Modeling maize canopy 3D architecture: application to reflectance simulation[END_REF].The results were based on a 4D real structural maize growth equation driven by cumulative degree days developed on the ELSYS platform. The reflectance simulations are carried out using RGM and DART and cross-validated, respectively. Section 3.1 illustrates the three-dimensional maize growth equation and the way to cooperate it with the static maize modelling model. Section 3.2 shows the maize canopy's continuous-time phase reflectance simulation results using the DART and RGM model based on this real structural maize canopy scene and the homogeneous assumption canopy scene. Differences in FC between

  of the top stem, and the radius of the bottom stem. Leaf width and length, leaf attachment point height, leaf breakpoint length and leaf inclination angle were all required. Individual maize plants were then modelled, and the entire maize scene was reproduced by replicating individual maize plants at specified locations according to plant density and row spacing (Figure 3.1).

Figure 3 . 1

 31 Figure 3.1 Flowchart of maize 3D scene modelling and reflectance simulation. Light ellipses indicate input parameters, dark ellipses indicate output parameters, and white ellipses indicate intermediate parameters. N max, Td, Smax, θmax, Hmax were entered into the maize growth equation to calculate the geometric parameters of a single maize plant at a given Td, including leaf order (N), base stem width (RN), leaf basal width (WN), leaf length (LN), leaf vertical height (HN), leaf inclination angle (θN) and leaf azimuth angle (φN). Leaf breakpoint location (Pbreak) was then calculated from N and LN. Pbreak, N, RN, WN, LN, HN, θN and φN were input into the ELSYS model to simulate the 3-dimensional structure of the single maize. Finally, the generated single maize structures, D, drows, SZA, Sun Azimuth Angle (SAA), View Zenith Angle (VZA) and View Azimuth

Figure 3 .

 3 2 shows the growth stages of maize at six equally spaced consecutive time series. The maize canopy was expressed in its true LAI and FC for a qualitative description of maize growth. Trends in FC and true LAI with Nmax, Td, Smax, θmax, Hmax 3D Maize Model D, drows, SZA, SAA, VZA, VAA were shown in Figure 3.2.g. True LAI and FC rose significantly between Td = 100 and Td = 1100, a period when maize was growing. After a long period of stability (Td = 1150 to Td = 1600), maize declined (old leaves began to fall off) after Td = 1600.

Figure 3 . 2

 32 Figure 3.2 Growth scenes of maize at equally spaced continuous time series. (a) Degree days (Td) = 300; (b) Td = 600; (c) Td = 900; (d) Td = 1200; (e) Td = 1500; (f) Td = 1800; (g) Variation of true LAI and Foliage Coverage (FC) with cumulative degree day. True LAI is plotted on the left y-axis, while FC is plotted on the right y-axis. FC shows a significant difference between the homogeneous and heterogeneous canopy, especially with the increase of LAI. Maize has gone through periods from emergence to male anthesis. Withering is simulated by the gradual loss of the bottom leaves. The DART model displays the 3D scenes.

Figure 3 . 3

 33 Figure 3.3 The simulated Bi-directional Reflectance Factors (BRFs) in (a) red and (b) NIR spectral bands in the hotspot (HS), the darkspot (DS), and the nadir for homogeneous and heterogeneous canopy as a function of Td. SAIL simulates the homogeneous canopy reflectance, and RGM simulates the heterogeneous canopy reflectance. DART simulated homogeneous (DART homo) and heterogeneous (DART hete) data are used to validate the results of SAIL and RGM separately. Overestimation of the homogeneous canopy's simulated reflectance in the nadir direction in the NIR band is observed.

Figure 3 .

 3 Figure 3.3 shows the comparison of the BRFs for the homogeneous and heterogeneous (or 3D

Chapter 4

 4 Inversion of leaf optical properties in urban areas at the sub-pixel scale The OPs of the leaves have a positive effect on the analysis of their biochemical parameters and the monitoring of vegetation growth status. DART calibration has been successfully used to simulate the short-wave net radiation in the cities. Here we show another potential usage of DART calibration to invert the OPs of leaves in the cities. Based on the DART calibration, a mixed pixel decomposition and physical inversion of leaf OPs are performed to invert the OPs of leaves from mixed pixels in the urban areas. The inversion results are validated using simulated cases (with or without noise) and satellite images, and accuracy assessment and sensitivity analyses are performed. Section 4.1 examine the DART calibration based on the DART simulated images. This first experiment is an ideal case where the only unknown parameters are the optical properties of landcover, and we assume that we know precisely other remote sensing parameters. A second analysis focuses on the common case for which certain acquisition parameters are known with a certain level of uncertainty. Sensitivity analysis allows us to identify the parameters whose uncertainty will greatly impact the accuracy of the optical properties inverted by the DART calibration. Section 4.2 examine the DART calibration based on one sentinel-2 satellite image of the Basel city, in which a detailed 3D database describing the main elements (buildings, road infrastructures, vegetation) exists. The image obtained after DART calibration is very similar to the real image. The optical properties inverted for the different landcover of the scene show a very realistic heterogeneity,

Figure 4 . 1

 41 Figure 4.1 Schematic urban scene: ground, water, shrubs, trees, and buildings with the flat, triangular prism and frustum of a square pyramid roof. (a) 3D view. (b) Top view.

Figure 4 . 2

 42 Figure 4.2 (a) DART pseudo satellite image of the 3-D scene. It has some spectral confusions (e.g., the different reflectance for nine roofs and the same reflectance of three roofs as ground reflectance). (b) DART image simulated with constant SS per urban component. (c) DART image after the calibration procedure. All images are RGB colour composites.

Figure 4 .

 4 Figure 4.2 shows the pseudo-satellite image and DART simulated image before and after calibration. There was spectral confusion in the pseudo-satellite image, including roofs with different OP and one type of roof with the same OP as the ground. In the initial image, all roofs shared the same reflectance. As expected, the calibrated images were visually very similar to the pseudo-satellite images.Figure 4.3 shows the evolution of the mean relative error for each component |ε| ��� component with iterations in the green band. The mean relative error at iteration 0 was

Figure 4 .

 4 Figure 4.4 shows the spatial image and violin plot of the relative error for each component in the calibrated image. The blank areas in the image indicated that there was no corresponding component in that pixel. The relative error of each component was calculated pixel-wise. Besides, the violin plot vertical axis coordinates range was restricted to the range of relative errors. Similar phenomena were observed for other bands. The calibration accuracy of the pure pixels appeared to

  satellite observation factors of the DART input parameters for the initial image were not precisely known. Their imprecision affected the inversion of the OPs. The five considered influencing factors correspond to two broad categories: (1) parameters at image acquisition: the SZA, the spatial image resolution and the MTF; (2) the accuracy of the image calibration (shifting along the x and y axes)

Figure 4 . 3

 43 Figure 4.3 Evolutions of the mean absolute value of the relative error of the component SS with the iterations in the green band. The mean absolute value of the relative error at iteration 0 is due to the selected initial SS value. Relative errors are calculated pixel-wise for each component, and a mean statistical value of the absolute value of the relative error is calculated. Similar trends are found in the other bands.

Figure 4 . 4

 44 Figure 4.4 Final calibration results of the error analysis for ground (a and b), roof (c and d), water (e and f), shrub (g and h), and trees (i and j) in the green band in the spatial and frequency domain.Relative errors are calculated pixel-wise. Blank in the image means there is no corresponding component in that pixel. Violin range is limited to within the range of the relative error. A similar phenomenon was also observed for other bands.

Figure 4 .

 4 3 shows that the calibration procedure had almost converged when the fifth iteration was performed. Therefore, a constant maximum iteration value of 8 was set for all sensitivity analyses. The eFAST sensitivity analysis was used to analyse SZA, satellite image resolution, pixel shifting, scene modelling inaccuracies and MTF. The eFAST methods were used for sensitivity analyses based on 420 cases. Sensitivity indices were calculated on a per-component per-band basis and superimposed over all bands (Figure 4.5).

Figure 4 .

 4 Figure 4.6 shows the mean relative error for each component for all 420 cases. The ground and roof performed well in all bands: the range was 0.015-1.572 for the ground and 0.024-1.7662 for the roof. The other components were less accurate than the ground and roof, especially those with low reflectance in the corresponding bands. If the OP of the components was low, significant errors tended to occur. For example, for vegetation (i.e., shrubs, trees), the range of the mean relative error in the green band was 0.159-2.709. Similarly, the mean relative error of the water increased sharply in SWIR1 and SWIR2, which had very low values of OPs. The range of the mean relative error for water in all bands was 0.053-15.683.

2 and Figure 4 . 3 )

 43 show that the DART calibration had an excellent overall accuracy under ideal experimental conditions. However, the inversion accuracy of the vegetation's OPs was significantly lower than the other urban components. Three explanations were proposed here. (1) As a voxel composed of leaf elements, vegetation could not be modelled as a simple surface like the other urban components.[START_REF] Li | Spatio-temporal Characteristics of Area Coverage and Observation Geometry of the MISR Land-surface BRF Product: A Case Study of the Central Part of Northeast Asia[END_REF] The OPs of the leaf include transmittance and reflectance. Here, they were represented by the single scattered albedo of the leaf with a pre-defined spectral ratio between its reflectance and transmittance.[START_REF] Zhen | Potentials and Limits of Vegetation Indices With BRDF Signatures for Soil-Noise Resistance and Estimation of Leaf Area Index[END_REF] The values of the OPs of the leaf were low in the spectral bands under consideration.

Figure 4 . 5

 45 Figure 4.5 Sensitivity analysis results. The eFAST method is used for sensitivity analysis based on 420 cases. Sensitivity indices are calculated per component per band and stacked over all bands. All subplots share the same x-axis. In descending order of magnitude, the effect magnitudes are SZA, satellite spatial resolution, pixel-shift, landscape model inaccuracy, and MTF.

Figure 4 . 6

 46 Figure 4.6 Mean absolute value of the relative errors on the retrieved SS of scene components for the blue, green, red, NIR, SWIR1, and SWIR2 bands. The relative errors were calculated for each component in terms of pixels. The absolute value of the relative error was the average of all 420 cases. The violin range was constrained to the boundary of the mean relative error.

  The location and geometry of all urban components, including the local digital elevation model (DEM), were taken from the local urban database. The location, height and canopy dimensions of the trees were set from the measurements. DART simulated the trunks as regular octahedrons and the crowns as ellipsoidal bodies filled with fragmented triangles. Only the necessary information about the urban database was given here. A more detailed description can be found in the literature[START_REF] Landier | Modélisation 3D du bilan radiatif des milieux urbains par inversion d'images satellites en cartes de réflectance et de température des matériaux urbains[END_REF].

Figure 4 . 7

 47 Figure 4.7 Basel 3-D scene. (a) Top view. (b) Side view of the red square in (a). Buildings (roofs: dark, walls: light grey), trees (green), river (blue), and ground (yellow).

Figure 4 . 8

 48 Figure 4.8 (a) PlanetScope image in the Basel area. (b) DART image simulated with constant SS per urban component. (c) DART image after the calibration procedure. All images are RGB colour composites.

Figure 4 . 9

 49 Figure 4.9 Retrieved SS of vegetation, river, roof, and ground by DART calibration. All images are RGB colour composites. The black pixels indicate no corresponding landcover in that pixel.

Figure 4 .

 4 Figure 4.8 shows the satellite image and DART simulated image before and after DART calibration. Before calibration, the difference between the DART simulated image and the satellite image was significant. However, after calibration, the differences became minimal. The distribution of OPs in space and frequency is depicted in Figure 4.9 and Figure 4.10, where the cross symbols indicate the initial OPs at iteration 0 of the DART calibration in Figure 4.8.b. Reasons for the bimodal problem of the water was due to the lack of bridges in the urban 3D database. Therefore,

Figure 4 . 10

 410 Figure 4.10 Frequency distribution of SS obtained from DART calibration in (a) blue, (b) green, (c) red, and (d) NIR bands. The cross symbol indicates the SS in iteration 0 of DART calibration in Fig. 8(b).

3 )

 3 Compared with NDVI, SAVI shifts the convergent point from the original point to the point (-XSAVI / 2, -XSAVI / 2), and the (1+ XSAVI) is multiplied to maintain the bounded conditions. The index isolines of SAVI is depicted in Figure5.2.b.

Figure 5 . 1

 51 Figure 5.1 The geometry explanation of vegetation indices. A vegetation pixel is depicted as (Rnir, Rred), and the isoline of TSAVI (upper) and RVI (lower) go through it. A soil pixel is depicted as (Rnir,soil, Rred,soil), and the soil line goes through it.

  (v) PVIPVI is the distance between the vegetation reflectance point and the soil line based on the distance equation from point to line. The index isolines of PVI is depicted in Figure5.2.d.

Figure 5 . 2

 52 Figure 5.2 Concepts of vegetation isolines for various vegetation indices. (a) NDVI; (b) SAVI; (c) TSAVI; (d) PVI. The convergence point of NDVI is located in the original point; the convergence point of SAVI is located in the angle dividing line of the first and third quadrants; the convergence point of TSAVI is located in the soil line; the isolines of PVI is parallel to the soil line.

Figure 5 . 3 -

 53 Figure 5.3 Canopy scene created by DART for reflectance simulation. (a) Homogeneous scene; (b) heterogeneous scene.

Figure 5 . 4

 54 Figure 5.4 DART simulated BRF in red and NIR spectral bands for homogeneous and heterogeneous canopy as a function of LAI with ALA equal to 40 • and soil reflectance in red band equal to 0.25 when SZA equal to (a) 30 • , (b) 45 • , and (c) 60 • . Similar trends are also found for other ALA and soil optical property values. The dynamic ranges of the homogeneous canopy are slightly larger than the heterogeneous canopy because the low-NIR-reflectance soil occupies a large fraction of the heterogeneous canopy.

Figure 5 .

 5 6 shows the scatter plots of the FC variation with the abscissas of the intersection of the vegetation isoline and the soil line for the corresponding LAI.

Figure 5 . 5

 55 Figure 5.5 Vegetation isolines with SZA equal to 45 • and ALA equal to 70 • for (a) homogeneous and (b) heterogeneous canopy from nadir observation. In the homogeneous canopy, the highest LAI vegetation isoline intersects with the soil line in the first quarter, whereas in the heterogeneous canopy, the lowest LAI vegetation isoline intersects with the soil line in the first quarter. Besides, the dynamic slope range in the heterogeneous canopy is far smaller than that in the homogeneous canopy with the increase of LAI. The vegetation reflectance points in the same isoline are due to the variation of soil reflectance. For homogeneous canopies, the slope of the vegetation isoline increased with increasing FC (Figure 5.5.a), and the abscissa of the intersection of the vegetation isoline and the soil line gradually increased and converged to an immobile point (Figure 5.6.a). This phenomenon may help explain

Figure 5 . 6

 56 Figure 5.6 Relationship between abscissae of crosspoints and foliage cover using (a) homogeneous and (b) heterogeneous simulated data set when SZA equal to 45 • . The symbols represent the ALA, and the colours represent the LAI. In the homogeneous scene, the right shift phenomenon of crosspoints is obvious. Besides, the abscissae of crosspoints in the homogeneous canopy are linearly positive with foliage cover, and the value of abscissae can vary from negative to positive. However, in the heterogeneous scenes, the right shift phenomenon of crosspoints is not obvious, and the value of abscissae can even be positive when LAI is low.

  decreasing function, where a, b, x and y represent the reflectance in the NIR band, the reflectance in the red band, the soil adjustment factor XSAVI and SAVI. Obviously, some preconditions can be derived, such as y > 0, 0 < b ≤ a < 1 and 0 ≤ x ≤ 1. In addition, usually for vegetation, the sum of the reflectance values in the NIR band (b) and the red band (a) in the nadir direction observation does not exceed 1. Therefore, a + b < 1. By assuming that x1 < x2, we have ( ) (

VIs calculated using S2 image from 4

 4 February 2019 to 10 November 2020 and L8 Tier 1 imagery from 16 April 2013 to 21 October 2020 for multispectral surface reflectance data. The use of multiple data sources helped to increase the sampling frequency of the time series. Differences in VIs due to different sensors were ignored because of the similarity in resolution between S2 (10 m) and L8 (30 m) satellites.

Figure 5 . 7

 57 Figure 5.7 Study area with the base map of the MODIS classification product (MCD12Q1V6). It has five types of land cover: water body (blue), evergreen broadleaf vegetation (dark green), deciduous broadleaf vegetation (yellow), annual grass vegetation (yellow-green), and urban and built-up lands (dark grey).

Figure 5 . 8

 58 Figure 5.8 Decomposition analysis for the long-term LAI in the Apiacás (2013-2020). (a) LAI is decomposed into (b) trend, (c) seasonal, and (d) residual. The full-year LAI is separated into four quarters. Then, the mean value of each quarter is calculated and is used for time series analysis.

Figure 5 . 9

 59 Figure 5.9 Stationarity analysis of MODIS LAI residual by rolling mean and standard deviation over the whole observation period.

Figure 5 . 10

 510 Figure 5.10 Time series of VIs (left axis) and LAI (right axis). The solid circles highlight a period when the LAI has a trend opposite to that of all VIs. The slash circles highlight a period when the trend of LAI is consistent with NDVI while being opposite to other VIs. The dashed circles highlight a period when the trend of LAI is consistent with all VIs.

Figure 5 . 11

 511 Figure 5.11 Scatter plot between VIs and (a) MODIS and (b) field-measured LAI. The solid lines represent the trendlines of the linear regression model.

2 (

 2 Figure 5.12.b), and R 2 and P fluctuate considerably in this region. The maximum R 2 (0.6417) and minimum P (0.0863) occur at XSAVI = -0.183. As XSAVI increases from -0.184 to 1, SAVI performance becomes progressively worse, with R 2 decreasing and P increasing.

Figure 5 . 12

 512 Figure 5.12 Linear regression between SAVI and (a) MODIS and (b) field-measured LAI with varying XSAVI. The optimal results are observed in the negative XSAVI region.

Figure 5 . 13

 513 Figure 5.13 Map of the Canadian boreal forest shows the BOREAS project's two specific study areas: SSA and NSA. Six sites are investigated, including old jack pine (SOJP), young jack pine (SYJP), old black spruce (SOBS), old aspen (SOA) in the SSA and young jack pine (NYJP), and old black spruce (NOBS) in the NSA.

Figure 5 .

 5 14 shows the variation curve of BRF values with LAI. From the figure, significant fluctuations in reflectance at different sampling sites exist because of the variations in the leaves OPs and the spatial distribution of the leaves of the tree species.

Figure 5 . 14

 514 Figure 5.14 Fitting reflectance in the red and NIR spectral bands from BOREAS study areas as a function of LAI when SZA is equal to (a) 30 • , (b) 45 • , and (c) 60 • .

  Simulated and measured data were used to test the ability of all VIs considered to invert the LAI. The effectiveness of all VIs for LAI inversion was assessed by regression analysis. Regression parameters for homogeneous scenes, heterogeneous scenes and field measurements are presented in

Figure 5 . 15 ,

 515 Figure 5.17. Four fit models, unsaturation (linear), saturation (exponential), and two common used

Figure 5 . 15 Figure 5 . 16

 515516 Figure 5.15 Standard goodness-of-fit statistics: R 2 and RMSE of regression between the ten VIs and LAI of the homogeneous simulated data, for four fitting models when SZA equal to (a) 30 • , (b)

Figure 5 . 17

 517 Figure 5.17 Standard goodness-of-fit statistics: R 2 and RMSE of regression between the ten VIs and LAI of field measurements, for four fitting models when SZA equal to (a) 30 • , (b) 45 • , and (c)

Figure 5 . 18

 518 Figure 5.18 Performances of VIs by considering 14 groups of random errors in the homogeneous simulated data using an exponential model when SZA is equal to (a) 30 • , (b) 45 • , and (c) 60 • . Multispectral VIs are drawing as dotted lines, VIs57.5 are drawing as dash-dot lines, and hotspotsignature VIs are drawing as solid lines. The good LAI estimation performance of multispectral VIs can be seen.

Figure 5 . 19

 519 Figure 5.19 Performances of VIs by considering 14 groups of random errors in heterogeneous simulated data using exponential model when SZA equal to (a) 30 • , (b) 45 • , and (c) 60 • . Multispectral VIs are drawing as dotted lines, VIs57.5 are drawing as dash-dot lines, and hotspot-signature VIs are drawing as solid lines. The good LAI estimation performance of hotspot-signature VIs can be seen.

Figure 5 . 20

 520 Figure 5.20 Performances of VIs by considering 14 groups of random errors in field measurements using exponential model when SZA equal to (a) 30 • , (b) 45 • , and (c) 60 • . Multispectral VIs are drawing as dotted lines, VIs57.5 are drawing as dash-dot lines, and hotspot-signature VIs are drawing as solid lines. The good LAI estimation performance of hotspot-signature VIs can be seen.

Figure 6 . 1

 61 Figure 6.1 Values of SNR for all VIs considered when SZA equal to (a) 30 • , (b) 45 • , and (c) 60 • .The SNR values corresponding to ALA in the nine groups are averaged and presented here. The four hotspot-signature VIs show good performance for the homogeneous and heterogeneous simulated data in most cases, followed by VIs57.5.The first evaluation indicator was defined by Leprieur and Roujean[START_REF] Leroy | Sun and view angle corrections on reflectances derived from NOAA/AVHRR data[END_REF] based on SNR as ( ) ( )

  [START_REF] Zhu | A regional maize yield hierarchical linear model combining landsat 8 vegetative indices and meteorological data: Case study in jilin province[END_REF] showed good resistance to soil noise for both canopies. Furthermore, in moderate LAI canopies with prevalent soil background and pronounced multiple scattering effects, it was observed that NDVI, NDVI57.5, NHVI and NHVI2 barely attenuated soil effects which presented relatively low SNR and the largest T values. The literature[START_REF] Rondeaux | Optimization of soil-adjusted vegetation indices[END_REF][START_REF] Huete | A soil-adjusted vegetation index (SAVI)[END_REF][START_REF] Qi | A modified soil adjusted vegetation index[END_REF][START_REF] Gilabert | A generalized soil-adjusted vegetation index[END_REF] also reported similar findings. For dense homogeneous canopies (where LAI was large), all VIs showed the lowest T values because the contribution of soil was less significant, as reported in the literature[START_REF] Gilabert | A generalized soil-adjusted vegetation index[END_REF]. Thus, HSVI emerged as the VI with the best resistance to soil noise.

Figure 6 . 2

 62 Figure 6.2 Efficiency of the different VIs as measured by means of T as a function of LAI for the homogeneous and the heterogeneous canopy when SZA equal to (a) 30 • , (b) 45 • , and (c) 60 • . The T values corresponding to ALA in the nine groups are averaged and presented here. The excellent performance of HSVI and HEVI2 can be observed overall.

Figure 6 . 3

 63 Figure 6.3 Relative error of atmospheric effects (Eatmosphere) for EVI, EVI2, NDVI, and SAVI. Surface reflectances calculated from L8 and S2 were considered true values, while the VIs calculated from TOA reflectance were considered measured values. The Eatmosphere,VI (x,y) was calculated pixel-wise, and the Eatmosphere,VI,i values for all vegetation pixels in a scene image i were averaged and plotted as a box plot. The outliers in the box plot were hidden for better visualization.

  L8 and S2 differed significantly from the Spectral Response Functions (SRF) of the Moderate Resolution Imaging Spectroradiometer (MODIS) satellites. Figure 6.4 shows the SRFs for the MODIS, S2 and L8 sensors in the blue, red and NIR bands. The differences in SRFs between sensors emphasised the importance of calibrating the correlation coefficient between the blue and red bands on sensors other than MODIS.

Figure 6 . 4

 64 Figure 6.4 The spectral response functions (SRFs) of MODIS, Sentinel-2, and Landsat-8 in blue, red, and NIR bands. Sentinel-2 and Landsat-8 SRFs are referred from the USGS Spectral Library Version 7 [234], and MODIS SRFs are referred from the European Organization for Meteorological Satellites (https://nwp-saf.eumetsat.int/downloads/rtcoef_rttov13/ir_srf/rtcoef_eos_1_modis-shifted_srf.html).

Figure 6 . 5

 65 Figure 6.5 The relative error of EVI2 (EEVI2) with EVI calculated from Landsat-8 and Sentinel-2 surface reflectance products. EVI is considered as "true" data, and EVI2 is considered as "measured" data. Values are calculated pixel-wise and aggregated to a mean value over the image. The figure is drawn using RainCloudPlots [237, 238], comprising violin (top), box (middle), and strip (bottom) plots.

( 3 )

 3 measurement and airborne spectrometers with wavelengths covered from the ultraviolet to the farinfrared (0.2 μm -200 μm). The spectres in the spectral library came from different components of plants or patches of vegetation, which included many types and different sources. Aerial spectrometer measurements were included in wooded vegetation plots where trees were too tall to be measured with a field spectrometer. For the spectra in the spectral library, four different spectrometer models were used. (1) Beckman 5270 (0.2 μm-3 μm), (2) standard, high resolution and high resolution next-generation models of the ASD field portable spectrometer (0.35 μm-2.5 μm), Nicolet Fourier Transform Infrared (FTIR) interferometer spectrometer (1.12 μm-216 μm), and (4) NASA Airborne Visible -Infrared Imaging Spectrometer Airborne Visible / InfraRed Imaging Spectrometer (AVIRIS) (0.37 μm-2.5 μm).

Figure 6 . 6

 66 Figure 6.6 The linear regression between reflectances in a blue and red band based on USGS spectral library. The shadow area indicates the 95% confidence interval for the regression estimate.

Figure 6 . 7

 67 Figure 6.7 Distribution of the linear regression coefficients (i.e. slope a and interception b) in Eq.(6.5) from (a) Sentinel-2; (b) Landsat-8. Coefficients are calculated pixel-wise, and the mean value is aggregated over the whole image. A white circle highlights the spectral correlation coefficient c derived from the field measured data, and a white star highlights the c of the original EVI2 [106]. The range of a and b is limited at [0.5, 2.4] (x-axis) and [-0.04 0.04] (y-axis) for a better visual effect.

Figure 6 .

 6 Figure 6.6 shows the linear regression coefficients a and b of equations (6.5) derived from the USGS Spectral Library. The overall regression was good, with correlation coefficients (R) equal to 0.765 (S2) and 0.747 (L8) and p-values (P) equal to 0.000 for both S2 and L8, indicating an excellent linear relationship. Additionally, the slope (a S2) and the intercept value (bS2) for S2 were 1.262 and 0.021; the slope (aL8) and the intercept value (bL8) for L8 were 1.269 and 0.028. These values were then used to validate the fit of the regression equation for cL8 and cS2.

Figure 6 .

 6 Figure 6.7 shows the distribution of linear regression coefficients a and b derived from remote

Figure 6 .

 6 Figure 6.8 shows the improvements in EVI2S2 and EVI2L8 after calibration, where EVI is considered the "true" value and EVI2S2 or EVI2L8 is considered the "measured" value. The relative error of EVI2 (EEVI2) was used to assess the robustness of EVI2 on the S2 and L8 sensors. Before calibration, the R 2 and RMSE of EVI2 were 0.87 and 0.08 (Figure6.8.a). After calibration, the R 2 and RMSE for EVI2S2 changed to 0.92 and 0.06 (Figure6.8.b). Similar excellent performance can also be observed for EVI2L8. Before calibration, the R 2 and RMSE of EVI2 were 0.92 and 0.06 (Figure6.8.c). After calibration, the R 2 and RMSE of EVI2S2 changed to 0.95 and 0.05 (Figure6.8.d).

Figure 6 . 8 .

 68 Figure 6.8.e shows the EEVI2 before and after calibration, and Figure 6.8 depicts the distribution of the statistics for EVI2. The error distribution for EVI2, EVI2L8 and EVI2S2 essentially satisfied a normal distribution. In addition, the expected EEVI2 values of EVI2L8 and EVI2S2 were closer to zero than those of EVI.

Figure 6 . 8

 68 Figure 6.8 Validation of EVI2S2 and EVI2L8 using the USGS Spectral Library. (a) EVI and EVI2 with S2 spectral respond function; (b) EVI and EVI2S2 with S2 spectral respond function; (c) EVI and EVI2 with L8 spectral respond function; (d) EVI and EVI2L8 with L8 spectral respond function; (e) Relative error of EVI2 (EEVI2) with EVI being considered as "true" value and EVI2S2 or EVI2L8 being considered as "measured" value. Values are calculated per spectral curve.

( 1 )

 1 Coupling maize growth equations adds the ability to simulate continuous-time phase data with RGM and DART models. The development of three-dimensional radiative transfer models and the popularity of highresolution imagery emphasises the importance of realistic scene modelling. Based on maize's threedimensional static modelling model, growth control factors with the cumulative degree day as a control variable are introduced. The period from seedling emergence to male flowering is simulated, development rules describing the structure of the entire plant growing season are added, and the time-series simulation of the three-dimensional radiative transfer model is achieved. The fulfilment of the growth model fills in the knowledge of the spatial and temporal variation of critical parameters of the ground, simulates changes in the spatial distribution of crops and enables continuous timephased data simulations.

( 1 )( 2 )( 3 )

 123 fixed observation direction. The high accuracy of LAI estimates can help increase the value of remote sensing technology for current hotspot issues such as large-scale agroforestry monitoring, global food security, global carbon cycle, global surface change monitoring, primary productivity analysis and global food crop security. This work introduced the plant growth model into the DART model. There are also interesting points where DART products can be introduced into the growth model, and they are listed below: Radiative budget (e.g., Absorbed Photosynthetically Active Radiation: APAR): The vegetation growth stage greatly depends on the time integral of APAR. On the other hand, APAR depends on the vegetation 3D architecture. Therefore, a DART simulation of the plant APAR, using meteorological data, indicates the energy received by a plant and consequently can influence the vegetation growth model. Remote sensing data (e.g., directional and albedo images): Comparing these data with actual remote sensing data can inform the vegetation radiative budget because DART offers the advantage of simulating simultaneously remote sensing and radiative budget products. Therefore, DART remote sensing data can be used to drive growth vegetation models. In addition, DART remote sensing data can be used to train machine learning algorithms to detect vegetation species (vegetation classification) and bio-physical parameters inversion (i.e., tree height), which are essential parameters in plant growth models. Fluorescence data: Sun-induced chlorophyll fluorescence (SIF) offers a physiologically based GPP proxy, potentially avoiding the need for LUE parameterizations. In fact, many studies have already shown an empirical linear scaling between SIF and GPP across various vegetation types under different environmental conditions, revealing the potential of SIF to monitor GPP across all spatial scales[239]. Plant SIF depends on the plant growth stage and consequently on the plant 3D architecture. It illustrates the link between DART SIF products, plant architecture and fluorescence parameters of the plant foliar elements.(4) Mock-up data. DART simulation of the 3D architecture of plants should accurately represent the 3D architecture of plants simulated by growth models. DART can work with different precision levels for representing this 3D architecture. characters of the L-system.

A

  simple example code of describing the several periods of a flowering plant development (Figure A.2) from the literature [159] is listed below: animation(1,1, 7,20, 1,7, 60,170, 3) /* axiom */ animation(kx,ky, n,m, h1,h2, delta1,delta2, nsteps) --> let(i:= 0) setname ('an') nsteps*( ky*( [kx*( [object(n,m, h1+(h2-h1)*real(i)/(nsteps-1), delta1+(delta2-delta1)*real(i)/(nsteps-1))] 90+ val(2*n) f 90-)] val(2*n) f) flush inc(i)) object(n,m,h,delta) --> trunk(n,m) layer(n,h,delta) trunk(n,m) --> m*(layer(n,1,90) 90^ f 90&) layer(n,h,delta) --> [n*(facet(h,delta,90-180./n)f val(360./n)+)] facet(h,delta,beta) --> [[beta+ delta^ { h f | h f ] f beta-| delta^ h f }]

Figure A. 2

 2 Figure A.2 Flower generated by using the example code from ELSYS

Table 1 . 1

 11 Summary of single angular indices .................................................................................

Table 2 . 1

 21 Parameters of some common used LAD. .......................................................................

Table 2 . 2

 22 Measurement parameters of the elm ...............................................................................

Table 2 . 3

 23 Components optical properties of the single tree scene ..................................................

Table 2 . 4

 24 Input parameters of the RGM model in the real structural single tree scene ..................

Table 2 .

 2

5 Components optical properties of the simplified heterogeneous canopy scene in the blue, green, red and NIR band. ................................................................................................................

Table 2.6 Input parameters of the simplified heterogeneous canopy scene. ...................................

Table 3 .

 3

1 Input parameters for maize scene modelling and canopy reflectance simulation. ..........

Table 4.1 Input optical property (OP) for the DART calibration procedure. ρ means reflectance and τ means transmissivity. SWIR means short wave infrared. ............................................................. Table 4.2 Input simulation parameters for the DART calibration procedure. ................................. Table 4.3 Satellite image parameters. .............................................................................................

Table 5 .

 5

1 Input parameters for the canopy reflectance simulation datasets. ................................... Table 5.2 Additional input parameters for the heterogeneous scenes. ............................................ Table 5.3 Dickey-Fuller test. ..........................................................................................................

Table 5 .

 5

4 Field measured LAI data from GLAIFM for evaluating the LAI estimation accuracy of VIs. ..................................................................................................................................................

Table 5 .

 5

5 Linear regression of the four VIs on LAI (slope, intercept, R 2 , P) .................................

Table 5 . 6

 56 Parameters of the linear and logarithmic fitted model in the homogeneous canopies. ...

Table 5 . 7

 57 Parameters of the exponential and power fit model in the homogeneous canopy. .........

Table 5 . 8

 58 Parameters of the linear and logarithmic fit model in the heterogeneous canopies. .......

Table 5 .

 5

9 Exponential and power fit model parameters in the heterogeneous canopies................. Table 5.10 Parameters of the linear and logarithmic fitted model in the field measured data. ....... Table 5.11 Exponential and power fit model parameters in field measured data. ..........................

Table 6 . 1

 61 Statistics of the relative error of atmospheric effects (Eatmosphere) of EVI, EVI2, NDVI, and SAVI of Landsat-8 (L8) and Sentinel-2 (S2).................................................................................

Table 6 . 2

 62 Statistics of relative error of EVI2 (EEVI2) for remote sensing data. ............................

Table 6 . 3

 63 Descriptive statistics for EVI2 relative error (EEVI2) in USGS Spectral Library data ...

  Vegetation is a general term for all kinds of plants that grow on Earth's surface. As one of the main components of the ecosystem, vegetation plays a vital role in human survival and development, including maintaining and improving ecological balance and the environment. Thus, it is inseparable from human survival and development. As a renewable energy source, vegetation plays an essential role in the earth's ecosystem. Vegetation is a major influencing factor and indicator for monitoring the dynamics of global change. Vegetation influences the energy balance of the earth's climate system and plays a vital role in climatic, hydrological and biochemical cycles. It is an indicator of climate and human activities in the natural world. Vegetation participates in and controls biogeochemical cycles, such as water, carbon and nitrogen. Documenting vegetation dynamics is essential for a better understanding of the terrestrial carbon cycle and improving forest management.

	GF1 NHVI SQL	GaoFen-1 Normalized Hotspot signature Vegetation Index Structured Query Language
	GLAIFM GO NIR NNT SR SRF	Global Leaf Area Index from Field Measurements Geometric Optical Near-InfraRed Neural Networks Tool Surface Reflectance Spectral Response Functions Chapter 1 Introduction
	GRVI NOBS SSA	simple ratio near-infrared/Green Ratio Vegetation Index Northern Old Black Spruce Southern Study Area
	GSD NSA SSO	Ground Sampling Distances Northern Study Area Sun Synchronous Orbit
	GUI NYJP STD	Graphical User Interface Northern Young Jack Pine STandard Deviation
	ALA HDS OLI/TIRS SUNY	Average Leaf inclination Angle Hot-Dark Spot vegetation index Operational Land Imager / Thermal InfraRed Sensor State University of New York
	ATBD HEVI2 ONERA SWIR	Algorithmic Theoretical Basis Document Hotspot-signature 2-band Enhanced Vegetation Index French Office National d'Études et de Recherches Aérospatiales Short Wave InfraRed
	AVHRR HSVI OP SYJP	Advanced Very High Resolution Radiometer Hotspot-signature Soil-adjusted Vegetation Index Optical Property Southern Young Jack Pine
	AVIRIS INRAE OPs SZA	Airborne Visible / InfraRed Imaging Spectrometer Institut National de Recherche pour l'Agriculture, l'alimentation et Optical Properties Sun Zenith Angle
	BOA OSAVI TIR	Bottom Of Atmosphere l'Environnement Optimized Soil Adjusted Vegetation Index Thermal InfraRed
	BOREAS ISS POLDER TOA	BOReal Ecosystem-Atmosphere Study International Space Station POLarization and Directionality of Earth's Reflectance Top Of Atmosphere
	BRDF KH PP TRGM	Bidirectional Reflectance Distribution Function KeyHole Principal Plane Thermal Radiosity-Graphics based Model
	BRF KM PROBA TRIM	Bidirectional Reflectance Factor Kubelka-Munk PRoject for On-Board Autonomy Three-dimensional Radiation Interaction Model
	CESBIO L8 PROSAIL TSAVI	Centre d'Etudes Spatiales de la BIOsphère Landsat-8 model of leaf optical PROperties SPECTra + Scattering by Arbitrary Transformed Soil Adjusted Vegetation Index
	CFMASK LAD URBANFLUXES URBan ANthropogenic heat FLUX from Earth observation Satellites C Function of MASK Leaf Angle Distribution Inclined Leaves
	CHRIS LAI PROSPECT USGS	Compact High Resolution Imaging Spectrometer Leaf Area Index model of leaf optical PROperties SPECTra United States Geological Survey
	CNES LDBM PVI UTC	Centre National d'Etudes Spatiales Leaf Dry Biomass Matter Perpendicular Vegetation Index Universal Time Coordinated
	CP LEDAPS QA UTM	Cross-principal Plane Landsat Ecosystem Disturbance Adaptive Processing System Quality Assurance Universal Transverse Mercator
	DART LIBERTY RAMI VAA	Discrete Anisotropic Radiative Transfer Leaf Incorporating Biochemistry Exhibiting Reflectance and RAdiative transfer Model Intercomparison View Azimuth Angle
	DEM RAPID VI	Digital Elevation Model Transmittance Yields Radiosity Applicable to Porous IndiviDual objects Vegetation Index
	DN LIDAR RGB VIs	Digital Number LIght Detection And Ranging Red-Green-Blue Vegetation Indices
	DOS LOPEX93 RGM VNIR	Disk Operating System Leaf Optical Properties EXperiment 93 Radiosity-Graphics based Model Visible and Near-InfraRed
	DSM LSMA RMSE VZA	Digital Surface Model Linear Spectral Mixture Analysis Root Mean Square Error View Zenith Angle
	DTM LUTs RPV WGS	Digital Terrain Model Look-Up Tables Rahman, Pinty Verstraete World Geodetic System
	DVI MASTER RT ZY3	Differential Vegetation Index Moderate-resolution imaging spectroradiometer / Advanced Spaceborne Radiative Transfer Ziyuan-3
	ECMWF RVI	European Centre for Medium-range Weather Forecasts Thermal Emission and reflection Radiometer airborne simulator Ratio Vegetation Index
	eFAST MESMA S2	extended Fourier Amplitude Sensitivity Test Multiple Endmember Spectral Mixture Analysis Sentinel-2
	ELSYS MISR SAA	Extended L-SYStem Multi-angle Imaging SpectroRadiometer Sun Azimuth Angle
	EMRPV MODIS SAIL	Exponential Modified Rahman, Pinty Verstraete MODerate-resolution Imaging Spectroradiometer Scattering by Arbitrary Inclined Leaves
	ESA MODTRAN SARVI	European Space Agency MODerate resolution atmospheric TRANsmission Soil and Atmosphere Resistant Vegetation Index
	ETM MRPV SAS	Enhanced Thematic Mapper Modified Rahman, Pinty Verstraete SubArctic Summer
	EUMETSAT MSAVI SAVI	EUropean organization for the exploitation of METeorological SATellites Modified Soil Adjusted Vegetation Index Soil Adjusted Vegetation Index
	EVI MSG/SEVIRI SCOPE	Enhanced Vegetation Index Meteosat Second Generation / Spinning Enhanced Visible and InfraRed Soil Canopy Observation, Photochemistry and Energy fluxes
	FAPAR SNR	Fraction of Absorbed Photosynthetically Active Radiation Imager Signal-to-Noise Ratio
	FC MSI SOA	Foliage Coverage MultiSpectral Instrument Southern Old Aspen
	FLIGHT MTF SOBS	3-D Forest LIGHT model Modulation Transfer Function Southern Old Black Spruce
	FPAR NASA SOJP	Fraction of Photosynthetically Active Radiation National Aeronautics and Space Administration Southern Old Jack Pine
	FTIR NCAR SOS	Fourier Transform InfraRed National Center for Atmospheric Research Successive Orders of Scattering
	GCRICAS NCEP SOSVRT	Global Change Research Institute of the Czech Academy of Sciences National Centers for Environmental Prediction Successive Orders of Scattering Vector Radiative Transfer
	GEE NDHD SPOT	Google Earth Engine Normalised Difference between Hotspot and Darkspot Satellite Pour l'Observation de la Terre
	GEMI NDVI SPRINT	Global Environmental Monitoring Index Normalized Difference Vegetation Index SPReading of photons for radiation INTerception

) Providing high-quality validation data for sensor design development

	the production of airborne simulators with the same capabilities before the launch of on-board
	sensors to assess their performance and potential applications, e.g., the National Aeronautics and
	Space Administration (NASA) conducts Moderate-resolution imaging spectroradiometer /
	Advanced Spaceborne Thermal Emission and Reflection Radiometer airborne simulator (MASTER)
	airborne simulation imaging experiments before the launch of MODerate-resolution Imaging
	Spectroradiometer (MODIS) to provide a reference value for its performance.

Simulated remote sensing data can provide a reference for performance tests before satellite launch. Remote sensor developers need to systematically use the simulated data to validate sensor parameters for future development. The sensor application department needs to use the simulated data to prioritize data processing and applicability modelling to use the relevant data rapidly. For some western countries with solid hardware manufacturing capabilities, it is common to prioritize

(2) Providing reliable support for quantitative remote sensing inversion modelling, spatial and temporal scale conversion and assimilation of multi-source remote sensing data

	platforms (airborne, satellite-based, multi-angle, hyperspectral). It can provide reliable support for
	research on quantitative remote sensing inversion modelling, spatial and temporal scale conversion and
	multi-source remote sensing data assimilation.
	Current research on the design of remote sensing sensors is difficult to combine high spatial
	resolution, high temporal resolution and high spectral resolution simultaneously, especially when
	interfered with by cloud cover, weather and other imaging conditions. Thus, it cannot obtain sufficient
	satellite data to meet the needs of scientific research applications such as disaster emergency monitoring
	and continuous time series analysis. The research on vegetation reflectance simulation and bio-optical
	property inversion based on a three-dimensional radiative transfer model helps to integrate remote
	sensing, computer, agronomy, mathematics and other methods by coupling soil-leaf-canopy-atmosphere
	radiative transfer model with optical remote sensor imaging models, constructing multi-dimensional
	remote sensing data simulation of continuous multi-temporal crop growth, and realizing multiple load

1.2 Research background 1.2.1 Canopy reflectance simulation study

Table 1 . 1

 11 Summary of single angular indices

	VI	Prerequisites	Strengths	Limitations
	RVI	Rred, Rnir	Good at high-density vegetation	Poor at low-density vegetation
			coverage	coverage
	NDVI	Rred, Rnir	Relate to canopy structure, LAI and	Sensitive to the soil brightness, soil
			canopy photosynthesis	colour, atmosphere, cloud and
				cloud shadow, and shadow
	SAVI	Rred, Rnir, XSAVI	Insensitive to soil noise	The determination of soil
				adjustment factor
	MSAVI Rred, Rnir	No soil adjustment factor	Sensitive to atmosphere
	PVI	Rred, Rnir, a, b
			Insensitive to soil noise	Require soil line parameters
	OSAVI Rred, Rnir, XOSAVI	Do not depend on soil line	Applications are not extensive
			parameters
	EVI	Rblue, Rred, Rnir XEVI, C1, C1 Good at soil and atmosphere noise	Require blue band reflectance
	EVI2	Rred, Rnir, XEVI, C1, C1, c	No requirement for blue band	Sensitive to sensor
			reflectance

Good at low-density vegetation coverage Poor at high-density vegetation coverage TSAVI Rred, Rnir, XTSAVI, a, b

  Three types of vegetation canopy scenes created by the discrete anisotropic radiative transfer (DART) model are presented, including a homogeneous assumption canopy scene, a simplified heterogeneous canopy scene and a real structural canopy scene, together with the modelling steps. The principles of radiative transfer modelling for RGM and DART are presented to provide the basis for subsequent chapters.

As an essential chapter, this chapter details the three 3D radiative transfer models mentioned in the introduction (section 1.2.1): SAIL, RGM and DART. Section 2.1 focuses on the vegetation remote sensing scene construction, together with their radiative transfer simulation for SAIL (section 2.1.1), DART (section 2.1.2) and RGM (section 2.1.3), respectively. One major limitation of RGM is the huge calculation load of view factors. Therefore, section 2.2 focuses on improving the RGM model: we speed up the RGM model using the multithreading technology based on the OpenACC. The results are compared and validated with the RGM and DART models.

Table 2 . 1

 21 Parameters of some common used LAD.

	Leaf angle distribution	Function forms	Mean leaf angle (°)
	Spherical	g l * (θ l )=sinθ l	57.30
	Uniform	g l * (θ l )=2/π	45.00
	Planophile	g l * (θ l )=2(1+cos2θ l )/π	26.76
	Erectophile	g l * (θ l )=2(1-cos2θ l )/π	63.24
	Plagiophilel	g l * (θ l )=2(1-cos4θ l )/π	45.00
	Extremophile	g l

* (θ l )=2(1+cos4θ l )/π 45.00

  DART Model (https://dart.omp.eu/#/) is a physically-based 3D radiative transfer model that simulates the radiative transfer of electromagnetic radiation from visible to thermal infrared wavelengths in a realistic 3D structural scene. It has been developed at the CESBIO in Toulouse,

	France, since 1992 and was patented in 2003. It provides radiation balance and measurements from
	satellite, airborne and ground-based passive (imaging spectrometer) to remote sensing active light
	detection and ranging (LiDAR). The DART [28, 88, 103-105, 129-133] modelled the surface
	explicitly, so the spatial extent and accuracy of the simulation are limited only by the input data and
	available computational power.
	(1) DART scene creation
	DART scenes can be built with both voxels or facets, and the voxels representing vegetation in

the DART model can be filled in two ways to represent the scattering properties, including turbid media or discrete triangular facets. A turbid medium is the light scattering of significant intensity due to irregular (randomly distributed) optical inhomogeneities in composition or generation.

Table 2 . 2

 22 Measurement parameters of the elm

	Parameters	Value	Parameters	Value
		Tree height	1.31		East	0.50
	Height (m)	Trunk height Canopy height	0.77 0.56	Crown width (m)	South West	0.65 0.53
		Height below canopy	0.69		North	0.53
		Base diameter	0.16	Average stem length (m)	0.0056
	Diameter at breast height (m)	Diameter at breast 1 Diameter at breast 2 Diameter at breast 3	0.15 0.14 0.12	Average leaf inclination Leaf area index Number of twig leaves	59.85° 5.27 12
		Diameter at breast 4	0.12	Average branch inclination	76.4° / 132.6°

Table 2 . 3

 23 Components optical properties of the single tree scene

	Grouping	Bands	Blue (440.0nm)	Green (550.0nm)	Red (660.0nm)	NIR (860.0nm)
	Leaf	Reflectance Transmittance	0.062 0.029	0.119 0.166	0.107 0.096	0.349 0.573

Table 2 . 5

 25 Components optical properties of the simplified heterogeneous canopy scene in the blue, green, red and NIR band.

Table 2 . 6

 26 Input parameters of the simplified heterogeneous canopy scene.

Table 3 . 1

 31 Input parameters for maize scene modelling and canopy reflectance simulation.

		Parameters	Value
		Spectral parameters *	Red	NIR
	Rsoil	Soil reflectance (%)	13.6	21.6
	Rleaf	Leaf reflectance (%)	8.48	46.58
	Tleaf	Leaf transmittance (%)	0.62	47.59
		Maize Scene
	Nmax	Maximum number of leaves produced	22
	Smax	Maximum total leaf area produced by a plant (m 2 )	0.5
	D	Plant density (number of plants / m 2 )	8.8
	drows	Row spacing (m)	0.8
	Hmax	Maximum canopy height (m)	2.5
	θmax	Maximum leaf inclination (°)	45
	Td	Cumulative degree day difference	100-1950 (step 50)
		Sun-target-observation geometry
	SZA	Sun zenith angle (°)	45
	SAA	Sun azimuth angle (°)	225
	VZA	View zenith angle (°)	0
	VAA	View azimuth angle (°)	0

* The optical properties (OPs) of the leaves were selected from Leaf Optical Properties Experiment 93 (LOPEX93) data.

Table 4 . 1

 41 Input optical property (OP) for the DART calibration procedure. ρ means reflectance and τ means transmissivity. SWIR means short wave infrared.

	Component OP	Ground (ρ)	Roof (ρ)	Shrubs (ρ)	Water (ρ)	Trees (ρ + τ)
		Blue (492.4 ± 66nm)	0.2169	0.0748	0.1173	0.0695	0.3578
		Green (559.8 ± 36nm)	0.2483	0.1122	0.2635	0.0649	0.4600
	Basic simulation	Red (664.6 ± 31nm)	0.2561	0.2501	0.0987	0.0543	0.6140
	initial OP	NIR (832.8 ± 106nm)	0.2366	0.3100	0.9289	0.0283	0.9289
		SWIR 1 (1613.7 ± 91nm)	0.2011	0.7165	0.6794	0.0053	0.8722
		SWIR 2 (2202.4 ± 175nm)	0.1693	0.7190	0.4178	0.0051	0.6497
		Blue (492.4 ± 66nm)	0.0574	0.1067 0.1131 0.0574	0.0927 0.0536	0.0277 0.0364 0.0326	0.1153 0.0829
		Green (559.8 ± 36nm)	0.0710	0.1456 0.1206 0.0710	0.1996 0.1329	0.0269 0.0250 0.0380	0.1733 0.1942
	Pseudo-satellite simulation of initial OP	Red (664.6 ± 31nm) NIR (832.8 ± 106nm)	0.1003 0.1681	0.3074 0.1308 0.1003 0.3257 0.1460 0.1681	0.0892 0.0405 0.9426 0.9013	0.0267 0.0206 0.0245 0.0262 0.0198 0.0198	0.0854 0.0820 0.9596 0.9013
		SWIR 1 (1613.7 ± 91nm)	0.3349	0.3623 0.1049 0.3349	0.5631 0.5551	0.0210 0.0186 0.0186	0.5199 0.6668
		SWIR 2 (2202.4 ± 175nm)	0.2635	0.3213 0.1009 0.2635	0.3091 0.2494	0.0189 0.0163 0.0163	0.2778 0.4911

  The calibration accuracy of the pure pixels appeared to be more accurate than the mixed pixels: |ε| ��� pure ground =0.0065 , |ε| ��� mixed ground =0.0693 ; |ε| ��� pure roof =0.0070 , |ε| ��� mixed roof =0.0089 ; |ε| ��� pure water =0.0015 , |ε| ��� mixed water =0.0268 ; |ε| ��� pure tree =0.2519 ,

Table 4 . 2

 42 Input simulation parameters for the DART calibration procedure.

	Parameters	No noise	With noise
	MTF	1	0.15 -0.3
	Solar zenith angle	30°	0 -60°	The eFAST sensitivity
	Pixel resolution	4 m	0.5 -15 m	analysis determines the
	X-axis pixel shifting	0 pixel	0 -1 pixel	step size for each input
	Y-axis pixel shifting	0 pixel	0 -1 pixel	parameter.
	Inaccurate city database	100%	95% -105%

Table 4 . 3

 43 Satellite image parameters.

	Parameters	Value
	Identity document	20190915_100328_1003
	Date	September 15, 2019, 10:03:28 UTC
	Source	Four-band PlanetScope
	Percentage of cloud computing	0%
	Spatial resolution	3 m
	Serial identifier	1003
	Off nadir angle	0.0°
	Sun altitude angle	42.2°
	Solar azimuth angle	152.0°
	Sensors	PS2
	Spectral band (nm)	Blue 455-515 500-590 590-670 780-860 Green Red NIR

  NDVI

	Based on Eq.(1.2), we can get
			nir	red	nir
	nir	red	red		red
	nir	red	nir	red	nir
			red		red

Table 5 . 3

 53 Dickey-Fuller test.

	Parameters	Value
	Test statistics	-6.979
	MacKinnon's approximate p-value 8.281 × 10 -10
	Lag time for use	1
	Number of observations used	26
	Threshold (1%)	-3.711
	Threshold (5%)	-2.981
	Threshold (10%)	-2.630

Table 5 . 4

 54 Field measured LAI data from GLAIFM for evaluating the LAI estimation accuracy of VIs.

	Name

of measuring station Latitude (°) Longitude (°) LAI Date

  VIs were compared using a linear regression model. The coefficients of the linear regression model are given in Table5.5. Figure5.11 and Table5.5 clearly show that all VIs have an intense saturation effect. The VIs show almost no upward trend with the increase of LAI. In addition, all VIs were weakly correlated with LAI. Overall, NDVI showed better performance in both cases, with the highest R 2 (0.1632 for MODIS LAI and 0.4313 for measured LAI), the lowest P (0.0173 for MODIS LAI and 0.2860 for measured LAI), and the second-highest is TSAVI. EVI and SAVI performed unsatisfactorily, especially with MODIS LAI data that were negatively correlated. In conclusion, using VIs to invert LAI in very dense vegetation areas may lead to significant errors.

	BOREAS NSA/OJP, Thompson	55.92	-98.62	4.38	July 1994
	BOREAS NSA/OBS, Thompson	55.91	-98.45	4.06	July 1994
	BOREAS NSA, Thompson	55.91	-98.52	8.41	July 1994
	BOREAS NSA, Thompson	55.80	-98.00	6.21	July 1994
	BOREAS NSA, Thompson	55.75	-97.80	5.44	July 1994
	BOREAS SSA, Prince Albert	54.06	-105.93	10.59 August 1994
	Arakawa River, Urawa	35.83	139.62	4.24 September 1985
	Westvaco, Summerville, SC	33.20	-80.25	10.4 February 1991
	The ability of the VIs to estimate LAI was examined, and the scatter plots of LAI (y-axis) and
	VIs (x-axis) are given in Figure 5.11.

Table 5 . 5

 55 Linear regression of the four VIs on LAI (slope, intercept, R 2 , P)

	VIs	LAI Type	Slope intercept R 2	P *
	NDVI	MODIS LAI In-situ measurement of LAI 6.2097 3.4104 0.4313 0.2860 2.4769 3.1819 0.1632 0.0173
	SAVI	MODIS LAI In-situ measurement of LAI 10.1760 3.9748 0.2915 0.4836 -2.4165 6.4499 -0.0904 0.1894
	TSAVI	MODIS LAI In-situ measurement of LAI 7.6553 4.9714 0.3386 0.4120 1.0364 4.7082 0.0454 0.5103
	EVI	MODIS LAI In-situ measurement of LAI 4.4827 5.3874 0.1017 0.8106 -3.0249 6.8936 -0.1504 0.0285
	* Hypothesis testing for two-sided P with a null hypothesis of a zero-slope using the Wald Test with a t-
	distribution test statistic.

Table 5 .

 5 6 and Table 5.7, Table 5.8 and Table 5.9, and Table 5.10 and Table 5.11, respectively. The corresponding standard fit statistics for the 10 VIs are compared in

Table 5 . 6

 56 Parameters of the linear and logarithmic fitted model in the homogeneous canopies.

	VIs	Model *	Linear SZA = 30° SZA = 45° SZA = 60° SZA = 30° SZA = 45° SZA = 60° Logarithmic
	NDVI	a b	6.4063 -2.0500	6.3285 -2.0923	6.2162 -2.1809	2.9696 3.6740	2.9756 3.5953	3.0051 3.4604
	SAVI	a b	8.7534 -1.8983	8.6769 -1.9523	8.5810 -2.0660	3.0747 4.8821	3.0868 4.8079	3.1301 4.6859
	EVI2	a b	7.7166 -1.5871	7.6389 -1.6390	7.5345 -1.7445	2.7929 4.5862	2.8011 4.5082	2.8339 4.3769
	NDVI57.5	a b	6.0725 -2.1270	5.9904 -2.2003	5.9404 -2.3170	2.9614 3.4162	2.9872 3.3092	3.0493 3.1917
	SAVI57.5	a b	8.6167 -2.0599	8.1519 -2.0912	7.6726 -2.1168	3.1268 4.7032	3.1109 4.4746	3.0978 4.2301
	EVI257.5	a b	7.5958 -1.7451	7.0715 -1.7582	6.5199 -1.7521	2.8353 4.4020	2.8038 4.1501	2.7687 3.8789
	NHVI	a b	9.9660 -0.5502	8.0588 -0.5706	6.7866 -0.5960	1.2336 4.1841	1.2596 3.9102	1.3257 3.7150
	NHVI2	a b	25.1249 -0.6349	21.2949 -0.6827	18.7819 -0.7489	1.3073 5.4268	1.3512 5.2435	1.4373 5.1713
	HSVI	a b	35.1257 -0.6228	29.8895 -0.6759	26.0295 -0.7161	1.3514 6.0509	1.4097 5.9153	1.4915 5.8491
	HEVI2	a b	32.3840 -0.5541	27.5053 -0.6035	23.8331 -0.6366	1.3075 5.9062	1.3616 5.7644	1.4341 5.6750

Table 5 . 7

 57 Parameters of the exponential and power fit model in the homogeneous canopy.

	VIs	Model *	Exponent SZA = 30° SZA = 45° SZA = 60°	SZA = 30°	Power SZA = 45°	SZA = 60°
	NDVI	a b	0.0001 11.6596	0.0001 12.5845	0.0001 13.4342	12.1888 10.9410	12.1704 11.7817	11.6137 12.3995
	SAVI	a b	0.0284 7.3482	0.0154 8.1770	0.0066 9.3108	25.8920 4.6696	32.6541 5.4084	43.0951 6.3311
	EVI2	a b	0.0633 5.6524	0.0412 6.1712	0.0234 6.8414	13.1243 3.7341	14.6867 4.2725	16.6426 4.9316
	NDVI57.5	a b	0.0001 11.2820	0.0001 13.1848	0.0001 17.5224	8.6426 10.5486	8.7964 12.4349	10.2528 16.6176
	SAVI57.5	a b	0.0033 10.3810	0.0148 7.8152	0.0381 6.1618	60.2814 7.1697	24.9076 5.4532	13.5794 4.2936
	EVI257.5	a b	0.0136 7.6138	0.0434 5.7140	0.0930 4.4265	20.8621 5.6084	11.0840 4.2816	7.0568 3.3343
	NHVI	a b	0.0835 8.9803	0.1660 5.9663	0.2500 4.3287	89.7483 3.6812	21.4376 2.6495	10.1741 2.0753
	NHVI2	a b	0.0441 25.9988	0.0956 18.2210	0.1510 13.9267	10014.4321 805.8288 194.1731 4.4615 3.3393 2.6764
	HSVI	a b	0.0563 32.9982	0.1130 23.4318	0.2070 16.8295	23131.5800 1719.9606 254.5588 4.2277 3.2084 2.3737
	HEVI2	a b	0.0898 26.9992	0.1464 19.9430	0.2451 14.5383	5258.0039 3.6517	771.4886 156.8637 2.9119 2.2019

Table 5 . 8

 58 Parameters of the linear and logarithmic fit model in the heterogeneous canopies.

	VIs	Model *	Linear SZA = 30° SZA = 45° SZA = 60° SZA = 30° SZA = 45° SZA = 60° Logarithmic
	NDVI	a b	7.8729 -2.3171	7.6749 -2.4487	7.3225 -2.6781	3.3747 4.3119	3.4223 4.1379	3.5240 3.8259
	SAVI	a b	13.7324 -2.6396	13.0677 -2.6028	12.0432 -2.5615	3.8605 6.4948	3.7568 6.2230	3.6283 5.8137
	EVI2	a b	13.0048 -2.3332	12.3225 -2.2946	11.2624 -2.2456	3.5739 6.2595	3.4674 5.9850	3.3312 5.5696
	NDVI57.5	a b	7.0462 -2.6869	6.8924 -2.8198	6.8876 -3.1464	3.4875 3.6592	3.5503 3.4872	3.7674 3.2728
	SAVI57.5	a b	11.2763 -2.4189	10.2241 -2.3001	8.8646 -2.1610	3.4383 5.5098	3.2923 5.1327	3.1229 4.6169
	EVI257.5	a b	10.5078 -2.1145	9.3752 -1.9827	7.8813 -1.8064	3.1496 5.2709	2.9955 4.8769	2.8038 4.3288
	NHVI	a b	9.1791 -0.8825	7.0631 -1.2617	5.9318 -1.8892	1.8945 4.7688	2.3130 4.2868	3.1206 3.7394
	NHVI2	a b	25.5830 -1.0611	21.7590 -1.5381	20.2524 -2.3213	2.0386 6.8410	2.5103 6.9912	3.3598 7.5477
	HSVI	a b	42.1506 -1.0973	36.1692 -1.5584	33.4403 -2.2164	2.2726 8.4914	2.8521 9.1066	3.6725 10.0329
	HEVI2	a b	40.9385 -0.9887	34.6456 -1.3899	31.2983 -1.9382	2.1854 8.3039	2.7020 8.8022	3.3898 9.4852

Table 5 . 9

 59 Exponential and power fit model parameters in the heterogeneous canopies.

	VIs	Model *	Exponent SZA = 30° SZA = 45° SZA = 60° SZA = 30° SZA = 45° SZA = 60° Power
	NDVI	a b	0.0964 4.7755	0.0536 5.3445	0.0055 7.6861	9.2658 3.2634	9.4684 3.9035	10.8040 6.3928
	SAVI	a b	0.0523 9.1001	0.0343 9.6392	0.0190 10.1963	70.5309 3.8354	86.0224 100.3810 4.3156 4.9717
	EVI2	a b	0.0856 7.9522	0.0634 8.2457	0.0424 8.4539	45.8887 3.3123	50.7466 3.6533	51.9962 4.1015
	NDVI57.5	a b	0.0028 8.2077	0.0002 10.8540	0.0001 16.3924	9.6727 7.1050	10.0167 9.8359	11.3722 15.2866
	SAVI57.5	a b	0.0089 11.3145	0.0251 8.7487	0.0697 6.2624	149.7779 52.6041 5.8634 4.7010	18.7216 3.5264
	EVI257.5	a b	0.0266 9.0266	0.0572 7.0199	0.1339 4.8947	63.0527 4.6664	27.5656 3.8121	11.3918 2.8415
	NHVI	a b	0.3119 4.5055	0.3208 3.1593	0.2559 2.5755	14.6215 2.1344	7.1111 2.1360	3.7794 2.3528
	NHVI2	a b	0.2407 13.6717	0.2021 11.3338	0.1252 10.7294	198.7988 145.3202 145.1828 2.4272 2.6415 3.1182
	HSVI	a b	0.2441 21.9301	0.2202 17.8326	0.1741 15.8897	539.9359 368.4561 340.4486 2.3469 2.4481 2.6907
	HEVI2	a b	0.2849 20.4285	0.2542 16.6526	0.2036 14.7390	363.7233 264.7138 243.8130 2.1671 2.2859 2.5203

Table 5 .

 5 10 Parameters of the linear and logarithmic fitted model in the field measured data.

	VIs	Model *	Linear SZA = 30° SZA = 45° SZA = 60° SZA = 30° SZA = 45° SZA = 60° Logarithmic
	NDVI	a b	4.6498 0.2837	4.7694 0.0426	3.2455 1.0764	3.4655 4.8081	3.8051 4.7407	2.8431 4.3225
	SAVI	a b	2.2827 2.7866	-0.6380 3.8819	-4.2412 5.1516	1.2828 4.9267	-0.0548 3.5933	-1.8177 1.7174
	EVI2	a b	1.7316 3.0108	-0.7760 3.9229	-3.7763 4.9369	1.0510 4.7286	-0.1245 3.5162	-1.6159 1.8574
	NDVI57.5	a b	9.4790 -4.4473	-6.0954 8.9971	-5.0803 8.3122	8.3484 4.9790	-5.2254 2.9469	-3.8769 3.2635
	SAVI57.5	a b	-6.5589 5.8180	-8.1280 6.4513	-6.8570 6.2299	-2.7449 0.5229	-3.3045 -0.0307	-1.9831 1.4127
	EVI257.5	a b	-5.7890 5.4716	-7.1026 6.0001	-6.1330 5.9238	-2.3988 0.7758	-2.8775 0.2885	-1.7936 1.5507
	NHVI	a b	4.5848 1.1239	2.7814 0.7873	0.1864 2.9006	3.0574 5.6888	3.0322 3.8314	1.1300 3.0587
	NHVI2	a b	22.1552 -0.6220	18.2641 -1.6869	7.2137 0.8512	4.4854 11.1883	5.2219 10.2372	3.2502 6.9991
	HSVI	a b	54.1838 -1.7004	53.6655 -3.6825	21.5021 0.1122	5.2466 15.9358	6.6833 17.0805	3.9000 10.8928
	HEVI2	a b	55.7362 -1.6911	53.2029 -3.3414	22.5904 0.1046	5.0053 15.5258	6.2135 16.3838	3.7857 10.8553

* Linear: y=ax+b; logarithmic: y=a•ln(x)+b; where a and b are the parameters of the fitted model, determined using the Marquardt-Levenberg algorithm, x represents VI and y represents LAI. SZA represents the solar zenith angle.

Table 5 .

 5 [START_REF] Yang | Review of Optical Remote Sensing Imaging Simulation of Farmland Radiation Transfer Process[END_REF] Exponential and power fit model parameters in field measured data.

	VIs	Model *	Exponent SZA = 30° SZA = 45° SZA = 60° SZA = 30° SZA = 45° SZA = 60° Power
	NDVI	a b	1.6311 1.1053	1.5656 1.1138	1.9869 0.7652	4.8122 0.8471	4.7060 0.9037	4.2722 0.6728
	SAVI	a b	3.0004 0.5162	3.8697 -0.1610	5.8925 -1.3725	4.8536 0.2878	3.5996 -0.0134	1.9396 -0.5868
	EVI2	a b	3.1536 0.3944	3.9132 -0.1984	5.5222 -1.2359	4.6422 0.2355	3.5307 -0.0309	2.0114 -0.5291
	NDVI57.5	a b	0.5389 2.2336	15.9786 -1.6926	10.4680 -1.1625	4.9744 1.9781	2.9834 -1.4344	3.2957 -0.8952
	SAVI57.5	a b	8.6646 -2.7139	11.2891 -3.5204	7.3022 -2.0090	0.9924 -1.1065	0.8880 -1.1965	2.0608 -0.4647
	EVI257.5	a b	7.5554 -2.4101	9.6873 -3.2003	6.9533 -1.9146	1.0837 -0.9783	0.9744 -1.0611	2.1116 -0.4266
	NHVI	a b	2.0857 0.9757	1.7859 0.6541	2.9769 0.0407	5.6254 0.6922	3.6295 0.7801	3.0649 0.2504
	NHVI2	a b	1.3496 5.0011	0.8046 4.9467	1.8607 1.6221	22.7432 1.1119	22.8017 1.5101	7.8036 0.7861
	HSVI	a b	0.9915 12.8240	0.3211 17.0814	1.5504 4.9349	94.1681 385.1687 21.1641 1.4088 2.3729 0.9732
	HEVI2	a b	0.8982 14.1849	0.2759 18.8194	1.5071 5.3470	111.8443 471.7868 22.5653 1.4661 2.4311 0.9839

Table 6 . 1

 61 Statistics of the relative error of atmospheric effects (Eatmosphere) of EVI, EVI2, NDVI, and SAVI of Landsat-8 (L8) and Sentinel-2 (S2).

	Statistics	L8	EVI	S2	L8	EVI2	S2	L8	NDVI	S2	L8	SAVI	S2
	25 percentile	0.252		0.022	-0.175	-0.297	-0.212	-0.290	-0.169	-0.278
	50 percentile	0.295		0.094	-0.111	-0.238	-0.137	-0.213	-0.105	-0.220
	75 percentile	0.324		0.158	-0.087	-0.189	-0.105	-0.162	-0.082	-0.173
	Max	7137.023		29965.000	9.900	3808.099 9.689 5133.701	8.934	4125.002
	Mean	-0.203		0.143	-0.376	-0.240	-0.240	-0.272	-0.369	-0.226
	Min	-10750.800	-29861.800	-97.576 -345.676 -9.557 -313.473 -97.324 -357.240
	STD	20.189		47.252	1.612	4.187	0.268		3.891	1.606	4.313

  investigate cross-sensor continuity of NDVI and EVI, investigations on EVI2 cross-sensor translation are quite few. It indeed induces uncertainty if the VIs from the different sensors are not calibrated [231]. Cross-sensor continuity studies are essential for vegetation because the vegetation monitor usually requires a high frequency of VIs observation in the growth period, while one single sensor is challenging to provide VIs with such a high temporal frequency [4]. Thus, it is necessary to verify the stability of EVI2 on other sensors except for MODIS, Remote sensing satellites commonly used for vegetation monitoring include Landsat and Sentinel. Landsat is a joint project of the US Geological Survey and NASA and has been continuously observing the Earth from 1972 to the present. Today, Landsat satellites image the Earth's entire surface at a resolution of 30 m approximately once every fortnight, including multispectral and thermal data. The data produced by the United States Geological Survey (USGS)

	evaluate the consistency between EVI2 and EVI on L8 and S2 sensors and calibrate the spectral
	correlation coefficient c of EVI2. The calibrated EVI2 should better apply to the L8 and S2 satellites
	and other sensors with similar spectral response function (SRF) by having the best similarity as the
	3-band EVI. First, we compared the atmosphere effect on EVI and EVI2 and evaluated their
	consistency on L8 and S2 sensors. Then, two new spectral correlation coefficients, c L8 and cS2, are
	derived by regressing the band coefficient between blue and red band reflectance using L8 and S2
	SR products. Finally, we use the United States Geological Survey (USGS) reflectance library
	version 7 [234] to validate the availability of cL8 and cS2.

such as the widely used Sentinel-2 (S2) and Landsat-8 (L8) sensors. A study has been conducted to provide a harmonized surface reflectance for these two satellites to provide a high time-frequency and spatial cover monitor

[START_REF] Claverie | The Harmonized Landsat and Sentinel-2 surface reflectance data set[END_REF]

.

Therefore, following the recommendation from Jiang et al.

[START_REF] Jiang | Development of a two band enhanced vegetation index without a blue band[END_REF]

, we aim to quantitatively

Table 6 . 2

 62 Statistics of relative error of EVI2 (EEVI2) for remote sensing data.

	Statistics	Landsat-8	Sentinel-2
	25 percentile	-0.044	-0.109
	50 percentile	-0.025	-0.063
	75 percentile	0.002	-0.021
	Max	1.322	1.2
	Mean	-0.025	-0.09
	Min	-6.054	-7.067
	STD	0.119	0.17

Table 6 . 3

 63 Descriptive statistics for EVI2 relative error (EEVI2) in USGS Spectral Library data

	Statistics	EVI2	Landsat-8 EVI2L8	EVI2	Sentinel-2 EVI2S2
	25 percentile	-0.060	-0.093	-0.054	-0.117
	50 percentile	-0.032	-0.065	-0.020	-0.083
	75 percentile	-0.006	-0.030	0.006	-0.047
	Max	0.831	0.641	1.075	0.653
	Mean	-0.028	-0.064	-0.017	-0.089
	Min	-0.587	-0.627	-0.577	-0.656
	STD	0.121	0.109	0.140	0.112

Table A . 3

 A3 Commonly used mathematical operators.

	Operator	Function
	+	Add real numbers and integers
	-	Subtract real numbers and integers
	*	Multiply real numbers and integers
	/	Divide real numbers and integers
	**	Power function
	^	Power function
	mod	Modulo function for integers
	and	Logical operators "and" on integers
	or	Logical operators "or" on integers
	not	Logical operators "not" on integers
	<	For comparing numbers using smaller
	>	For comparing numbers using larger
	<=	For comparing numbers using smaller or equal
	>=	For comparing numbers using larger or equal
	=	Equality test, both for numbers and vectors.
	:=	Set the value of a variable.

* Linear: y=ax+b; logarithmic: y=a•ln(x)+b; where a and b are the parameters of the fitted model, determined using the Marquardt-Levenberg algorithm, x represents VI and y represents LAI. SZA represents the solar zenith angle.

all, I would like to thank my supervisor, Professor Jean-Philippe Gastellu-Etchegorry, from the Centre d'Etudes Spatiales and a National Aeronautics and Space Administration (NASA) researcher in the USA. The original code of the Radiosity-Graphics based Model (RGMSincere thanks to Prof. Binet Renaud from Centre National d'Etudes Spatiales Sincere thanks to Prof. Elise Koeniguer from the French Office National d'

geometry registration errors and the accuracy of the 3D scene were often not well known. Therefore, they differed in satellite and base simulations. MTF was modelled with a Gaussian filter but was only used in satellite simulations because the current DART calibration methods could not consider it. Similarly, the geometric calibration of satellite images and the geometry of the 3D city scene were not perfectly accurate. Therefore, they were fixed in the satellite simulation and varied in the base simulation to simulate geometric inaccuracies. Each influence factor (Table 4.2) produced 70 samples, which resulted in 420 samples in total. The range of MTF was set to 0.15 -0.3, corresponding to the 10 m spatial resolution band of the Sentinel-2 (S2) image [START_REF] Sentinel | SENTINEL-2 products respect stringent radiometric and geometric image quality requirements[END_REF]. The range of SZA was set to 0 -60°. The range of spatial resolution was set to 0.5 -15 m, corresponding to most high-resolution satellites, including Satellite Pour l'Observation de la Terre (SPOT), RapidEye, Advanced Land Observing Satellite (Alos), Ziyuan-3 (ZY3), Gaofen-1 (GF1), KeyHole (KH), IKONOS, QuickBird, WorldView, GeoEye, and the Pleiades. According to the literature [START_REF] Plyer | A new coregistration algorithm for recent applications on urban SAR images[END_REF][START_REF] Brigot | Adaptation and evaluation of an optical flow method applied to coregistration of forest remote sensing images[END_REF],

the co-registration accuracy of satellite images can be subpixels. Therefore, the maximum range of the pixel shifting was set to one pixel. Furthermore, the input parameters used to construct the 3D scene (e.g., the geometry of buildings, trees, water and leaf area) were not precisely known. As a result, the coefficient of variation of these parameters in the base simulation with noisy interference scenes was between 95% and 105% compared to the pseudo-satellite simulation.

The global sensitivity analysis focused on the effect of differences in the values of the DART input parameters on the accuracy of the inversion-derived OP. Global sensitivity analyses were performed only in the presence of noisy interference cases. The global sensitivity analysis model of the extended Fourier Amplitude Sensitivity Test (eFAST) [START_REF] Saltelli | A quantitative model-independent method for global sensitivity analysis of model output[END_REF] was used to assess the effects of these factors, including the MTF, SZA, satellite spatial resolution, pixel-shifting errors in the x and y axes, and inaccuracies of the 3D scene.

The eFAST is a sensitivity analysis method that calculates a model's predicted values and variance. The calculation is based on turning a multi-dimensional problem into a one-dimensional integration over all unknown model inputs. A search curve through the entire parameter space is created to avoid multidimensional integration. The Fourier series representation decomposition obtains each input variable's fractional contribution to the model's predicted variance. The analysis consists of four stages: (1) establishing the range and distribution of the input parameters and specifying the predicted values and output variances in the form of integrals over the input parameter space; (2) transforming the given multidimensional integrals into single-dimensional integrals over the input parameter space; (3) estimating the predicted values and performance variances; and (4)

estimating the sensitivity indices.

Summary

The accuracy of the DART calibration to invert the OPs of each component of the urban scene in the visible, NIR and two short-wave infrared bands is assessed in noise-free ideal experiments and experiments with artificially added noise. Sensitivity analyses are also carried out for a number of parameters. In the ideal case without noise, the average relative errors for all bands for ground, roof, water, trees and shrubs are 0.013, 0.005, 0.027, 0.297 and 0.250, respectively. In noise interference case, the corresponding values are 0.233, 0.507, 3.088, 0.834 and 1.256.

(1) Calibration accuracy is influenced by SZA, satellite spatial resolution, pixel shifting, 3D city database inaccuracy and MTF in descending order of magnitude.

(

The calibration is insensitive to the number of bands and is highly resistant to shadow interference.

(3) Some uncontrollable external factors may significantly affect inversion accuracies, such as co-registration errors or user-defined leaf reflectance and transmittance ratios.

Chapter 5 Using vegetation indices to estimate leaf area index

Leaf area index (LAI), defined as the one-sided green leaf area per unit ground area in broadleaf canopies and one-half the total needle surface area per unit ground area in coniferous canopies, is an essential indicator for describing the canopy structure of forest ecosystems. LAI estimates the stand's primary productivity and evaluates forest condition over large areas. A number of studies link the LAI and radiometric measurements of vegetation. Despite their wide variety and is designed to correlate with LAI [START_REF] Baret | Potentials and limits of vegetation indices for LAI and APAR assessment[END_REF], all VIs inevitably suffer from saturation effects [START_REF] Huete | A comparison of vegetation indices over a global set of TM images for EOS-MODIS[END_REF]. For LAI values more than a certain threshold, the derivative of VIs is relative to LAI decreases. It is the so-called saturation effect [START_REF] Baret | Potentials and limits of vegetation indices for LAI and APAR assessment[END_REF][START_REF] Huete | The use of vegetation indices in forested regions: issues of linearity and saturation[END_REF]. Saturation is primarily due to the low sensitivity of reflectance in the red band [START_REF] Huete | The use of vegetation indices in forested regions: issues of linearity and saturation[END_REF][START_REF] Zhen | Potentials and Limits of Vegetation Indices With BRDF Signatures for Soil-Noise Resistance and Estimation of Leaf Area Index[END_REF]. It directly limits the application of VIs in areas of dense vegetation areas.

Therefore, this chapter aims to alleviate the saturation of vegetation indices for LAI estimation. The chapter is presented in the papers:

(1) Z. Zhen, S. Chen, T. Yin, E. Chavanon, N. Lauret, J. Guilleux, M.

Summary

The robustness of LAI estimation using VIs is compared. We recommend using the negative soil adjustment factor of SAVI in dense vegetation cover area for the single observation sensors; we recommend using hotspot-signature vegetation indices for the multi-angle observation sensors.

(1) The right shift phenomenon indicates that as the FC increases, the intersection points between the vegetation isoline and the soil line gradually move towards the positive red band axis.

If the intercept of the vegetation isoline is smaller than that of the soil line, the final intersections can reach the positive area of the red band axis. In areas with low LAI, vegetation isolines tend to be parallel. In areas with high LAI, vegetation isolines tend to converge at a point.

(2) The vegetation isolines behaviours may help explain why VIs such as PVI, which assume parallel vegetation isolines, perform better in areas with low FC. In contrast, VIs such as RVI, NDVI, SAVI, and TSAVI assume a convergence of vegetation isolines to a common point perform better in areas with high FC.

(3) Therefore, based on the intersection right shift phenomenon, we recommend using PVI for extremely low vegetation cover, using SAVI and TSAVI for relatively low vegetation cover, using NDVI for relatively high vegetation cover, using SAVI with a negative soil adjustment factor for extremely high vegetation cover.

(4) Since the optimal soil adjustment factor is the negative value of the abscissa of the intersection of the vegetation isoline and the soil line, the optimal soil adjustment factor for SAVI can take a negative value in high FC area according to the right shift phenomenon. In the experiment, the optimal soil adjustment factor is approximately -0.148 when the mean LAI is 5.35 and -0.183 when the mean LAI is 6.72.

(5) Hotspot signature vegetation indices show better performance for LAI estimation than single angular indices in the heterogeneous. However, this advantage cannot be observed in the pure homogeneous canopy.

Chapter 6 Analysis of the influence of environmental noise on vegetation indices

It is generally believed that some algebraic combination of remote sensing spectral bands can reflect some helpful information about vegetation structure and vegetation cover states, such as leaf density and distribution, leaf moisture content, tree age, lack of minerals and pests and diseases. A good vegetation index should be sensitive to these factors. On the other hand, other factors also affect spectral reflectances, such as soil properties, sunshine and atmospheric conditions, and sensor geometry parameters. A good vegetation index should be insensitive to these factors. Therefore, studies should be carried out to evaluate the effect of environmental noise on VIs. Section 6.1 compares the performances of VIs for soil noise influence. An assumption that hotspot vegetation indices can better enhance the vegetation signal and reduce soil noise is proposed and validated using simulated reflectance data. Section 6.2 studies the atmosphere effect on EVI, EVI2, NDVI and SAVI. Atmospheric absorption and scattering can lead to variations in wavelength radiance, thus affecting VIs.

Therefore, it is essential to understand the behaviours of VIs to atmosphere effect, Section 6.3 studies the effect of spectral response functions on EVI2. EVI2 is a 2-band alternative to 3-band EVI using a spectral correlation coefficient between red and blue from MODIS.

Therefore, it is necessary to evaluate the performance of EVI2 on different sensors with different spectral response functions and calibrate the spectral correlation coefficients of EVI2 for other sensors.

The chapter is partly presented in the paper:

(1) Z.

Comparison of vegetation indices for soil noise resistance

The spectral signal of the soil affects the mixed spectra of vegetation and soil, which can reduce the accuracy of monitoring vegetation LAI using remote sensing. Therefore, VIs should be designed to reduce soil response and improve vegetation response. The DART simulated data were used to analyse the resistance of the VI to soil noise based on two indicators, the Signal-to-Noise Ratio (SNR) and the dependent coefficient (T).

Two assessment criteria, the signal-to-noise ratio [START_REF] Leroy | Sun and view angle corrections on reflectances derived from NOAA/AVHRR data[END_REF] and the LAI-dependent parameter [START_REF] Gilabert | Designing a generalized soiladjusted vegetation index (GESAVI)[END_REF],

were used to assess the soil noise tolerance of VIs (multispectral VIs, VIs57.5 and hotspot signature VIs). products also include land/water, cloud mask and European Centre for Medium-range Weather Forecasts (ECMWF) data (total ozone, total water vapour and mean sea level pressure) [236].

The Level-2A products provide BOA reflectance images of the associated Level-1C products.

Each Level-2A product, therefore, also consists of a 100 km 2 strip of mapping geometry (UTM / WGS84 projection). Level-2A products from 2017 onwards can be downloaded directly from European Space Agency (ESA). Users can also use the relevant Level-1C products as input for Level-2A generation via the S2 toolbox.

Chapter 7 Conclusion and perspectives

Recently, many significant achievements have been achieved for canopy reflectance simulation and its bio-optical parameter inversion based on 3D radiative transfer models. However, many important issues remain unsolved, such as continuous-time phase simulation, inversion of vegetation leaf OPs in urban areas, and VIs saturation effects. We present a coupled growth model to achieve a 3D radiative transfer simulation capability in continuous time phase; an accurate inversion of the OPs of leaves in urban areas based on DART calibration, and an analysis of the vegetation isolines behaviours and a proposed intersection right shift phenomenon to mitigate the saturation effect of the SAVI in high FC areas. The first half of the thesis focuses on modelling and reflectance simulation of 3D vegetation scenes. The second half focuses on the inversion of vegetation bio-optical property using the 3D radiative transfer models. The main work and conclusions are summarised, the innovations of the thesis are distilled, and prospects for further research are presented.

Major conclusions

Key conclusions and related discussions include:

(1) The maize growth equation is coupled to compare the differences in simulated reflectance between homogeneous, simplified heterogeneous and real structural remote sensing scenes between the three radiative transfer models SAIL, DART and RGM. The results show that the assumption of homogeneous canopy leads to a non-negligible overestimation of nadir-observed reflectance in the NIR band. It is due to multiple scattering are mainly from vegetation. However, the proportion of vegetation seen at the nadir is minimal, and the homogeneity assumption does not well consider it.

The NIR is a vital band, and the nadir is a vital viewing angle for vegetation monitoring using remote sensing data. Therefore, the importance of considering the three-dimensional structure of the tree canopy is emphasised when using nadir observations.

(2) Inversion of leaf OPs from the mixed pixel in urban areas is achieved, simulation accuracy is assessed, and sensitivity analyses are performed. If only OP is unknown, the mean relative errors for all bands for ground, roof, water, tree and shrub are 0.013, 0.005, 0.027, 0.297 and 0.250, respectively. In the case of inaccuracy input parameters, it is 0.233, 0.507, 3.088, 0.834 and 1.256.

In addition, the inversion accuracy is affected by SZA, satellite spatial resolution, pixel shifting, scene model inaccuracy and MTF in decreasing order of magnitude. The calibration without noise interference is insensitive to the number of bands and is highly resistant to shadow interference.

Some uncontrollable external factors may significantly affect inversion accuracies, such as coregistration errors or user-defined ratios between leaf reflectance and transmittance.

( The main difference between ELSYS and other L-system includes that: Compared with the traditional L-system, the ELSYS program, the extension of ELSYS includes: (1) Using variable name rather than the single letter for predecessor and successors in the iteration rules; (2) Data structure such as integer, real, and vector are supported; (3) Supporting the parameter passing among iteration rules; (4) Supporting logical control statement (i.e., conditional and loop statement); [START_REF] Zhu | A regional maize yield hierarchical linear model combining landsat 8 vegetative indices and meteorological data: Case study in jilin province[END_REF] built-in functions and operators. These extensions make ELSYS much power than the traditional Lsystem.

All variables in ELSYS are case-sensitive, and there are two types of variables: user-defined variables and system variables. ELSYS has three types of variables and data types: integers, real numbers and 3D real vectors. Vectors are used to express position and direction in 3D space or to describe colours using Red-Green-Blue (RGB), and they are similar to arrays and can be accessed via subscripts. User-defined variables must start with a letter, and the type of variables is defined by the initial letter of the variable: i, j, k, l, m, n, I, J, K, L, M, N are integers, u, v, w, U, V, W are vectors and other variables starting with an initial letter are real numbers. The ELSYS syntax rules are described here only in a preliminary way; for system variables and more detailed syntax rules, see the literature [START_REF] Goel | A High-Level Language for L-systems and Its Applications[END_REF].

Some commonly used built-in operators are listed in Table A Input source code is written using the same process, and the code consists of four main sections of content:

(1) Global variable declaration (optional): This section consists of multiple global (variable 1, variable 2, ......, variable n) statements that declare global variables. The scope of a global variable is the entire source file.

(2) Function declaration (optional): This section consists of multiple def(function name (variable 1, variable 2, ......, variable n):= function body) statements that declare global functions.

(3) Determination of the axioms of the L-system: This section determines the axiomatic