* Devroye and Gyorfi [1985] L. Devroye, L. Gyorfi, Nonparametric Density Estimation: The L1 View, Wiley, New York, 1985\.
  * Sriperumbudur et al. [2010] B. K. Sriperumbudur, A. Gretton, K. Fukumizu, B. Schölkopf, G. R. Lanckriet,  Hilbert space embeddings and metrics on probability measures,  J. Mach. Learn. Res. 11 (2010) 1517–1561.
  * Arthur and Vassilvitskii [2007] D. Arthur, S. Vassilvitskii,  k-means++: The advantages of careful seeding,  in: Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’07, Society for Industrial and Applied Mathematics, New Orleans, LA, USA, 2007, pp. 1027–1035. URL: https://dl.acm.org/doi/10.5555/1283383.1283494.
  * McLachlan and Basford [1988] G. J. McLachlan, K. E. Basford, Mixture Models: Inference and Applications to Clustering, volume 84 of Statistics: Textbooks and Monographs, M. Dekker, New York, NY, USA; Basel, Switzerland, 1988. ISBN-13: 978-0-8247-7691-6.
  * Ng et al. [2002] A. Y. Ng, M. I. Jordan, Y. Weiss,  On spectral clustering: Analysis and an algorithm,  in: T. G. Dietterich, S. Becker, Z. Ghahramani (Eds.), Advances in Neural Information Processing Systems 14, MIT Press, 2002, pp. 849–856. URL: https://papers.nips.cc/paper/2092-on-spectral-clustering-analysis-and-an-algorithm.
  * von Luxburg [2007] U. von Luxburg,  A tutorial on spectral clustering,  Statistics and Computing 17 (2007) 395–416. doi:10.1007/s11222-007-9033-z.
  * Chang and Lin [2011] C.-C. Chang, C.-J. Lin,  LIBSVM: A library for support vector machines,  ACM Trans. Intell. Syst. Technol. 2 (2011). URL: https://doi.org/10.1145/1961189.1961199. doi:10.1145/1961189.1961199.
  * Schmidhuber [2015] J. Schmidhuber,  Deep learning in neural networks: An overview,  Neural Networks 61 (2015) 85–117. URL: https://www.sciencedirect.com/science/article/pii/S0893608014002135. doi:https://doi.org/10.1016/j.neunet.2014.09.003.
  * Bach [2017] F. Bach,  Breaking the curse of dimensionality with convex neural networks,  Journal of Machine Learning Research 18 (2017) 1–53. URL: http://jmlr.org/papers/v18/14-546.html.
  * Branch et al. [1999] M. A. Branch, T. F. Coleman, Y. Li,  A subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems,  SIAM Journal on Scientific Computing 21 (1999) 1–23. URL: https://doi.org/10.1137/S1064827595289108. doi:10.1137/S1064827595289108. arXiv:https://doi.org/10.1137/S1064827595289108.
  * Zhang et al. [2007] J. Zhang, M. Marszałek, S. Lazebnik, C. Schmid,  Local features and kernels for classification of texture and object categories: A comprehensive study,  International Journal Of Computer Vision 73 (2007) 213–238. doi:10.1007/s11263-006-9794-4.
  * Smola et al. [1998] A. J. Smola, B. Schölkopf, K.-R. Müller,  The connection between regularization operators and support vector kernels,  Neural Networks 11 (1998) 637–649. URL: https://www.sciencedirect.com/science/article/pii/S089360809800032X. doi:https://doi.org/10.1016/S0893-6080(98)00032-X.
  * Rudi and Rosasco [2017] A. Rudi, L. Rosasco,  Generalization properties of learning with random features,  in: I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, R. Garnett (Eds.), Advances in Neural Information Processing Systems, volume 30, Curran Associates, Inc., 2017. URL: https://proceedings.neurips.cc/paper_files/paper/2017/file/61b1fb3f59e28c67f3925f3c79be81a1-Paper.pdf.
  * Do Carmo [1976] M. Do Carmo, Differential Geometry of Curves and Surfaces, Prentice Hall, 1976\.
  * Pitkäranta et al. [2001] J. Pitkäranta, A.-M. Matache, C. Schwab,  Fourier mode analysis of layers in shallow shell deformations.,  Computer Methods in Applied Mechanics and Engineering 190 (2001) 2943–2975.
  * Niemi and Pitkäranta [2008] A. H. Niemi, J. Pitkäranta,  Bilinear finite elements for shells: Isoparametric quadrilaterals,  International Journal for Numerical Methods in Engineering 75 (2008) 212–240. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.2252. doi:https://doi.org/10.1002/nme.2252. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.2252.
  * Malinen [2002] M. Malinen,  On the classical shell model underlying bilinear degenerated shell finite elements: general shell geometry,  International Journal for Numerical Methods in Engineering 55 (2002) 629–652.
  * Malinen and Pitkäranta [2000] M. Malinen, J. Pitkäranta,  A benchmark study of reduced-strain shell finite elements: quadratic schemes,  International Journal for Numerical Methods in Engineering 48 (2000) 1637–1671.

##  Appendix A Software and Reproducibility of Results

Version 2.0.0 of the Python package LUQ [47] was developed as part of this
work and is available at https://github.com/CU-Denver-UQ/LUQ. The software
requires numpy ($\geq$1.17.0 and $<$1.23.0), scipy ($\geq$1.7.0 and $<$1.8),
matplotlib ($\geq$3.1.0), and scikit-learn ($>=$1.0.2). The example data,
numerics, and figures presented in this paper can be reproduced using the
Jupyter notebooks and/or Python scripts found within the LUQ GitHub
repository. The illustrative example is found in the example folder named
“wave” and is separated by parts to align with the presentation in this
manuscript. The shell examples are found in the example folder named “shells”
and require the package mat73 for reading the data into memory. See the README
file on github for more information. We note that the examples are
reproducible machine-to-machine, but due to back-end multi-threading
mechanisms within the used packages, there may exist minor deviations between
machines and in comparison to those reported in the manuscript.

##  Appendix B Shell Models

###  B.1 Shell Geometry

Our reference structures are thin shells of revolution. They can formally be
characterized as domains in $\mathbb{R}^{3}$ of type

$\Omega=\left\\{\mathbf{x}+z\mathbf{n(x)}\ | \ \mathbf{x}\in\Gamma,-d/2<z<d/2\right\\},$ (13)

where $d$ is the (constant) thickness of the shell, $\Gamma$ is a (mid)surface
of revolution, and $\mathbf{n(x)}$ is the unit normal to $\Gamma$. The
geometry model is based on principal curvature coordinates, where only four
parameters are needed to specify the curvature and metric on $\Gamma$, namely
the radii of principal curvature $R_{1}$, $R_{2}$, and the Lamé parameters,
$A_{1}$, $A_{2}$, which relate coordinates changes to arc lengths. The
displacement vector field of the midsurface $\mathbf{u}=\\{u,v,w\\}$ is the
set of projections to directions

$\mathbf{e}_{1}=\frac{1}{A_{1}}\frac{\partial\Psi}{\partial
x_{1}},\quad\mathbf{e}_{2}=\frac{1}{A_{2}}\frac{\partial\Psi}{\partial
x_{2}},\quad\mathbf{e}_{3}=\mathbf{e}_{1}\times\mathbf{e}_{2},$ (14)

where $\Psi(x_{1},x_{2})$ is a suitable parametrization of the surface of
revolution, $\mathbf{e}_{1},\mathbf{e}_{2}$ are the unit tangent vectors along
the principal curvature lines, and $\mathbf{e}_{3}$ is the unit normal. In
other words

$\mathbf{u}=u\,\mathbf{e}_{1}+v\,\mathbf{e}_{2}+w\,\mathbf{e}_{3}.$

####  B.1.1 Profile functions and parametrization

When a plane curve is rotated (in three dimensions) about a line in the plane
of the curve, it sweeps out a surface of revolution. Consider a plane curve
rotating about the $x$-axis, the so-called profile function in the $xy$-plane,
$y=f(x)$, which induces a surface $\Gamma_{f}$.

Let $I=[\alpha,\beta]\subset\mathbb{R}$ be a bounded closed interval, and let
$f(x):\>I\rightarrow\mathbb{R}^{+}$ be a regular function. The shell
midsurface $\Gamma_{f}$ is parameterized by means of the mapping

$\displaystyle\bm{\Psi}_{f}\ :\ I\times[0,2\pi]\ \longrightarrow\
\mathbb{R}^{3}$ (15)
$\displaystyle\bm{\Psi}_{f}(x_{1},x_{2})=(x_{1},f(x_{1})\cos{x_{2}},f(x_{1})\sin{x_{2}})\
.$

For $\Gamma_{f}$ we have

$A_{1}(x)=\sqrt{1+[f^{\prime}(x)]^{2}},\quad A_{2}(x)=f(x),\quad
R_{1}(x)=-\frac{A_{1}(x)^{3}}{f^{\prime\prime}(x)},\quad
R_{2}(x)=A_{1}(x)A_{2}(x).$ (16)

###  B.2 Two-Dimensional Model

The shell defined as $\Omega\subset\mathbb{R}^{3}$ above is a three-
dimensional body for which linear elasticity could be considered “exact” for
small deformations. Classical dimension reduction models are reasonably
accurate for these shells [66]. We consider a model allowing for continuous
solutions, namely the Reissner-Mindlin model (often attributed to Naghdi). The
displacement field $\mathbf{u}$ has two additional components $\theta$,
$\psi$, that is, it has five components $u$, $v$, $w$, $\theta$, $\psi$, each
of which is a function of two variables on the mid-surface of the shell. The
first two components represent the tangential displacements of the mid-
surface, $w$ is the transverse deflection, and whereas $\theta$, $\psi$ are
dimensionless rotations. In this section we assume that the shell consists of
homogeneous isotropic material with Young modulus $E$ and Poisson ratio $\nu$.

The total energy of the shell has the following expression in our dimension
reduction model:

$\mathcal{F}(\mathbf{u})=\frac{1}{2}D(a(\mathbf{u},\mathbf{u})+d^{2}\,b(\mathbf{u},\mathbf{u}))-q(\mathbf{u}),$
(17)

where $D=E\,d/(12(1-\nu^{2}))$ is a scaling factor (see for instance [67]),
$q$ is the external load potential, and $a(\mathbf{u},\mathbf{u})$ and
$b(\mathbf{u},\mathbf{u})$ represent the portions of total deformation energy
that are stored in membrane and transverse shear deformations and bending
deformations, respectively. The latter quadratic forms are independent of $d$
and defined as

$\displaystyle a(\mathbf{u},\mathbf{u})$ $\displaystyle=$ $\displaystyle
a_{m}(\mathbf{u},\mathbf{u})+a_{s}(\mathbf{u},\mathbf{u})$ (18)
$\displaystyle=$ $\displaystyle
12\int_{\omega}\Bigl{[}\nu(\beta_{11}(\mathbf{u})+\beta_{22}(\mathbf{u}))^{2}+(1-\nu)\sum_{i,j=1}^{2}\beta_{ij}(\mathbf{u})^{2}\Bigr{]}\
A_{1}A_{2}\ d\gamma+$ $\displaystyle
6(1-\nu)\int_{\omega}\Bigl{[}(\rho_{1}(\mathbf{u})^{2}+\rho_{2}(\mathbf{u}))^{2}\Bigr{]}\
A_{1}A_{2}\ d\gamma,$ $\displaystyle b(\mathbf{u},\mathbf{u})$
$\displaystyle=$
$\displaystyle\int_{\omega}\Bigl{[}\nu(\kappa_{11}(\mathbf{u})+\kappa_{22}(\mathbf{u}))^{2}+(1-\nu)\sum_{i,j=1}^{2}\kappa_{ij}(\mathbf{u})^{2}\Bigr{]}\
A_{1}A_{2}\ d\gamma,$ (19)

where $\beta_{ij}$, $\rho_{i}$, and $\kappa_{ij}$ stand for the membrane,
transverse shear, and bending strains, respectively. The strain-displacement
relations are linear and involve at most first derivatives of the displacement
components and hence are suitable for standard finite element methods.

Using the identities above, the bending, membrane, and shear strains [68],
$\kappa_{ij},\beta_{ij},$ and $\rho_{i}$, respectively, can be written as

$\displaystyle\kappa_{11}$ $\displaystyle=$
$\displaystyle\frac{1}{A_{1}}\frac{\partial\theta}{\partial
x}+\frac{\psi}{A_{1}A_{2}}\frac{\partial A_{1}}{\partial y},$
$\displaystyle\kappa_{22}$ $\displaystyle=$
$\displaystyle\frac{1}{A_{2}}\frac{\partial\psi}{\partial
y}+\frac{\theta}{A_{1}A_{2}}\frac{\partial A_{2}}{\partial x},$
$\displaystyle\kappa_{12}$ $\displaystyle=$
$\displaystyle\kappa_{21}=\frac{1}{2}\left[\frac{1}{A_{1}}\frac{\partial\psi}{\partial
x}+\frac{1}{A_{2}}\frac{\partial\theta}{\partial
y}-\frac{\theta}{A_{1}A_{2}}\frac{\partial A_{1}}{\partial
y}-\frac{\psi}{A_{1}A_{2}}\frac{\partial A_{2}}{\partial x}\right.$
$\displaystyle-\frac{1}{R_{1}}\left(\frac{1}{A_{2}}\frac{\partial u}{\partial
y}-\frac{v}{A_{1}A_{2}}\frac{\partial A_{2}}{\partial x}\right)$
$\displaystyle-\left.\frac{1}{R_{2}}\left(\frac{1}{A_{1}}\frac{\partial
v}{\partial x}-\frac{u}{A_{1}A_{2}}\frac{\partial A_{1}}{\partial
y}\right)\right],$

$\displaystyle\beta_{11}$ $\displaystyle=$
$\displaystyle\frac{1}{A_{1}}\frac{\partial u}{\partial
x}+\frac{v}{A_{1}A_{2}}\frac{\partial A_{1}}{\partial y}+\frac{w}{R_{1}},$
$\displaystyle\beta_{22}$ $\displaystyle=$
$\displaystyle\frac{1}{A_{2}}\frac{\partial v}{\partial
y}+\frac{u}{A_{1}A_{2}}\frac{\partial A_{2}}{\partial x}+\frac{w}{R_{2}},$
$\displaystyle\beta_{12}$ $\displaystyle=$
$\displaystyle\beta_{21}=\frac{1}{2}\left(\frac{1}{A_{1}}\frac{\partial
v}{\partial x}+\frac{1}{A_{2}}\frac{\partial u}{\partial
y}-\frac{u}{A_{1}A_{2}}\frac{\partial A_{1}}{\partial
y}-\frac{v}{A_{1}A_{2}}\frac{\partial A_{2}}{\partial x}\right),$
$\displaystyle\rho_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{A_{1}}\frac{\partial w}{\partial
x}-\frac{u}{R_{1}}-\theta,$ $\displaystyle\rho_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{A_{2}}\frac{\partial w}{\partial
y}-\frac{v}{R_{2}}-\psi.$

######  Remark 1.

When the shell parametrisations defined above are used, all terms of the form
$\partial A_{i}/\partial y$ are identically zero.

The energy norm $|||\cdot|||$ is defined in a natural way in terms of the
deformation energy:

$\mathcal{E}(\mathbf{u}):=|||\mathbf{u}|||^{2}=a(\mathbf{u},\mathbf{u})+d^{2}\,b(\mathbf{u},\mathbf{u}).$
(20)

Similarly for bending, membrane, and shear energies:

$B(\mathbf{u}):=d^{2}\,b(\mathbf{u},\mathbf{u}),\quad
M(\mathbf{u}):=a_{m}(\mathbf{u},\mathbf{u}),\quad
S(\mathbf{u}):=a_{s}(\mathbf{u},\mathbf{u}).$ (21)

The load potential has the form

$q(\mathbf{v})=\int_{\omega}\mathbf{f}(x,y)\,\cdot\mathbf{v}\,A_{1}A_{2}\
dx\,dy.$

We are interested in problems where the load acts in the transverse direction
of the shell surface, i.e., $\mathbf{f}(x,y)=[0,0,f_{w}(x,y),0,0]^{T}$. It can
be shown that if for the load $\mathbf{f}\in[L^{2}(\omega)]^{5}$ holds, the
variational problem has a unique weak solution
$\mathbf{u}\in[H^{1}(\omega)]^{5}$. The corresponding result is true in the
finite dimensional case, when the finite element method is employed.

The existence of the asymptotic displacement field corresponding to the limit
$t\to 0$ is possible if a suitable scaling of the load amplitude is done in
the limit process. For discussion on details, see Malinen and Pitkäranta [69].

###  B.3 One-dimensional Model

The shell model can be further reduced to a one-dimensional one, if the load
is periodic in the angular direction and the Young modulus is constant in the
angular direction, that is, $E=E(x)$. Due to the periodic structure one can
choose between two periodic transverse loadings acting on the $w$-component,
either $(f_{1})_{w}(x,y)=F(x)\cos(k\,y)$ or $(f_{2})_{w}(x,y)=F(x)\sin(k\,y)$,
$k=0,1,2,\ldots$, and then the corresponding solution fields $\mathbf{u}(x,y)$
are either

$\mathbf{u}_{1}(x,y)=\left(\begin{array}[]{c}u(x)\cos(k\,y)\\\
v(x)\sin(k\,y)\\\ w(x)\cos(k\,y)\\\ \theta(x)\cos(k\,y)\\\
\psi(x)\sin(k\,y)\\\ \end{array}\right)\text{ or
}\mathbf{u}_{2}(x,y)=\left(\begin{array}[]{c}u(x)\sin(k\,y)\\\
v(x)\cos(k\,y)\\\ w(x)\sin(k\,y)\\\ \theta(x)\sin(k\,y)\\\
\psi(x)\cos(k\,y)\\\ \end{array}\right).$ (22)

######  Remark 2.

The two ansatz simply represent the trigonometric basis in the angular
direction.

Using the ansatz $\mathbf{u}_{1}(x,y)$ above the energies are defined over an
interval $I$. Here the three energies are defined separately, with $B$, $M$,
$S$, denoting bending, membrane, and shear, respectively.

$\displaystyle d^{2}\mathcal{A}_{B}^{1D}(\mathbf{u},\mathbf{u})$
$\displaystyle=$ $\displaystyle
d^{2}\int_{I}E(x)\Bigl{[}\nu(\kappa_{11}(\mathbf{u})+\kappa_{22}(\mathbf{u}))^{2}+(1-\nu)\sum_{i,j=1}^{2}\kappa_{ij}(\mathbf{u})^{2}\Bigr{]}\
A_{1}A_{2}\ dx,$ (23)
$\displaystyle\mathcal{A}_{M}^{1D}(\mathbf{u},\mathbf{u})$ $\displaystyle=$
$\displaystyle
12\int_{I}E(x)\Bigl{[}\nu(\beta_{11}(\mathbf{u})+\beta_{22}(\mathbf{u}))^{2}+(1-\nu)\sum_{i,j=1}^{2}\beta_{ij}(\mathbf{u})^{2}\Bigr{]}\
A_{1}A_{2}\ dx,$ (24)
$\displaystyle\mathcal{A}_{S}^{1D}(\mathbf{u},\mathbf{u})$ $\displaystyle=$
$\displaystyle
6(1-\nu)\int_{I}E(x)\Bigl{[}(\rho_{1}(\mathbf{u})^{2}+\rho_{2}(\mathbf{u}))^{2}\Bigr{]}\
A_{1}A_{2}\ dx.$ (25)

The strains can be written in terms of the wave number $k$:

$\displaystyle\kappa_{11}$ $\displaystyle=$
$\displaystyle\frac{1}{A_{1}}\frac{\partial\theta}{\partial x},$
$\displaystyle\kappa_{22}$ $\displaystyle=$
$\displaystyle\frac{k\,\psi}{A_{2}}+\frac{\theta}{A_{1}A_{2}}\frac{\partial
A_{2}}{\partial x},$ $\displaystyle\kappa_{12}$ $\displaystyle=$
$\displaystyle\frac{1}{2}[\frac{1}{A_{1}}\frac{\partial\psi}{\partial
x}-\frac{k\,\theta}{A_{2}}-\frac{\psi}{A_{1}A_{2}}\frac{\partial
A_{2}}{\partial x}$
$\displaystyle+\frac{1}{R_{1}}(\frac{k\,u}{A_{2}}+\frac{v}{A_{1}A_{2}}\frac{\partial
A_{2}}{\partial x})-\frac{1}{R_{2}}\frac{1}{A_{1}}\frac{\partial v}{\partial
x}]=\kappa_{21},$ $\displaystyle\beta_{11}$ $\displaystyle=$
$\displaystyle\frac{1}{A_{1}}\frac{\partial u}{\partial x}+\frac{w}{R_{1}},$
$\displaystyle\beta_{22}$ $\displaystyle=$
$\displaystyle\frac{k\,v}{A_{2}}+\frac{u}{A_{1}A_{2}}\frac{\partial
A_{2}}{\partial x}+\frac{w}{R_{2}},$ $\displaystyle\beta_{12}$
$\displaystyle=$ $\displaystyle\frac{1}{2}(\frac{1}{A_{1}}\frac{\partial
v}{\partial x}-\frac{k\,u}{A_{2}}-\frac{v}{A_{1}A_{2}}\frac{\partial
A_{2}}{\partial x})=\beta_{21},$ $\displaystyle\rho_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{A_{1}}\frac{\partial w}{\partial
x}-\frac{u}{R_{1}}-\theta,$ $\displaystyle\rho_{2}$ $\displaystyle=$
$\displaystyle-\frac{k\,w}{A_{2}}-\frac{v}{R_{2}}-\psi.$

The load potential becomes

$\mathcal{Q}^{1D}(\mathbf{v})=\int_{\alpha}^{\beta}F(x)\,w\,A_{1}A_{2}\,dx.$

Notice that we have omitted the coefficient $\pi$ ($2\pi$, when $k=0$) from
the load potential and the energy definitions (23)–(25). This has to be taken
into account whenever one- and two-dimensional energies are compared.


# Using Neural Networks for Fast SAR Roughness Estimation of High Resolution
Images

Li Fan, and Jeova Farias Sales Rocha Neto J. F. S. Rocha Neto is with the
Department of Computer Science, Haverford College, Haverford, PA, 19041 USA
e-mail<EMAIL_ADDRESS>Fan is with Department of Statistical
Science, Duke University, Durham, NC, 27708, USA. e-mail<EMAIL_ADDRESS>
###### Abstract

The analysis of Synthetic Aperture Radar (SAR) imagery is an important step in
remote sensing applications, and it is a challenging problem due to its
inherent speckle noise. One typical solution is to model the data using the
$G_{I}^{0}$ distribution and extract its roughness information, which in turn
can be used in posterior imaging tasks, such as segmentation, classification
and interpretation. This leads to the need of quick and reliable estimation of
the roughness parameter from SAR data, especially with high resolution images.
Unfortunately, traditional parameter estimation procedures are slow and prone
to estimation failures. In this work, we proposed a neural network-based
estimation framework that first learns how to predict underlying parameters of
$G_{I}^{0}$ samples and then can be used to estimate the roughness of unseen
data. We show that this approach leads to an estimator that is quicker, yields
less estimation error and is less prone to failures than the traditional
estimation procedures for this problem, even when we use a simple network.
More importantly, we show that this same methodology can be generalized to
handle image inputs and, even if trained on purely synthetic data for a few
seconds, is able to perform real time pixel-wise roughness estimation for high
resolution real SAR imagery.

###### Index Terms:

Synthetic Aperture Radar Images, Neural Networks, Image Analysis, Statistical
Modeling.

##  I Introduction

Synthetic Aperture Radar (SAR) data is crucial to remote sensing and earth
monitoring. Its ability to capture high resolution snapshots of targets and
landscapes independently of the weather conditions and sunlight has opened the
way to important advancements in environmental monitoring, emergency response,
evaluation of damages in natural catastrophes, urban planning and ecology, to
mention a few. Because of its widespread use, there is a considerable
industrial appeal for SAR image understanding algorithms. This task, however,
is challenging due to the degrading speckle noise inherent to such data, which
prevents the application of image processing techniques that are common in
other image domains.

Thankfully, this same appeal led to the developments of statistical models
that describe the SAR data despite its noise pattern. In particular, the
$G_{I}^{0}$ [1] distribution for intensity data, highly successful in both
theory and practice [2], arouse to prominence in the remote sensing community
because of its capacity to model well textured, extremely textured and
textureless terrain data. To do so, this model relies on three parameters, one
of which is roughness. It directly corresponds to the captured texture and has
been extensively used for SAR data understanding [3, 4, 5, 6, 7]. Estimating
such parameter became crucial in many SAR imaging techniques [7, 4]. However,
the best current available algorithms for such a task, i.e. Maximum Likelihood
[2], Log-cumulant [8] and Minimum Distance Estimators [5], are too
computationally expensive to be applied to high resolution images within
reasonable time.

Taking a different approach, deep learning has been successfully applied to
SAR [9] and many other imaging domains [10]. In these techniques, one
typically relies on training a large neural network on a substantial amount of
labeled data, or augmentations thereof. Practically speaking, this
considerably sized process can be accomplished by implementing parallelizable
learning techniques in powerful parallel processing machines, such as Graphics
Processing Units (GPUs). Nowadays, many efficient and simple-to-use libraries
are available for manipulating data, and training deep networks on GPUs, which
prompts researchers and practicians alike to design their methodologies in a
way to leverage these powerful tools.

In this work, we propose using neural network-based learning algorithm for
estimating roughness maps in SAR images. In our methodology, we generate sets
of $G_{I}^{0}$ samples of varying sizes, compute their sample moments and use
them to predict the parameters that generated the data. Once the network is
trained on this fully synthetic data, we proceed with network inference on
real data. To the best of our knowledge, our work is the first of its kind to
train a neural network to estimate parameters of distributions based on sample
moments from synthetic data. As we shall show, in fact, a very small network
composed, trained for a few seconds on a small dataset, is sufficient to
produce better estimates on unseen data than maximum likelihood and log-
cumulant-based estimators in terms of mean squared error, while being also
faster and less prone to failures. Finally, this same process can be fully
implemented, from sample moment computation and training to network inference,
on a GPU for high resolution image-sized input, with pixel-wise roughness
estimation accomplished in a few milliseconds on average. Overall, we propose
a simple, principled and efficient method for training neural networks for SAR
image understanding, using the statistical baggage this type of data carries.

##  II Preliminaries

###  II-A SAR Statistical Modeling and Parameter Estimation

The intensity return $Z$ in monopolarized SAR image can be effectively modeled
by the product of two independent random variables: $X$, the backscatter data,
and $Y$, the inherent speckle noise [2]. Typically, $Y$ is modeled as a
unitary-mean Gamma distributed random variable with shape parameter
corresponding to the number of looks $L\in\mathbb{N}^{*}$ used to capture the
data, and considered known during estimation [2]. Assuming that $X$ obeys the
reciprocal of Gamma law, Frery et Al. [1] showed that $Z\sim
G_{I}^{0}(\alpha,\gamma,L)$ with density:

$f_{\mathcal{G}^{0}_{I}}(z,\theta)=\frac{L^{L}\Gamma(L-\alpha)}{\gamma^{\alpha}\Gamma(-\alpha)\Gamma(L)}z^{L-1}(\gamma+Lz)^{\alpha-L},$
(1)

where $-\alpha,\gamma>0$ are called roughness and scale, respectively.

The $G_{I}^{0}$ model has become a popular choice for SAR data due to its
expressiveness [1], mathematical tractability [2], and effectiveness when
applied to various imaging tasks [11, 4]. In fact, it was shown that the
roughness parameter $\alpha$ plays an important role in SAR understanding [5].
With roughness close to zero, typically $\alpha>-3$, it suggests the targeted
region is highly textured, therefore urban. For $\alpha\in[-3,-6]$, we have
evidence for moderately textured regions, which can correspond to forests.
Finally, $\alpha<-6$ implies textureless areas, such as seas and pastures.
Hence, using $\alpha$ estimates instead of pixel intensities on SAR images led
to improvements in image segmentation [7, 3] and region discrimination [12].

This direct application of SAR statistical modelling fueled the study and
improvement of estimators for the $G_{I}^{0}$ distribution [12, 5]. While
early work employed a Maximum Likelihood Estimator (MLE) [2] for parameter
estimation, later developments led to the application of Method of Moments [2]
, Method of Log Cumulants (LCUM) [8, 6, 7] and Minimum-Distance Estimators [5,
4] to this problem. Despite some of their successes, all these methods,
especially LCUM, heavily rely on slow optimization procedures and are prone to
high estimation failures [5]. These issues hinder their application in image-
sized data understanding, as the current praxis relies on generating roughness
maps by sweeping the image with a window that collects the intensity data
centered at each pixel and then proceeds with the parameter estimation
algorithm [7], a process that can be excessively slow.

###  II-B Neural Networks

We focus on supervised learning algorithms that, given a dataset
$\mathcal{D}=\\{(\mathbf{x}_{i},\mathbf{y}_{i})\\}_{i=1}^{N}$, aim to compute
a model $f_{\theta}(\cdot)$ parameterized by $\theta$ that best predicts the
output $\mathbf{y}_{i}$’s when given the input $\mathbf{x}_{i}$’s. The
prediction error, to be minimized, is measured via a loss, typically the Mean
Squared Error (MSE):

$\operatorname{MSE}=\frac{1}{N}\sum_{(\mathbf{x},\mathbf{y})\in\mathcal{D}}\lVert\mathbf{y}-f_{\theta}(\mathbf{x}))\rVert^{2}_{2},$
(2)

where $\lVert\cdot\rVert_{2}$ is the Euclidean vector norm. In Neural
Networks, we assume that $f_{\theta}$ is composed of units of a non-linear
activations $\sigma(\cdot)$ and parameters
$\theta=\\{\mathcal{W},\mathcal{B}\\}$, where
$\mathcal{W}=\\{W_{i}\\}_{i=1}^{Q}$ are weight matrices and
$\mathcal{B}=\\{\mathbf{b}_{i}\\}_{i=1}^{Q}$ a set of bias vectors. The
Multilayer Perceptron (MLP), a popular network, uses a model
$NN_{\theta}(\mathbf{x})$ formed by the composition of affine transformations
of the inputs followed by the element-wise application of $\sigma(\cdot)$:

$NN_{\theta}(\mathbf{x})=\sigma(W_{Q}\ldots\sigma(W_{2}\sigma(W_{1}\mathbf{x}+\mathbf{b}_{1})+\mathbf{b}_{2})\ldots+\mathbf{b}_{Q}).$
(3)

Each transformation is typically referred to as a layer of the MLP, where the
first and last are the input and output layers, and others are called hidden
layers. A network with multiple hidden layers is commonly known as a “deep”
neural network.

Another, arguably more popular, neural network model is the Convolutional
Neural Network (CNN). Different from the MLP model, in CNN, the weights in
each layer are shared so that they effectively perform a convolution operation
on the input. For example, consider an image $I$ of shape $c\times w\times l$,
where $w$ and $l$ are 2D spatial dimensions of $I$, and $c$ is its number of
channels (three for RGB data). One can design a matrix of weights whose
application on a vectorized/flattened version of $I$ corresponds to a
convolution of $I$ on its original shape and a kernel matrix of shape $c\times
k\times k$, where $k\ll w,l$. Using such convolutional layers brings several
advantages for learning on image-sized data. First, the number of weights to
be learned is drastically reduced, from $O(cwl)$ to $O(ck^{2})$. Second, the
number of weights to be learned no longer depends on the size of the input.
Third, one can show that this process learns separate feature detectors, which
improves the network performance on important imaging tasks such as detection,
denoising and segmentation [10].

More importantly to us is the setting where MLPs and CNNs are equivalent.
Suppose we have the input vectors of the dataset $\mathcal{D}$ in
$\mathbb{R}^{d}$ and output vectors in $\mathbb{R}^{d^{\prime}}$. Learning a
dense layer for $\mathcal{D}$ would require finding a matrix
$W\in\mathbb{R}^{d\times d^{\prime}}$. Now, consider the tensors
$X\in\mathbb{R}^{d\times w\times l}$ and $Y\in\mathbb{R}^{d^{\prime}\times
w\times l}$, where $w$ and $l$ are chosen such that $wl=N$, and each
$\mathbf{x}_{n}$ (resp. $\mathbf{y}_{n}$) is placed in a transversal tube
$X[:,i,j]$ (resp. $Y[:,i,j]$)111We use a notation similar to Numpy’s for array
indexing (for example, $X[i,j]$ is the element in $X$ at indices $(i,j)$) and
slicing (or example, $X[i,:]$ is all the elements in the $i$-th row of
$X$).[13]. It is easy to see that training $d^{\prime}$ kernel matrices of
shape $d\times 1\times 1$ is the same as learning the values in $W$ [14]. As
discussed in the next section, using these $1\times 1$ convolutional layers
will enable us to see our parameter estimation procedure as an inference on a
pre-trained CNN, leveraging the speed of highly parallel computing engines
during estimation and preventing us from using time-consuming loops for
generating roughness maps.

##  III Proposed Methodology

###  III-A Estimating $G_{I}^{0}$ Using Neural Networks

Our approach for estimation relies on neural network-based learning procedure.
Let $L\in\mathbb{N}^{*}$ be known [2] and
$\mathcal{H}=\\{(\mathcal{Z}_{i},\alpha_{i})\\}_{i=1}^{N}$ be a dataset
composed of samples $\mathcal{Z}_{i}$ from a
$G_{I}^{0}(\alpha_{i},\gamma_{i},L)$ distribution of varying sizes, where we
assume $\gamma_{i}=-\alpha_{i}-1$ 222We follow the praxis of assuming that the
samples have unit mean [2].. We also assume that $\alpha_{i}$ is always above
a given lower bound (set to $-15$ in our experiments) [5].

Our goal is to train a neural network to predict $\alpha_{i}$ from
$\mathcal{Z}_{i}$. However, the data in $\mathcal{Z}_{i}$ cannot be used
directly, as the input vector of our network needs to be of fixed size. The
final estimator also needs to be invariant to sample permutations in each
input set. To solve these issues, the authors in [15] proved that that all
permutation invariant functions $t(\cdot)$ on a set $\mathcal{S}$ can be
represented as
$t(\mathcal{S})=\rho\left(\sum_{s\in\mathcal{S}}\phi(s)\right)$, where
$\rho(\cdot)$ and $\phi(\cdot)$ are continuous functions. Setting
$\rho(x)=\frac{1}{|\mathcal{Z}_{i}|}x$ and $\phi(x)=(\log x)^{m}$ and applying
the resulting $t(\cdot)$ on a sample $\mathcal{Z}_{i}$ from $\mathcal{D}$, we
get:

$\mu_{m}(\mathcal{Z}_{i})=\frac{1}{|\mathcal{Z}_{i}|}\sum_{z\in\mathcal{Z}_{i}}(\log
z)^{m},$ (4)

where $\mu_{m}(\mathcal{Z}_{i})$ conveniently represents the sample log-moment
of order $m$, a statistic commonly used in SAR literature [6, 7]. In practice,
we create a vector of $N_{m}$ moments
$\boldsymbol{\mu}_{i}=[\mu_{1}(\mathcal{Z}_{i})\,\,\mu_{2}(\mathcal{Z}_{i})\,\,\ldots\,\,\mu_{N_{m}}(\mathcal{Z}_{i})]$
and train our network on the dataset
$\mathcal{D}=\\{(\boldsymbol{\mu}_{i},\alpha_{i})\\}_{i=0}^{N}$. We use an MLP
as our network architecture, since its structure can be adapted to infer full
sized roughness maps, as explained in the next section.

Input: $L$: # looks. $R$: dataset size. $T$: # epochs. $N_{m}$: # log-mom.

$\mathcal{A}$: set of $\alpha$ values, $\mathcal{S}$: set of sample sizes.

Function  _Main(__)_ :

$\mathcal{D}\leftarrow$
SyntheticDataset(_$L,w,h,\mathcal{A},\mathcal{K},N_{m},R$_)

$NN_{\theta}\leftarrow$ TrainNeuralNet(_$\mathcal{D},T$_)

Function  _$\mathcal{D}=$
SyntheticDataset(_$L,\mathcal{A},\mathcal{S},N_{m},R$_)_:

$\mathcal{D}\leftarrow$ Empty set

repeat  _R_ times

$\alpha,s\leftarrow$ Sample uniformly from $\mathcal{A}$ and $\mathcal{S}$,
resp.

$\gamma\leftarrow-\alpha-1$

$\mathcal{Z}\leftarrow\\{z|z\sim
G_{I}^{0}(\alpha,\gamma,L),|\mathcal{Z}|=s\\}$,

$\boldsymbol{\mu}\leftarrow$ ComputeMoments(_$\mathcal{Z},N_{m}$_)

$\mathcal{D}\leftarrow\mathcal{D}\cup\\{(\boldsymbol{\mu},\alpha)\\}$

Function  _$\boldsymbol{\mu}=$ ComputeMoments(_$\mathcal{Z},N_{m}$_)_:

foreach  _$m\in\\{1,2,\ldots,N_{m}\\}$_ do

$\boldsymbol{\mu}[m]\leftarrow\frac{1}{|\mathcal{Z}|}\sum_{z\in\mathcal{Z}}(\log
z)^{m}$

Function  _$NN_{\theta}=$ TrainNeuralNet(_$\mathcal{D}$ , $T$_)_:

$NN_{\theta}\leftarrow$ Multilayer Perceptron with parameters $\theta$.

repeat  _T_ times

foreach  _batch $(\mathbf{x},\mathbf{y})\in\mathcal{D}$_ do

Optimize $\theta$ to minimize
$\lVert\mathbf{y}-NN_{\theta}(\mathbf{x})\rVert_{2}^{2}$

Algorithm 1 Training Roughness Estimator

Input: $L$: # looks. $w$ and $h$: image dimensions.

$R$: dataset size, $T$: # epochs. $N_{m}$: # log-moments.

$\mathcal{A}$: set of $\alpha$ values, $\mathcal{K}$: set of kernel sizes,

Function  _Main(__)_ :

$\mathcal{D}\leftarrow$
SyntheticDataset(_$L,w,h,\mathcal{A},\mathcal{K},N_{m},R$_)

$NN_{\theta}\leftarrow$ TrainNeuralNet(_$\mathcal{D},T$_)

Function  _$\mathcal{D}=$
SyntheticDataset(_$L,w,h,\mathcal{A},\mathcal{K},N_{m},R$_)_:

$\mathcal{D}\leftarrow$ Empty set

repeat  _R_ times

$\alpha,k\leftarrow$ Sample uniformly from $\mathcal{A}$ and $\mathcal{K}$,
resp.

$\gamma\leftarrow-\alpha-1$

$I\leftarrow\text{Matrix in }\mathbb{R}^{w\times h},\text{with }I[i,j]\sim
G_{I}^{0}(\alpha,\gamma,L)$

$A\leftarrow\text{Matrix in }\mathbb{R}^{w\times h},\text{with }A[i,j]=\alpha$

$M\leftarrow$ ComputeImgMoments(_$I,N_{m},k,\alpha,\gamma,L$_)

$\mathcal{D}\leftarrow\mathcal{D}\cup\\{(M,A)\\}$

Function  _$M=$ ComputeImgMoments(_$I,N_{m},k,\alpha,\gamma,L$_)_:

$I\leftarrow I+$padding of $G_{I}^{0}(\alpha,\gamma,L)$ samples.

foreach  _$m\in\\{1,2,\ldots,N_{m}\\}$_ do

$M[i,:,:]\leftarrow\operatorname{AvgPool}(\log(I)^{m},k\times k\text{
kernel})$

Function  _$NN_{\theta}=$ TrainNeuralNet(_$\mathcal{D}$ , $T$_)_:

$NN_{\theta}\leftarrow$ FCNN with only $1\times 1$ convolutions.

repeat  _T_ times

foreach  _batch $(X,Y)\in\mathcal{D}$_ do

Optimize $\theta$ to minimize $\lVert Y-NN_{\theta}(X)\rVert_{2}^{2}$

Algorithm 2 Training Roughness Map Estimator

Algorithm 1 gives an overview of this training procedure. Following Main(), a
synthetic $G_{I}^{0}$ dataset is created by drawing samples of various
predefined sample sizes and $\alpha$ values. For each set, ComputeMoments()
runs on each sample to find its first $N_{m}$ log-moments and concatenate the
resulting moment vector to the roughness value that generated it. Finally, an
MLP is trained on these pairs of moment vectors and the $\alpha$ parameter
over a given amount of epochs.

After the network is trained, one can estimate the parameters of a given
unseen sample set by (1) computing its moments via ComputeMoments() and (2)
feeding them through the trained model. We hope that, although training
requires extra computational time, the inference step in this process is fast.

###  III-B Adaptation to Roughness Map Estimation

In many practical scenarios involving parameter estimation in SAR, one wishes
to estimate the roughness of all pixel locations in a potentially high
resolution image $I$ [7]. While we could apply the algorithm described
previously to all $k\times k$ sized windows in $I$ for a given $k>0$, in this
section we show the functions described in Algorithm 1 can be easily adapted
to image size data and be fully implemented on a GPU. This transition will
enable us to fully parallelize our estimation algorithm in both training and
inference phases, which will consequently highly reduce our computation time.

In this new setting, let $I\in\mathbb{R}^{w\times h}$. For our training phase,
each input in our training set will be composed of samples from the
$G_{I}^{0}$ distribution for unique values of $\alpha$, $\gamma$ and $L$
disposed in a $w\times h$ matrix. In the inference phase, $I$ is a real SAR
image of any size, meaning that our method will be able to estimate the
roughness of an image with a network trained on purely synthetic and random
data.

Our goals are (1) to efficiently compute all the desired $N_{m}$ moments of
the data surrounding each pixel in $I$ on a $k\times k$ window, composing a
tensor $M\in\mathbb{R}^{N_{n}\times w\times h}$, and (2) to feed this data to
an appropriate network that estimates the roughness parameters of each pixel
location on the image grid. For step (1), one can use an Average Pooling
Layer, $\operatorname{AvgPool}(\cdot)$, where a convolution kernel sweeps the
input data, averaging the pixel values within that window [16]. Here we add an
appropriate padding composed of $G_{I}^{0}$ samples to
$\operatorname{AvgPool}$’s input, so its output size remains the same. Now, if
the input data is $\log(I)^{m}$, where the exponentiation is computed pixel-
wise, one can estimate each pixel’s log-moments of order $m$ via
$\operatorname{AvgPool}(\log(I)^{m})$. Each channel in $M$ is computed by
applying this procedure to all $m\in\\{1,2,\cdots,N_{m}\\}$.

For step (2), we can use the connection between MLPs and $1\times 1$
convolutional networks described in Section II-B and turn the MLP architecture
used in Algorithm 1 in to a Fully Convolutional Neural Network (FCNN). In that
network, the input of its convolutions will be the transversal tubes on $M$,
which correspond to the sample log-moments of each pixel in $I$. During
training, its output will be a matrix $A\in\mathbb{R}^{w\times h}$, such that
$A[i,j]=\alpha,\forall i,j$, where $\alpha$ is the roughness parameter that
generated $I$ along with $\gamma=-\alpha-1$ and a predetermined $L$. Algorithm
2 explains this training algorithm in more detail and is analogous to
Algorithm 1 in its execution.

##  IV Numerical Experiments and Discussion

Figure 1: MSE values for the synthetic data. Missing values on graphs
correspond to the method failing in all simulations.

###  IV-A Data and Algorithmic Setup and Assessment Methodology

We assess the performance of our proposed methodology on two settings. The
first is a synthetic one where we aim at estimating the roughness parameter
from samples of the $G_{I}^{0}$ distribution. We qualitatively compare the
resulting network from Algorithm 1 to the standard estimators for this
problem. In our second experiment, we qualitatively evaluate a model trained
using Algorithm 2 on a real SAR image.

To sample synthetic SAR data from the $G_{I}^{0}$ model, we make use of its
multiplicative nature and sample $X\sim\Gamma(1,L)$ and
$Y^{\prime}\sim\Gamma(-\alpha,\gamma)$ to get $Z=X/Y^{\prime}\sim
G_{I}^{0}(\alpha,\gamma,L)$ [2]. In our real experiments, we use a $1500\times
1500$ SAR image acquired by Airborne SAR (AIRSAR) in HH polarization with
$L=1$ in C-band over Lake Superior, near Sault Saint Marie, MI.

We choose simple and analogous network architectures to be trained in both
Algorithms 1 and 2. The MLP in Algorithm 1 consists of two hidden layers, one
with eight and the other with four units, followed by $\tanh()$ activation
functions. The FCNN in Algorithm 2 is composed of two $1\times 1$
convolutional layers with eight and four filters, respectively, followed again
by $\tanh()$ activations. These activations were chosen as they performed
better than the traditional ReLU activation. For experimental simplicity, no
batch normalization, dropout or data preprocessing were used during training.

(a) Original Image

(b) $k=3$ (0.097 s)

(c) $k=11$ (0.637 s)

(d) $k=20$ (2.206 s)

(e) $k=45$ (11.887 s)

Figure 2: Qualitative results on a real SAR image. (a) The original intensity image. (b)-(e) The estimated roughness maps for different kernel sizes along with their respective elapsed computation time. The last color bar is shared by all maps. The images are better visualized in the digital version of this paper. TABLE I: Failure rates (%) on synthetic data. $L$ | MLE | LCUM | Neural Network
---|---|---|---
2 Log-moments | 4 Log-moments
$1$ | 71.84 | 37.44 | 0.14 | 3.24
$3$ | 59.84 | 25.55 | 1.46 | 6.22
$8$ | 47.99 | 18.59 | 6.04 | 6.96

Our networks were trained following the predicament in the TrainNeuralNet()
functions in both Algorithms 1 and 2 using Adaptive Moment Estimation (ADAM)
optimizer [17] with learning rate of 0.001 and batch size of 32. For Algorithm
1, we set $\mathcal{S}=\\{100,1000,10000\\}$ and tested $N_{m}=2$ and $4$, and
for Algorithm 2, we have $(w,h)=(10,10)$ and $\mathcal{K}=\\{2,5,\ldots,11\\}$
and only used $N_{m}=2$. Both algorithms use
$\mathcal{A}=\\{-15,-13.5,\cdots,-1.5\\}$, $R=1000$ datapoints and $T=300$
epochs, enough for their training convergence.

We compare our algorithms to LCUM and MLE results. To optimize both methods,
we use Scipy’s fsolve, a Python wrapper for an implementation of Powell’s dog
leg method [18]. We start each optimization with $\alpha=-1.0001$.

In our synthetic experiments, the quantitative assessment of our methods’
estimation performance compared the estimated roughness with their ground
truth values via Mean Square Error (MSE), as depicted in 2. For each
experiment, we also compared the failure rates for each algorithm, where we
consider an $\alpha$ estimate to fail if (1) it is not in the interval
$[-1.5,-15]$ [5], or (2) if its optimization procedure, in the case of MLE and
LCUM, did not converge.

We used a Tesla T4(R) GPU with 16Gb of RAM for our neural network training.
The other methods were run on an Intel(R) Xeon(R) CPU at 2.30GHz with 26Gb of
RAM. All methods were implemented in Python 3 and Pytorch 333The code used to
generate the results from this paper can be found at
https://github.com/jeovafarias/SAR-Roughness-Estimation-Neural-Nets.git..
Training each network took around 30 seconds each.

###  IV-B Results on Synthetic Data

We start by evaluating our estimation methods on synthetic samples of size
$s\in\\{9,25,49,121,1000\\}$, generated using roughness
$\alpha\in\\{-1.5,-7,-15\\}$ and number of looks $L\in\\{1,3,8\\}$. For each
setting, we performed a Monte Carlo experiment with 1000 $G_{I}^{0}$ samples.
Figure 1 compares MLE, LCUM and the networks trained using Algorithm 1 with 2
and 4 log-moments when estimating $\alpha$ for each sample. We only considered
the samples whose estimation did not fail to compute the MSE values. Note that
we run our networks on sample sizes that they were not trained on and, despite
that, both of them outperform its counterparts in most scenarios. The one
trained on fewer log-moments generally performed better, potentially
corroborating the premise in LCUM that also only uses two log-moments in its
algorithm. This also means that our neural network-based approach more
competently utilizes the information contained in those log-moments for
estimation than LCUM, especially for lower $\alpha$ values.

Table I shows the failure rates of each method for various $L$ values using
the same data used in Figure 1. Our networks also overperform both MLE and
LCUM on this domain for all depicted scenarios, especially for lower $L$,
where the estimation is typically more challenging [2]. Counterintuitively,
more failures are detected in our methods as $L$ increases. Further
experimentation is necessary to better understand this phenomenon, and it is
left to future work. We also found that estimation on our networks is
generally 3 times faster than the other methods.

###  IV-C Results on Real Data

Figure 2 shows the roughness maps estimated by networks trained according to
Algorithm 2 and using kernel sizes $k=\\{3,11,20,45\\}$ during inference. Note
that the $\alpha$ values in all maps follow their expected understanding:
higher values are related to urban areas; moderate values suggest forest zones
and low values correspond to lake regions [5]. This demonstrates our
methodology’s ability to correctly estimate the desired pixel roughness,
despite the networks having only been trained on purely synthetic data.
Furthermore, we see the blurring effect one expects from applying larger
kernels in their convolutions. More interestingly, however, is the observation
that the network is still able to predict an expected map using kernel sizes
it was not trained one (such as $k=\\{20,45\\}$). Finally and more
importantly, Figure 2 also provides the timings to the computation of each
map. Here, we note that we are able to estimate the roughness of each
individual pixel of an image as large as $1500\times 1500$ pixel in less than
a tenth of a second for $k=3$. Processing an image at this rate is usually
considered real time [19], meaning that this estimation can be performed as
the images are being acquired. As $k$ increases, the map computation becomes
slower, but it is still quickly accomplished for low $k$.

##  V Conclusion

We proposed a neural-network based algorithm for roughness estimation in
$G_{I}^{0}$ modeled SAR data, a task that is increasingly crucial in SAR image
understanding. Practically, it consists of using log-moments from the samples
as the input of a network that, when trained, outputs the underlying roughness
value. This same network can be easily adapted to process SAR images, from
moment computation to parameter estimation, leading to quick estimation of
pixel-level roughness. We empirically demonstrate that such networks trained
on purely synthetic data are able to outperform traditional estimation methods
and return reliable roughness maps from high resolution SAR images in real
time. Overall, this result shows that one can use cheap synthetic data to
train performant networks for SAR imaging tasks, where labeled data
acquisition is expensive. SAR data is in fact ideal for this approach, since
there are good statistical models for them that allow the generation of
synthetic training data, and they are structured such that important
operations on them can be implemented on GPUs. Future work will consist of
applying similar techniques to other data domains, where these aspects can
also be found.

## References

  * [1] A. C. Frery, H.-J. Muller, C. d. C. F. Yanasse, and S. J. S. Sant’Anna, “A model for extremely heterogeneous clutter,” _IEEE transactions on geoscience and remote sensing_ , vol. 35, no. 3, pp. 648–659, 1997.
  * [2] A. C. Frery, J. Wu, and L. Gomez, _SAR Image Analysis-A Computational Statistics Approach: With R Code, Data, and Applications_. John Wiley & Sons, 2022.
  * [3] J. F. R. Neto, A. M. Braga, R. C. Marques, and F. N. de Medeiros, “Level-set formulation based on an infinite series of sample moments for sar image segmentation,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 17, no. 5, pp. 908–911, 2019.
  * [4] J. Cassetti and A. C. Frery, “An improved minimum-distance texture estimator for speckled data under the g 0 model,” _Journal of Mathematical Imaging and Vision_ , vol. 64, no. 6, pp. 609–624, 2022.
  * [5] J. Gambini, J. Cassetti, M. M. Lucini, and A. C. Frery, “Parameter estimation in sar imagery using stochastic distances and asymmetric kernels,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 8, no. 1, pp. 365–375, 2015.
  * [6] S. Khan and R. Guida, “Application of mellin-kind statistics to polarimetric g distribution for sar data,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 52, no. 6, pp. 3513–3528, 2013.
  * [7] F. Á. Rodrigues, J. R. Neto, R. P. Marques, F. S. de Medeiros, and J. S. Nobre, “Sar image segmentation using the roughness information,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 13, no. 2, pp. 132–136, 2016.
  * [8] J. M. Nicolas, “Introduction aux statistiques de deuxième espèce: Applications des logs-moments et des logs-cumulants à l’analyse des lois d’images radar,” _TS. Traitement du signal_ , vol. 19, no. 3, pp. 139–167, 2002.
  * [9] X. X. Zhu, S. Montazeri, M. Ali, Y. Hua, Y. Wang, L. Mou, Y. Shi, F. Xu, and R. Bamler, “Deep learning meets sar: Concepts, models, pitfalls, and perspectives,” _IEEE Geoscience and Remote Sensing Magazine_ , vol. 9, no. 4, pp. 143–172, 2021.
  * [10] S. Minaee, Y. Boykov, F. Porikli, A. Plaza, N. Kehtarnavaz, and D. Terzopoulos, “Image segmentation using deep learning: A survey,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 44, no. 7, pp. 3523–3542, 2021.
  * [11] R. H. Nobre, F. A. Rodrigues, R. C. Marques, J. S. Nobre, J. F. Neto, and F. N. Medeiros, “Sar image segmentation with rényi’s entropy,” _IEEE Signal Processing Letters_ , vol. 23, no. 11, pp. 1551–1555, 2016.
  * [12] A. D. Nascimento, R. J. Cintra, and A. C. Frery, “Hypothesis testing in speckled data with stochastic distances,” _IEEE Transactions on geoscience and remote sensing_ , vol. 48, no. 1, pp. 373–385, 2009.
  * [13] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” _SIAM review_ , vol. 51, no. 3, pp. 455–500, 2009.
  * [14] M. Lin, Q. Chen, and S. Yan, “Network in network,” _Preprint arXiv:1312.4400_ , 2013.
  * [15] M. Zaheer, S. Kottur, S. Ravanbakhsh, B. Poczos, R. R. Salakhutdinov, and A. J. Smola, “Deep sets,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [16] F. Bieder, R. Sandkühler, and P. C. Cattin, “Comparison of methods generalizing max-and average-pooling,” _Preprint arXiv:2103.01746_ , 2021\.
  * [17] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [18] J. J. Moré, B. S. Garbow, and K. E. Hillstrom, “User guide for minpack-1,” CM-P00068642, Tech. Rep., 1980.
  * [19] G. Bradski, “Real-time computer vision with opencv,” _Communications of the ACM_ , vol. 55, no. 6, pp. 61–70, 2012.


# Phononic real Chern insulator with protected corner modes in graphynes

Jiaojiao Zhu  Research Laboratory for Quantum Materials, Singapore University
of Technology and Design, Singapore 487372, Singapore  Weikang Wu
<EMAIL_ADDRESS>Division of Physics and Applied Physics, School of
Physical and Mathematical Sciences, Nanyang Technological University,
Singapore, 637371, Singapore  Research Laboratory for Quantum Materials,
Singapore University of Technology and Design, Singapore 487372, Singapore
Jianzhou Zhao<EMAIL_ADDRESS>Co-Innovation Center for New Energetic
Materials, Southwest University of Science and Technology, Mianyang 621010,
China  Research Laboratory for Quantum Materials, Singapore University of
Technology and Design, Singapore 487372, Singapore  Cong Chen  School of
Physics, and Key Laboratory of Micro-nano Measurement-Manipulation and
Physics, Beihang University, Beijing 100191, China  Qianqian Wang  Research
Laboratory for Quantum Materials, Singapore University of Technology and
Design, Singapore 487372, Singapore  Xian-Lei Sheng  School of Physics, and
Key Laboratory of Micro-nano Measurement-Manipulation and Physics, Beihang
University, Beijing 100191, China  Lifa Zhang  NNU-SULI Thermal Energy
Research Center (NSTER) & Center for Quantum Transport and Thermal Energy
Science (CQTES), School of Physics and Technology, Nanjing Normal University,
Nanjing 210023, China  Y. X. Zhao<EMAIL_ADDRESS>National Laboratory of
Solid State Microstructures and Department of Physics, Nanjing University,
Nanjing 210093, China  Collaborative Innovation Center of Advanced
Microstructures, Nanjing University, Nanjing 210093, China  Shengyuan A. Yang
Research Laboratory for Quantum Materials, Singapore University of Technology
and Design, Singapore 487372, Singapore

###### Abstract

Higher-order topological insulators have attracted great research interest
recently. Different from conventional topological insulators, higher-order
topological insulators do not necessarily require spin-orbit coupling, which
makes it possible to realize them in spinless systems. Here, we study phonons
in 2D graphyne family materials. By using first-principle calculations and
topology/symmetry analysis, we find that phonons in both graphdiyne and
$\gamma$-graphyne exhibit a second-order topology, which belongs to the
specific case known as real Chern insulator. We identify the nontrivial
phononic band gaps, which are characterized by nontrivial real Chern numbers
enabled by the spacetime inversion symmetry. The protected phonon corner modes
are verified by the calculation on a finite-size nanodisk. Our study extends
the scope of higher-order topology to phonons in real materials. The spatially
localized phonon modes could be useful for novel phononic applications.

##  I Introduction

The study of topological insulators has developed into a vast research field
Hasan and Kane (2010); Qi and Zhang (2011); Shen (2012); Bansil  _et al._
(2016). A conventional $d$-dimensional topological insulator features a
nontrivial bulk topological invariant defined for the valence states below the
bulk band gap, and it has protected gapless states at its $(d-1)$-dimensional
boundaries. Later on, it was realized that with certain spatial symmetries,
there could exist a class of higher-order topological insulators (HOTIs), in
which the gapless states appear not at $(d-1)$\- but at $(d-n)$-dimensional
boundaries with $n>1$ Zhang  _et al._ (2013); Benalcazar  _et al._ (2017a);
Langbehn  _et al._ (2017); Song  _et al._ (2017); Benalcazar  _et al._
(2017b); Schindler  _et al._ (2018a). For example, a two-dimensional (2D)
second-order topological insulator would have gapped bulk and edge spectra,
meanwhile, gapless excitations occur at its 0D corners. The concept of HOTI
has attracted great interest Schindler  _et al._ (2018b); Wang  _et al._
(2019); Yue  _et al._ (2019); Xu  _et al._ (2019); Sheng  _et al._ (2019);
Zhang  _et al._ (2020). Unlike the conventional topological insulator, a HOTI
does not necessarily require the spin-orbit coupling. In other words, HOTI can
be realized in spinless systems. Thus, its impact is not limited to electronic
systems, but also spreads into bosonic and even classical systems, such as
photonic/acoustic metamaterials Mittal  _et al._ (2019); El Hassan  _et al._
(2019); Xue  _et al._ (2019); Ni  _et al._ (2019); Qi  _et al._ (2020),
electric circuit arrays Imhof  _et al._ (2018); Ezawa (2018), and mechanical
networks Serra-Garcia  _et al._ (2018); Ma  _et al._ (2019).

Phonons are the fundamental bosonic excitations in a crystal material,
describing the collective vibrations of atomic lattice. Phonons play a crucial
role in the material’s thermal properties, interact closely with other
quasiparticles or collective excitations, such as electrons, photons, magnons,
etc., and underly remarkable effects such as superconductivity. Recently,
there is a surge of interest in transferring the topological physics to
phonons Zhang  _et al._ (2010); Liu  _et al._ (2017); Zhang  _et al._ (2018);
Miao  _et al._ (2018); Zhang  _et al._ (2019); Xia  _et al._ (2019); Long  _et
al._ (2020); Liu  _et al._ (2020); Wang  _et al._ (2020); Li  _et al._ (2020,
2021); Zhu  _et al._ (2021). However, so far, most studies are focusing on the
gapless phonon states, such as nodal points and nodal lines in the phonon
spectrum, which are analogues of topological semimetals in the electronic
context. Meanwhile, gapped topological phonon states receive less attention.
Particularly, to our knowledge, there is no proposal of a realistic material
which exhibits a HOTI state in phonons.

In this work, we fill this gap by showing that phonons in the 2D graphyne
family materials exhibit the second-order topology. This study is motivated by
our recent proposal of these materials as the first example of 2D electronic
HOTIs Sheng  _et al._ (2019); Chen  _et al._ (2021). We have shown that their
HOTI state belongs to the specific case known as real Chern insulator Zhao and
Lu (2017), which is characterized by a nontrivial real Chern number $\nu_{R}$
and enabled by the spacetime inversion symmetry $\mathcal{PT}$. Here, we show
that the discussion can be naturally extended to phonons. Phonons are
spinless, obey bosonic statistics, and respect the same crystal symmetries
including $\mathcal{PT}$. Hence, each phononic band gap in the spectrum can be
characterized by a real Chern number $\nu_{R}\in\mathbb{Z}_{2}$. When
$\nu_{R}=1$, it indicates that the corresponding phononic gap is topologically
nontrivial, and there will be protected phonon modes localized at corners for
a finite sample that preserves $\mathcal{PT}$. Taking graphdiyne and graphyne
from the family as examples and using first-principles calculations, we
compute the real Chern numbers for their phononic band gaps via the parity
eigenvalue approach. The phonon corner modes are explicitly verified by the
calculation on a disk geometry. Our approach can be readily applied to other
materials to reveal possible higher-order topology in phonons. In addition to
the fundamental interest, the results here also open a new route to realize
localized phonons, which might be useful for novel device applications.

##  II Calculation Method

Density functional theory (DFT) calculations were conducted by using the
Vienna ab initio simulation package (VASP) Kresse and Furthmüller (1996a, b).
The projector augmented wave (PAW) pseudopotentials were adopted in the
calculation Blöchl (1994); Kresse and Joubert (1999). Generalized gradient
approximation (GGA) in the form of Perdew-Burke-Ernzerhof (PBE) realization
was adopted for the exchange-correlation potential Perdew  _et al._ (1996).
The valence electrons treated in the calculations include C (2s22p2). The
kinetic energy cutoff was fixed to 520 eV. $\Gamma$-centered $10\times
10\times 1$ $k$ point mesh was adopted for the self-consistent calculations.
The energy and force convergence criteria were set to be 10-7 eV and 0.001
eV/Å, respectively. We used density functional perturbation theory (DFPT)
Gonze and Lee (1997) in combination with the Phonopy package Togo and Tanaka
(2015) to obtain the force constants and phonon spectra.A supercell of
$3\times 3\times 1$ is adopted for the calculation of force constants. For
computing the phonons for the nanodisk geometry, we first calculated the
second rank tensor of force constants in Cartesian coordinates from DFPT, from
which we extracted the parameters for constructing the phononic tight-binding
model. The spatial distribution of phonon modes was calculated by using
Pybinding package Mol .

##  III Structure and phonon spectrum

Graphynes represent a family of 2D carbon allotropes composed of $sp$ and
$sp^{2}$ hybridized carbon atoms Baughman  _et al._ (1987); Diederich (1994);
Huang  _et al._ (2018); Gao  _et al._ (2019). They can be derived from the
graphene structure by inserting acetylenic linkages in different manners.
Their structural models were first proposed by Baughman _et al._ in 1987
Baughman  _et al._ (1987). Some members of this family have been successfully
realized in experiment. The most prominent example is graphdiyne, which was
first synthesized in 2010 by Li _et al._ via a cross-coupling reaction method
Huang  _et al._ (2018). Meanwhile, as the representative of the family,
$\gamma$-graphyne (often simply referred to as graphyne) has been realized in
experiment in the form of small fragments Diederich (1994); Bunz  _et al._
(1999). Hence, in this work, we will focus on these two examples. We will
first discuss the results of graphdiyne in the following sections. The results
for $\gamma$-graphyne are similar and will be briefly summarized in Sec. V.

The structure of graphdiyne is illustrated in Fig. 1(a). The lattice is
completely flat with single-atom thickness. It can be viewed as formed by
inserting one diacetylenic linkage between two neighboring benzene rings in
the graphene structure while maintaining its $p6m$ symmetry. Graphdiyne
exhibits a high $\pi$ conjugation, which helps to stabilize the 2D planar
structure and lower the energy. In fact, it was shown that graphdiyne is the
most stable among 2D carbon allotropes that contain diacetylenic linkages
Haley  _et al._ (1997). Each unit cell of the structure contains 18 carbon
atoms. The optimized lattice constant from our first-principles calculations
is 9.460 Å, which is consistent with previous results.

Figure 1: (a) Crystal structure of graphdiyne. The shaded region indicates the
unit cell. (b) Brillouin zone with high symmetry points labeled. (c)
Calculated phonon spectrum of graphdiyne. The arrows indicate three of the
five global phononic band gaps. Other two gaps are small hence are not
indicated here.

The calculated phonon spectrum of graphdiyne is plotted in Fig. 1(c). One
observes that the spectrum spans a range up to about 70 THz, which is
comparable to that of graphene ($\sim 50$ THz), implying its relatively strong
bonding. There are totally 54 phonon branches (3 acoustic branches plus 51
optical branches), corresponding to the 18 atoms in a unit cell. Among the
three acoustic branches, two have a linear dispersion approaching the $\Gamma$
point, whereas the remaining out-of-plane (ZA) acoustic branch has a quadratic
dispersion, which is a characteristic feature of 2D materials Zabel (2001);
Zhu  _et al._ (2014). The sound speed for the longitudinal acoustic phonons is
about 18.2 km/s, which is much larger than that of the MoS2 ($\sim$ 6.5 km/s)
Zhu  _et al._ (2015), and close to graphene ($\sim$ 21.2 km/s) Kaasbjerg  _et
al._ (2012).

Notably, we can observe several phononic band gaps in the spectrum, which are
labeled from Gap 1 to Gap 5 in Fig. 1(c) with increasing energy. Each gap
separates the whole phonon spectrum into two groups: the phonon bands below
the gap and the bands above, analogous to the valence and conduction bands in
the electronic context. Thus, for each gap, we can define a real Chern number
based on the bands below the gap (i.e., the “valence” bands) to characterize
the topology of the gap.

Before proceeding, we should comment on two important differences between
phonons and electrons. First, electrons are fermions, whereas phonons are
bosons. They follow distinct statistics. It follows that for electronic
systems, the only meaningful band gap to analyze is the one that near the
Fermi energy; but for phonons, all the band gaps are meaningful to discuss.
Second, the phonon energy scale (few tens of THz) is much smaller than the
electron’s. This makes the whole phonon spectrum detectable in experiment,
e.g., by optical or neutron scattering. For example, it is demonstrated in
recent experiment with inelastic x-ray scattering technique, which can achieve
the meV-resolution Baron  _et al._ (2001); Mohr  _et al._ (2007); Baron
(2015). This means that all the phononic band gaps are also experimentally
accessible.

The graphdiyne system preserves the $\mathcal{PT}$ symmetry. $\mathcal{PT}$
operates locally at every $k$ point of the Brillouin zone (BZ). For a spinless
system like phonons, we have $(\mathcal{PT})^{2}=1$. Up to a possible unitary
transformation, it can always be represented as $\mathcal{PT}=\mathcal{K}$,
where $\mathcal{K}$ is the complex conjugation. This means that with
$\mathcal{PT}$ symmetry, all the phonon eigenmodes $u_{n}(\bm{k})$ can be made
real, i.e., $u_{n}(\bm{k})=u^{*}_{n}(\bm{k})$. In Ref. Zhao and Lu (2017), it
has been shown that the real bands in 2D are characterized by a
$\mathbb{Z}_{2}$ invariant, the real Chern number $\nu_{R}$, which is
analogous to the conventional Chern number for complex bands. When there are
only two bands below the gap, $\nu_{R}$ has a integral expression in terms of
the real Berry curvature of these bands, as shown in Ref. Zhao and Lu (2017).
However, when there are more than two bands, an integral expression for
$\nu_{R}$ is not available. Nevertheless, for systems that also preserve the
inversion symmetry $\mathcal{P}$ (which is the case here), $\nu_{R}$ can be
evaluated by using the parity analysis approach Ahn  _et al._ (2018).

In this approach, one needs to analyze the parity eigenvalues of phonon modes
at the four inversion-invariant momentum points $\Gamma_{i}$ $(i=1,\cdots,4)$
of the BZ, which include $\Gamma$ and the three $M$ points [see Fig. 1(b)].
For each phononic band gap, at point $\Gamma_{i}$, one counts the number of
phonon bands below the gap that have the negative parity eigenvalue. Denote
this number by $n^{\Gamma_{i}}_{-}$. Then, $\nu_{R}$ can be obtained as Ahn
_et al._ (2018)

$(-1)^{\nu_{R}}=\prod_{i=1}^{4}(-1)^{\left\lfloor(n_{-}^{\Gamma_{i}}/2)\right\rfloor}.$
(1)

where $\left\lfloor\cdots\right\rfloor$ is the floor function.

We apply this approach to analyze the 5 phononic band gaps in graphdiyne. The
obtained results are listed in Table 1. One can see that Gap 1 and Gap 2 have
a nontrivial real Chern number $\nu_{R}=1$, so these gaps correspond to a
phononic real Chern insulator state.

Table 1:  Topological characterization of the phononic band gaps in graphdiyne. Here, $N$ is the number of bands below the gap. $n_{-}^{\Gamma}$ and $n_{-}^{M}$ are the number of bands (out of the $N$ bands) with negative parity eigenvalues at $\Gamma$ and $M$, respectively. $\nu_{R}$ gives the corresponding real Chern number. Gap No. |  Frequency range (THz) |  Bands below the gap ($N$) |  $n_{-}^{\Gamma}$ |  $n_{-}^{M}$ |  $\nu_{R}$
---|---|---|---|---|---
1 |  $15.0-15.2$ |  28 |  16 |  14 |  1
2 |  $15.6-17.6$ |  30 |  16 |  14 |  1
3 |  $23.2-24.4$ |  36 |  18 |  18 |  0
4 |  $39.8-40.1$ |  43 |  22 |  22 |  0
5 |  $44.2-64.4$ |  48 |  24 |  24 |  0

Let’s consider Gap 2, which is the gap between 15.6 THz and 17.6 THz. Below
this gap, there are 30 phonon branches. At $\Gamma$, there are 14 branches
with positive parity eigenvalues and 16 branches with negative parity
eigenvalues, i.e., $n^{\Gamma}_{+}=14$ and $n^{\Gamma}_{-}=16$. In comparison,
at the three $M$ points, we find $n^{M}_{+}=16$ and $n^{M}_{-}=14$. This
mismatch signals a double band inversion between $\Gamma$ and $M$, which is
captured by the nontrivial real Chern number. Meanwhile, Gap 1 has a smaller
size from 15.00 THz to 15.25 THz. The analysis is similar, so we will not
elaborate here.

##  IV Corner phonon mode

It has been shown that a real Chern insulator has a second-order topology Chen
_et al._ (2021), i.e., there must exist 0D corner modes for any sample that
preserves the $\mathcal{PT}$ symmetry. Since we have identified phononic band
gaps in graphdiyne that feature nontrivial real Chern numbers, it means that
in a graphdiyne nanodisk, we should have localized corner phonon modes in
these gaps.

To confirm this prediction, we calculate the phonon spectrum for the
graphdiyne nanodisk geometry. The calculation is based on an _ab initio_ tight
binding model extracted from the first-principles calculations. In Fig. 2(a),
we plot the phonon spectrum close to Gap 2. One observes that at 17.02 THz
inside Gap 2, there are 6 degenerate phonon modes. We have checked the spatial
distribution of these modes and found that they are localized at the six
corners of the disk, as illustrated in Fig. 2(c). As a comparison, in Fig.
2(d), we also plot the spatial distribution of the nearby modes at 15.59 THz,
and one can easily see that the conventional modes is not localized but spread
throughout the sample. In Fig. 2(b), we plot the local density of states at
one corner (taking the mode weight within one unit cell of the corner). One
can clearly observe the sharp peak due to the localized corner modes.

In the above calculation, we have taken a hexagonal shaped disk geometry. This
is the natural sample geometry that is obtained from the bottom-up synthesis
approach Matsuoka  _et al._ (2017). Hence, our results here can be directly
tested in experiment. Moreover, it must be emphasized that the existence of
the corner modes is general and not depend on the specific sample geometry. As
demonstrated in Ref. Chen  _et al._ (2021), such modes must exist when the
sample geometry preserves $\mathcal{PT}$. In Fig. 3, we show the result for a
diamond shaped sample. Again, localized corner modes in Gap 2 can be observed.

Figure 2: (a) Phonon spectrum for the graphdiyne nanodisk in (c). Here, we
focus on the part of frequency range close to the bulk Gap 2. (b) Corner
projection spectrum, which is the phonon local density of states at one corner
of the disk. The arrow indicates the peak due to the corner modes. (c) Spatial
distribution of the phonon modes at 17.02 THz, i.e. the red colored modes
inside Gap 2 in (a), which are localized at the corners of the disk. For
comparison, (d) shows the distribution of phonon modes at a nearby frequencies
(15.59 THz), which is extended in the bulk of the disk. Figure 3: (a) Phonon
spectrum for a diamond shaped graphdiyne nanodisk in (b). Again, we focus on
the frequency range close to Gap 2. (b) shows the distribution of the two red
colored modes at 16.54 THz, which are localized at corners.

##  V Discussion and Conclusion

We have demonstrated the phononic HOTI state in graphdiyne. More specifically,
it is the real Chern insulator state. In the bulk, there exist phononic band
gaps that are characterized by nontrivial real Chern numbers. At the boundary,
the nontrivial invariant dictates the existence of corner phonon modes.

Figure 4: (a) Crystal structure of $\gamma$-graphyne. The unit cell is marked
by the shaded region. (b) Phonon spectrum of $\gamma$-graphyne. The arrows
indicate three global gaps.

The results for $\gamma$-graphyne is similar. The structure of
$\gamma$-graphyne is shown in Fig. 4(a), where the benzene rings are connected
by acetylenic rather than diacetylenic bond. There are 5 phononic band gaps in
its phonon spectrum, of which Gap 1 and Gap 5 are identified to be nontrivial
and feature $\nu_{R}=1$ (see Supplemental Material sup  for more details).
Here, let’s consider Gap 1. In Fig. 5(c), we show the corner phonon modes
correspond to the gap around 20.42 THz for a graphyne nanodisk. Clearly, our
approach can be applied to study other members of the graphyne family as well,
as all typical members of the family possess both the $\mathcal{P}$ and
$\mathcal{T}$ symmetries. Particularly, we expect that the phononic real Chern
insulator states can appear in the graphyne-$n$ structures Gao  _et al._
(2019); Huang  _et al._ (2018), which are generalizations of $\gamma$-graphyne
and graphdiyne with extended linkage containing $n$ $-$C$\equiv$C$-$ bonds.

Figure 5: (a) Phonon spectrum for the $\gamma$-graphyne nanodisk shown in (c).
Here, we focus on the frequency range close to Gap 1 in Fig. 4(b). (b)
Projection spectra at one corner of the disk. The arrow indicate the peak due
to the corner modes. (c) Distribution for the red-colored phonon modes at
20.42 THz in (a). For comparison, (d) shows the distribution for bulk modes at
20.05 THz.

Experimentally, graphdiyne has been successfully synthesized in experiment Li
_et al._ (2010); Huang  _et al._ (2018); Gao  _et al._ (2019) and samples with
the hexagonal disk geometry can be naturally realized via the bottom-up growth
method Matsuoka  _et al._ (2017). To probe the corner phonon modes, one can
excite them with an infrared light in resonance with their frequency. As the
corner modes are spatially localized at the corners, they can manifest in
spatially resolved fluorescence pattern, or they can be detected via local
probes such as scanning tunneling microscopy, by comparing the measurements at
the center and at the corner of the disk.

The localized phonon mode in a spectral gap could be useful for studying many
interesting effects. For example, it may offer a new platform to realize the
phonon laser effect, which was previously demonstrated in trapped
ions/nanosphere and micro-cavities Vahala  _et al._ (2009); Huang and Jing
(2019); Pettit  _et al._ (2019); Ohtani  _et al._ (2019). The effect could be
useful for constructing novel devices for ultra-sensitive force sensing and
high-precision imaging applications.

In conclusion, we have reported the phononic HOTI state, more specifically,
the real Chern insulator state, in the graphyne family materials. Taking
graphdiyne and $\gamma$-graphyne as examples, we show that their phonon
spectra have nontrivial phononic band gaps characterized by the real Chern
numbers enabled by the $\mathcal{PT}$ symmetry. We evaluate the real Chern
numbers by the parity analysis approach. The topological corner phonon modes
associated with the nontrivial gaps are explicitly demonstrated via
calculation on a disk geometry. Our study extends the scope of higher-order
topology to phonons in real material systems. The spatially localized corner
phonon modes could be useful for novel phononic applications and phonon-photon
coupling at nanoscale.

###### Acknowledgements.

The authors thank D. L. Deng for helpful discussions. This work is supported
by Singapore Ministry of Education AcRF Tier 2 (MOE2019-T2-1-001). We
acknowledge computational support from the Texas Advanced Computing Center and
the National Supercomputing Centre Singapore.

## References

  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Shen (2012) S.-Q. Shen,  _Topological insulators_ , Vol. 174 (Springer, 2012).
  * Bansil  _et al._ (2016) A. Bansil, H. Lin, and T. Das, Rev. Mod. Phys. 88, 021004 (2016).
  * Zhang  _et al._ (2013) F. Zhang, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 110, 046404 (2013).
  * Benalcazar  _et al._ (2017a) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science 357, 61 (2017a).
  * Langbehn  _et al._ (2017) J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Phys. Rev. Lett. 119, 246401 (2017).
  * Song  _et al._ (2017) Z. Song, Z. Fang, and C. Fang, Phys. Rev. Lett. 119, 246402 (2017).
  * Benalcazar  _et al._ (2017b) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Phys. Rev. B 96, 245115 (2017b).
  * Schindler  _et al._ (2018a) F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. Parkin, B. A. Bernevig, and T. Neupert, Sci. Adv. 4, eaat0346 (2018a).
  * Schindler  _et al._ (2018b) F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, _et al._ , Nat. Phys. 14, 918 (2018b).
  * Wang  _et al._ (2019) Z. Wang, B. J. Wieder, J. Li, B. Yan, and B. A. Bernevig, Phys. Rev. Lett. 123, 186401 (2019).
  * Yue  _et al._ (2019) C. Yue, Y. Xu, Z. Song, H. Weng, Y.-M. Lu, C. Fang, and X. Dai, Nat. Phys. 15, 577 (2019).
  * Xu  _et al._ (2019) Y. Xu, Z. Song, Z. Wang, H. Weng, and X. Dai, Phys. Rev. Lett. 122, 256402 (2019).
  * Sheng  _et al._ (2019) X.-L. Sheng, C. Chen, H. Liu, Z. Chen, Z.-M. Yu, Y. Zhao, and S. A. Yang, Phys. Rev. Lett. 123, 256402 (2019).
  * Zhang  _et al._ (2020) R.-X. Zhang, F. Wu, and S. D. Sarma, Phys. Rev. Lett. 124, 136407 (2020).
  * Mittal  _et al._ (2019) S. Mittal, V. V. Orre, G. Zhu, M. A. Gorlach, A. Poddubny, and M. Hafezi, Nat. Photon. 13, 692 (2019).
  * El Hassan  _et al._ (2019) A. El Hassan, F. K. Kunst, A. Moritz, G. Andler, E. J. Bergholtz, and M. Bourennane, Nat. Photon. 13, 697 (2019).
  * Xue  _et al._ (2019) H. Xue, Y. Yang, F. Gao, Y. Chong, and B. Zhang, Nat. Mater. 18, 108 (2019).
  * Ni  _et al._ (2019) X. Ni, M. Weiner, A. Alu, and A. B. Khanikaev, Nat. Mater. 18, 113 (2019).
  * Qi  _et al._ (2020) Y. Qi, C. Qiu, M. Xiao, H. He, M. Ke, and Z. Liu, Phys. Rev. Lett. 124, 206601 (2020).
  * Imhof  _et al._ (2018) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, _et al._ , Nat. Phys. 14, 925 (2018).
  * Ezawa (2018) M. Ezawa, Phys. Rev. B 98, 201402 (2018).
  * Serra-Garcia  _et al._ (2018) M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, Nature 555, 342 (2018).
  * Ma  _et al._ (2019) G. Ma, M. Xiao, and C. T. Chan, Nat. Rev. Phys. 1, 281 (2019).
  * Zhang  _et al._ (2010) L. Zhang, J. Ren, J.-S. Wang, and B. Li, Phys. Rev. Lett. 105, 225901 (2010).
  * Liu  _et al._ (2017) Y. Liu, C.-S. Lian, Y. Li, Y. Xu, and W. Duan, Phys. Rev. Lett. 119, 255901 (2017).
  * Zhang  _et al._ (2018) T. Zhang, Z. Song, A. Alexandradinata, H. Weng, C. Fang, L. Lu, and Z. Fang, Phys. Rev. Lett. 120, 016401 (2018).
  * Miao  _et al._ (2018) H. Miao, T. Zhang, L. Wang, D. Meyers, A. Said, Y. Wang, Y. Shi, H. Weng, Z. Fang, and M. Dean, Phys. Rev. Lett. 121, 035302 (2018).
  * Zhang  _et al._ (2019) T. Zhang, H. Miao, Q. Wang, J. Lin, Y. Cao, G. Fabbris, A. Said, X. Liu, H. Lei, Z. Fang, _et al._ , Phys. Rev. Lett. 123, 245302 (2019).
  * Xia  _et al._ (2019) B. Xia, R. Wang, Z. Chen, Y. Zhao, and H. Xu, Phys. Rev. Lett. 123, 065501 (2019).
  * Long  _et al._ (2020) Y. Long, J. Ren, and H. Chen, Phys. Rev. Lett. 124, 185501 (2020).
  * Liu  _et al._ (2020) Y. Liu, X. Chen, and Y. Xu, Adv. Funct. Mater. 30, 1904784 (2020).
  * Wang  _et al._ (2020) R. Wang, B. Xia, Z. Chen, B. Zheng, Y. Zhao, and H. Xu, Phys. Rev. Lett. 124, 105303 (2020).
  * Li  _et al._ (2020) J. Li, L. Wang, J. Liu, R. Li, Z. Zhang, and X.-Q. Chen, Phys. Rev. B 101, 081403 (2020).
  * Li  _et al._ (2021) J. Li, J. Liu, S. A. Baronett, M. Liu, L. Wang, R. Li, Y. Chen, D. Li, Q. Zhu, and X.-Q. Chen, Nat. Commun. 12, 1 (2021).
  * Zhu  _et al._ (2021) J. Zhu, W. Wu, J. Zhao, H. Chen, L. Zhang, and S. A. Yang, arXiv preprint arXiv:2104.14816 (2021).
  * Chen  _et al._ (2021) C. Chen, W. Wu, Z.-M. Yu, Z. Chen, Y. Zhao, X.-L. Sheng, and S. A. Yang, Phys. Rev. B 104, 085205 (2021).
  * Zhao and Lu (2017) Y. Zhao and Y. Lu, Phys. Rev. Lett. 118, 056401 (2017).
  * Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996a).
  * Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996b).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Gonze and Lee (1997) X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997).
  * Togo and Tanaka (2015) A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * (47) Dean Moldovan, Miša Anđelković, and Francois Peeters. (2020). pybinding v0.9.5: a Python package for tight-binding calculations (v0.9.5). Zenodo. https://doi.org/10.5281/zenodo.4010216.
  * Baughman  _et al._ (1987) R. Baughman, H. Eckhardt, and M. Kertesz, J. Chem. Phys. 87, 6687 (1987).
  * Diederich (1994) F. Diederich, Nature 369, 199 (1994).
  * Huang  _et al._ (2018) C. Huang, Y. Li, N. Wang, Y. Xue, Z. Zuo, H. Liu, and Y. Li, Chem. Rev. 118, 7744 (2018).
  * Gao  _et al._ (2019) X. Gao, H. Liu, D. Wang, and J. Zhang, Chem. Soc. Rev. 48, 908 (2019).
  * Bunz  _et al._ (1999) U. H. Bunz, Y. Rubin, and Y. Tobe, Chem. Soc. Rev. 28, 107 (1999).
  * Haley  _et al._ (1997) M. M. Haley, S. C. Brand, and J. J. Pak, Angew. Chem. Int. Ed. Engl. 36, 836 (1997).
  * Zabel (2001) H. Zabel, J. Phys. Condens. Matter 13, 7679 (2001).
  * Zhu  _et al._ (2014) L. Zhu, G. Zhang, and B. Li, Phys. Rev. B 90, 214302 (2014).
  * Zhu  _et al._ (2015) L. Zhu, T. Zhang, Z. Sun, J. Li, G. Chen, and S. A. Yang, Nanotechnology 26, 465707 (2015).
  * Kaasbjerg  _et al._ (2012) K. Kaasbjerg, K. S. Thygesen, and K. W. Jacobsen, Phys. Rev. B 85, 165440 (2012).
  * Baron  _et al._ (2001) A. Baron, Y. Tanaka, D. Miwa, D. Ishikawa, T. Mochizuki, K. Takeshita, S. Goto, T. Matsushita, H. Kimura, F. Yamamoto, _et al._ , Nucl. Instrum. Meth. A 467, 627 (2001).
  * Mohr  _et al._ (2007) M. Mohr, J. Maultzsch, E. Dobardžić, S. Reich, I. Milošević, M. Damnjanović, A. Bosak, M. Krisch, and C. Thomsen, Phys. Rev. B 76, 035439 (2007).
  * Baron (2015) A. Q. Baron, arXiv preprint arXiv:1504.01098 (2015).
  * Ahn  _et al._ (2018) J. Ahn, D. Kim, Y. Kim, and B.-J. Yang, Phys. Rev. Lett. 121, 106403 (2018).
  * Matsuoka  _et al._ (2017) R. Matsuoka, R. Sakamoto, K. Hoshiko, S. Sasaki, H. Masunaga, K. Nagashio, and H. Nishihara, J. Am. Chem. Soc. 139, 3145 (2017).
  * (63) See Supplemental Material.
  * Li  _et al._ (2010) G. Li, Y. Li, H. Liu, Y. Guo, Y. Li, and D. Zhu, Chem. Commun. 46, 3256 (2010).
  * Vahala  _et al._ (2009) K. Vahala, M. Herrmann, S. Knünz, V. Batteiger, G. Saathoff, T. Hänsch, and T. Udem, Nat. Phys. 5, 682 (2009).
  * Huang and Jing (2019) R. Huang and H. Jing, Nat. Photon. 13, 372 (2019).
  * Pettit  _et al._ (2019) R. M. Pettit, W. Ge, P. Kumar, D. R. Luntz-Martin, J. T. Schultz, L. P. Neukirch, M. Bhattacharya, and A. N. Vamivakas, Nat. Photon. 13, 402 (2019).
  * Ohtani  _et al._ (2019) K. Ohtani, B. Meng, M. Franckié, L. Bosco, C. Ndebeka-Bandou, M. Beck, and J. Faist, Sci. Adv. 5, eaau1632 (2019).