# Effects of varying initial conditions of ship encountering
wave groups in computing extreme motion statistics

Xianliang Gong, Yulin Pan
(Department of Naval Architecture and Marine Engineering,
University of Michigan, USA)

## ABSTRACT

In computing ship motion statistics (e.g., exceeding probability) in an
irregular wave field, it is a common practice to represent the irregular waves
by a large number of wave groups and compute the motion statistics from the
distribution of these wave groups. While this procedure significantly reduces
the computational cost, the uncertainties introduced in this reduced-order
computation have not been quantified. In general, the representation of a
continuous wave field by separated groups loses information about wave phases,
frequency modulation and initial conditions of ship when encountering the wave
groups, among which the last one is arguably the most influential factor for
the ship motion statistics. In this paper, we test three methods to
incorporate the ship initial conditions into the computation of roll motion
exceeding probability, namely the methods of natural initial condition, pre-
computed initial condition and constant initial condition. For different input
wave spectra and ship motion dynamics (in terms of the parametric excitation
level modeled by a nonlinear roll equation), the performances of the three
methods are carefully benchmarked and the superiority of the natural initial
condition method is demonstrated. We finally show that the computation using
the natural initial condition method can be greatly accelerated through a
sequential sampling algorithm making use of the variational heteroscedastic
Gaussian process regression.

## INTRODUCTION

Extreme ship motions can be triggered from various mechanisms (e.g., large
waves, parametric rolls) for ships travelling in an irregular wave field.
Although these motions occur with a low probability, they may cause severe
damage to ships especially at high sea states. A reliable quantification of
the extreme motion statistics, e.g., probability of ship motion exceeding a
certain threshold, is of vital importance for assessment of the system
reliability and reduction of the failure probability in the ship design
process.

The computation of the motion exceeding probability, however, is a non-trivial
task. Considering the high dimensions of the wavefield (as input to the ship-
motion computation), the rareness of the extreme motions, and the high
computational cost in CFD simulations of ship motion in waves, the direct
Monte Carlo sampling in the whole input space (or a long simulation covering
all wave conditions) is computationally prohibitive. Many methods have been
developed to reduce the computational cost of the problem, with one critical
effort to reduce the dimension of the input irregular wave field through a
wave group representation (e.g. Anastopoulos and Spyrou, 2017; Cousins and
Sapsis, 2016). Under such a representation (specifically the second reference)
only two parameters of the group height $A$ and group length $L$ (with a given
probability distribution) are needed to describe the wave field. Although this
procedure allows the computation to be done with a much lower cost, certain
information has been lost in using separate wave groups to represent a
continuous wave field, namely the wave phases, the frequency modulation (of
waves in a group) and the initial condition of the ship when encountering a
wave group. The uncertainties associated with this lost information are
usually ignored (or insufficiently quantified) in most existing works adopting
the method of wave group parameterization.

Another critical development in accelerating the computation of motion
exceeding probability is to use sequential sampling (or sequential Bayesian
experimental design) method (e.g. Echard et al., 2011; Hu and Mahadevan, 2016;
Mohamad and Sapsis, 2018; Gong et al., 2021) to reduce the Monte-Carlo
sampling numbers required for the computation. The principle of these methods
is to establish a Gaussian process as a surrogate model through the existing
samples (and their responses), and optimize for the location of the next
sample in the $(A,L)$ group parameter space. In applying these methods to
naval platforms, the problem associated with the initial conditions (as well
as other lost information) was usually ignored or somehow avoided by
considering a static platform (Mohamad and Sapsis, 2018). In order to
incorporate the effect of varying initial conditions into the sequential
sampling framework, the authors of this paper have developed a new approach
(VHGP-Seq) in Gong and Pan (2022), which uses the variational heteroscedastic
Gaussian process as a surrogate model (thus quantifying the varying
uncertainty in the responses from different groups) and solves a new
optimization problem to locate the next-best sample. The VHGP-Seq method has
been validated successfully in one case of computing the exceeding probability
of roll motions, but this is not sufficient to guarantee the effectiveness of
the method for wide ranges of the input wave spectra and ship roll dynamics.

The purposes of the present paper are twofold. In the first part, we test the
performance of three existing methods to handle the initial conditions of ship
encountering wave groups in computing the motion exceeding probability.
Specifically, the first method (Anastopoulos and Spyrou, 2017; Gong et al.,
2021) uses constant (usually zero displacement and velocity) initial condition
of the ship motion for each wave group and therefore neglect the variation of
the response induced from the variation of initial conditions (hereafter named
constant initial condition). The second method (Anastopoulos and Spyrou, 2019;
Silva and Maki, 2021) generates a distribution of initial conditions prior to
the computation of statistics through simulation of ship motion in the same
wave field (hereafter named pre-computed initial condition). The third method
(Gong and Pan, 2022) simulates the ship motion from several groups ahead of
the objective group and thus incorporates the initial condition in a natural
way (hereafter named natural initial condition). The three methods are applied
to cases with different input wave spectra, motion thresholds and ship motion
dynamics (in terms of the level of parametric excitation simulated by a
nonlinear roll equation), and are compared to the true solutions generated as
the ship goes through a sufficiently long continuous wave field (or time
series). We demonstrate the superiority of the method of natural initial
condition in term of its agreement to the true solution in all tested cases.

In the second part of the paper, we test the performance of the VHGP-Seq
sampling method in reproducing the true solution of the exceeding probability
(or solution from the method of natural initial condition) for a number of
cases (as a representative subset of all cases used to test the methods of
initial conditions). For all these cases, we show that the exceeding
probability from the VHGP-Seq method converges to the true solution in
$O(100)$ samples, instead of more than $O(1000)$ samples required in the
random sampling approach.

## METHOD

Figure 1: (a) Wave field containing successive wave groups with surface
elevation $\eta(t)$ ( ) and envelop $\rho(t)$ ( ). (b) $\rho(t)$ ( ) fitted by
Gaussian wave groups $\rho_{i}(t)$ ( ). (c) The joint PDF of $L$ and $A$
computed from a large number of detected groups.

We start from a narrow-band wave field with a sufficiently long time series
$\eta(t)$, which contains successive wave groups
$\\{g_{1},g_{2},\cdots,g_{n}\\}$ (see illustration from figure 1(a), and
mathematical definitions shortly). Our objective is to compute the exceeding
probability (probability of motion larger than a threshold, say $\theta_{t}$)
when a ship goes through the wave field described by $\eta(t)$. In particular,
we define the ground truth of the group-based exceeding probability as

$P_{true}=\sum_{i=1}^{n}I(r_{i})/n,$ (1)

where $r_{i}$ is the maximum ship motion within a segment of $\eta(t)$
corresponding to the group $g_{i}$, simulated as the ship continuously goes
through $\eta(t)$, $n$ is the total number of wave groups in $\eta(t)$, and
$I$ is an indicator function:

$I(r_{i})=\left\\{\begin{aligned} &1&&if\;r_{i}>\theta_{t}\\\
&0&&if\;r_{i}\leq\theta_{t}\end{aligned}\right..$

The computation of $P_{true}$, as defined in (1), involves a simulation of the
ship response from $\eta(t)$, which needs to be extremely long to cover (many
times) all wave conditions associated with a given wave spectrum. When coupled
with sophisticated motion simulations such as CFD, the computational cost can
become prohibitively high.

In order to reduce the overall computational cost, one may seek to
parameterize the time series $\eta(t)$ by features of the wave groups and
their associated probability distribution. There are two approaches to achieve
this goal, the first by assuming a Markov-chain property of
$\\{g_{i}\\}_{i=1}^{n}$ and computing the group distribution from the
transitional probability following Kimura (1980) and Longuet-Higgins (1984),
and the second by constructing the groups directly from the envelope process
of $\eta(t)$ (Cousins and Sapsis, 2016). We follow the second approach in the
current paper, where we compute the envelope process $\rho(t)$ from $\eta(t)$
through the Hilbert transform (Shum and Melville, 1984), and then construct
Gaussian-like wave groups $\rho_{i}(t)$ which best fit $\rho(t)$ locally:

$\rho_{i}(t)\sim a_{i}\exp\frac{-(t-t_{i})^{2}}{2l_{i}^{2}},$ (2)

where $t_{i}$, $a_{i}$, and $l_{i}$ are respectively the temporal location,
amplitude and length of group $g_{i}$ (see figure 1(b)). This construction
relies on a group detection algorithm proposed by Cousins and Sapsis (2016)
and later improved by Gong et al. (2021). Based on the detected groups, we can
construct a two-dimensional description of the wave field with probability
distribution $p_{A,L}(a,l)$ (figure 1(c)), and then compute an estimation of
the group-based exceeding probability as

$P_{est}=\int I(r(a,l))p_{A,L}(a,l)\mathrm{d}a\mathrm{d}l,$ (3)

where $r(a,l)$ is the maximum ship response in a group with parameter $(a,l)$
(with choice of initial conditions that will be described later). While (3)
provides a reduced-order computation relative to (1), the uncertainties
associated with this simplification have not been fully quantified.
Specifically, in using separated $(a,l)$ groups to represent $\eta(t)$, the
uncertainties can come from certain lost information, including the wave
phases in the group, frequency modulation in the group and the initial
condition of the ship encountering a group. Among these factors, the ship
initial condition is arguably (as will also be confirmed in our results) the
most influential one for the motion exceeding probability. In the next three
sections, we review different existing approaches to handle the encountering
initial conditions, followed by a brief review of the VHGP-Seq sampling method
and a presentation of our test framework.

Figure 2: Flowchart of the three methods to incorporate the effects of initial
conditions.

### Approaches for initial conditions

We summarize three existing approaches to handle the initial conditions
associated with $r(a,l)$.

  1. 1.

constant initial condition (CIC): We compute $r(a,l)$ for all groups using a
constant initial condition (usually with zero displacement and velocity), and
compute the exceeding probability directly by (3).

  2. 2.

pre-computed initial condition (PIC): We generate a distribution of ship
motions (say, roll motion $\theta$ and velocity $\dot{\theta}$) from a
simulation covering part of $\eta(t)$, which is then used as an ensemble of
initial conditions to compute the response from an $(a,l)$ group. In this
case, the response becomes a random variable denoted by $R(a,l)$, and
accordingly we compute the exceeding probability by

$P_{est}=\int\mathbb{P}(R(a,l)>\theta_{t})p_{A,L}(a,l)\mathrm{d}a\mathrm{d}l.$
(4)

  3. 3.

natural initial condition (NIC): For a selected pair of parameters $(a,l)$, we
simulate the responses from a number of detected $(a,l)$ groups. For each of
these $(a,l)$ groups, we start the simulation from two or three groups ahead
of the objective group to naturally capture the initial condition. The
quantity $\mathbb{P}(R(a,l)>\theta_{t})$ can then be approximated by the
fraction of the responses that exceed $\theta_{t}$. Finally the exceeding
probability is computed from (4). This approach is in essence the same as a
Monte-Carlo sampling of random groups from $\eta(t)$ with applications of NIC
to simulate the response from each of them.

A sketch of these three approaches are further shown in figure 2. In
application of all these methods, a necessary procedure is to recover a time
series of surface elevation $\eta_{i}(t)$ from a given $(a,l)$ group
$g_{i}(t)$ (as the input to simulate the response). In this work, the time
series $\eta_{i}(t)$ is obtained by

$\eta_{i}(t)=\rho_{i}(t)\cos(\omega_{p}t),$ (5)

where $\omega_{p}$ is the peak angular frequency of the wave spectrum under
consideration, $\rho_{i}(t)$ is connected to $a$ and $l$ from (2). It is clear
from (5) that the information about wave phases (and frequency modulation) of
the original series $\eta(t)$ is lost in all three methods. Other than this
lost information, the NIC method seems to be the most sophisticated one to
incorporate the effect of initial conditions. Therefore, it is not
unreasonable to expect that the method of NIC may provide the closest result
to the true solution (1). In fact, A direct comparison between the result from
NIC and the true solution (1) will reveal whether the ship initial condition
is the dominant factor of uncertainty and whether the proposed method is
sufficient to capture its effect on the exceeding probability.

We also remark that direct applications of the PIC and NIC methods (i.e., with
random sampling in the $(A,L)$ parameter space) to compute the exceeding
probability are still relatively expensive. Taking the NIC method as an
example, while it is possible to explore the reduction of the computational
cost by simulating only a subset of wave groups out of those in (1), this
reduction may not be sufficient. It should be noted here that our first focus
of the paper is to understand the accuracy of the three methods in
incorporating the effect of initial conditions, irrespective of the
computational cost. The reduction of the computational cost will further be
explored in the context of VHGP-Seq method that is next discussed.

### VHGP-Seq sampling method

We next briefly review the VHGP-Seq sampling method, that can be used to
significantly reduce the number of samples required in the computation using
the NIC method. The principle of the VHGP-Seq is to first establish a
surrogate model of the stochastic response $R(a,l)$ using the variational
heteroscedastic Gaussian process (VHGP) from existing samples. In this
surrogate model, the random variable $R(a,l)$ (for given $a$ and $l$) is
approximated by a normal distribution, and two Gaussian processes are used to
model the mean and variance of $R(a,l)$. Then the next-best sample location is
chosen based on an acquisition optimization problem which maximizes the
contribution of the sample to the computation in (4).

Since the variational algorithm is used in VHGP, the surrogate model can be
established in an efficient manner, with a computational cost that is only
twice of a standard Gaussian process regression. We remark that the VHGP-Seq
method is so far the only available method for sequential sampling in problems
with a stochastic response function whose standard deviation varies with the
input. The NIC problem is clearly an example of this type of problems that are
generally encountered in many other contexts. For more details of VHGP-Seq, we
refer the readers to our previous paper (Gong and Pan, 2022).

### Test framework

Our purpose is to benchmark the performances of CIC, PIC and NIC in computing
the exceeding probability, as well as the performance of VHGP-Seq in
accelerating the computation. Since the true solution (1) is needed in the
comparisons, computation of the ship motion using CFD becomes impractical
because of the prohibitively high computational cost. In this work, we instead
employ a nonlinear roll equation (Umeda et al., 2004) to calculate the time
series of ship roll $\theta(t)$ from the wave input $\eta(t)$:

$\displaystyle\ddot{\theta}+\alpha_{1}\dot{\theta}+\alpha_{2}\dot{\theta}^{3}+$
$\displaystyle(\beta_{1}+\epsilon_{1}\eta(t))\theta+\beta_{2}\theta^{3}$
$\displaystyle=\epsilon_{2}\eta(t).$ (6)

The equation (6) phenomenologically models the ship roll motion in irregular
waves incorporating nonlinear resonance and parametric roll mechanisms. We use
empirical parameters $\alpha_{1}=0.01$, $\alpha_{2}=0.052$,
$\beta_{1}=0.1175$, $\beta_{2}=-0.09$, and $\epsilon_{2}=0.004$. The parameter
$\epsilon_{1}$ is varied to model different levels of parametric excitation,
as one important mechanism leading to large responses, with its value chosen
from $\epsilon_{1}=\\{0,0.02,0.04\\}$. Given input $\eta(t)$ and initial
condition, (6) is numerically integrated with a 4th-order Runge-Kutta method
to calculate the roll motion $\theta(t)$.

In addition, we test all methods for four different input wave spectra. The
first three are drawn from spectra of a Gaussian form

$F(k)\sim\exp\frac{-(k-k_{0})^{2}}{2\mathcal{K}^{2}},$ (7)

with significant wave height $H_{s}=12m$, peak (carrier) wavenumber
$k_{0}=0.018m^{-1}$ (corresponding to peak period $T_{p}=15s$). The parameter
$\mathcal{K}$, which models the bandwidth of the spectrum, is chosen from
$\mathcal{K}=\\{0.02,0.03,0.04\\}m^{-1}$ representing spectra of increasing
bandwidth (although all of them are relatively narrow-banded). We term the
three spectra with increasing $\mathcal{K}$ as S1, S2, and S3 hereafter. The
fourth spectrum, termed S4, is in the JONSWAP form with parameters
$H_{s}=12m$, $T_{p}=15s$ and $\gamma=6$. Compared to S3, S4 has similar peak-
mode region but a longer tail in power-law form (instead of exponential tail
in S3).

The true solutions $P_{true}$ in all cases are computed by (1) from a time
series $\eta(t)$ of eight million $T_{p}$. To compute $P_{est}$ using CIC, we
use constant initial conditions of $\theta=\dot{\theta}=0$ for all wave groups
in simulating (6). For PIC, we establish the distribution of initial
conditions using simulation results for thirty thousands $T_{p}$. For NIC, we
simulate (6) starting from two wave groups ahead of the objective $(a,l)$ wave
group. We make sure that sufficient number of random samples are used for
convergent solutions in the computations applying CIC, PIC and NIC.

## RESULTS

###  CIC, PIC, and NIC for initial conditions

We first benchmark the performances of CIC, PIC, and NIC methods in computing
$P_{est}$. In figure 3 we plot the true solution $P_{true}$, as well as
$P_{est}$ from CIC, PIC, and NIC methods as functions of the threshold
$\theta_{t}$, for different spectra S1$\sim$S4 and parametric excitation level
$\epsilon_{1}$. For the same spectrum and a given $\theta_{t}$, $P_{true}$
increases with the increase of $\epsilon_{1}$ (by examining figures in the
same row). This behavior signifies the effect of parametric roll in generating
large responses.

(a) Gaussian Spectrum S1 with $\mathcal{K}=0.02m^{-1}$

(b) Gaussian Spectrum S2 with $\mathcal{K}=0.03m^{-1}$

(c) Gaussian Spectrum S3 with $\mathcal{K}=0.04m^{-1}$

(d) Jonswap Spectrum S4 with $\gamma=6$

Figure 3: First column: the joint PDFs of $(A,L)$ obtained from wave fields
described by different wave spectra. Second to forth columns: $P_{est}$ as a
function of threshold $\theta_{t}$ (in radians) computed by PIC method ( ),
NIC method ( ) and CIC method ( ), as well as true solution $P_{true}$ ( ) for
$\epsilon_{1}$= 0, 0.02 and 0.04 from left to right, and different spectra
from top to bottom.

Comparing to $P_{true}$, the CIC method significantly underestimates the
exceeding probability for large motion thresholds (i.e., $P_{est}\ll P_{true}$
for large $\theta_{t}$) for all cases. The reason for these deviations lies in
the fact that some large motions are induced by non-zero initial conditions
when the ship encounters a wave group, and these scenarios are completely
ignored when applying the CIC method. For example, the parametric roll
phenomenon may take a few wave groups to develop. If the initial condition for
the target group (where the large response from parametric roll is supposed to
happen) is artificially set to zero, the large response is then erroneously
missed in the CIC method.

Relative to CIC, the PIC method shows a better performance in terms of its
better agreement with $P_{true}$ in most cases. We also see that the
performance of PIC is better for larger values of $\epsilon_{1}$, i.e., higher
level of parametric excitation. This behavior can be argued from the following
reasoning: In principle, the PIC method tends to provide an overestimation of
initial conditions because its distribution of initial conditions is drawn
from a time series $\theta(t)$ which includes both those at the starting point
of and within the group (with the latter resulting in an overestimation). This
tendency of overestimation somehow accounts for the increased probability of
parametric roll with the increase of $\epsilon_{1}$, resulting in a favorable
estimation seen in figure 3.

Finally, the results from the NIC method closely follow $P_{true}$ in all
cases (hardly distinguishable with a log axis). We recall that potential
deviation of the NIC method from the computation (1) can be caused by the
neglect of wave phases and frequency modulation of $\eta(t)$ in NIC, and that
the NIC attempts to trigger an initial condition by simulating the previous
two groups (which cannot be validated beforehand). The favorable agreement
achieved in figure 3 indicates important understanding on both aspects: (i)
the effects of wave phases and frequency modulation in the wave group are
indeed not significant to the exceeding probability (or they are of secondary
importance compared to that of varying initial conditions). (ii) the procedure
in the NIC method is sufficient to capture the correct initial condition for
the ship encountering a given wave group.

Figure 4: (upper panel) mean and (lower panel) standard deviation of $R(a,l)$
for three choices of $(a,l)$ parameters for the S1 spectrum, computed by PIC (
), NIC ( ), CIC ( ) and the ground truth ( ).

The performances of the three methods can be further elucidated by examining
their computations regarding wave groups sharing the same $(a,l)$ parameter.
In figure 4, the responses $R(a,l)$ from three selected $(a,l)$ parameters are
plotted, in terms of the mean and standard deviation of $R(a,l)$ from the
three methods and the true solutions. It can be seen that the CIC method, by
applying zero initial conditions, well captures the mean of $R(a,l)$ but
ignores the standard deviation, i.e., the variation in the ensemble of
$R(a,l)$. Under this situation, it can be expected that some initial
conditions leading to exceeding motions are neglected in CIC, resulting in its
underestimation of the exceeding probability. The PIC method, on the other
hand, provides overestimation of both the mean and standard deviation of
$R(a,l)$ for smaller values of $\epsilon_{1}$, but relatively accurate
estimation at larger $\epsilon_{1}$. This explains the performance of PIC
observed in figure 3. Finally, the NIC method captures accurately both the
mean and standard deviation of $R(a,l)$ for chosen $(a,l)$ parameters.

We remark that the above conclusions are drawn based on (6) with a set of
selected parameters. To generalize the conclusion with high confidence, it is
desired to test the results for more parameter values in (6). However,
considering that our current set of parameters (and the equation) are not
chosen in a specific way, one can anticipate that the conclusion holds for
much broader situations of sea conditions and ship motion dynamics.

### VHGP-Seq to accelerate the NIC method

While the NIC method provides very accurate estimation of the exceeding
probability, a direct application with random sampling to compute (4) can be
prohibitively expensive if CFD is used to obtain the motion for given $a$ and
$l$. For example, the results in figure 3 are obtained with 200,000 random
samples of wave groups.

To accelerate the NIC method, we apply our recently developed VHGP-Seq
sampling method to reduce the number of required samples (for an accurate
estimation of the exceeding probability). Figure 5 shows a typical example of
$P_{est}$ computed by the NIC method with random and VHGP-Seq samplings for
the S2 spectrum, $\theta_{t}=0.9$ and $\epsilon_{1}=0.02$. In particular, we
include the mean and standard deviation of $P_{est}$ computed by 100
realizations of both methods. For NIC with random sampling, the mean of
$P_{est}$ always lies in the vicinity of $P_{true}$ due to the basic
properties of Monte-Carlo random sampling (i.e., the Monte-Carlo random
sampling is an unbiased estimator from law of large numbers). However, to
obtain meaningful solution in one realization we require the standard
deviation of $P_{est}$ to be reduced to a low level, i.e., the result
converges to the true solution in a single experiment. From figure 5 it is
clear that it takes random sampling method more than O(1000) samples for the
standard deviation of $P_{est}$ to reach a satisfactory state for estimation.
In contrast, when VHGP-Seq sampling is applied, less than O(100) samples are
required for the convergence of the estimated result to the true solution,
i.e., the mean of $P_{est}$ agreeing with $P_{true}$ and the standard
deviation of $P_{est}$ reduced to a low level.

Figure 5: The exceeding probability for the S2 spectrum, $\theta_{t}=0.9$ and
$\epsilon_{1}=0.02$, with $P_{est}$ computed by the NIC method with VHGP-Seq
sampling ( ) and random sampling ( ). The solid line and shaded region
represent the mean and one-sided standard deviation of $P_{est}$ respectively,
computed from 100 realizations of both methods. The true solution $P_{true}$
is indicated by ( ).

Furthermore, we apply the VHGP-Seq sampling method to more cases selected from
figure (3). In Table 1 we show the results of the NIC method with VHGP-Seq and
random samplings for $\epsilon_{1}=0.02$ and different other selected
parameters. Specifically, we present the mean and standard deviation of
$P_{est}$ computed by 100 realizations of the random and VHGP-Seq samplings.
In all cases, with 120 samples it is found that the VHGP-Seq sampling provides
results with the means closer to $P_{true}$ and standard deviations much
smaller than those from random sampling. With these evidences we conclude that
O(100) samples are sufficient for the VHGP-Seq method to provide satisfactory
estimation of exceeding probability of $O(10^{-3})$ when coupled with the NIC
method.

Table 1: Mean and standard deviation of $P_{est}$ at 120 samples of the VHGP-
Seq and random sampling methods compared to the true solution $P_{true}$ for
$\epsilon_{1}=0.02$ and other selected parameters. The standard deviation is
expressed by its percentage of the mean.

S | $\theta_{t}$ | $P_{true}$ | mean (std) at 120 samples
---|---|---|---
sequential | random
S1 | 0.7 | 0.0421 | 0.0416(5%) | 0.0420(44%)
S2 | 0.9 | 0.0172 | 0.0174(8%) | 0.0156(67%)
S3 | 1 | 0.0232 | 0.0233(8%) | 0.0242(61%)
S4 | 1.05 | 0.0076 | 0.0082(12%) | 0.0105(101%)

## CONCLUSIONS AND DISCUSSION

Figure 6: $\eta(t)$ ( ) generated by JONSWAP spectrum and the detected
Gaussian wave groups ( ).

In this paper, we investigate the uncertainties associated with the wave group
representation of surface elevation series $\eta(t)$ in estimating the motion
exceeding probability. We find that the uncertainties are mainly caused by the
loss of initial conditions of the ship encountering wave groups when the
reduced-order group representation is used. Three methods to handle the
initial condition problem are tested, namely the methods of constant initial
condition (CIC), pre-computed initial condition (PIC) and natural initial
condition (NIC). For a broad range of cases with different input wave spectra,
levels of parametric excitation and motion thresholds, NIC outperforms PIC and
CIC, and provides accurate estimation of the exceeding probability in
remarkable agreement with the true solution. Furthermore, we explore the
application of the VHGP-Seq sampling method to accelerate the computation
using the NIC method, and show that only O(100) samples are needed to obtain
converged result of $P_{est}$ of $O(10^{-3})$ instead of more than O(1000)
samples in random sampling.

We finally remark on the group-based probability used in this paper. The
exceeding probability (1) is always well defined given a group detection from
$\eta(t)$, and can be computed with the methods presented in this paper.
However, for certain cases the practical relevancy of (1) needs to be further
discussed. This is especially for cases with an input wave spectrum consisting
of many short wave components so that the wave group structure cannot be
uniquely defined. An example is the JONSWAP spectrum used in this paper, which
has a realistic power-law tail that decays much slower than the exponential
tail in the Gaussian-form spectrum. Figure 6 shows a segment of $\eta(t)$
described by the JONSWAP spectrum together with the detected groups. While
these groups provide the basis for applying all calculations presented in this
paper, the group-based exceeding probability may vary depending on which and
how many groups (as the number of groups serves as the denominator in
computing the probability) are detected from $\eta(t)$. This information can
vary with the settings in the group detection algorithm, so that the group-
based probability may not be intuitive to interpret in practice.

In general, a more robust measure that we can apply for this type of problem
is the temporal exceeding probability, i.e., ratio of the motion exceeding
time and the total time. This measure remains (largely) an invariant even if
the group-parameterization framework is applied for the computation. More
specifically, as long as the groups that lead to large exceeding motions are
detected (which is practically true), the results of temporal exceeding
probability will remain a constant. However, the computation of temporal
exceeding probability involves other issues to be resolved. We will provide
detailed algorithm and results in an upcoming paper.

## ACKNOWLEDGEMENT

This research is supported by Office of Naval Research grant N00014-20-1-2096.
We thank the program manager Dr. Woei-Min Lin for several helpful discussions
on the research.

## References

  * Anastopoulos and Spyrou (2017) Panayiotis A Anastopoulos and Kostas J Spyrou.  Evaluation of the critical wave groups method for calculating the probability of extreme ship responses in beam seas.  In _Proceedings of ISSW_ , pages 131–138, 2017.
  * Anastopoulos and Spyrou (2019) Panayiotis A Anastopoulos and Kostas J Spyrou.  Evaluation of the critical wave groups method in calculating the probability of ship capsize in beam seas.  _Ocean Engineering_ , 187:106213, 2019.
  * Cousins and Sapsis (2016) Will Cousins and Themistoklis P Sapsis.  Reduced-order precursors of rare events in unidirectional nonlinear water waves.  _Journal of Fluid Mechanics_ , 790:368–388, 2016.
  * Echard et al. (2011) Benjamin Echard, Nicolas Gayton, and Maurice Lemaire.  Ak-mcs: an active learning reliability method combining kriging and monte carlo simulation.  _Structural Safety_ , 33(2):145–154, 2011.
  * Gong and Pan (2022) Xianliang Gong and Yulin Pan.  Sequential bayesian experimental design for estimation of extreme-event probability in stochastic input-to-response systems.  _Computer Methods in Applied Mechanics and Engineering_ , 395:114979, 2022.
  * Gong et al. (2021) Xianliang Gong, Zhou Zhang, Kevin J Maki, and Yulin Pan.  Full resolution of extreme ship response statistics.  _arXiv preprint arXiv:2108.03636, 33rd Symposium on Marine Hydrodynamics_ , 2021.
  * Hu and Mahadevan (2016) Zhen Hu and Sankaran Mahadevan.  Global sensitivity analysis-enhanced surrogate (gsas) modeling for reliability analysis.  _Structural and Multidisciplinary Optimization_ , 53(3):501–521, 2016.
  * Kimura (1980) Akira Kimura.  Statistical properties of random wave groups.  In _Coastal Engineering 1980_ , pages 2955–2973. 1980.
  * Longuet-Higgins (1984) Michael Selwyn Longuet-Higgins.  Statistical properties of wave groups in a random sea state.  _Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences_ , 312(1521):219–250, 1984.
  * Mohamad and Sapsis (2018) Mustafa A Mohamad and Themistoklis P Sapsis.  Sequential sampling strategy for extreme event statistics in nonlinear dynamical systems.  _Proceedings of the National Academy of Sciences_ , 115(44):11138–11143, 2018.
  * Shum and Melville (1984) KT Shum and W Kendall Melville.  Estimates of the joint statistics of amplitudes and periods of ocean waves using an integral transform technique.  _Journal of Geophysical Research: Oceans_ , 89(C4):6467–6476, 1984.
  * Silva and Maki (2021) Kevin M Silva and Kevin J Maki.  Towards a computational fluid dynamics implementation of the critical wave groups method.  _Ocean Engineering_ , 235:109451, 2021.
  * Umeda et al. (2004) Naoya Umeda, Hirotada Hashimoto, Dracos Vassalos, Shinichi Urano, and Kenji Okou.  Nonlinear dynamics on parametric roll resonance with realistic numerical modelling.  _International shipbuilding progress_ , 51(2):205–220, 2004.


# PatchNet: A Simple Face Anti-Spoofing Framework
via Fine-Grained Patch Recognition

Chien-Yi Wang1 Yu-Ding Lu2 Shang-Ta Yang1 Shang-Hong Lai1
1Microsoft AI R&D Center, Taiwan 2HTC
{chiwa, shanya<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Face anti-spoofing (FAS) plays a critical role in securing face recognition
systems from different presentation attacks. Previous works leverage auxiliary
pixel-level supervision and domain generalization approaches to address unseen
spoof types. However, the local characteristics of image captures, i.e.,
capturing devices and presenting materials, are ignored in existing works and
we argue that such information is required for networks to discriminate
between live and spoof images. In this work, we propose PatchNet which
reformulates face anti-spoofing as a fine-grained patch-type recognition
problem. To be specific, our framework recognizes the combination of capturing
devices and presenting materials based on the patches cropped from non-
distorted face images. This reformulation can largely improve the data
variation and enforce the network to learn discriminative feature from local
capture patterns. In addition, to further improve the generalization ability
of the spoof feature, we propose the novel Asymmetric Margin-based
Classification Loss and Self-supervised Similarity Loss to regularize the
patch embedding space. Our experimental results verify our assumption and show
that the model is capable of recognizing unseen spoof types robustly by only
looking at local regions. Moreover, the fine-grained and patch-level
reformulation of FAS outperforms the existing approaches on intra-dataset,
cross-dataset, and domain generalization benchmarks. Furthermore, our PatchNet
framework can enable practical applications like Few-Shot Reference-based FAS
and facilitate future exploration of spoof-related intrinsic cues.

##  1 Introduction

Face anti-spoofing (FAS) is a crucial technique to prevent face recognition
systems from security attacks. With the advance of deep neural network,
several learning-based approaches were proposed to discriminate live faces
from physical presentation attacks.

Figure 1: The face recognition model learns a face embedding space to
discriminate between identities. Our fine-grained patch-type recognition model
learns a patch embedding space to discriminate between patches with different
capture characteristics.

Previous face anti-spoofing methods are highly limited by the scale and
variation of the datasets. Commonly used datasets [1, 18, 27, 5, 34] contain
less than 100 identities during training, and the spoof images are captured
under limited variation. Based on our observation, training on such datasets
with a binary classification model is prone to overfit to the biases
introduced by the data collection, and the learned features are vulnerable in
the unseen testing scenarios. Therefore, previous face anti-spoofing works
[18, 21, 29, 30, 15, 14] leverage auxiliary pixel-wise supervision (e.g., the
facial depth map and reflection map) as a strong prior knowledge to achieve
better generalization ability under testing scenarios with unseen illumination
or spoof types. The other FAS works [19, 33] propose to adopt Generative
Adversarial Network (GAN) to disentangle the feature maps of live faces and
spoof images by reconstructing new live and spoof facial images. Despite the
effectiveness of these spoof-detecting techniques, it is still remained as an
open question to describe the intrinsic cues learned from networks. Yu  _et
al_. [29] rephrase FAS as a structural material recognition problem, which
assumes that the discrimination of the structural materials between human
facial skin and physical spoofing carriers is the essence for FAS tasks.
Following the similar motivation, we believe that the capability of
recognizing and comparing different fine-grained material types is the key to
learn robust intrinsic cues for FAS.

Figure 2: Overview of our proposed PatchNet framework. We address the face
anti-spoofing with a fine-grained patch-type recognition model. The patch-type
classes are pre-defined by the combination of the capture device and the
presenting material, and the patch inputs are extracted from the face captures
by non-distorted augmentation operations. Asymmetric Angular Margin Softmax
Loss is employed in the last classification layer to impose larger angular
margin on live classes. Self-Supervised Similarity Loss is applied to enforce
the patch feature invariance within a single capture.

In this paper, we propose PatchNet which learns discriminative features based
on patches cropped from the entire face regions. Inspired by previous works
[28, 2, 6], the patch-level inputs can enhance the data variation and enforce
the network to learn spoof-specific features in the local region, and thus
prevent the network from overfitting to the biases introduced by datasets.
Instead of resizing the input face images into the same size as adopted by
recent FAS works, we directly crop the fixed-size patches from raw facial
captures to avoid the distortion of discriminative FAS cues. With the patch-
level inputs, our PatchNet aims at classifying the corresponding fine-grained
categories, i.e., the capturing devices and presenting materials, and we
denote each category as a specific “patch-type”. To enforce the network to
learn robust spoof-related feature to recognize unseen patch types during
testing, we adopt the angular margin-based softmax loss that is commonly used
in face recognition tasks [23, 8, 25], which aims to optimize the face
embedding on the normalized hypersphere (Figure 1). Moreover, since the patch
type classes are not symmetric between live and spoof faces, we propose
“asymmetric angular margin loss” and impose a larger margin on live type
classes. Inspired by the recent works on self-supervised learning [3, 4, 10],
and the fact that material patterns are presented spatially in the entire face
region, we also propose “self-supervised similarity loss” to regularize the
features with location and rotation invariance.

To demonstrate the effectiveness of PatchNet, we conduct extensive experiments
on intra-dataset, cross-dataset, and domain generalization benchmark datasets,
and PatchNet achieves the state-of-the-art performance under most testing
scenarios. Moreover, we also conduct the ablation study to further investigate
the proposed components.

Our contributions are summarized as follows:

  * •

We reformulate face anti-spoofing as a fine-grained patch recognition problem,
and design a simple framework called PatchNet to learn an embedding space to
encode intrinsic cues from local patches to represent captures’
characteristics.

  * •

We propose novel Asymmetric Margin-based Softmax Loss and Self-supervised
Similarity Loss to supervise the PatchNet training. While the former helps to
learn a more generalized patch type embedding space to address the asymmetry
between live and spoof, the latter can enforce the patch feature to be
invariant within a single capture.

  * •

The proposed framework could achieve state-of-the-art performance on intra-
dataset, cross-dataset, and domain generalization benchmarks simultaneously
without auxiliary pixel-wise supervision and domain generalization techniques.
Moreover, the learned patch embedding space can enable applications like Few-
Shot reference FAS and patch type retrieval, which can boost the FAS
performance in certain deployment scenarios.

##  2 Related Works

Auxiliary-based Methods. Most of the recent works leverage auxiliary tasks as
the prior knowledge to guide the feature learning toward more generalizable
cues. Liu _et al_. [18] proposed to employ the depth map and rPPG as strong
supervision signals for live samples to regularize the features. Kim _et al_.
[13] further leveraged the reflection maps as the supervision signal for spoof
samples. Many other FAS methods [21, 29, 9, 20, 30, 15, 14] also heavily rely
on similar auxiliary pixel-wise supervision to improve their FAS model
performance. Even though the feature learning can benefit from such
supervision, the pseudo ground truths for those tasks are not accurate, and
the generation of those supervision signals takes high computation resources.

Domain Generalization FAS Methods. In the face anti-spoofing community, domain
generalization techniques are developed to address the domain shift between
different anti-spoofing datasets. Shao _et al_. [21] employed meta-learning
techniques to simulate the target domain shift during the training process to
regularize the feature learning directions. Wang _et al_. [24] proposed to
learn domain-independent features via a disentangled representation learning
framework. The most related work to ours is [12], which treats live and spoof
samples asymmetrically and applies adversarial loss and triplet loss to
regularize the features in the normalized space. Actually, domains are hard to
define in FAS tasks, as even within the same dataset, there are captures with
very different capture devices. While people are using generalization methods
to find common features across collections and spoof types, we aim to break
the concept of domain and propose to learn a generic embedding space that
encodes capture characteristics explicitly.

##  3 Proposed Method

###  3.1 Overview

As illustrated in Fig. 2, we reformulate face anti-spoofing as a fine-grained
patch-type recognition problem and propose a simple training framework to
learn the patch features efficiently. First, we apply certain transform on the
original image to obtain the patch inputs, and the patch features are
extracted by an encoder and then normalized in the feature space. Based on the
meta-info from the training dataset, we split the categories finely based on
the presenting materials and capture devices. For example, in CASIA-FASD,
there are two different spoof mediums and three different capture resolutions,
so there are nine fine-grained patch types (three live and six spoof types).

Inspired by the latest face recognition approaches, during training we employ
angular margin-based softmax loss, which can force feature cluster for each
category to be compactly distributed and enable better generalization ability.
Furthermore, as the distribution discrepancies between spoof samples are
larger than live samples, we treat live and spoof samples asymmetrically:
force the model to learn a more compact cluster within live samples while
leaving spoof samples more dispersed in the feature space. We modify the
angular margin-based softmax loss and apply asymmetric margin onto live and
spoof patch types: imposing larger angular margin on live types to push more
compact boundaries. Finally, the self-supervised similarity loss further
regularizes the patch features by applying the positive part of the
contrastive loss on two transformed patch views from a single whole face
image. Given that spoof-specific discriminative information is present
spatially in the entire face region, the features between two different patch
views from the same face capture should be similar.

###  3.2 Patch Features Extraction

We want to avoid any transform which can lead to image distortion or the
reduction of the important spoof-related information. Given the cropped face
region $x_{i}$ from the raw capture, the two augmented patch views from
$x_{i}$ are $x_{i}^{t_{1}}=t_{1}(x_{i})$ and $x_{i}^{t_{2}}=t_{2}(x_{i})$,
where $t_{1},t_{2}\sim\mathcal{T}$. $\mathcal{T}$ is the sequence of non-
distorted augmentation operations, which only have random horizontal flip,
random rotation, and fixed size cropping. The two input patches are then
passed into the encoder $E_{\theta}$ and normalization layer to get the final
features:
$f_{i}^{t_{1}}=Normalize(E_{\theta}(x_{i}^{t_{1}})),f_{i}^{t_{2}}=Normalize(E_{\theta}(x_{i}^{t_{2}}))$.

###  3.3 Fine-Grained Patch Recognition

Assuming we have $N$ patch type classes in the training dataset, which
consists of $k$ live and $N-k$ spoof classes. Each input patch $t(x_{i})$
belongs to one fine-grained ground truth class
$y_{i}\in\\{L_{1},L_{2},...L_{k},S_{1},S_{2},...S_{N-k}\\}$, and the Angular-
Margin Softmax Loss is applied to regularize the patch features. The Angular-
Margin Softmax Loss has many variants [23, 8, 17, 25, 16] and is commonly used
in face recognition to improve the generalization ability to open-set
identities. In this work, we employ AM-Softmax [23] loss to optimize the fine-
grained patch recognition model and modify it to address the asymmetric nature
in face anti-spoofing.

####  3.3.1 Preliminaries

The formulation of the original Softmax loss is given by

$\displaystyle\mathcal{L}_{S}$
$\displaystyle=-\frac{1}{n}\sum_{i=1}^{n}{log\frac{e^{W_{y_{i}}^{T}\bm{f}_{i}}}{\sum_{j=1}^{c}{e^{W_{j}^{T}\bm{f}_{i}}}}}$
(1)
$\displaystyle=-\frac{1}{n}\sum_{i=1}^{n}{log\frac{e^{\|W_{y_{i}}\|\|\bm{f}_{i}\|cos(\theta_{y_{i}})}}{\sum_{j=1}^{c}{e^{\|W_{j}\|\|\bm{f}_{i}\|cos(\theta_{j})}}}},$

where $\bm{f}$ is the input of the fully connected layer for classification
($\bm{f}_{i}$ denotes the $i$-th sample), $W_{j}$ is the $j$-th column of the
fully connected layer, and $y_{i}$ is the ground truth label of the $i$-th
sample. The term $W_{y_{i}}^{T}\bm{f}_{i}$ is also called the target logit of
the $i$-th sample.

The large-margin property is introduced by Sphereface [16], which defines a
general function $\psi(\theta)$ to impose the angular margin between feature
and weight vectors. After applying feature and weight normalization
($\|W_{y_{i}}\|=\|\bm{f}_{i}\|=1$), the loss function becomes

$\mathcal{L}_{S}=-\frac{1}{n}\sum_{i=1}^{n}{log\frac{e^{\psi(\theta_{y_{i}})}}{e^{\psi(\theta_{y_{i}})}+\sum_{j=1,j\neq
y_{i}}^{c}{e^{cos(\theta_{j})}}}},$ (2)

where in AM-Softmax [23] the function $\psi(\theta)$ is defined as

$\psi(\theta)=cos\theta-m$ (3)

During implementation, the input after normalizing both the feature and the
weight is actually
$x=cos\theta_{y_{i}}=\frac{W_{y_{i}}^{T}f_{i}}{\|W_{y_{i}}\|\|f_{i}\|}$, so in
the forward propagation it only needs to compute

$\displaystyle\Psi(x)$ $\displaystyle=x-m$ (4)

Then it scales the cosine values using a hyper-parameter $s$ and the final AM-
Softmax loss function becomes

$\displaystyle\mathcal{L}_{AMS}$
$\displaystyle=-\frac{1}{n}\sum_{i=1}^{n}{log\frac{e^{s\cdot\left(cos\theta_{y_{i}}-m\right)}}{e^{s\cdot\left(cos\theta_{y_{i}}-m\right)}+\sum_{j=1,j\neq
y_{i}}^{c}{e^{s\cdot cos\theta_{j}}}}}$ (5)
$\displaystyle=-\frac{1}{n}\sum_{i=1}^{n}{log\frac{e^{s\cdot\left(W_{y_{i}}^{T}\bm{f}_{i}-m\right)}}{e^{s\cdot\left(W_{y_{i}}^{T}\bm{f}_{i}-m\right)}+\sum_{j=1,j\neq
y_{i}}^{c}{e^{sW_{j}^{T}\bm{f}_{i}}}}}.$

####  3.3.2 Asymmetric AM-Softmax Loss

We impose different angular margin $m_{l}$ and $m_{s}$ on live and spoof
categories, respectively. Denote live category set as
$L=\\{L_{1},L_{2},...L_{k}\\}$ and spoof category set as
$S=\\{S_{1},S_{2},...S_{N-k}\\}$. The modified AM-Softmax Loss of one feature
sample $f_{i}$ becomes

$\displaystyle\mathcal{L}_{AAMS}(\bm{f}_{i})$
$\displaystyle=\begin{cases}-log\frac{e^{s\cdot\left(W_{y_{i}}^{T}\bm{f}_{i}-m_{l}\right)}}{e^{s\cdot\left(W_{y_{i}}^{T}\bm{f}_{i}-m_{l}\right)}+\sum\limits_{j=1,j\neq
y_{i}}^{N}{e^{s\cdot W_{j}^{T}\bm{f}_{i}}}}&\text{${y_{i}}\in{L}$}\\\
-log\frac{e^{s\cdot\left(W_{y_{i}}^{T}\bm{f}_{i}-m_{s}\right)}}{e^{s\cdot\left(W_{y_{i}}^{T}\bm{f}_{i}-m_{s}\right)}+\sum\limits_{j=1,j\neq
y_{i}}^{N}{e^{s\cdot W_{j}^{T}\bm{f}_{i}}}}&\text{${y_{i}}\in{S}$}\end{cases}$
(6)

The final Asymmetric Recognition Loss on two augmented patch views from the
image is formulated as

$\displaystyle\mathcal{L}_{Asym}=-\frac{1}{n}\sum_{i=1}^{n}(\mathcal{L}_{AAMS}(\bm{f}_{i}^{t_{1}})+\mathcal{L}_{AAMS}(\bm{f}_{i}^{t_{2}}))$
(7)

###  3.4 Self-Supervised Similarity Loss

Given two different patch views from the same face image, the self-supervised
similarity constraint is applied to enforce the features to be similar.
Therefore, the spoof-related feature can be learned with patch location and
rotation invariance.

$\displaystyle\mathcal{L}_{Sim}(\bm{f}_{i}^{t_{1}},\bm{f}_{i}^{t_{2}})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left\|\bm{f}_{i}^{t_{1}}-\bm{f}_{i}^{t_{2}}\right\|_{2}$
(8)

###  3.5 Training and Testing

####  3.5.1 Total Loss

The total loss $L$ of the proposed framework during training is

$\displaystyle\mathcal{L}=\alpha_{1}\mathcal{L}_{Asym}+\alpha_{2}\mathcal{L}_{Sim}$
(9)

where $\alpha_{1}$ and $\alpha_{2}$ are the weights to balance the influence
of loss components. In all of the experiments, we set
$\alpha_{1}=\alpha_{2}=1.0$.

####  3.5.2 Testing Strategy

Given a test face image, we uniformly crop patches from the whole image for
the network inference, with the patch size the same as the one in the training
process. Assuming we have $P$ cropped patches features
($f^{1},f^{2},...f^{P}$) from one face image, then the average live
probability can be obtained by the sum of live class probabilities in the last
fully connected layer:

$\displaystyle LiveProb=\frac{1}{P}\sum_{i=1}^{P}\sum_{y\in L}Softmax(s\cdot
W_{y}^{T}\bm{f^{i}})$ (10)

##  4 Experiments

###  4.1 Datasets and Protocols

Databases. Five databases OULU-NPU [1] (denoted as O), SiW [18] (denoted as
S), CASIA-FASD [34] (denoted as C), Replay-Attack [5] (denoted as I), MSU-MFSD
[27] (denoted as M) are used in the testing protocols. OULU-NPU and SiW are
large-scale high-resolution databases containing four and three protocols to
validate the generalization (e.g., unseen environment and spoof mediums) of
models, respectively, which are utilized for intra-dataset testing. CASIA-
MFSD, Replay-Attack, and MSU-MFSD are databases that contain low-resolution
videos with much fewer video clips and are used for cross-dataset testing to
validate the generalization ability to testing data with large distribution
shift. There are three capture devices with quality ranging from low to high
in CASIA-FASD, two devices in SiW and MSU-MFSD, and only one device in the
other datasets. The fine-grained class number and the other statistics of
databases are shown in Tab. 1. Note that in Oulu-NPU, even the collections are
captured by six different types of phones, the quality and fine details are
pretty similar, so we only split the patch type into five classes in total.
More details and sample images can be found in the supplementary material.

Dataset | # Subjects | # Clips | # Classes
---|---|---|---
Train | Test | Train | Test | Live | Spoof
OULU-NPU (O) | 20 | 20 | 1800 | 1800 | 1 | 4
SiW (S) | 90 | 75 | 2442 | 2036 | 2 | 12
CASIA-FASD (C) | 20 | 30 | 480 | 720 | 3 | 6
ReplayAttack (I) | 30 | 20 | 360 | 240 | 1 | 3
MSU-MFSD (M) | 15 | 20 | 120 | 160 | 2 | 6

Table 1: Statistics of the face anti-spoofing datasets.

Performance Metrics. In intra-dataset testing on OULU-NPU and SiW, we follow
the original protocols and metrics, i.e., Attack Presentation Classification
Error Rate (APCER), Bona Fide Presentation Classification Error Rate (BPCER),
and Average Classification Error Rate (ACER) for a fair comparison. Half Total
Error Rate (HTER) and Area Under Curve (AUC) are adopted in the cross-dataset
testing between OULU-NPU, CASIA-MFSD, Replay-Attack, and MSU-MFSD.

###  4.2 Implementation Details

All face anti-spoofing datasets above are stored in video format originally.
We randomly select three frames from each video clip and use the state-of-the-
art face detector RetinaFace [7] to crop the face for training. We set the
fixed patch crop size as 160, and set the hyperparameter
$s=30.0,m_{l}=0.4,m_{s}=0.1$ in all protocols. We use ResNet18 [11] as the
patch feature encoder, and we did not see much performance difference while
using an encoder with larger capacity (as shown in the supplementary
materials). Models are trained with SGD optimizer and the initial learning
rate is 0.002. We train models with maximum 200 epochs while the learning rate
halves every 90 epochs. During testing, we uniformly crop the fixed-size
patches from the face input image: the minimum x and y coordinates are
$size/2.0$, and the maximum x and y coordinates are $width-(size/2.0)$ and
$height-(size/2.0)$, respectively. In all of the experiments during testing,
we uniformly sample 3 patch anchors on each side, which results in $P=9$
patches for score averaging.

###  4.3 Intra-Dataset Testing

We conduct experiments on Oulu-NPU [1] and SiW [18] for intra-dataset testing
results. We compare the results with the most recent face anti-spoofing
methods in the following.

####  4.3.1 Results on Oulu-NPU

Oulu-NPU [1] has four challenging protocols, which evaluate the model
robustness against the unseen environment, unseen spoof mediums, unseen
capture devices, and all of the above, respectively. The number of classes
during training are 5, 3, 5, and 3, respectively. As shown in Tab. 2, our
simple patch-based recognition approach achieves the best performance in all
protocols. It clearly verifies the better generalization ability of the
features learned through the patch recognition proxy tasks.

Prot. | Method | APCER(%) | BPCER(%) | ACER(%)
---|---|---|---|---
1 | Disentangle [33] | 1.7 | 0.8 | 1.3
SpoofTrace [19] | 0.8 | 1.3 | 1.1
BCN [29] | 0.0 | 1.6 | 0.8
CDCN [32] | 0.4 | 1.7 | 1.0
NAS-FAS [31] | 0.4 | 0.0 | 0.2
PatchNet (Ours) | 0.0 | 0.0 | 0.0
2 | Disentangle [33] | 1.1 | 3.6 | 2.4
SpoofTrace [19] | 2.3 | 1.6 | 1.9
BCN [29] | 2.6 | 0.8 | 1.7
CDCN [32] | 1.5 | 1.4 | 1.5
NAS-FAS [31] | 1.5 | 0.8 | 1.2
PatchNet (Ours) | 1.1 | 1.2 | 1.2
3 | Disentangle [33] | 2.8$\pm$2.2 | 1.7$\pm$2.6 | 2.2$\pm$2.2
SpoofTrace [19] | 1.6$\pm$1.6 | 4.0$\pm$5.4 | 2.8$\pm$3.3
BCN [29] | 2.8$\pm$2.4 | 2.3$\pm$2.8 | 2.5$\pm$1.1
CDCN [32] | 2.4$\pm$1.3 | 2.2$\pm$2.0 | 2.3$\pm$1.4
NAS-FAS [31] | 2.1$\pm$1.3 | 1.4$\pm$1.1 | 1.7$\pm$0.6
PatchNet (Ours) | 1.8$\pm$1.47 | 0.56$\pm$1.24 | 1.18$\pm$1.26
4 | Disentangle [33] | 5.4$\pm$2.9 | 3.3$\pm$6.0 | 4.4$\pm$3.0
SpoofTrace [19] | 2.3$\pm$3.6 | 5.2$\pm$5.4 | 3.8$\pm$4.2
BCN [29] | 2.9$\pm$4.0 | 7.5$\pm$6.9 | 5.2$\pm$3.7
CDCN [32] | 4.6$\pm$4.6 | 9.2$\pm$8.0 | 6.9$\pm$2.9
NAS-FAS [31] | 4.2$\pm$5.3 | 1.7$\pm$2.6 | 2.9$\pm$2.8
PatchNet (Ours) | 2.5$\pm$3.81 | 3.33$\pm$3.73 | 2.9$\pm$3.0

Table 2: The results of testing on OULU-NPU protocols.

Prot. | Method | APCER(%) | BPCER(%) | ACER(%)
---|---|---|---|---
1 | Disentangle [33] | 0.07 | 0.50 | 0.28
SpoofTrace [19] | 0.00 | 0.00 | 0.00
BCN [29] | 0.55 | 0.17 | 0.36
CDCN [32] | 0.07 | 0.17 | 0.12
DualStage [26] | 0.00 | 0.00 | 0.00
| NAS-FAS [31] | 0.07 | 0.17 | 0.12
| PatchNet (Ours) | 0.00 | 0.00 | 0.00
2 | Disentangle [33] | 0.08$\pm$0.17 | 0.13$\pm$0.09 | 0.10$\pm$0.04
SpoofTrace [19] | 0.00$\pm$0.00 | 0.00$\pm$0.00 | 0.00$\pm$0.00
BCN [29] | 0.08$\pm$0.17 | 0.15$\pm$0.00 | 0.11$\pm$0.08
CDCN [32] | 0.00$\pm$0.00 | 0.13$\pm$0.09 | 0.06$\pm$0.04
DualStage [26] | 0.00$\pm$0.00 | 0.00$\pm$0.00 | 0.00$\pm$0.00
| NAS-FAS [31] | 0.00$\pm$0.00 | 0.09$\pm$0.10 | 0.04$\pm$0.05
| PatchNet (Ours) | 0.00$\pm$0.00 | 0.00$\pm$0.00 | 0.00$\pm$0.00
3 | Disentangle [33] | 9.35$\pm$6.14 | 1.84$\pm$2.60 | 5.59$\pm$4.37
SpoofTrace [19] | 8.3$\pm$3.3 | 7.5$\pm$3.3 | 7.9$\pm$3.3
BCN [29] | 2.55$\pm$0.89 | 2.34$\pm$0.47 | 2.45$\pm$0.68
CDCN [32] | 1.67$\pm$0.11 | 1.76$\pm$0.12 | 1.71$\pm$0.11
DualStage [26] | 4.77$\pm$5.04 | 2.44$\pm$2.74 | 3.58$\pm$3.93
| NAS-FAS [31] | 1.58$\pm$0.23 | 1.46$\pm$0.08 | 1.52$\pm$0.13
| PatchNet (Ours) | 3.06$\pm$1.1 | 1.83$\pm$0.83 | 2.45$\pm$0.45

Table 3: The results of testing on SiW protocols.

####  4.3.2 Results on SiW

SiW [18] is another commonly used high-quality dataset with more identities.
The collection is captured by two different quality devices: Canon EOS T6 and
Logitech C920. Compared to Oulu-NPU, it includes more environment variations
and spoof mediums. The numbers of fine-grained patch type classes during
training in protocol 1, 2, 3-1, and 3-2 are 14, 8, 6, and 10, respectively. As
shown in Tab. 3, our method performs the best for the first two protocols and
achieves competitive results in protocol 3.

###  4.4 Ablation Study

In this subsection, all ablation studies are conducted on Protocol 1
(different illumination conditions and location between the train and test
sets) of OULU-NPU [1] to explore the details of our patch-based recognition
framework.

Output Class | Input Extraction | Loss Functions | ACER(%)
---|---|---|---
Binary | Fine | Resize | PatchCrop | $\mathcal{L}_{Asym}$ | $\mathcal{L}_{Sim}$
✓ |  | ✓ |  |  |  | 6.25
✓ |  |  | ✓ |  |  | 3.54
| ✓ | ✓ |  |  |  | 5.63
| ✓ |  | ✓ |  |  | 1.88
| ✓ |  | ✓ |  | ✓ | 1.46
| ✓ |  | ✓ | ✓ |  | 0.63
| ✓ |  | ✓ | ✓ | ✓ | 0.0

Table 4: Ablation study of each component in PatchNet on OULU-NPU protocol 1.

Figure 3: Comparison between choices of patch crop size.

$m_{l}$ | $m_{s}$ | ACER(%)
---|---|---
0.0 | 0.0 | 1.46
0.2 | 0.2 | 0.83
0.2 | 0.1 | 0.41
0.4 | 0.4 | 0.63
0.4 | 0.3 | 0.2
0.4 | 0.2 | 0.2
0.4 | 0.1 | 0.0
0.4 | 0.0 | 0.41

Table 5: Ablation study of margin choices in $L_{Asym}$.

Impact of Each Component. Tab. 4 shows the ablation study of each component in
our proposed framework. The first row is the naive baseline (ACER: 6.25%)
which formulates the face anti-spoofing as a binary classification problem
(trained with the standard Cross Entropy loss) with the resized 256x256 face
input. Surprisingly, by only adopting the fine-grained classes and raw frame
cropping strategies, we can improve the performance significantly to 1.88%
ACER. It shows that the naive baseline model could overfit to the high-level
biases in the public FAS dataset, which only contains limited background and
identities. Moreover, the fine details in cropped patches from raw frames are
very critical to discriminate between different patch types in high-quality
datasets like OULU-NPU. From the lower part of the table, we can observe the
advantage of the proposed margin-based classification loss and the self-
supervised similarity loss. It is clear that both regularization techniques
can facilitate the encoder to learn more intrinsic features related to the
capture device’s characteristics and presenting materials.

Impact of Patch Crop Size. Fig. 3 demonstrates the ACER(%) on OULU-NPU
protocol 1 between different crop sizes. We can observe that larger patch
sizes during training might be prone to overfit to biases from the face
capture. With the regularization of patch recognition loss and patch-based
augmentation to increase training data variance, the overall performance does
not differ much when enlarging the patch size. However, when the patch size is
too small (e.g., 64), the performance degrades significantly as the capture
characteristics can not be learned with very limited information.

Asymmetric Margin Choices. We conduct ablation experiments to verify the
effectiveness of our asymmetric margin design in the angular margin softmax
loss. From Tab. 5, we can observe that adding angular margins can
significantly improve the generalization capability and outperform the model
without any margin. Furthermore, adding a very large margin to spoof patch
types which have more diverse appearances would hurt the discrimination power
of learned features. We find that PatchNet works well on all testing protocols
with the margin choice $m_{l}=0.4$, $m_{s}=0.1$.

###  4.5 Cross-Dataset Testing

####  4.5.1 Experiments between C and I

First, following the related works, CASIA-MFSD (C) [34] and ReplayAttack (I)
[5] are used for cross-dataset experiments, and the results are measured in
HTER. During training, the numbers of fine-grained patch-type classes are 9
and 4 in $C\rightarrow I$ and $I\rightarrow C$ protocols, respectively. The
results are shown in Tab. 6. Given the limited number of clips and low-quality
videos from the ReplayAttack dataset, it is hard to learn generalizable
features which can perform very well on other datasets, so the error rate in
protocol $I\rightarrow C$ is still high compared with $C\rightarrow I$. Our
proposed framework can achieve competitive performance compared to previous
works in both protocols.

Method | Train | Test | Train | Test
---|---|---|---|---
| CASIA-
---
MFSD(C)
| Replay-
---
Attack(I)
| Replay-
---
Attack(I)
| CASIA-
---
MFSD(C)
STASN [28] | 31.5 | 30.9
Disentangle [33] | 22.4 | 30.3
BCN [29] | 16.6 | 36.4
CDCN [32] | 15.5 | 32.6
DC-CDN [30] | 6.0 | 30.1
PatchNet (Ours) | 9.9 | 26.2

Table 6: The results of cross-dataset testing between CASIA-MFSD and Replay-
Attack. The evaluation metric is HTER(%).

Method | O&C&I to M | O&M&I to C | O&C&M to I | I&C&M to O
---|---|---|---|---
HTER(%) | AUC(%) | HTER(%) | AUC(%) | HTER(%) | AUC(%) | HTER(%) | AUC(%)
PatchNet w/ coarse cls [12] | 10.24 | 96.45 | 15.67 | 92.47 | 21.65 | 91.08 | 16.26 | 91.33
PatchNet w/o margin | 10.0 | 96.61 | 18.0 | 91.57 | 17.25 | 90.47 | 15.04 | 92.42
PatchNet w/o $\mathcal{L}_{Sim}$ | 8.9 | 97.42 | 13.44 | 93.99 | 15.1 | 92.10 | 14.24 | 92.93
PatchNet (Ours) | 7.10 | 98.46 | 11.33 | 94.58 | 14.6 | 92.51 | 11.82 | 95.07

Table 7: Evaluations of different components of the proposed method on four
cross-dataset protocols.

Method | O&C&I to M | O&M&I to C | O&C&M to I | I&C&M to O
---|---|---|---|---
HTER(%) | AUC(%) | HTER(%) | AUC(%) | HTER(%) | AUC(%) | HTER(%) | AUC(%)
Auxiliary [18] | 22.72 | 85.88 | 33.52 | 73.15 | 29.14 | 71.69 | 30.17 | 77.61
MADDG [20] | 17.69 | 88.06 | 24.50 | 84.51 | 22.19 | 84.99 | 27.89 | 80.02
PAD-GAN [24] | 17.02 | 90.10 | 19.68 | 87.43 | 20.87 | 86.72 | 25.02 | 81.47
RFM [21] | 13.89 | 93.98 | 20.27 | 88.16 | 17.30 | 90.48 | 16.45 | 91.16
NAS-FAS [31] | 16.85 | 90.42 | 15.21 | 92.64 | 11.63 | 96.98 | 13.16 | 94.18
SSDG-R [12] | 7.38 | 97.17 | 10.44 | 95.94 | 11.71 | 96.59 | 15.61 | 91.54
ANRL [14] | 16.03 | 91.04 | 10.83 | 96.75 | 17.85 | 89.26 | 15.67 | 91.90
DRDG [15] | 15.56 | 91.79 | 12.43 | 95.81 | 19.05 | 88.79 | 15.63 | 91.75
PatchNet (Ours) | 7.10 | 98.46 | 11.33 | 94.58 | 13.4 | 95.67 | 11.82 | 95.07

Table 8: Comparison results between the proposed PatchNet and state-of-the-art
methods on four domain generalization protocols.

####  4.5.2 Domain Generalization Experiments

Some recent FAS works [24, 12, 21] consider each dataset as one domain and
promote the domain generalization benchmark in FAS, which utilizes three
datasets for training, and the remaining one as testing. As we aim to
distinguish the patch type in the fine-grained manner, our proposed framework
can be directly used to evaluate such benchmark without employing further
generalization techniques (e.g., adversarial training or meta-learning). With
access to more different patch types with more diverse capture devices, our
framework is capable of learning more discriminative features through the
patch recognition proxy task. There are four protocols in this benchmark:
O&C&I to M, M&I&O to C, M&C&O to I, and M&C&I to O. During training, we
directly combine the fine-grained patch classes from the three training
datasets, which results in 18, 17, 22, and 21 classes, respectively.

The testing results are shown in Tab. 8. The proposed PatchNet achieves
competitive results on all protocols. Due to the high variance of capture
types in C dataset, it is hard to learn robust local features to address both
high and low resolution scenarios. We also conduct an ablation study in this
benchmark to explore the influence of each component in our framework and show
the results in Tab. 7. In the first ablation experiment, we split the patch
classes using the strategy proposed by SSDG [12]: It aggregates the live
samples as one class and treats spoof samples from each other dataset as one
class, which results in 4 classes. The results verify that the proposed fine-
grained class split, $L_{Asym}$, and $L_{Sim}$ are all important to regularize
the network to generalize better in challenging FAS tasks.

###  4.6 Visualizations

Patch Feature Distribution. In Fig. 4, we visualize the patch features by
t-SNE [22] from OULU-NPU protocol 1, which both the training and testing sets
consist of 5 patch types (live, print1, print2, screen1, screen2). We can
observe in (b) that the model trained without margin cannot distinguish print1
and print2 types very well. The distribution for the live samples is more
compact, and clusters are separated better for the model trained with the
margin. The training and testing feature sets are aligned well in the feature
space.

Figure 4: The t-SNE features visualizations on Oulu-NPU Protocol 1. (a)
training with asymmetric margin (b) training without margin. Figure 5: The
t-SNE visualizations of normalized patch features in cross-dataset protocol
M&I&O to C. The fine-grained patch type class is denoted by
(Dataset)(SensorID)_(Liveness)(MediumID). Figure 6: Patch score map of 5
different types in Oulu-NPU Protocol 1. From left to right: Live, Print1,
Print2, Screen1, and Screen2. The number represents the live probability.

In Fig. 5, we visualize the feature distribution in the protocol M&I&O to C.
We can observe that the features from the training set are separated well in
the fine-grained manner. As dataset C contains three different quality capture
devices, the corresponding three live classes are located separately in the
embedding space. Aligning with the device characteristics, we observe that
features from C3_L are close to O_L, which both types are captured from high-
quality capture devices. Moreover, low-quality devices like C1_L exhibit
similar appearance with that from similar quality device M1_L. It implies that
the patch embedding space encodes the capture characteristic well.

Patch Score Map. In Fig. 6, we compute the scores from patches across the
whole image and visualize the live probability in the overlapped heat map. The
samples are from the testing set of OULU-NPU protocol 1. The patch scores
across the whole face input image are primarily consistent, except for some
background and boundary parts, which is expected as the spoof cues should be
learned from the face region but not the background’s biases.

###  4.7 Applications

Testing on each Capture Device. Following the observation in Sec. 4.6, we
argue that instead of testing on the whole dataset, it is more feasible to
test the anti-spoofing performance on each device separately (e.g., in the
CASIA-FASD dataset). We re-organize the testing dataset in O&C&I to M and
O&M&I to C protocols, which is splitting the testing data by their capture
device ID into 2 (M1, M2) and 3 (C1, C2, C3) groups, respectively. The cross-
dataset anti-spoofing performance on the new split dataset of these two
protocols is shown in Tab. 9. We observe that anti-spoofing performance on
both devices in the M dataset are better than the average performance using
the whole dataset, which is reasonable as the testing sets are smaller.
However, only one device (the high-quality one C3) in the C dataset is better
than the average performance, which aligns with the t-SNE visualization in
Fig. 5. The capture images of C2 have many noises and image compression
effects, which can lead to significant degradation of the discriminative power
of features. With the detailed performance report on each device, we can
further improve the anti-spoofing system or improve the quality of the
problematic capture devices as well.

Few-Shot Reference Anti-Spoofing. With the learned patch embedding space which
encodes intrinsic patch features to discriminate between patch types, the
distance between features in the embedding space can be used to measure the
similarity between patch types. As the features are already normalized in the
space, we can compute the cosine distance between features to enable a new
application: few-shot reference anti-spoofing. While there is a new capture
device, it is easier to acquire some live face image samples, and those sample
features can be used as the reference in the embedding space. We compute the
distances between other testing features and live reference to obtain the
similarity scores. Higher similarity scores mean the samples are more likely
to be live. The 5-shot and 10-shot live reference anti-spoofing performance
are reported in Tab. 9. We can conclude that with the live reference features
and the learned patch embedding space, the performance can be boosted a lot in
practical scenarios.

Method | O&C&I to M | O&M&I to C
---|---|---
(AUC: 98.46%) | (AUC: 94.58%)
M1 | M2 | C1 | C2 | C3
PatchNet | 99.54 | 98.63 | 94.29 | 88.49 | 98.13
PatchNet w/ 5-shot | 99.7 | 99.6 | 94.3 | 89.8 | 98.6
PatchNet w/ 10-shot | 99.8 | 99.6 | 95.3 | 90.7 | 99.2

Table 9: Split testing and few-shot live reference testing on each capture
device. The few-shot testing score is averaged by 10 experiment runs. AUC(%)
score is reported. Figure 7: Patch type retrieval using a live sample from O
dataset as the query. Top retrieval results are shown, and the numbers below
are the cosine similarity.

Patch Type Retrieval. The normalized patch embedding space can be used for the
patch type retrieval application, and can boost the FAS performance on certain
capture devices. For example, after training the feature space with M&C&I
datasets to recognize 21 patch types, we can measure the similarity between
testing patch types and training patch types. As shown in Fig. 7, we take a
live sample from O dataset as the query, and retrieve the training patch types
by computing the cosine distance with each type weight vector. The top-7
retrieval patch types are (M2_L, C3_L, M1_L, M1_S2, C2_L, M2_S2, I_L). Some
spoof types from M dataset have higher ranking than live types from C2 and I
datasets, which have low quality captures. To optimize the performance while
testing on O dataset, we can 1) remove ambiguous live types (C2_L, C1_L, I_L)
during testing, or 2) re-define the class for (C2_L, C1_L, I_L) as “spoof” for
training. Both strategies can boost the performance to 95.27% and 95.87% AUC,
respectively.

##  5 Conclusions and Future Work

In this paper, we reformulate face anti-spoofing as a fine-grained patch type
recognition task and present a simple training framework called PatchNet to
efficiently learn the patch embedding space which encodes the spoof-related
capture characteristics. The novel loss functions are designed to enhance the
feature discrimination power. Extensive experiments on challenging FAS
protocols verify the effectiveness of the proposed method. We note that
exploration of a generic embedding space to discriminate different captures is
still at an early stage. Future directions include: 1) Learning a more
generalized embedding space by datasets with more variations or transferred
from the material perception task, and 2) Investigation into the Few-Shot FAS
protocols which have more practical values.

## References

  * [1] Z. Boulkenafet, J. Komulainen, L. Li, X. Feng, and A. Hadid.  Oulu-npu: A mobile face presentation attack database with real-world variations.  In IEEE International Conference on Automatic Face Gesture Recognition (FG), 2017.
  * [2] R. Cai, H. Li, S. Wang, C. Chen, and A. C. Kot.  Drl-fas: A novel framework based on deep reinforcement learning for face anti-spoofing.  IEEE Transactions on Information Forensics and Security, 2021.
  * [3] Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton.  A simple framework for contrastive learning of visual representations.  In ICML, 2020.
  * [4] Xinlei Chen and Kaiming He.  Exploring simple siamese representation learning.  arXiv preprint arXiv:2011.10566, 2020.
  * [5] Ivana Chingovska, André Anjos, and Sébastien Marcel.  On the effectiveness of local binary patterns in face anti-spoofing.  In BIOSIG-proceedings of the international conference of biometrics special interest group (BIOSIG). IEEE, 2012.
  * [6] Debayan Deb and Anil K Jain.  Look locally infer globally: A generalizable face anti-spoofing approach.  IEEE Transactions on Information Forensics and Security, 2020.
  * [7] Jiankang Deng, Jia Guo, Evangelos Ververas, Irene Kotsia, and Stefanos Zafeiriou.  Retinaface: Single-shot multi-level face localisation in the wild.  In CVPR, 2020.
  * [8] Jiankang Deng, Jia Guo, Niannan Xue, and Stefanos Zafeiriou.  Arcface: Additive angular margin loss for deep face recognition.  In CVPR, 2019.
  * [9] Haocheng Feng, Zhibin Hong, Haixiao Yue, Yang Chen, Keyao Wang, Junyu Han, Jingtuo Liu, and Errui Ding.  Learning generalized spoof cues for face anti-spoofing.  arXiv preprint arXiv:2005.03922, 2020.
  * [10] Jean-Bastien Grill, Florian Strub, Florent Altché, Corentin Tallec, Pierre H Richemond, Elena Buchatskaya, Carl Doersch, Bernardo Avila Pires, Zhaohan Daniel Guo, Mohammad Gheshlaghi Azar, et al.  Bootstrap your own latent: A new approach to self-supervised learning.  arXiv preprint arXiv:2006.07733, 2020.
  * [11] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In CVPR, 2016.
  * [12] Yunpei Jia, Jie Zhang, Shiguang Shan, and Xilin Chen.  Single-side domain generalization for face anti-spoofing.  In CVPR, 2020.
  * [13] Taewook Kim, YongHyun Kim, Inhan Kim, and Daijin Kim.  Basn: Enriching feature representation using bipartite auxiliary supervisions for face anti-spoofing.  In CVPR Workshops, 2019.
  * [14] Shubao Liu, Ke-Yue Zhang, Taiping Yao, Mingwei Bi, Shouhong Ding, Jilin Li, Feiyue Huang, and Lizhuang Ma.  Adaptive normalized representation learning for generalizable face anti-spoofing.  In Proceedings of the 29th ACM International Conference on Multimedia, 2021.
  * [15] Shubao Liu, Ke-Yue Zhang, Taiping Yao, Kekai Sheng, Shouhong Ding, Ying Tai, Jilin Li, Yuan Xie, and Lizhuang Ma.  Dual reweighting domain generalization for face presentation attack detection.  arXiv preprint arXiv:2106.16128, 2021.
  * [16] Weiyang Liu, Yandong Wen, Zhiding Yu, Ming Li, Bhiksha Raj, and Le Song.  Sphereface: Deep hypersphere embedding for face recognition.  In CVPR, 2017.
  * [17] Weiyang Liu, Yandong Wen, Zhiding Yu, and Meng Yang.  Large-margin softmax loss for convolutional neural networks.  In ICML, 2016.
  * [18] Yaojie Liu, Amin Jourabloo, and Xiaoming Liu.  Learning deep models for face anti-spoofing: Binary or auxiliary supervision.  In CVPR, 2018.
  * [19] Yaojie Liu, Joel Stehouwer, and Xiaoming Liu.  On disentangling spoof trace for generic face anti-spoofing.  In ECCV, 2020.
  * [20] Rui Shao, Xiangyuan Lan, Jiawei Li, and Pong C Yuen.  Multi-adversarial discriminative deep domain generalization for face presentation attack detection.  In CVPR, 2019.
  * [21] Rui Shao, Xiangyuan Lan, and Pong C Yuen.  Regularized fine-grained meta face anti-spoofing.  In AAAI, 2020.
  * [22] Laurens Van der Maaten and Geoffrey Hinton.  Visualizing data using t-sne.  Journal of machine learning research, 2008.
  * [23] Feng Wang, Jian Cheng, Weiyang Liu, and Haijun Liu.  Additive margin softmax for face verification.  IEEE Signal Processing Letters, 2018.
  * [24] Guoqing Wang, Hu Han, Shiguang Shan, and Xilin Chen.  Cross-domain face presentation attack detection via multi-domain disentangled representation learning.  In CVPR, 2020.
  * [25] Hao Wang, Yitong Wang, Zheng Zhou, Xing Ji, Dihong Gong, Jingchao Zhou, Zhifeng Li, and Wei Liu.  Cosface: Large margin cosine loss for deep face recognition.  In CVPR, 2018.
  * [26] Yu-Chun Wang, Chien-Yi Wang, and Shang-Hong Lai.  Disentangled representation with dual-stage feature learning for face anti-spoofing.  In WACV, 2022.
  * [27] D. Wen, H. Han, and A. K. Jain.  Face spoof detection with image distortion analysis.  IEEE Transactions on Information Forensics and Security, 2015.
  * [28] Xiao Yang, Wenhan Luo, Linchao Bao, Yuan Gao, Dihong Gong, Shibao Zheng, Zhifeng Li, and Wei Liu.  Face anti-spoofing: Model matters, so does data.  In CVPR, 2019.
  * [29] Zitong Yu, Xiaobai Li, Xuesong Niu, Jingang Shi, and Guoying Zhao.  Face anti-spoofing with human material perception.  In ECCV, 2020.
  * [30] Zitong Yu, Yunxiao Qin, Hengshuang Zhao, Xiaobai Li, and Guoying Zhao.  Dual-cross central difference network for face anti-spoofing.  arXiv preprint arXiv:2105.01290, 2021.
  * [31] Z Yu, J Wan, Y Qin, X Li, SZ Li, and G Zhao.  Nas-fas: Static-dynamic central difference network search for face anti-spoofing.  IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI), 2020.
  * [32] Zitong Yu, Chenxu Zhao, Zezheng Wang, Yunxiao Qin, Zhuo Su, Xiaobai Li, Feng Zhou, and Guoying Zhao.  Searching central difference convolutional networks for face anti-spoofing.  In CVPR, 2020.
  * [33] Ke-Yue Zhang, Taiping Yao, Jian Zhang, Ying Tai, Shouhong Ding, Jilin Li, Feiyue Huang, Haichuan Song, and Lizhuang Ma.  Face anti-spoofing via disentangled representation learning.  In ECCV, 2020.
  * [34] Z. Zhang, J. Yan, S. Liu, Z. Lei, D. Yi, and S. Z. Li.  A face antispoofing database with diverse attacks.  In IAPR International Conference on Biometrics (ICB), 2012.