# Evidence for the Direct Two-Photon Transition from $\psi(3686)$ to $J/\psi$

M. Ablikim1, M. N. Achasov5, D. J. Ambrose40, F. F. An1, Q. An41, Z. H. An1,
J. Z. Bai1, R. B. Ferroli18, Y. Ban27, J. Becker2, N. Berger1, M. B.
Bertani18, J. M. Bian39, E. Boger20,a, O. Bondarenko21, I. Boyko20, R. A.
Briere3, V. Bytev20, X. Cai1, A. C. Calcaterra18, G. F. Cao1, J. F. Chang1, G.
Chelkov20,a, G. Chen1, H. S. Chen1, J. C. Chen1, M. L. Chen1, S. J. Chen25, Y.
Chen1, Y. B. Chen1, H. P. Cheng14, Y. P. Chu1, D. Cronin-Hennessy39, H. L.
Dai1, J. P. Dai1, D. Dedovich20, Z. Y. Deng1, A. Denig19, I. Denysenko20,b, M.
Destefanis44, W. M. Ding29, Y. Ding23, L. Y. Dong1, M. Y. Dong1, S. X. Du47,
J. Fang1, S. S. Fang1, L. Fava44,c, F. Feldbauer2, C. Q. Feng41, C. D. Fu1, J.
L. Fu25, Y. Gao36, C. Geng41, K. Goetzen7, W. X. Gong1, W. Gradl19, M.
Greco44, M. H. Gu1, Y. T. Gu9, Y. H. Guan6, A. Q. Guo26, L. B. Guo24, Y. P.
Guo26, Y. L. Han1, X. Q. Hao1, F. A. Harris38, K. L. He1, M. He1, Z. Y. He26,
T. Held2, Y. K. Heng1, Z. L. Hou1, H. M. Hu1, J. F. Hu6, T. Hu1, B. Huang1, G.
M. Huang15, J. S. Huang12, X. T. Huang29, Y. P. Huang1, T. Hussain43, C. S.
Ji41, Q. Ji1, X. B. Ji1, X. L. Ji1, L. K. Jia1, L. L. Jiang1, X. S. Jiang1, J.
B. Jiao29, Z. Jiao14, D. P. Jin1, S. Jin1, F. F. Jing36, N. Kalantar-
Nayestanaki21, M. Kavatsyuk21, W. Kuehn37, W. Lai1, J. S. Lange37, J. K. C.
Leung35, C. H. Li1, Cheng Li41, Cui Li41, D. M. Li47, F. Li1, G. Li1, H. B.
Li1, J. C. Li1, K. Li10, Lei Li1, N. B. Li24, Q. J. Li1, S. L. Li1, W. D. Li1,
W. G. Li1, X. L. Li29, X. N. Li1, X. Q. Li26, X. R. Li28, Z. B. Li33, H.
Liang41, Y. F. Liang31, Y. T. Liang37, G. R. Liao36, X. T. Liao1, B. J. Liu1,
B. J. Liu34, C. L. Liu3, C. X. Liu1, C. Y. Liu1, F. H. Liu30, Fang Liu1, Feng
Liu15, H. Liu1, H. B. Liu6, H. H. Liu13, H. M. Liu1, H. W. Liu1, J. P. Liu45,
Kun Liu27, Kai Liu6, K. Y. Liu23, P. L. Liu29, S. B. Liu41, X. Liu22, X. H.
Liu1, Y. B. Liu26, Y. Liu1, Z. A. Liu1, Zhiqiang Liu1, Zhiqing Liu1, H.
Loehner21, G. R. Lu12, H. J. Lu14, J. G. Lu1, Q. W. Lu30, X. R. Lu6, Y. P.
Lu1, C. L. Luo24, M. X. Luo46, T. Luo38, X. L. Luo1, M. Lv1, C. L. Ma6, F. C.
Ma23, H. L. Ma1, Q. M. Ma1, S. Ma1, T. Ma1, X. Y. Ma1, Y. Ma11, F. E. Maas11,
M. Maggiora44, Q. A. Malik43, H. Mao1, Y. J. Mao27, Z. P. Mao1, J. G.
Messchendorp21, J. Min1, T. J. Min1, R. E. Mitchell17, X. H. Mo1, C. Morales
Morales11, C. Motzko2, N. Yu. Muchnoi5, Y. Nefedov20, C. Nicholson6, I. B..
Nikolaev5, Z. Ning1, S. L. Olsen28, Q. Ouyang1, S. P. Pacetti18,d, J. W.
Park28, M. Pelizaeus38, K. Peters7, J. L. Ping24, R. G. Ping1, R. Poling39, E.
Prencipe19, C. S. J. Pun35, M. Qi25, S. Qian1, C. F. Qiao6, X. S. Qin1, Y.
Qin27, Z. H. Qin1, J. F. Qiu1, K. H. Rashid43, G. Rong1, X. D. Ruan9, A.
Sarantsev20,e, J. Schulze2, M. Shao41, C. P. Shen38,f, X. Y. Shen1, H. Y.
Sheng1, M. R. Shepherd17, X. Y. Song1, S. Spataro44, B. Spruck37, D. H. Sun1,
G. X. Sun1, J. F. Sun12, S. S. Sun1, X. D. Sun1, Y. J. Sun41, Y. Z. Sun1, Z.
J. Sun1, Z. T. Sun41, C. J. Tang31, X. Tang1, E. H. Thorndike40, H. L. Tian1,
D. Toth39, M. U. Ulrich37, G. S. Varner38, B. Wang9, B. Q. Wang27, K. Wang1,
L. L. Wang4, L. S. Wang1, M. Wang29, P. Wang1, P. L. Wang1, Q. Wang1, Q. J.
Wang1, S. G. Wang27, X. F. Wang12, X. L. Wang41, Y. D. Wang41, Y. F. Wang1, Y.
Q. Wang29, Z. Wang1, Z. G. Wang1, Z. Y. Wang1, D. H. Wei8, P. Weidenkaff19, Q.
G. Wen41, S. P. Wen1, M. W. Werner37, U. Wiedner2, L. H. Wu1, N. Wu1, S. X.
Wu41, W. Wu26, Z. Wu1, L. G. Xia36, Z. J. Xiao24, Y. G. Xie1, Q. L. Xiu1, G.
F. Xu1, G. M. Xu27, H. Xu1, Q. J. Xu10, X. P. Xu32, Y. Xu26, Z. R. Xu41, F.
Xue15, Z. Xue1, L. Yan41, W. B. Yan41, Y. H. Yan16, H. X. Yang1, T. Yang9, Y.
Yang15, Y. X. Yang8, H. Ye1, M. Ye1, M. H. Ye4, B. X. Yu1, C. X. Yu26, J. S.
Yu22, S. P. Yu29, C. Z. Yuan1, W. L. Yuan24, Y. Yuan1, A. A. Zafar43, A. Z.
Zallo18, Y. Zeng16, B. X. Zhang1, B. Y. Zhang1, C. C. Zhang1, D. H. Zhang1, H.
H. Zhang33, H. Y. Zhang1, J. Zhang24, J. G. Zhang12, J. Q. Zhang1, J. W.
Zhang1, J. Y. Zhang1, J. Z. Zhang1, L. Zhang25, S. H. Zhang1, T. R. Zhang24,
X. J. Zhang1, X. Y. Zhang29, Y. Zhang1, Y. H. Zhang1, Y. S. Zhang9, Z. P.
Zhang41, Z. Y. Zhang45, G. Zhao1, H. S. Zhao1, J. W. Zhao1, K. X. Zhao24, Lei
Zhao41, Ling Zhao1, M. G. Zhao26, Q. Zhao1, S. J. Zhao47, T. C. Zhao1, X. H.
Zhao25, Y. B. Zhao1, Z. G. Zhao41, A. Zhemchugov20,a, B. Zheng42, J. P.
Zheng1, Y. H. Zheng6, Z. P. Zheng1, B. Zhong1, J. Zhong2, L. Zhou1, X. K.
Zhou6, X. R. Zhou41, C. Zhu1, K. Zhu1, K. J. Zhu1, S. H. Zhu1, X. L. Zhu36, X.
W. Zhu1, Y. M. Zhu26, Y. S. Zhu1, Z. A. Zhu1, J. Zhuang1, B. S. Zou1, J. H.
Zou1, J. X. Zuo1 (BESIII Collaboration) 1 Institute of High Energy Physics,
Beijing 100049, P. R. China
2 Bochum Ruhr-University, 44780 Bochum, Germany
3 Carnegie Mellon University, Pittsburgh, PA 15213, USA
4 China Center of Advanced Science and Technology, Beijing 100190, P. R. China
5 G.I. Budker Institute of Nuclear Physics SB RAS (BINP), Novosibirsk 630090,
Russia
6 Graduate University of Chinese Academy of Sciences, Beijing 100049, P. R.
China
7 GSI Helmholtzcentre for Heavy Ion Research GmbH, D-64291 Darmstadt, Germany
8 Guangxi Normal University, Guilin 541004, P. R. China
9 GuangXi University, Nanning 530004,P.R.China
10 Hangzhou Normal University, Hangzhou 310036, P. R. China
11 Helmholtz Institute Mainz, J.J. Becherweg 45,D 55099 Mainz,Germany
12 Henan Normal University, Xinxiang 453007, P. R. China
13 Henan University of Science and Technology, Luoyang 471003, P. R. China
14 Huangshan College, Huangshan 245000, P. R. China
15 Huazhong Normal University, Wuhan 430079, P. R. China
16 Hunan University, Changsha 410082, P. R. China
17 Indiana University, Bloomington, Indiana 47405, USA
18 INFN Laboratori Nazionali di Frascati , Frascati, Italy
19 Johannes Gutenberg University of Mainz, Johann-Joachim-Becher-Weg 45, 55099
Mainz, Germany
20 Joint Institute for Nuclear Research, 141980 Dubna, Russia
21 KVI/University of Groningen, 9747 AA Groningen, The Netherlands
22 Lanzhou University, Lanzhou 730000, P. R. China
23 Liaoning University, Shenyang 110036, P. R. China
24 Nanjing Normal University, Nanjing 210046, P. R. China
25 Nanjing University, Nanjing 210093, P. R. China
26 Nankai University, Tianjin 300071, P. R. China
27 Peking University, Beijing 100871, P. R. China
28 Seoul National University, Seoul, 151-747 Korea
29 Shandong University, Jinan 250100, P. R. China
30 Shanxi University, Taiyuan 030006, P. R. China
31 Sichuan University, Chengdu 610064, P. R. China
32 Soochow University, Suzhou 215006, China
33 Sun Yat-Sen University, Guangzhou 510275, P. R. China
34 The Chinese University of Hong Kong, Shatin, N.T., Hong Kong.
35 The University of Hong Kong, Pokfulam, Hong Kong
36 Tsinghua University, Beijing 100084, P. R. China
37 Universitaet Giessen, 35392 Giessen, Germany
38 University of Hawaii, Honolulu, Hawaii 96822, USA
39 University of Minnesota, Minneapolis, MN 55455, USA
40 University of Rochester, Rochester, New York 14627, USA
41 University of Science and Technology of China, Hefei 230026, P. R. China
42 University of South China, Hengyang 421001, P. R. China
43 University of the Punjab, Lahore-54590, Pakistan
44 University of Turin and INFN, Turin, Italy
45 Wuhan University, Wuhan 430072, P. R. China
46 Zhejiang University, Hangzhou 310027, P. R. China
47 Zhengzhou University, Zhengzhou 450001, P. R. China
a also at the Moscow Institute of Physics and Technology, Moscow, Russia
b on leave from the Bogolyubov Institute for Theoretical Physics, Kiev,
Ukraine
c University of Piemonte Orientale and INFN (Turin)
d Currently at INFN and University of Perugia, I-06100 Perugia, Italy
e also at the PNPI, Gatchina, Russia
f now at Nagoya University, Nagoya, Japan

###### Abstract

The two-photon transition $\psi(3686)\to\gamma\gamma J/\psi$ is studied in a
sample of 106 million $\psi(3686)$ decays collected by the BESIII detector.
The branching fraction is measured to be $(3.1\pm
0.6(\mathrm{stat})^{+0.8}_{-1.0}(\mathrm{syst}))\times 10^{-4}$ using
$J/\psi\to e^{+}e^{-}$ and $J/\psi\to\mu^{+}\mu^{-}$ decays, and its upper
limit is estimated to be $4.5\times 10^{-4}$ at the 90% conference level. This
work represents the first measurement of a two-photon transition among
charmonium states. The orientation of the $\psi(3686)$ decay plane and the
$J/\psi$ polarization in this decay are also studied. In addition, the product
branching fractions of sequential $E1$ transitions
$\psi(3686)\to\gamma\chi_{cJ},\chi_{cJ}\to\gamma J/\psi(J=0,1,2)$ are
reported.

###### pacs:

14.40.Pq, 13.20.Gd, 14.40.-n

The XYZ Brambilla:2010cs  particles, which do not fit potential model
expectations in QCD theory, have been a key challenge to the QCD description
of charmonium-like states Asner:2008nq . To fully understand those states, it
is necessary to consider the coupling of a charmonium state to a $D\bar{D}$
meson pair. These coupled-channel effects, which also play an important role
in the charmonium transitions of low lying states (_i_._e_., from $\psi(3686)$
to $J/\psi$), are especially relevant for the radiative transition processes
Eichten:2004uh . In the well-known electric dipole transitions, the strength
of coupled-channel effects will likely be hard to establish, since the
accompanying relativistic corrections may be more important Li:2009zu .
However, the two-photon transition $\psi(3686)\to\gamma\gamma J/\psi$ is more
sensitive to the coupled-channel effect and thus provides a unique opportunity
to investigate these issues He:2010pb .

Two-photon spectroscopy has been a very powerful tool for the study of the
excitation spectra of a variety of systems with a wide range of sizes, such as
molecules, atomic hydrogen and positronium Pachucki:1996jw . Studying the
analogous process in quarkonium states is a natural extension of this work, in
order to gain insight into non-perturbative QCD phenomena. But so far, two-
photon transitions in quarkonia have eluded experimental observation
Bai:2004cg ; Adam:2005uh ; :2008kb . For example, in a study of
$\psi(3686)\to\gamma\chi_{cJ}(J=0,1,2)$ reported by CLEO-c :2008kb , the upper
limit for $\mathcal{B}(\psi(3686)\to\gamma\gamma J/\psi)$ was estimated to be
$1\times 10^{-3}$.

This Letter presents the first evidence for the two-photon transition
$\psi(3686)\to\gamma\gamma J/\psi$, as well as studies of the orientation of
the $\psi(3686)$ decay plane and the $J/\psi$ polarization in the decay. The
branching fractions of double $E1$ transitions $\psi(3686)\to\gamma(\gamma
J/\psi)_{\chi_{cJ}}$ through $\chi_{cJ}$ intermediate states are also
reported. The data analyzed were obtained by the BESIII experiment :2009vd
viewing electron-positron collisions at the BEPCII collider. An integrated
luminosity of 156.4 $\rm{pb}^{-1}$ was obtained at a center-of-mass energy
$\sqrt{s}=M(\psi(3686))=3.686\,\mathrm{GeV}$. The number of $\psi(3686)$
decays in this sample is estimated to be $(1.06\pm 0.04)\times 10^{8}$
Ablikim:2010zn . In addition, 42.6 $\rm{pb}^{-1}$ of continuum data were taken
below the $\psi(3686)$, at $\sqrt{s}=3.65\,\mathrm{GeV}$, to evaluate the
potential backgrounds from non-resonant events.

The upgraded BEPCII Zhang:2010zz  at Beijing is a two-ring electron-positron
collider. The BESIII detector :2009vd  is an approximately cylindrically
symmetric detector which covers $93\%$ of the solid angle around the collision
point. In order of increasing distance from the interaction point, the sub-
detectors include a 43-layer main wire drift chamber (MDC), a time-of-flight
(TOF) system with two layers in the barrel region and one layer for each end-
cap, and a 6240 cell CsI(Tl) crystal electro-magnetic calorimeter (EMC) with
both barrel and endcap sections. The barrel components reside within a
superconducting solenoid magnet providing a 1.0 T magnetic field aligned with
the beam axis. Finally, there is a muon chamber consisting of nine layers of
resistive plate chambers within the return yoke of the magnet. The momentum
resolution for charged tracks in the MDC is $0.5\%$ for transverse momenta of
$1\,\mathrm{GeV}/c$. The energy resolution for showers in the EMC is $2.5\%$
for $1\,\mathrm{GeV}$ photons.

Figure 1: Up (a): distributions of $M_{\gamma\gamma-\mathrm{recoil}}$ in data
(points) and in the combined dataset (solid line) of MC simulation of
$\psi(3686)$ decays (shaded histogram) and continuum backgrounds (dashed
line), before the KF is applied. The arrows indicate the window to select a
$J/\psi$ candidate; Down: scatter plots of
$M_{\gamma_{\mathrm{sm}}-\mathrm{recoil}}$ versus $M_{\gamma\gamma}$ for the
$\gamma\gamma e^{+}e^{-}$ channel, in data (b), continuum data (c), MC
simulated signal (d), after applying the KF constrain and the
$M_{\gamma\gamma-\mathrm{recoil}}$ window. The corresponding plots for the
$\gamma\gamma\mu^{+}\mu^{-}$ channel are very similar.

This work studies $\psi(3686)\to\gamma\gamma J/\psi$ followed by
$J/\psi\to\ell^{+}\ell^{-}$ ($\ell$ denotes $e$ or $\mu$), which is referred
to as the signal process. Events selected contain exactly two oppositely
charged good tracks in the MDC tracking system, corresponding to the dilepton
from $J/\psi$ decay. The requirements to judge a track as good include
$|\cos\theta|<0.93$ ($\theta$ is the polar angle with respect to the beam
direction), and the minimum distance of approach between the track and the
production vertex less than 10 cm along the beam axis and less than 1 cm
projected in the perpendicular plane. The lepton is identified with the ratio
of EMC shower energy to MDC track momentum, $E/p$, which must be larger than
0.7 for an electron, or smaller than 0.6 for a muon. To suppress non-$J/\psi$
decay leptons, we require the momentum of each lepton to be larger than
$0.8\,\mathrm{GeV}/c$. A vertex fit (VF) constrains the production vertex,
which is updated run-by-run, and the tracks of the dilepton candidates to a
common vertex; only events with $\chi^{2}_{\mathrm{VF}}/\mathrm{d.o.f.}<20$
are accepted.

Reconstructed EMC showers unmatched to either charged track and with an energy
larger than $25\,\mathrm{MeV}$ in the barrel region ($|\cos\theta|<0.80$) or
larger than $50\,\mathrm{MeV}$ in the end-caps ($0.86<|\cos\theta|<0.92$) are
used as photon candidates. To reject bremsstrahlung photons, showers matching
the initial momentum of either lepton within $10^{\circ}$ are also discarded.
Showers from noise, not originating from the beam collision, are suppressed by
requiring the EMC cluster time to lie within a 700 ns window near the event
start time.

Events are required to have only two photon candidates. A kinematic fit (KF)
constrains the vertexed dilepton to the nominal mass of the intermediate
$J/\psi$, and the resulting $J/\psi$ and photon candidates to the known
initial four-momentum of the $\psi(3686)$. The KF fit quality
$\chi^{2}_{\mathrm{KF}}$ is required to be
$\chi^{2}_{\mathrm{KF}}/\mathrm{d.o.f.}<12$. For convenience, we use
$\gamma_{\mathrm{lg}}$ ($\gamma_{\mathrm{sm}}$) to denote the larger (smaller)
energy photon. As indicated in Fig. 1(a), $J/\psi$ candidates are identified
with the requirement that the recoil mass of the two photons,
$M_{\gamma\gamma-\mathrm{recoil}}$, is within
$(3.08,3.14)\,\mathrm{GeV}/c^{2}$.

Scatter plots of recoiling mass $M_{\gamma_{\mathrm{sm}}-\mathrm{recoil}}$
from the lower energy photon $\gamma_{\mathrm{sm}}$ versus invariant mass of
two photons $M_{\gamma\gamma}$ are shown in Fig. 1, where clear resonance
bands are seen from the decays $\psi(3686)\to\gamma\chi_{cJ}(J=0,1,2)$ (three
horizontal bands) and $\psi(3686)\to\pi^{0}(\eta)J/\psi$ (two vertical bands).
As indicated in Fig. 1(c), the continuum backgrounds are most dominant at the
top of the plots, of which the primary sources include the bhabha scattering,
the dimuon process and the ISR production of $J/\psi$. These backgrounds are
excluded by discarding events with
$M_{\gamma_{\mathrm{sm}}-\mathrm{recoil}}>3.6\,\mathrm{GeV}/c^{2}$. To
suppress backgrounds from $\psi(3686)\rightarrow\pi^{0}(\eta)J/\psi$, the
diphoton invariant mass $M_{\gamma\gamma}$ is required to be larger than
$0.15\,\mathrm{GeV}/c^{2}$ and the recoil momentum of the diphoton must be
larger than $0.25\,\mathrm{GeV}/c$.

Figure 2: (color online) Plot a: unbinned maximum likelihood fit to the
distribution of $M_{\gamma_{\mathrm{sm}}-\mathrm{recoil}}$ in data with
combination of the two $J/\psi$-decay modes. Thick lines are the sum of the
fitting models and long-dashed lines are the $\chi_{cJ}$ shapes. Short-dashed
lines represent the two-photon signal processes. Shaded histograms are
$\psi(3686)$-decay backgrounds (yellow) and non-$\psi(3686)$ backgrounds
(green), with the fixed amplitude and shape taken from MC simulation and
continuum data. Plot b: the number of standard deviations, $n_{\sigma}$, of
data points from the fitted curves in plot a. The rates of the signal process
and sequential $\chi_{cJ}$ processes are derived from these fits. Plot c:
distributions of $M_{\gamma\gamma-\mathrm{recoil}}$ in data (signals and known
backgrounds) with the kinematic requirement
$3.44\,\mathrm{GeV}/c^{2}<M_{\gamma_{\mathrm{sm}}-\mathrm{recoil}}<3.48\,\mathrm{GeV}/c^{2}$
and with the removal of $\chi^{2}_{\mathrm{KF}}$ and
$M_{\gamma\gamma-\mathrm{recoil}}$ restrictions. Plot d: stacked histograms of
the three $\chi_{cJ}$ components in plot c.

Monte Carlo (MC) simulations of $\psi(3686)$ decays are used to understand the
backgrounds and also to estimate the detection efficiency. At BESIII, the
simulation includes the beam energy spread and treats the initial-state
radiation with KKMC Jadach:2000ir . Specific decay modes from the PDG pdg2012
are modeled with EVTGEN ref:bes3gen , and the unknown decay modes with
Lundcharm Chen:2000tv . The detector response is described using GEANT4
ref:geant4 . For the $\psi(3686)\to\gamma\gamma J/\psi$ channel, the momenta
of decay particles are simulated according to the measured polarization
structure in this work. Generic $\psi(3686)$ decay samples serve for
understanding the background channels; dominant backgrounds were generated
with high statistics. Angular distributions of the cascade $E1$ transitions
$\psi(3686)\to\gamma\chi_{cJ}\to\gamma\gamma J/\psi$ are assumed to follow the
formulae in Ref. Karl:1975jp . Note that the $\chi_{cJ}$ line shapes were
simulated with the Breit-Wigner distributions weighted with
$E^{3}_{\gamma^{*}_{1}}E^{3}_{\gamma^{*}_{2}}$ to account for the double $E1$
transitions, and extended out to $\pm 200\,\mathrm{MeV}/c^{2}$ away from the
nominal masses, using masses and widths in PDG pdg2012 . Here,
$E_{\gamma^{*}_{1}}$($E_{\gamma^{*}_{2}}$) is the energy of the radiative
photon $\gamma^{*}_{1}$($\gamma^{*}_{2}$) in the rest frame of the mother
particle $\psi(3686)$($\chi_{cJ}$).

The yield of the signal process $\psi(3686)\to\gamma\gamma J/\psi$, together
with those of the cascade $E1$ transition processes, is estimated by a global
fit to the spectrum of $M_{\gamma_{\mathrm{sm}}-\mathrm{recoil}}$. The fit
results are shown in Fig. 2. The shape and magnitude of $\psi(3686)$-decay
backgrounds were fixed based on MC simulation. Non-$\psi(3686)$ decay
backgrounds are estimated in continuum data, scaling by luminosity and the
$1/s$ dependence of the cross sections. This scaling is verified by the good
description of the $J/\psi$ backgrounds in the
$M_{\gamma\gamma-\mathrm{recoil}}$ distribution shown in Fig. 1(a). The
distributions of the signal process and the cascade $E1$ process are taken
from the reconstructed shapes in MC simulation of the modes and smeared with
an asymmetric Gaussian with free parameters, which is used to compensate for
the difference in line shape between MC and data. By taking the MC shape,
detector resolution and wrong assignment of the $E1$ photon are taken into
account. The quality of goodness-of-fit test,
$\chi^{2}/$d.o.f.$=108.0/94=1.15$ in the $\gamma\gamma e^{+}e^{-}$ mode and
$124.8/94=1.33$ in the $\gamma\gamma\mu^{+}\mu^{-}$ mode. The observed signal
yields are given in Table 1. The $\psi(3686)\to\gamma\gamma J/\psi$ transition
is observed with a statistical significance of 6.6$\sigma$, as determined by
the ratio of the maximum likelihood value and the likelihood value for a fit
with null-signal hypothesis. When the systematic uncertainties are taken into
account with the assumption of Gaussian distributions, the significance is
evaluated to be 3.8$\sigma$, which corresponds to a probability of a
background fluctuation to the observed signal yield of $7.2\times 10^{-5}$.
The upper limit for $\mathcal{B}(\psi(3686)\to\gamma\gamma J/\psi)$ is
estimated to be $4.5\times 10^{-4}$ at the 90% confidence level including
systematic uncertainties.

In calculating $\mathcal{B}(\gamma\gamma J/\psi)$, a correction factor is
included due to the interferences among $\chi_{cJ}$ states. This effect was
checked by the variations of the observed signals in the global fit with
inclusion of a floating interference component, which is modeled by the
detector-smeared shape of a theoretical calculation He:2010pb . It is found
that relative changes on the signal yields are negative with lower bound of
$-10\%$. Hence, a correction factor $0.95$ is assigned and $5\%$ is taken as
systematic uncertainty.

A cross-check on our procedures is performed with the
$M_{\gamma\gamma-\mathrm{recoil}}$ spectrum for the events in the region
$3.44\,\mathrm{GeV}/c^{2}<M_{\gamma_{\mathrm{sm}}-\mathrm{recoil}}<3.48\,\mathrm{GeV}/c^{2}$
without restrictions on $\chi^{2}_{\mathrm{KF}}$ and
$M_{\gamma\gamma-\mathrm{recoil}}$, as shown in Fig. 2(c). An excess of data
above known backgrounds can be seen around the $J/\psi$ nominal mass, which is
expected from the sought-after two-photon process. With the inclusion of the
estimated yields of the signal process, the excess is well understood. The
high-mass peak above the $J/\psi$ peak comes from the backgrounds of
$\psi(3686)\to\pi^{0}\pi^{0}J/\psi$ decays. This satellite peak can be well
described in MC simulation. In Fig. 2(d), the three $\chi_{cJ}$ tails show
distinguishable distributions; the small left bump is from the $\chi_{c1}$
tail, while the $\chi_{c0}$ tail is dominant at the right side. The
distribution in data in Fig. 2(c) can only be well described by the simulated
$\chi_{cJ}$ shapes.

Table 1: For different channels: the number of observed signals $n_{e}$ ($n_{\mu}$) and detection efficiency $\epsilon_{e}$ ($\epsilon_{\mu}$) in $\gamma\gamma e^{+}e^{-}$ ($\gamma\gamma\mu^{+}\mu^{-}$) mode; the absolute branching fractions. On the bottom, the relative branching fractions $R_{MN}\equiv\mathcal{B}_{\chi_{cM}}/\mathcal{B}_{\chi_{cN}}$, where $\mathcal{B}_{\chi_{cJ}}\equiv\mathcal{B}(\psi(3686)\to\gamma(\gamma J/\psi)_{\chi_{cJ}})$ are listed. Here the first errors are statistical and the second are systematic. Channels | $n_{e}$ |  $\epsilon_{e}$(%) | $n_{\mu}$ |  $\epsilon_{\mu}$(%) |  $\mathcal{B}$($\times 10^{-4}$)
---|---|---|---|---|---
$\gamma\gamma J/\psi$ |  564$\pm$116 | 22.4 |  536$\pm$128 | 30.0 | $3.1\pm 0.6^{+0.8}_{-1.0}$
$\gamma(\gamma J/\psi)_{\chi_{c0}}$ |  1801$\pm$60 | 19.3 |  2491$\pm$69 | 26.0 | $15.1\pm 0.3\pm 1.0$
$\gamma(\gamma J/\psi)_{\chi_{c1}}$ |  59953$\pm$253 | 28.5 |  81922$\pm$295 | 38.2 | $337.7\pm 0.9\pm 18.3$
$\gamma(\gamma J/\psi)_{\chi_{c2}}$ |  32171$\pm$187 | 27.5 |  44136$\pm$219 | 37.1 | $187.4\pm 0.7\pm 10.2$
$R_{21}\equiv\frac{\mathcal{B}_{\chi_{c2}}}{\mathcal{B}_{\chi_{c1}}}(\%)$ | $R_{01}\equiv\frac{\mathcal{B}_{\chi_{c0}}}{\mathcal{B}_{\chi_{c1}}}(\%)$ | $R_{02}\equiv\frac{\mathcal{B}_{\chi_{c0}}}{\mathcal{B}_{\chi_{c2}}}(\%)$
$55.47\pm 0.26\pm 0.11$ | $4.45\pm 0.09\pm 0.18$ | $8.03\pm 0.17\pm 0.33$

The angle of the normal axis of the $\psi(3686)$ decay plane with respect to
the $\psi(3686)$ polarization vector (aligned to the beam axis), $\beta$, can
be determined in our data. The event rate may be expressed, to leading order,
as $\frac{d\,N}{d\,\cos\beta}\propto 1+a\cos^{2}\beta$. The measurement was
carried out in the rest frame of the $\psi(3686)$ and the decay plane of the
$\psi(3686)$ was determined with the momenta of the two decay particles
$J/\psi$ and $\gamma_{\mathrm{lg}}$. The signal yields in each angular bin
were extracted by the global fit to the corresponding dataset following the
aforementioned procedure. After correction of the extracted signal yields with
the detection efficiency, Fig. 3(a) shows the fit to the distribution of
$|\cos\beta|$ for the sum of the two dilepton modes; we obtain $a=0.53\pm
0.68$.

Table 2: Summary of the systematic uncertainties on the measurement of $\mathcal{B}_{\mathrm{sig}}$ of $\gamma\gamma J/\psi$ signal process, $\mathcal{B}_{\chi_{cJ}}$ for $\chi_{cJ}$ intermediate processes and the relative branching fractions $R_{MN}$, following the notation convention in Table 1. The tot systematic uncertainty is the square root of the sum. A dash (–) means the uncertainty is negligible. Values inside the parentheses are for the $\gamma\gamma\mu^{+}\mu^{-}$ mode, while values outside are for the $\gamma\gamma e^{+}e^{-}$ mode. Numbers without brackets represent uncertainties that are common to both modes. systematic uncertainty(%) | ${\footnotesize\mathcal{B}_{\mathrm{sig}}}$ | $\mathcal{B}_{\chi_{c0}}$ | $\mathcal{B}_{\chi_{c1}}$ | $\mathcal{B}_{\chi_{c2}}$ | ${\footnotesize R_{01}}$ | ${\footnotesize R_{02}}$ |  ${\footnotesize R_{21}}$
---|---|---|---|---|---|---|---
lepton track | $2(2)$ | $2(2)$ | $2(2)$ |  $2(2)$ |  |  |
photon shower | 2 | 2 | 2 |  2  |  |  |
number of photons | 10(3) | 1(1) | 1(1) | 1(1) | 2(–) | 2(–) |  –(–)
KF, $\chi^{2}_{\mathrm{KF}}$ requirement | 2(2) | 2(2) | 2(2) |  2(2)  |  |  |
$\chi_{cJ}$ widths | ${}^{+15}_{-25}$ | 3 | – | – | 4 | 4 |  0.2
$M_{\gamma_{\mathrm{sm}}-\mathrm{recoil}}$ resolution | 4(5) | –(–) | –(–) | –(–) | –(–) | –(–) |  –(–)
other background | 4(2) | 1(1) | –(–) | –(–) | –(–) | –(–) |  –(–)
$\chi_{cJ}$ interference | 5 | 1 | – | – | 1 | 1 |  –
fitting | 8(5) | 1(1) | –(–) | –(–) | 1(1) | 1(1) |  –(–)
spin-structure | 20 | 1 | – | – | 1 | 1 |  –
number of $\psi(3686)$ | 4 | 4 | 4 |  4  |  |  |
$\mathcal{B}({\tiny J/\psi\to\ell^{+}\ell^{-}})$ | 1 | 1 | 1 |  1  |  |  |
total | correlated | 14(8) | 3(3) | 3(3) | 3(3) | 2(1) | 2(1) |  –(–)
uncorrelated | ${}^{+25}_{-33}$ | 6 | 5 | 5 | 4 | 4 | 0.2

Figure 3: (a) The corrected distribution of the normal angle $\beta$ of the
$\psi(3686)$ decay plane, and (b) the helicity angle $\theta_{\ell}$ of
$J/\psi$ decays. The curves in (a) and (b) present the fits of functions
$P_{0}(1+a\cos^{2}\beta)$ and $P_{0}(1+\alpha\cos^{2}\theta_{\ell})$,
respectively.

The polarization of $J/\psi$ should be helpful in understanding the mechanism
of the transition process Artoisenet:2007xi . The polarization parameter
$\alpha$ can be evaluated from the angular distribution of the decay rate,
expressed as $\frac{d\,N}{d\,\cos\theta_{\ell}}\propto
1+\alpha\cos^{2}\theta_{\ell}$. Here,
$\alpha=\frac{\Gamma_{\mathrm{T}}-2\Gamma_{\mathrm{L}}}{\Gamma_{\mathrm{T}}+2\Gamma_{\mathrm{L}}}$
(with $\Gamma_{\mathrm{T}}$ and $\Gamma_{\mathrm{L}}$ being the transversely
and longitudinally polarized decay widths, respectively) and the helicity
angle $\theta_{\ell}$ is defined as the angle of the lepton in the $J/\psi$
rest frame with respect to the $J/\psi$ boost direction in the laboratory
frame. For fully transverse (longitudinal) polarization, $\alpha=+1(-1)$.
Figure 3(b) shows the distribution of $|\cos\theta_{\ell}|$ for the sum of the
two dilepton modes, after correcting the signal yields for the detection
efficiency and the lepton final state radiation effect. Our fit result is
$\alpha=0.08\pm 0.42$.

Sources of systematic errors on the measurement of branching fractions are
listed in Table 2. Uncertainties associated with the efficiency of the lepton
tracking and identification were studied with a selected control sample of
$\psi(3686)\to\pi^{+}\pi^{-}(\ell^{+}\ell^{-})_{J/\psi}$. The potential bias
due to limiting the maximum number of photon candidates was studied by varying
the limit. Throughout the photon energy region in this work, detection and
energy resolution of photon are well-modeled within a 1% uncertainty
Ablikim:2010zn ; Collaboration:2010rc . Detector resolution of the $\chi_{cJ}$
tails is taken into account up to the accuracy of the MC simulation. The
corresponding systematic uncertainty is evaluated by scanning the sizes of
smearing parameters within their errors. For the signal process, the dominant
uncertainties are from the description of $\chi_{cJ}$ line shapes, _e_._g_.,
$\chi_{cJ}$ widths. The sensitivity to the $\chi_{cJ}$ widths is studied by a
comparison of the signal yields based on different settings of the $\chi_{cJ}$
widths in modeling the $\chi_{cJ}$ resonances within the current world-average
uncertainties. Relative changes of the signal detection efficiencies are
assigned as 20%, by varying the input spin-structure within the measurement
uncertainties and weighting the efficiencies in the Dalitz-like plot of Fig.
1(d).

Many sources of systematic uncertainties in Table 2 cancel out when extracting
the $\psi(3686)$ decay plane parameter $a$ and the $J/\psi$ polarization
parameter $\alpha$. The quadrature sum of the remaining systematic
uncertainties are ${}^{+0.68}_{-0.27}$ and ${}^{+0.07}_{-0.14}$ for $a$ and
$\alpha$, respectively.

To summarize, the first measurement of the two-photon transition
$\psi(3686)\to\gamma\gamma J/\psi$ was carried out at the BESIII experiment.
The branching fraction is given in Table 1, as well as those of the cascade
$E1$ transitions. The measurement of the two-photon process is consistent with
the upper limit obtained in Ref. :2008kb . The results for the signal process
are presented without considering the possible interferences between the
direct transition and the $\chi_{cJ}$ states, due to a lack of theoretical
guidance. The distribution of the normal angle of the $\psi(3686)$ decay plane
is characterized by the parameter $a=0.53\pm
0.68(\mathrm{stat})^{+0.68}_{-0.27}(\mathrm{syst})$, indicating a preference
for a positive value. The $J/\psi$ polarization parameter $\alpha$ was
evaluated as $0.08\pm 0.42(\mathrm{stat})^{+0.07}_{-0.14}(\mathrm{syst})$,
demonstrating a competitive mixing of the longitudinal and transverse
components. These results will help constrain the strength of the coupled-
channel effect in future theoretical calculation. The reported branching
fractions $\mathcal{B}(\psi(3686)\to\gamma(\gamma J/\psi)_{\chi_{cJ}})$ are
consistent with the world average results pdg2012 . The reported relative
branching fractions of $\mathcal{B}_{\chi_{cJ}}$ are obtained with the world’s
best precision.

X.-R. Lu thanks Zhi-Guo He and De-Shan Yang for useful suggestions. We thank
the staff of BEPCII and the computing center for their hard efforts. We are
grateful for support from our institutes and universities and from the
agencies: Ministry of Science and Technology of China, National Natural
Science Foundation of China, Chinese Academy of Sciences, Istituto Nazionale
di Fisica Nucleare, U.S. Department of Energy, U.S. National Science
Foundation, and National Research Foundation of Korea.

## References

  * (1) N. Brambilla et al., Eur. Phys. J. C 71, 1534 (2011).
  * (2) D. M. Asner et al., Int. J. Mod. Phys. A 24, Suppl. 1, (2009).
  * (3) E. J. Eichten, K. Lane and C. Quigg, Phys. Rev. D 69, 094019 (2004); X. Liu, Phys. Lett. B 680, 137 (2009); T. Barnes, AIP Conf. Proc. 1257, 11 (2010).
  * (4) B. Q. Li and K. T. Chao, Phys. Rev. D 79, 094004 (2009).
  * (5) Z. G. He et al., Phys. Rev. D 83, 054028 (2011).
  * (6) K. Pachucki et al., J. Phys. B 29, 177 (1996); A. Quattropani and F. Bassani, Phys. Rev. Lett. 50, 1258 (1983).
  * (7) J. Z. Bai et al. [BES Collaboration], Phys. Rev. D 70, 012006 (2004); M. Ablikim et al. [BES Collaboration], Phys. Rev. D 71, 092002 (2005).
  * (8) N. E. Adam et al. [CLEO Collaboration], Phys. Rev. Lett. 94, 232002 (2005).
  * (9) H. Mendez et al. [CLEO Collaboration], Phys. Rev. D 78, 011102(R) (2008).
  * (10) M. Ablikim et al. [BESIII Collaboration], Nucl. Instrum. Meth. A 614, 345 (2010).
  * (11) M. Ablikim et al. [BESIII Collaboration], Phys. Rev. D 81, 052005 (2010).
  * (12) C. Zhang, Sci. China G 53, 2084 (2010).
  * (13) S. Jadach et al., Phys. Rev. D 63, 113009 (2001).
  * (14) J. Beringer et al. [Particle Data Group], Phys. Rev. D86, 010001 (2012).
  * (15) R. G. Ping, Chinese Phys. C 32, 599 (2008).
  * (16) J. C. Chen et al., Phys. Rev. D 62, 034003 (2000).
  * (17) S. Agostinelli et al. [GEANT Collaboration], Nucl. Instrum. Meth. A 506, 250 (2003); J. Allison et al., IEEE Trans. Nucl. Sci. 53, 270 (2006).
  * (18) G. Karl, S. Meshkov and J. L. Rosner, Phys. Rev. D 13, 1203 (1976); M. A. Doncheski, H. Grotch and K. J. Sebastian, Phys. Rev. D 42, 2293 (1990).
  * (19) P. Artoisenet, J. P. Lansberg and F. Maltoni, Phys. Lett. B 653, 60 (2007).
  * (20) M. Ablikim et al. [BESIII Collaboration], Phys. Rev. Lett. 104, 132002 (2010).


# Comparing remote sensing-based forest biomass mapping approaches using new
forest inventory plots in contrasting forests in northeastern and southwestern
China

Wenquan Dong  Edward T.A. Mitchard  Yuwei Chen  Man Chen  Congfeng Cao  Peilun
Hu  Cong Xu  Steven Hancock

###### Abstract

Large-scale high spatial resolution aboveground biomass (AGB) maps play a
crucial role in determining forest carbon stocks and how they are changing,
which is instrumental in understanding the global carbon cycle, and
implementing policy to mitigate climate change. The advent of the new space-
borne Light Detection and Ranging (LiDAR) sensor, NASA’s Global Ecosystem
Dynamics Investigation (GEDI) instrument, provides unparalleled possibilities
for the accurate and unbiased estimation of forest AGB at high resolution,
particularly in dense and tall forests, where Synthetic Aperture Radar (SAR)
and passive optical data exhibit saturation. However, GEDI is a sampling
instrument, collecting disperate footprints, and its data must be combined
with that from other continuous cover satellites to create high-resolution
maps, using local machine learning methods. In this study, we developed local
models to estimate forest AGB from GEDI L2A data, as the models used to create
GEDI L4 AGB data incorporated minimal field data from China. We then applied
an ensemble model of decision trees (Light Gradient Boosting Machine,
LightGBM), and random forest regression to generate wall-to-wall AGB maps at
25 m resolution, using extensive GEDI footprints as well as Sentinel-1 C-band
data, Advanced Land Observing Satellite (ALOS) -2 Phased Array L-band SAR
(PALSAR) -2 data and Sentinel-2 optical data. We found many GEDI footprints
had errors not detected by their quality flags; filtering using remote sensing
data to find realistic possible values improved results. Through a 5-fold
cross-validation, LightGBM demonstrated a slightly better performance than
Random Forest across two contrasting regions. In the northeastern region, the
coefficient of determination (R2) of the two models is both 0.70 and the root-
mean-square-error (RMSE) ranges between 28.44 and 28.77 Mg ha-1. Meanwhile, in
the southwestern region, the R2 values range from 0.60 to 0.62, and the RMSE
values range between 37.71 to 38.29 Mg ha-1. However, in both regions, the
computation speed of LightGBM is substantially faster than that of the random
forest model, requiring roughly one-third of the time to compute on the same
hardware. Through the validation against field data, the 25 m resolution AGB
maps generated using the local models developed in this study exhibited higher
accuracy compared to the GEDI L4B AGB data. We found in both regions an
increase in error as slope increased. In the northeastern region, the R2
values decreased from 0.93 in flat areas to 0.42 in areas with >30 degrees
slope, while the RMSE increased from 14 to 38 Mg ha-1. In the southwestern
region, the R2 values decreased from 0.75 in flat areas to 0.50 in areas with
>30 degrees slope, while the RMSE increased from 20 to 44 Mg ha-1. The trained
models were tested on nearby but different regions and exhibited good
performance. This study not only proposes an approach to map forest AGB using
extensive GEDI and continuous SAR, optical data in China but also enables the
estimation of AGB changes using the trained model, even in the absence of GEDI
but presence of SAR and optical data.

###### keywords:

Forest aboveground biomass , LightGBM , Random forest , GEDI , Contrasting
regions

[inst1]organization=School of GeoSciences, University of
Edinburgh,addressline=King’s Buildings, city=Edinburgh , postcode=EH9 3FF,
country=UK

[inst2]organization=Space Intelligence Ltd,addressline=93 George Street,
city=Edinburgh , postcode=EH2 3ES, country=UK

[inst3]organization=Department of Remote Sensing and Photogrammetry, Finnish
Geospatial Research Institute, addressline=Geodeetinrinne 2, city=Kirkkonummi,
postcode=02430, country=Finland

[inst5]organization=Institute for Logic, Language and Computation, University
of Amsterdam, addressline=Science Park 900, city=Amsterdam, postcode=1098 XH,
NL, country=Netherlands

[inst6]organization=Department of Forest Science, University of Helsinki,
addressline=Yliopistonkatu 4, city=Helsinki, postcode=00100, country=Finland

[inst7]organization=Jihua Lab, addressline=No 28 huandao south road,
city=Foshan, postcode=528200, country=China

Unique field data collected co-located with GEDI footprints to estimate forest
aboveground biomass from GEDI L2A with locally tuned models.

New filters using SAR data proved to be highly beneficial for extracting high-
quality GEDI footprints.

Compared gradient boosting ensemble machine learning algorithm and random
forest

The trained models performed well in other areas of China.

##  1 Introduction

Forests play an essential role in the health and sustainability of the
planet’s ecosystems (Watson et al., 2018). Accurate information on forest
aboveground biomass (AGB) is critical for effective forest management, carbon
accounting, and climate change mitigation efforts (Xu et al., 2021; Houghton
et al., 2009). Traditional forest inventory methods for estimating forest
biomass are labor-intensive and time-consuming, making them impractical for
large-scale monitoring and assessment. Approaches based on earth observations
have emerged as a promising alternative for estimating forest AGB (Mitchard,
2018; Tucker et al., 2023), offering the potential to rapidly and consistently
map AGB over large areas. Optical, SAR, and LiDAR are the three main types of
earth observation data.

Passive optical data measures light reflected from the Earth’s surface, and is
widely used due to the ability to capture detailed spectral information, which
can be used to assess canopy cover and potentially estimates of AGB. Optical
data provides the longest temporal record, with the series of Landsat that
commenced data collection in 1972. It has traditionally been an efficient tool
to map wall-to-wall forest AGB by combining with field data through regression
models (McNicol et al., 2018; Lu et al., 2004; Zhu and Liu, 2015; H. Nguyen et
al., 2019; Zheng et al., 2004). However, optical data can only see the top of
the forest canopy, meaning similar canopy with different heights or trunk
widths, are not be distinguished well. Additionally, optical data saturates at
low AGB values and its utility can be significantly hampered by atmospheric
conditions such as clouds and haze.

SAR is a type of active sensor that transmits its own energy and subsequently
records the quantity of that energy reflected back (Meyer, 2019). The long
wavelengths of the SAR sensors used for biomass mapping (C, L, and longer
wavelengths if available) enables it to see through clouds, which ensures
consistent data acquisition. The continuous, reliable, all-weather, day-and-
night imaging capabilities of SAR make it particularly useful for forest AGB
estimation (Mitchard et al., 2009). SAR can further penetrate through the
canopy and interact with leaves, branches and trunks, and thus SAR is
sensitive to forest AGB (Balzter et al., 2007). SAR data which has longer
wavelength tends to have a higher saturation point than those with shorter
wavelengths (Woodhouse et al., 2012; Sinha et al., 2015). However, even the
L-band SAR, which has the longest wavelength among currently available space-
borne SAR data, has a relatively low saturation point for AGB at around
100-150 Mg ha-1 (Mitchard et al., 2009; Englhart et al., 2011; Mitchard et
al., 2012).

LiDAR data employs laser light to estimate tree height and structure, thereby
facilitating the construction of a detailed representation of a forest. This
technology allows for the direct estimation of 3D structural parameters such
as height, stem density, and canopy cover, providing invaluable insights into
forest characteristics and dynamics. Therefore, LiDAR suffers little from
saturation effects (Wulder et al., 2012; Duncanson et al., 2022), although it
cannot penetrate through clouds or haze. With the launch of GEDI, the high-
resolution LiDAR instrument which is designed to measure vegetation structure
(Dubayah et al., 2020), it is possible to quantify forest AGB in distributed
25m diameter footprints, a footprint size chosen to minimize slope effects
compared to the previous spaceborne LiDAR the Ice Cloud and Land Elevation
Satellite (ICESat) with its 65m footprints, over very large areas.

However, useful instruments do not necessarily produce accurate forest AGB
maps. China encompasses diverse forest cover types and complex terrain, making
it a challenging and highly variable environment for AGB estimation (Dai et
al., 2011). China is among the countries with the highest degree of inaccuracy
when it comes to estimating AGB using GEDI data. The AGB estimated by GEDI for
China is more than double the estimates provided by the Food and Agriculture
Organization (FAO) (Dubayah et al., 2022) and other studies (Su et al., 2016;
Huang et al., 2019; Chang et al., 2021). One potential reason for this is the
lack of training data in China (Duncanson et al., 2022). Forest biomass
estimation using remote sensing relies on the correlation between remotely
sensed data and ground-based biomass measurements, which vary greatly
depending on the forest type and terrain conditions such as elevation and
slope. In this study, we collected new field plots co-located with GEDI
footprints to develop relationships to convert GEDI relative height (RH)
metrics into forest AGB.

However, unlike SAR and optical data, GEDI data are not continuous. The
discontinuous LiDAR data is often combined with optical and SAR data to
generate wall-to-wall estimates (Saatchi et al., 2011). Their combined use can
provide a more comprehensive and accurate estimation of forest carbon. Random
forest (Breiman, 2001) is one of the most popular algorithm to extrapolate
LiDAR footprints to wall-to-wall map using optical and SAR imagery (Dong et
al., 2024; Baccini et al., 2012; Huang et al., 2019; Narine et al., 2020).

In this study, we used random forest as well as a more advanced machine
learning algorithm, LightGBM (Ke et al., 2017), to generate wall-to-wall AGB
maps with a resolution of 25 m. LightGBM is a gradient boosting ensemble
method that is known for its efficiency and ability to handle large-scale data
(Yan et al., 2021). However, the universality of the two models across
different different regions is not yet fully understood. We tested the
performance of the two models in two contrasting regions, and validated the
results using the field data. In addition, the trained models were applied to
respective areas. The remote sensing data used in the models were from GEDI,
Sentinel-1, Sentinel-2, and ALOS-2, all of which are open-source. We aim to
propose a universally applicable method for estimating AGB in China, or if
necessary, demonstrate the need for region-specific models across different
regions.

##  2 Data and methods

###  2.1 Study area

This study was conducted in two contrasting areas of China: one in
northeastern China, and the other one in southwestern China (Fig. 1). Each of
the study areas covers an area of 500 x 500 km, but the types of forests and
terrain are completely different. The forest types differ in terms of their
structure, species composition, and environmental conditions, providing an
ideal opportunity to test the universality of remote sensing-based forest
biomass mapping approaches which is of great significance for achieving large-
scale forest AGB predictions in areas with complex forest compositions.

Fig. 1: Location of the study area and field sites. The field plots are shown
on the composited ALOS PALSAR-2 mosaic (HH in red, HV in green, and HV/HH in
blue).

The northeastern study area, primarily located in Jilin Province as well as
the southern part of Heilongjiang Province, is characterized by the temperate
forest. The area experiences a continental monsoon climate, with long, cold
winters and short, warm summers (Liu et al., 2018). The average temperatures
in January and July are -20 and 21 ∘C, respectively. The average annual
precipitation ranges from 600 to 1100 mm (Dong et al., 2005). The forests are
dominated by Korean pine (Pinus koraiensis) mixed with deciduous species (e.g.
Betula platyphylla) (Wang, 2006) (Fig. 2). A portion of the primary forest has
been destroyed by large-scale industrial logging, and replaced by secondary
forests and plantations (Yu et al., 2011), which mainly consist of Betula
platyphylla, Larix gmelinii, and Pinus koraiensis (Shi et al., 2015).

The southwestern study area is mainly situated in Yunnan Province. Forests in
this region predominantly blanket the mountainous terrain, as this area is
dominated by plateau and mountainous areas, with elevations ranging from 285
to 5350 m, and an average elevation of 2039 m (Digital Elevation Model). The
average temperature fluctuates between roughly 9–11 ∘C in January and climbs
to around 22 ∘C in July, while the average annual precipitation varies from
1100 to 1600 mm (Shi and Chen, 2018). The natural conditions described above
contribute to the richness and diversity of Yunnan’s forest resources,
including needleleaf forest, subtropical evergreen broad-leaved forests and
deciduous broadleaf forest (Zhang et al., 2019; Zhu et al., 2022; Frayer et
al., 2014).

Fig. 2: Field photos showing different forest types and terrains in the two
study areas.

###  2.2 Field data of AGB

The field measurements were carried out during October and December 2021. We
measured 26 plots in northeastern region with 24 plots were under GEDI
footprints, and 16 plots in southwestern region with 12 plots were under GEDI
footprints. For the collection of field data, we used the geolocation of the
GEDI Version 1 footprints.The geolocation uncertainty associated with these
footprints is reported to be 23.8 meters (Roy et al., 2021), which may have
introduced some degree of uncertainty in the co-localization of GEDI
footprints and field plots. In order to align with the GEDI data, all field
plots were circular with a diameter of 25 m. We endeavored to maintain a
substantial distance between field points (typically greater than 2 km) and
ensured they were at least 50 meters away from non-forest areas. This strategy
was adopted to avoid issues related to mixed pixels.

After identifying the GEDI footprints to be measured (Fig. 1), we navigated to
the center of each footprint using a Garmin eTrex 20x handheld GPS device. To
improve the accuracy of our measurements, we averaged the GPS readings over a
period of 5 minutes. A tape was then used to delineate a circular plot,
centered on the core of the footprint, with a diameter of 25 meters.
Subsequently, we measured the height and Diameter at Breast Height (DBH) of
trees within the plot that have a DBH $\geq$ 5 cm. This methodology ensures a
comprehensive and accurate assessment of the field plot. DBH was measured at
1.3 m up the trunk from the ground using a diameter tape. For the majority of
the field plots, tree height was estimated using a Vertex IV with a
Transponder T3. However, for two plots in the northeastern region, the tree
height was measured using a clinometer. The field AGB of individual tree was
then obtained using tree height and DBH, calculated via the following
allometric equation (Luo et al., 2018).

$W=0.1355*(D^{2}*H)^{0.817}$ (1)

where, $W$ is the AGB of each tree (kg), $D$ and $H$ are DBH (cm) and height
(m) of the tree. We sum the AGB of all trees within the plot that have a DBH
$\geq$ 5cm. This total AGB was then converted to a per hectare basis for
standardization, expressed in Megagrams per hectare (Mg ha-1).

###  2.3 LiDAR data

####  2.3.1 Canopy height metrics from GEDI

The GEDI instrument captures the waveform of the light that has bounced back
from the Earth’s surface to the GEDI sensor on the International Space
Station, with a 25 m spatial footprint every 60 m along orbits. The GEDI
instrument is equipped with three lasers. Out of these, two lasers operate at
full power while the third one is divided into two coverage beams. These four
beams are dithered across-track, creating a total of eight ground tracks of
data, comprising four full-power beams and four coverage beams (Dubayah et
al., 2020).

The dataset essentially presents the distribution of Earth’s surface features,
like vegetation and terrain, within the footprint of a laser shot.
Specifically, it includes details on the vertical structure of the features,
which is of great importance when assessing characteristics like canopy height
or the vertical distribution of foliage in forests (Dubayah et al., 2020). In
this study, L2A dataset acquired from 2019 to 2022 was used to produce
continuous forest AGB map.

####  2.3.2 Airborne laser scanning data

Airborne laser scanning (ALS) are capable of collecting data on vegetation
structure, which can potentially serve as a proxy for traditional field
measurements (Potapov et al., 2021). We employed ALS point clouds, collected
in the summer of 2022, with a point density of 18 points$/m^{2}$ in our
research. The ALS data consists two areas, the first area of 18 ha and the
other area of 22 ha, co-located with 11 GEDI footprints respectively.

###  2.4 SAR data

To obtain as much information as possible on forest structure, we utilized
both C-band and L-band SAR in this study.

####  2.4.1 Sentinel-1 C-band data

Sentinel-1 Interferometric Wide swath (IW) product acquired during 2020 and
2022 was used in the estimation. Sentinel-1 was accessed and processed in
Google Earth Engine (GEE). The Sentinel-1 Ground Range Detected (GRD) data had
been calibrated and terrain-corrected, and the values were converted to
decibels.

The backscatter coefficient of VV (vertical-vertical) and VH (vertical-
horizontal) polarizations was exploited, as different polarizations carries
information about the structure of the forest (Flores-Anderson et al., 2019).
Their cross ratio (VV/VH) was also included, as the cross ratio of Sentinel-1
is likely to mitigate errors associated with the acquisition system (Veloso et
al., 2017) and has been demonstrated to effectively differentiate between
vegetation densities (Vreugdenhil et al., 2018). To reduce the impact of
temporal variations and anomalies that might occur in individual images, we
calculated temporal means of the data. To further reduce the speckle noise, we
applied a focal mean speckle filter with a radius of 50 m to the data.

####  2.4.2 ALOS-2 PALSAR-2 L-band data

The annual mosaics of L-band SAR data from ALOS-2 PALSAR-2 was utilized in
this study. The 25 m resolution yearly mosaic dataset was created by the Japan
Aerospace Exploration Agency (JAXA), using data acquired in the target year.
The ortho-rectified and slope corrected data was accessed in GEE, and we also
applied a focal mean speckle filter with a radius of 50 m to the means of
annual mosaic of 2020 to 2022 to reduce the impact of speckle noise. The
digital numbers (DN) of HV and HH were then converted to gamma naught values
in decibel using the the following equation (Shimada et al., 2009).

$\gamma_{0}=10\log_{10}(\text{DN}^{2})-83.0\,\text{dB}$ (2)

where, $\gamma_{0}$ is the backscatter coefficient in dB. In addition to the
both polarizations, we also prepared the ratio (HV/HH) data, which has been
identified as effective strategies for enhancing the saturation point (Hayashi
et al., 2019; Sarker et al., 2012).

###  2.5 Optical data

Optical data were also used to extrapolate LiDAR footprints to wall-to-wall
image. The MultiSpectral Instrument (MSI) on board the Sentinel-2 satellite
measures reflected radiance across 13 spectral bands. These bands span a range
from the Visible and Near Infrared (VNIR) to the Short Wave Infrared (SWIR)
spectra. This broad coverage enables detailed observations of various surface
conditions and phenomena. The ortho-rectified and atmospherically corrected
L2A surface reflectance (SR) dataset, excluding band 10, is publicly
accessible in GEE. In addition to utilizing the mean values of each band of
L2A imagery acquired between 2020 and 2022, we also computed the mean and
maximum values of the Normalized Difference Vegetation Index (NDVI), which
provide valuable information regarding the average and peak level of
vegetation growth.

###  2.6 Topographic data

The NASA Digital Elevation Model (NASADEM) is a global digital elevation model
with a resolution of 30 m. NASADEM combines data from different satellite
missions, including the Shuttle Radar Topography Mission (SRTM), the Advanced
Spaceborne Thermal Emission and Reflection Radiometer (ASTER), and ICESat, to
generate a seamless global elevation model (JPL, 2020). In this study, we
utilized the terrain slope derived from NASADEM as an additional variable for
analysis. The terrain slope was computed using a 3 × 3 moving window approach.
The integration of terrain slope information provided valuable insights into
the topographic characteristics of the study areas, enhancing the impact of
topographic variations on the AGB predictions.

###  2.7 Forest AGB estimation

####  2.7.1 GEDI data filtering

Before proceeding with model fitting, we implemented multiple filters to
identify and remove poor-quality footprints. L2A footprints collected from
2019 to 2022 were filtered using the following criteria (Dubayah et al.,
2021).

(a) Only shots that have a quality_flag of 1 were used. This flag signifies
that the shot meets a summation of energy level, sensitivity, amplitude,
quality of real-time surface tracking, and the difference from DEM.

(b) Shots acquired during degrade period were excluded.

(c) Only shots with a sensitivity $>$ 0.98 were used. A higher sensitivity
indicates a greater likelihood of penetrating to the ground beneath a specific
level of canopy cover.

(d) Shots were included if they were acquired at night, when the negative
impact of solar illumination on the quality of GEDI waveform is eliminated,
thus yielding better tree height prediction results (Liu et al., 2021).

After we filtered the GEDI footprints, we then compared GEDI L1B data with ALS
data to verify the quality of GEDI data, and to eliminate poor-quality
footprints in GEDI. In order to compare with GEDI footprints, we used ALS data
to correct geolocation errors in GEDI observations and to simulate GEDI data
using three ground-finding algorithms: the Gaussian fitting, the lowest
inflection point and the lowest maximum (Hancock et al., 2019).The simulated
GEDI observations from ALS data was then compared with actual GEDI footprints
(fig. 3).

Fig. 3: RH98 of actual GEDI against GEDI simulator. (a)-(c) is the scatter
plot for the first area, all the footprints are comprised within full-power
beam. (d)-(f) is the scatter plot for the second area, all the footprints are
comprised within coverage power beam. The dotted line represents 1:1 line. The
solid line represent the fitted line, and is extended to cover the entire
range of the x-axis. This was done to provide a comprehensive visual
representation of the relationships between GEDI RH metrics and ALS simulated
RH metrics.

As expected, footprints of GEDI power beams have a higher relationship with
ALS simulated footprints, with a higher R2 and a smaller RMSE. Since this also
had been proved in other study areas (Duncanson et al., 2020; Li et al., 2023;
Liu et al., 2021), we only included power beams in this study.

By applying the approach proposed by Hancock et al. (2019) to compare GEDI and
ALS data, we found that the horizontal geolocation errors of GEDI data in
China are substantial, with a 17.9 m error for the first study area and a 9.8
m error for the second area. We noticed a significant difference in the
geolocation error between the two areas, indicating that the geolocation error
of some footprints is larger than others. This means that the RH metrics at
the actual geolocation of some GEDI footprints may differ significantly from
the RH metrics at the geolocation it shows, which could introduce large error
into AGB predictions.

Because of this finding, we proceeded with an additional filtering process on
the GEDI footprints using ALOS-2 PALSAR-2 HV, HH polarizations, and Sentinel-1
VH, VV polarizations, which have been proven to have a close relationship with
forest height and AGB (Mitchard et al., 2012; Duncanson et al., 2020; Cartus
et al., 2012; Ge et al., 2022; Santi et al., 2017). We fitted the relationship
between GEDI RH98 and SAR imagery in the Northeast region. The reliability of
the fitting relationships were assessed using R2 and RMSE, calculated using
the following equations:

$R^{2}=1-\frac{\sum_{i=1}^{n}(SAR_{i}-f(RH_{i}))^{2}}{\sum_{i=1}^{n}(SAR_{i}-\overline{SAR})^{2}}$
(3) $RMSE=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(SAR_{i}-f(RH_{i}))^{2}}$ (4)

where $SAR_{i}$ is the value of SAR imagery, $f(RH_{i})$ is the fitted value,
$overline{SAR}$ is the average value of $SAR_{i}$. The R2 and RMSE here can be
influenced by the fitted function.

In northeastern region, cross-polarised SAR showed strong relationships with
GEDI RH98 (Fig. 4a, c), while the relationships between co-polarised SAR and
GEDI RH98 were not as expected (Fig. 4b, d). In contrast, in southwestern
region, neither polarization showed a significant relationship with GEDI RH
metrics, which might be substantially influenced by the steep slopes. There
remains a substantial amount of GEDI data that does not meet the required
quality standards. Generally, the backscatter of SAR from shorter trees is
expected to be lower than that from taller trees. We therefore filtered the
GEDI footprints in both regions using the relationships developed in
northeastern regions, and only retained the GEDI data that fell within a range
of plus and minus 3 dB from the fitted curve.

Fig. 4: ALOS PALSAR-2 backscatter in (a) HV, (b) HH polarizations, and
Sentinel-1 backscatter in (c) VH, (d) VV polarizations were plotted against
GEDI RH98 for 1,405,451 GEDI footprints in the northeastern region. The black
solid line represents the fitted line. The red solid line above the black
solid line represents the fitted line plus 3 dB, and the red solid line below
the black solid line represents the fitted line minus 3 dB. The color of each
point in the scatter plot represents the estimated density of points at that
location, calculated using a Gaussian kernel density estimate.

####  2.7.2 AGB of GEDI footprints

The GEDI L4B provides highly accurate biomass maps; however, there remains
considerable uncertainty in biomass estimation for China (Dubayah et al.,
2022). Furthermore, the reliability of the models used in GEDI L4A to convert
RH metrics to AGB for Asian is relatively low, as evidenced by small R2 values
and large %RMSE (Duncanson et al., 2022). The GEDI L4 biomass products heavily
rely on field data, and little data from China is included in the calibration
process. Consequently, there is a pressing need to establish local models for
calculating AGB from RH metrics, using field data collected within China.

In order to obtain AGB values for the GEDI footprints, we analysed the
correlation between field AGB and RH values ranging from RH30 to RH98. In the
northeastern region, the closest relationship was observed between field AGB
and RH98, while in the mountainous southwestern region, field AGB exhibited a
stronger correlation with RH80 (Fig. 5). We converted RH98 and RH80 to AGB
using the fitted relationships for both regions, respectively.

Fig. 5: Field measured AGB are plotted against (a) RH98 in the northeastern
region, and (b) RH80 in the southwestern region.

####  2.7.3 Accuracy assessment

We evaluated the model performance using a 5-fold cross-validation approach.
The dataset was divided into five distinct subsets. Four subsets were employed
for training, while the remaining one was used to test the model’s accuracy.
This procedure was repeated five times, ensuring each subset served as the
test set once. Subsequent to these iterations, the results from all five test
sets were aggregated to provide a comprehensive assessment of the model’s
accuracy.

The model performance and the accuracy of the results were evaluated using R2
and RMSE, computed using the following equations:

$R^{2}=1-\frac{\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}}$
(5) $RMSE=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}$ (6)
$\mathrm{Bias}=\frac{1}{K}\sum_{i=1}^{K}(\hat{y}_{i}-y_{i})$ (7)

where $y_{i}$ is the field measured AGB or the AGB derived from GEDI RH
metrics, $\hat{y}_{i}$ is the model predicted value, $\bar{y}$ is the average
value of $y_{i}$, and $K$ is the cumulative number of observations.

####  2.7.4 Wall-to-wall AGB map

In this study, we extrapolated the extensive GEDI footprints into continuous
AGB maps using LightGBM and random forest in Python 3.8. A comprehensive set
of 25 remote sensing layers were employed in the generation of wall-to-wall
maps. These layers comprised a wide array of data, including bands and ratios
from SAR, optical bands, the maximum and average NDVI values, data from the
DEM and slope, along with gridded latitude and longitude information.

Using the 5-fold cross-validation, each trained model from the five iterations
was employed to predict the AGB maps across the study area, resulting in five
separate AGB maps for each of the two regions. To derive a robust and
consolidated AGB map, we computed the mean of these five predictions for each
pixel. Additionally, to assess the variability and reliability of our
predictions, we calculated the standard deviation across the five results for
each pixel, which was then utilized to generate an uncertainty map (Fig. 6).

Fig. 6: Framework for estimating AGB and its uncertainty.

To show the spatial distribution of forest AGB in both regions, the AGB maps
were masked by a forest cover map in a 2021. The 2021 forest cover was derived
from the Hansen et al forest map (Hansen et al., 2013), based on the 2000
forest cover map as a baseline. To generate the 2021 map, the degraded areas
between 2000 and 2021 were excluded, and the forest areas that experienced an
increase from 2000 to 2012 were added to the map (there is no information of
forest gain after 2012).

##  3 Results

####  3.0.1 Model performance

Filtering GEDI data using parameters such as quality_flag, degrade,
sensitivity, solar elevation (acquisition time), and beam intensity had been
demonstrated to be an effective method for enhancing data quality. 1,405,451
and 2,042,330 GEDI footprints, filtered by the parameters above, were used to
generate AGB maps for the northeastern and southwestern regions, respectively.
The GEDI data was divided into training and test datasets at a ratio of 6:4.
Scatter plots of the estimates compared to the test data are shown in Fig. 7,
the variation of test data density is also illustrated. This initial analysis,
which utilized GEDI footprints irrespective of GEDI RH metrics against SAR
backscatter, indicated that both models exhibited strong performance in the
northeastern region. However, neither model performed well in the southwestern
region. There are many outliers in Fig. 7, which may be caused by geolocation
uncertainties in the GEDI data (Fig. 3 and section 2.7.1), or steep slopes.

Fig. 7: The predicted AGB values against AGB derived from GEDI RH metrics
before filtering using earth observation data. (a) LightGBM in northeastern
region, (b)random forest in northeastern region, (c) LightGBM in southwestern
region, (d) random forest in southwestern region. The color of each point in
the scatter plot represents the estimated density of points at that location.

We then utilized the GEDI footprints filtered using SAR data as described in
section 2.7.1, with 689,013 and 673,744 footprints remaining for the
northeastern and southwestern regions. We employed two different models to
predict forest AGB. Predicted AGB from both models against test dataset of AGB
derived from GEDI RH metrics are plotted in Fig. 8. The performance of both
models in the northeastern region did not show significant improvement, except
for a slight reduction in the RMSE and a notable reduction in the number of
outliers. However, both models showed significant improvement in the
southwestern region. The R2 and RMSE of the LightGBM model improved from 0.35
and 65.22 Mg ha-1 to 0.62 and 37.71 Mg ha-1, respectively. Similarly, the
Random Forest model saw its R2 and RMSE improve from 0.31 and 67.46 Mg ha-1 to
0.60 and 38.29 Mg ha-1, respectively.

Fig. 8: The predicted AGB values against AGB derived from GEDI RH metrics
after filtering using earth observation data. (a) LightGBM in northeastern
region, (b)random forest in northeastern region, (c) LightGBM in southwestern
region, (d) random forest in southwestern region. The color of each point in
the scatter plot represents the estimated density of points at that location.

Both models demonstrated similar accuracy in the two regions. However,
LightGBM outperforms Random Forest in terms of computational speed. Table 1
presents the computational time required for training and prediction using the
LightGBM and Random Forest models. The reduced computational time of LightGBM
can be attributed to its optimized algorithm and efficient implementation. The
notable difference in computational speed between LightGBM and Random Forest
makes LightGBM a favorable choice for applications where efficient processing
is crucial, especially when dealing with large-scale datasets.

Table 1: Computation time for LightGBM and random forest. NE represents the northeastern region, and SW represents the southwestern region. Study area | Model | Number of | Computation
---|---|---|---
GEDI footprints | time (min)
NE | LightGBM | 1,405,451 | 89
NE | random forest | 1,405,451 | 233
NE | LightGBM | 689,013 | 51
NE | random forest | 689,013 | 145
SW | LightGBM | 2,042,330 | 97
SW | random forest | 2,042,330 | 329
SW | LightGBM | 673,744 | 69
SW | random forest | 673,744 | 187

####  3.0.2 Accuracy of the AGB maps

The accuracy of the AGB maps produced by both models was evaluated using 42
field plots. The estimated AGB from LightGBM and random forest correspond well
with the field AGB in both regions (Fig. 9). Overall, the accuracy of LightGBM
is slightly higher than that of Random Forest. The AGB map generated by
LightGBM exhibits an R2 of 0.59, an RMSE of 55.36 Mg ha-1, and a bias of 17.61
Mg ha-1. In contrast, the AGB map produced by the random forest model yielded
an R2 of 0.56, an RMSE of 56.72 Mg ha-1, and a bias of 18.33 Mg ha-1. Compared
to the field data, both models have overestimated biomass in low biomass areas
and underestimated biomass in high biomass areas, as is common in biomass
mapping using machine learning (Shendryk, 2022; Huang et al., 2019). This
trend is also evident when evaluating model performance using AGB derived from
the GEDI footprints (Fig. 8).

For comparison, we also conducted an analysis using collected field data
against the 1 km resolution GEDI L4B Version 2 AGB map. Due to gaps in the
GEDI L4B data product, three field plots did not overlap with the GEDI L4B
data. Consequently, we conducted the accuracy assessment using the remaining
39 field plots. The results indicated an R2 of 0.36, an RMSE (Root Mean Square
Error) of 68.77 Mg ha-1, and a bias of -30.34 Mg ha-1.

Fig. 9: Scatter plot between predicted AGB and field measured AGB. (a)
LightGBM, (b)random forest, (c) GEDI L4B (Version 2). The solid line represent
the fitted line, and the dotted line represents 1:1 line.

####  3.0.3 Spatial distribution of AGB in both regions

The wall-to-wall forest AGB maps by LightGBM and random forest in both regions
are shown in Fig. 10. In the northeastern region, high AGB values are
predominantly concentrated in the eastern part of the study area. Conversely,
in the southwestern region, the distribution of AGB exhibits a more dispersed
pattern. In both the northeast and southwest regions, the spatial distribution
patterns of AGB predicted by the LightGBM and Random Forest models show strong
agreement, perhaps unsurprisingly given they are based on the same Earth
Observation data inputs. This is further evidenced in the histograms, where
both models generate similar overall AGB distributions 11. In the northeastern
region, the mean difference in AGB estimates between LightGBM and Random
Forest was -1.59 Mg ha-1. Similarly, in the southwestern region, the mean AGB
estimate difference was -1.44 Mg ha-1. Although the overall spatial
distribution patterns show strong agreement, difference maps in Fig. 10 reveal
there are still discernible differences between the LightGBM and Random Forest
AGB predictions. To better illustrate these spatial differences at a finer
scale, zoomed-in 2500 m x 2500 m maps are presented in Fig. 12. Focusing on
this local scale enables better visualization and analysis of the details and
spatial heterogeneity in the AGB distribution. The AGB maps generated by the
Random Forest model exhibit a deeper characterization of extreme values (Fig.
12b, e), both at the high and low ends of the AGB range. In contrast, the AGB
maps generated by the LightGBM model display a relatively smoother transition
between different AGB values (Fig. 12a, d), with more gradual change.

Fig. 10: Spatial distribution of AGB and the differences between the AGB maps
from the two models at 25 m resolution. (a) LightGBM in the northeastern
region, (b) Random forest in the northeastern region, (c) Difference Map
between LightGBM and RF Estimated AGB in the northeastern region (d) LightGBM
in the southwestern region, (e) Random forest in the southwestern region, (f)
Difference Map between LightGBM and RF Estimated AGB in the southwestern
region. (a) Fig. 11: Frequency distribution histograms of AGB, in 10 Mg ha-1,
for (a)The northeastern region, (b)The southwestern region. Fig. 12: Spatial
distribution of AGB for two 2500 m x 2500 m areas. (a) (d) LightGBM, (b) (e)
Random forest, (c) (f) difference maps.

We also generated uncertainty maps based on the standard deviation of results
from 5-fold cross-validation (Fig. 13). Both LightGBM and Random Forest showed
relatively low uncertainty in the two study regions. For LightGBM, the
uncertainty in the northeastern region was 3.28 Mg ha-1, while in the
southwestern region it was 4.22 Mg ha-1. Meanwhile, Random Forest uncertainty
was 2.96 Mg ha-1 in the northeast and 3.44 Mg ha-1 in the southwest. However,
LightGBM exhibited slightly higher uncertainty, indicating it may be more
sensitive to differences across folds in the training data. Further analysis
is warranted to determine the factors contributing to the uncertainty
patterns.

Fig. 13: Spatial distribution of uncertainty and the differences between the
uncertainty maps from the two models at 25 m resolution. (a) LightGBM in the
northeastern region, (b) Random forest in the northeastern region, (c)
Difference Map between LightGBM and RF uncertainty in the northeastern region,
(d) LightGBM in the southwestern region, (e) Random forest in the southwestern
region, (f) Difference Map between LightGBM and RF uncertainty in the
southwestern region.

##  4 Discussion

###  4.1 Influence of steep slopes

Both regions are characterized by a significant amount of mountainous terrain,
resulting in many GEDI footprints falling in areas with substantial slopes,
especially in the southwestern region (Fig. 14). This could potentially impact
the results, as steep slopes could lead to overestimation or underestimation
of canopy height by GEDI data, depending on the steepness of the slope and
spatial distribution of trees within footprints (Chen, 2010; Schneider et al.,
2020). Slopes also impact the backscatter of the SAR data (Shimada et al.,
2011). Ideally the SAR data should be corrected for increased or decreased
backscatter caused by terrain slopes. However, the relationship between slope
and backscatter is complex and depends on forest structure and properties, so
slope corrections are not always feasible (Ramachandran et al., 2021).
Additionally, some steep slopes can cause SAR layover and shadow effects,
resulting in no useful SAR data for those areas (Samuele et al., 2021). These
terrain factors introduce uncertainties in both the GEDI lidar and SAR data
which should be considered, especially in regions with substantial topography.

As the slope of the GEDI footprints increases, there is a noticeable decrease
in the R2 values, and the RMSE gradually increases (Fig. 15) in both regions.
Overall, the two models demonstrated similar accuracy in flat areas. However,
in areas with steep slopes, LightGBM outperformed the random forest model,
with a larger R2 and smaller RMSE. We chose not to remove GEDI footprints from
areas with steep slopes, because their removal could lead to unreliable AGB
predictions in these steep areas, as much remaining forest is genuinely
located on steep slopes in these regions. In the northeastern region, the
average slope of GEDI is 10.23 degrees, and for the GEDI footprints within the
Hansen forest mask, it is slightly higher at 11.97 degrees. The average slope
of GEDI in the northeastern region measures 20.38 degrees, whereas the average
slope of GEDI footprints within the Hansen forest mask stands at 21.18
degrees. More field data collected in areas with steep slopes are needed to
evaluate the impact of removing GEDI footprints in steep areas on the accuracy
of AGB estimation in these areas.

Fig. 14: GEDI RH metrics are plotted against slope. The slope bin of 5
indicates slope from 0 to 5, and similar with all other bins. The colors of
the boxes represent the quantity of GEDI footprints in different slope bins.
(a) RH98 in northeastern region, (b) RH80 in southwestern region, (c) RH98 in
northeastern region after filtering using SAR data, (d) RH80 in southwestern
region after filtering using SAR data. Fig. 15: Model performance in both
regions are plotted against slope.

###  4.2 Scalability of the approach

The models trained with GEDI and remote sensing imagery in the northeastern
and southwestern regions were applied directly to two 100 x 100 km areas (Fig.
1). These test areas are located 90 km from the northeastern region and 450 km
from the southwestern region, respectively. In the northern test area, the
results of both models closely aligned with the AGB derived from GEDI data
(Fig. 16). The LightGBM and random forest models achieved R2 values of 0.76
and 0.77, with respective RMSE values of 31.01 and 29.44 Mg ha-1. The fitted
line and the 1:1 line between the model predictions and the AGB derived from
GEDI were very close, indicating no significant saturation point within the
AGB range below 300 Mg ha-1. In the middle southern test area, the R2 values
were slightly lower, with LightGBM and random forest models achieving R2
values of 0.55 and 0.56, respectively. However, the RMSE values remained
relatively low, at 48.79 and 46.80 Mg ha-1, respectively, indicating high
accuracy in this area as well.

Fig. 16: Performance of trained models in other area. (a) Trained LightGBM in
the northern test region, (b) Trained random forest in the northern test
region, (c) Trained LightGBM in the middle southern region, (d) Trained random
forest in the middle southern region. The color of each point in the scatter
plot represents the estimated density of points at that location.

###  4.3 Linking field AGB and GEDI RH metrics

In this study, we measured canopy height and DBH of trees in 42 field plots,
with 36 of them co-located with GEDI footprints. There are various allometric
equations to estimate forest AGB using height and DBH, depending on tree
species and regions. The AGB estimates derived from different allometric
equations can vary significantly. In this study, we used a general allometric
equation developed for the trees in China (Luo et al., 2018). The selection of
field plots significantly influence the conversion of GEDI RH metrics into
AGB. It is crucial that the selection of field plots is representative. To
this end, we have endeavored to ensure that our field points are evenly
distributed spatially, thereby enhancing their representativeness.

###  4.4 Geolocation uncertainty of GEDI footprint

In this study, we used GEDI L2A version 2 dataset, which improved geolocation
accuracy compared to version 1. However, upon comparison with ALS data, we
found that the geolocation uncertainty in the GEDI data remains substantial.
This may introduce errors into the AGB estimation of this study from two
perspectives. Firstly, our field plots and GEDI footprints are at the same
location and of the same size. However, there may be a discrepancy of about
ten meters with the actual location of the GEDI footprints. This discrepancy
could potentially result in the relationship we’ve established between field
AGB and GEDI RH not being as accurate as intended. The effect had been
somewhat mitigated by our strategy for selecting field points, which ensures
that our field points were at least 50 meters away from non-forest areas. In
addition, the geolocation uncertainty may introduce errors when matching GEDI
footprints with continuous remote sensing images in the models. The
geolocation uncertainty of GEDI footprints can be greatly improved through the
use of ALS data. However, ALS data is scarce and expensive to acquire. Our
research provides a reference for regions and countries where ALS data is
lacking, especially for large regions. Our method of using known relationships
between AGB and key remote sensing metrics to filter out GEDI footprints that
are far from the expected values is shown here to be a reliable way to remove
low reliability GEDI footprints, whatever the cause of error.

##  5 Conclusions

In this study, we presented an experiment designed to determine the
universality of remote sensing-based forest biomass mapping approaches at 25 m
resolution, using new field plots co-located with GEDI footprints in two
contrasting regions in northeastern and southwestern China. Both LightGBM and
random forest performed well in the two regions compared to the GEDI
validation data and field AGB. We also applied the models trained in these two
contrasting regions to a northern area and a central-southern area,
respectively, and achieved good results.Our preliminary test of the
scalability of the models across different areas seems to demonstrate their
robustness and ability to capture the relationships between the GEDI data and
continuous remote sensing data in different areas with similar forest types.
The promising outcomes can be attributed to several key factors. Firstly, the
availability of a large volume of GEDI footprints and multiple sources of
remote sensing imagery, including SAR and optical data, contributes to the
accuracy and comprehensiveness of the prediction. In addition, the stringent
removal of poor-quality data from the GEDI dataset enhances the reliability of
the results. Our study demonstrates the significant potential of using GEDI
data to produce high resolution AGB maps of China in the GEDI era (2019 and
onwards).

LightGBM and random forest models exhibit similar performance in terms of
accuracy. However, LightGBM outperforms random forest in terms of
computational speed. With the growing number of space-borne sensors acquiring
imagery and LiDAR footprints at high resolutions and frequencies, LightGBM’s
computational speed advantage empowers researchers to efficiently process
large-scale and high-frequency AGB estimation (i.e. national scale annual
estimates).

##  6 Acknowledgements

Edward Mitchard’s time was part funded by NERC grant SEOSAW (Grant number
NE/P008755/1). The field work was supported by Davis Expedition Fund,
Elizabeth Sinclair Irvine Bequest and Centenary Agroforestry 89 Fund, Moray
Endowment Fund and Meiklejohn fund. The authors would like to thank ESA and
JAXA for providing the Sentinel-1, Sentinel-2 and ALOS-2 PALSAR-2 mosaic data.

## References

  * Baccini et al. (2012) Baccini, A., Goetz, S., Walker, W., Laporte, N., Sun, M., Sulla-Menashe, D., Hackler, J., Beck, P., Dubayah, R., Friedl, M., et al., 2012\.  Estimated carbon dioxide emissions from tropical deforestation improved by carbon-density maps.  Nature climate change 2, 182–185.
  * Balzter et al. (2007) Balzter, H., Rowland, C., Saich, P., 2007.  Forest canopy height and carbon estimation at monks wood national nature reserve, uk, using dual-wavelength sar interferometry.  Remote Sensing of Environment 108, 224–239.
  * Breiman (2001) Breiman, L., 2001.  Random forests.  Machine learning 45, 5–32.
  * Cartus et al. (2012) Cartus, O., Santoro, M., Kellndorfer, J., 2012.  Mapping forest aboveground biomass in the northeastern united states with alos palsar dual-polarization l-band.  Remote Sensing of Environment 124, 466–478.
  * Chang et al. (2021) Chang, Z., Hobeichi, S., Wang, Y.P., Tang, X., Abramowitz, G., Chen, Y., Cao, N., Yu, M., Huang, H., Zhou, G., et al., 2021.  New forest aboveground biomass maps of china integrating multiple datasets.  Remote Sensing 13, 2892\.
  * Chen (2010) Chen, Q., 2010.  Retrieving vegetation height of forests and woodlands over mountainous areas in the pacific coast region using satellite laser altimetry.  Remote Sensing of Environment 114, 1610–1627.
  * Dai et al. (2011) Dai, L., Wang, Y., Su, D., Zhou, L., Yu, D., Lewis, B.J., Qi, L., 2011\.  Major forest types and the evolution of sustainable forestry in china.  Environmental Management 48, 1066–1078.
  * Dong et al. (2024) Dong, W., Mitchard, E.T., Santoro, M., Chen, M., Wheeler, C.E., 2024.  A new circa 2007 biomass map for china differs significantly from existing maps.  Scientific Data 11, 287\.
  * Dong et al. (2005) Dong, X., Li-min, D., Guo-fan, S., Lei, T., Hui, W., 2005.  Forest fire risk zone mapping from satellite images and gis for baihe forestry bureau, jilin, china.  Journal of forestry research 16, 169–174.
  * Dubayah et al. (2022) Dubayah, R., Armston, J., Healey, S.P., Bruening, J.M., Patterson, P.L., Kellner, J.R., Duncanson, L., Saarela, S., Ståhl, G., Yang, Z., et al., 2022\.  Gedi launches a new era of biomass inference from space.  Environmental Research Letters 17, 095001.
  * Dubayah et al. (2020) Dubayah, R., Blair, J.B., Goetz, S., Fatoyinbo, L., Hansen, M., Healey, S., Hofton, M., Hurtt, G., Kellner, J., Luthcke, S., et al., 2020\.  The global ecosystem dynamics investigation: High-resolution laser ranging of the earth’s forests and topography.  Science of remote sensing 1, 100002\.
  * Dubayah et al. (2021) Dubayah, R., Hofton, M., Blair, J., Armston, J., Tang, H., Luthcke, S., 2021\.  Gedi l2a elevation and height metrics data global footprint level v002. nasa eosdis land processes daac.
  * Duncanson et al. (2022) Duncanson, L., Kellner, J.R., Armston, J., Dubayah, R., Minor, D.M., Hancock, S., Healey, S.P., Patterson, P.L., Saarela, S., Marselis, S., et al., 2022\.  Aboveground biomass density models for nasa’s global ecosystem dynamics investigation (gedi) lidar mission.  Remote Sensing of Environment 270, 112845.
  * Duncanson et al. (2020) Duncanson, L., Neuenschwander, A., Hancock, S., Thomas, N., Fatoyinbo, T., Simard, M., Silva, C.A., Armston, J., Luthcke, S.B., Hofton, M., et al., 2020\.  Biomass estimation from simulated gedi, icesat-2 and nisar across environmental gradients in sonoma county, california.  Remote Sensing of Environment 242, 111779.
  * Englhart et al. (2011) Englhart, S., Keuck, V., Siegert, F., 2011.  Aboveground biomass retrieval in tropical forests—the potential of combined x-and l-band sar data use.  Remote sensing of environment 115, 1260–1271.
  * Flores-Anderson et al. (2019) Flores-Anderson, A.I., Herndon, K.E., Thapa, R.B., Cherrington, E., 2019\.  Sar handbook: Comprehensive methodologies for forest monitoring and biomass estimation. nasa.
  * Frayer et al. (2014) Frayer, J., Müller, D., Sun, Z., Munroe, D.K., Xu, J., 2014.  Processes underlying 50 years of local forest-cover change in yunnan, china.  Forests 5, 3257–3273.
  * Ge et al. (2022) Ge, S., Su, W., Gu, H., Rauste, Y., Praks, J., Antropov, O., 2022.  Improved lstm model for boreal forest height mapping using sentinel-1 time series.  Remote Sensing 14, 5560\.
  * H. Nguyen et al. (2019) H. Nguyen, T., Jones, S., Soto-Berelov, M., Haywood, A., Hislop, S., 2019.  Landsat time-series for estimating forest aboveground biomass and its dynamics across space and time: A review.  Remote Sensing 12, 98\.
  * Hancock et al. (2019) Hancock, S., Armston, J., Hofton, M., Sun, X., Tang, H., Duncanson, L.I., Kellner, J.R., Dubayah, R., 2019\.  The gedi simulator: A large-footprint waveform lidar simulator for calibration and validation of spaceborne missions.  Earth and Space Science 6, 294–310.
  * Hansen et al. (2013) Hansen, M.C., Potapov, P.V., Moore, R., Hancher, M., Turubanova, S.A., Tyukavina, A., Thau, D., Stehman, S.V., Goetz, S.J., Loveland, T.R., et al., 2013\.  High-resolution global maps of 21st-century forest cover change.  science 342, 850–853.
  * Hayashi et al. (2019) Hayashi, M., Motohka, T., Sawada, Y., 2019.  Aboveground biomass mapping using alos-2/palsar-2 time-series images for borneo’s forest.  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 12, 5167–5177.
  * Houghton et al. (2009) Houghton, R., Hall, F., Goetz, S.J., 2009.  Importance of biomass in the global carbon cycle.  Journal of Geophysical Research: Biogeosciences 114\.
  * Huang et al. (2019) Huang, H., Liu, C., Wang, X., Zhou, X., Gong, P., 2019\.  Integration of multi-resource remotely sensed data and allometric models for forest aboveground biomass estimation in china.  Remote Sensing of Environment 221, 225–234.
  * JPL (2020) JPL, N., 2020.  Nasadem merged dem global 1 arc second v001 [data set]. nasa eosdis land processes daac.  2020 .
  * Ke et al. (2017) Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., Liu, T.Y., 2017.  Lightgbm: A highly efficient gradient boosting decision tree.  Advances in neural information processing systems 30\.
  * Li et al. (2023) Li, X., Wessels, K., Armston, J., Hancock, S., Mathieu, R., Main, R., Naidoo, L., Erasmus, B., Scholes, R., 2023.  First validation of gedi canopy heights in african savannas.  Remote Sensing of Environment 285, 113402.
  * Liu et al. (2021) Liu, A., Cheng, X., Chen, Z., 2021.  Performance evaluation of gedi and icesat-2 laser altimeter data for terrain and canopy height retrievals.  Remote Sensing of Environment 264, 112571.
  * Liu et al. (2018) Liu, Q., Shan, Y., Shu, L., Sun, P., Du, S., 2018\.  Spatial and temporal distribution of forest fire frequency and forest area burnt in jilin province, northeast china.  Journal of Forestry Research 29, 1233–1239.
  * Lu et al. (2004) Lu, D., Mausel, P., Brondızio, E., Moran, E., 2004\.  Relationships between forest stand parameters and landsat tm spectral responses in the brazilian amazon basin.  Forest ecology and management 198, 149–167.
  * Luo et al. (2018) Luo, Y., Wang, X., Ouyang, Z., 2018.  A china’s normalized tree biomass equation dataset.
  * McNicol et al. (2018) McNicol, I.M., Ryan, C.M., Mitchard, E.T., 2018.  Carbon losses from deforestation and widespread degradation offset by extensive growth in african woodlands.  Nature communications 9, 3045\.
  * Meyer (2019) Meyer, F., 2019.  Spaceborne synthetic aperture radar: Principles, data access, and basic processing techniques.  Synthetic Aperture Radar (SAR) Handbook: Comprehensive Methodologies for Forest Monitoring and Biomass Estimation , 21–64.
  * Mitchard (2018) Mitchard, E.T., 2018.  The tropical forest carbon cycle and climate change.  Nature 559, 527–534.
  * Mitchard et al. (2012) Mitchard, E.T., Saatchi, S.S., White, L.J., Abernethy, K.A., Jeffery, K.J., Lewis, S.L., Collins, M., Lefsky, M.A., Leal, M.E., Woodhouse, I.H., et al., 2012\.  Mapping tropical forest biomass with radar and spaceborne lidar in lopé national park, gabon: overcoming problems of high biomass and persistent cloud.  Biogeosciences 9, 179–191.
  * Mitchard et al. (2009) Mitchard, E.T., Saatchi, S.S., Woodhouse, I.H., Nangendo, G., Ribeiro, N., Williams, M., Ryan, C.M., Lewis, S.L., Feldpausch, T., Meir, P., 2009\.  Using satellite radar backscatter to predict above-ground woody biomass: A consistent relationship across four different african landscapes.  Geophysical Research Letters 36.
  * Narine et al. (2020) Narine, L.L., Popescu, S.C., Malambo, L., 2020.  Using icesat-2 to estimate and map forest aboveground biomass: A first example.  Remote Sensing 12, 1824\.
  * Potapov et al. (2021) Potapov, P., Li, X., Hernandez-Serna, A., Tyukavina, A., Hansen, M.C., Kommareddy, A., Pickens, A., Turubanova, S., Tang, H., Silva, C.E., et al., 2021\.  Mapping global forest canopy height through integration of gedi and landsat data.  Remote Sensing of Environment 253, 112165.
  * Ramachandran et al. (2021) Ramachandran, N., Saatchi, S., Tebaldini, S., d’Alessandro, M.M., Dikshit, O., 2021.  Evaluation of p-band sar tomography for mapping tropical forest vertical backscatter and tree height.  Remote Sensing 13, 1485\.
  * Roy et al. (2021) Roy, D.P., Kashongwe, H.B., Armston, J., 2021.  The impact of geolocation uncertainty on gedi tropical forest canopy height estimation and change monitoring.  Science of Remote Sensing 4, 100024\.
  * Saatchi et al. (2011) Saatchi, S.S., Harris, N.L., Brown, S., Lefsky, M., Mitchard, E.T., Salas, W., Zutta, B.R., Buermann, W., Lewis, S.L., Hagen, S., et al., 2011\.  Benchmark map of forest carbon stocks in tropical regions across three continents.  Proceedings of the national academy of sciences 108, 9899–9904.
  * Samuele et al. (2021) Samuele, D.P., Filippo, S., Orusa, T., Enrico, B.M., 2021\.  Mapping sar geometric distortions and their stability along time: A new tool in google earth engine based on sentinel-1 image time series.  International Journal of Remote Sensing 42, 9135–9154.
  * Santi et al. (2017) Santi, E., Paloscia, S., Pettinato, S., Fontanelli, G., Mura, M., Zolli, C., Maselli, F., Chiesi, M., Bottai, L., Chirici, G., 2017\.  The potential of multifrequency sar images for estimating forest biomass in mediterranean areas.  Remote Sensing of Environment 200, 63–73.
  * Sarker et al. (2012) Sarker, M.L.R., Nichol, J., Ahmad, B., Busu, I., Rahman, A.A., 2012.  Potential of texture measurements of two-date dual polarization palsar data for the improvement of forest biomass estimation.  ISPRS Journal of Photogrammetry and Remote Sensing 69, 146–166.
  * Schneider et al. (2020) Schneider, F.D., Ferraz, A., Hancock, S., Duncanson, L.I., Dubayah, R.O., Pavlick, R.P., Schimel, D.S., 2020.  Towards mapping the diversity of canopy structure from space with gedi.  Environmental Research Letters 15, 115006.
  * Shendryk (2022) Shendryk, Y., 2022.  Fusing gedi with earth observation data for large area aboveground biomass mapping.  International Journal of Applied Earth Observation and Geoinformation 115, 103108.
  * Shi et al. (2015) Shi, B., Gao, W., Jin, G., 2015.  Effects on rhizospheric and heterotrophic respiration of conversion from primary forest to secondary forest and plantations in northeast china.  European journal of soil biology 66, 11–18.
  * Shi and Chen (2018) Shi, H., Chen, J., 2018.  Characteristics of climate change and its relationship with land use/cover change in yunnan province, china.  International journal of climatology 38, 2520–2537.
  * Shimada et al. (2011) Shimada, M., Isoguchi, O., Motooka, T., Shiraishi, T., Mukaida, A., Okumura, H., Otaki, T., Itoh, T., 2011\.  Generation of 10m resolution palsar and jers-sar mosaic and forest/non-forest maps for forest carbon tracking, in: 2011 IEEE International Geoscience and Remote Sensing Symposium, IEEE. pp. 3510–3513.
  * Shimada et al. (2009) Shimada, M., Isoguchi, O., Tadono, T., Isono, K., 2009\.  Palsar radiometric and geometric calibration.  IEEE Transactions on Geoscience and Remote Sensing 47, 3915–3932.
  * Sinha et al. (2015) Sinha, S., Jeganathan, C., Sharma, L.K., Nathawat, M.S., 2015\.  A review of radar remote sensing for biomass estimation.  International Journal of Environmental Science and Technology 12, 1779–1792.
  * Su et al. (2016) Su, Y., Guo, Q., Xue, B., Hu, T., Alvarez, O., Tao, S., Fang, J., 2016.  Spatial distribution of forest aboveground biomass in china: Estimation through combination of spaceborne lidar, optical imagery, and forest inventory data.  Remote Sensing of Environment 173, 187–199.
  * Tucker et al. (2023) Tucker, C., Brandt, M., Hiernaux, P., Kariryaa, A., Rasmussen, K., Small, J., Igel, C., Reiner, F., Melocik, K., Meyer, J., et al., 2023\.  Sub-continental-scale carbon stocks of individual trees in african drylands.  Nature 615, 80–86.
  * Veloso et al. (2017) Veloso, A., Mermoz, S., Bouvet, A., Le Toan, T., Planells, M., Dejoux, J.F., Ceschia, E., 2017.  Understanding the temporal behavior of crops using sentinel-1 and sentinel-2-like data for agricultural applications.  Remote sensing of environment 199, 415–426.
  * Vreugdenhil et al. (2018) Vreugdenhil, M., Wagner, W., Bauer-Marschallinger, B., Pfeil, I., Teubner, I., Rüdiger, C., Strauss, P., 2018.  Sensitivity of sentinel-1 backscatter to vegetation dynamics: An austrian case study.  Remote Sensing 10, 1396\.
  * Wang (2006) Wang, C., 2006.  Biomass allometric equations for 10 co-occurring tree species in chinese temperate forests.  Forest Ecology and Management 222, 9–16.
  * Watson et al. (2018) Watson, J.E., Evans, T., Venter, O., Williams, B., Tulloch, A., Stewart, C., Thompson, I., Ray, J.C., Murray, K., Salazar, A., et al., 2018\.  The exceptional value of intact forest ecosystems.  Nature ecology & evolution 2, 599–610.
  * Woodhouse et al. (2012) Woodhouse, I.H., Mitchard, E.T.A., Brolly, M., Maniatis, D., Ryan, C.M., 2012.  Radar backscatter is not a’direct measure’of forest biomass.  Nature climate change 2, 556–557.
  * Wulder et al. (2012) Wulder, M.A., White, J.C., Nelson, R.F., Næsset, E., Ørka, H.O., Coops, N.C., Hilker, T., Bater, C.W., Gobakken, T., 2012.  Lidar sampling for large-area forest characterization: A review.  Remote sensing of environment 121, 196–209.
  * Xu et al. (2021) Xu, L., Saatchi, S.S., Yang, Y., Yu, Y., Pongratz, J., Bloom, A.A., Bowman, K., Worden, J., Liu, J., Yin, Y., et al., 2021\.  Changes in global terrestrial live biomass over the 21st century.  Science Advances 7, eabe9829.
  * Yan et al. (2021) Yan, J., Xu, Y., Cheng, Q., Jiang, S., Wang, Q., Xiao, Y., Ma, C., Yan, J., Wang, X., 2021.  Lightgbm: accelerated genomically designed crop breeding through ensemble learning.  Genome Biology 22, 1–24.
  * Yu et al. (2011) Yu, D., Zhou, L., Zhou, W., Ding, H., Wang, Q., Wang, Y., Wu, X., Dai, L., 2011.  Forest management in northeast china: history, problems, and challenges.  Environmental Management 48, 1122–1135.
  * Zhang et al. (2019) Zhang, G., Wang, M., Liu, K., 2019.  Forest fire susceptibility modeling using a convolutional neural network for yunnan province of china.  International Journal of Disaster Risk Science 10, 386–403.
  * Zheng et al. (2004) Zheng, D., Rademacher, J., Chen, J., Crow, T., Bresee, M., Le Moine, J., Ryu, S.R., 2004.  Estimating aboveground biomass using landsat 7 etm+ data across a managed landscape in northern wisconsin, usa.  Remote sensing of environment 93, 402–411.
  * Zhu and Liu (2015) Zhu, X., Liu, D., 2015.  Improving forest aboveground biomass estimation using seasonal landsat ndvi time-series.  ISPRS Journal of Photogrammetry and Remote Sensing 102, 222–231.
  * Zhu et al. (2022) Zhu, Z., Deng, X., Zhao, F., Li, S., Wang, L., 2022\.  How environmental factors affect forest fire occurrence in yunnan forest region.  Forests 13, 1392\.