  We develop a new method based on Gaussian process to reconstruct the mass
distribution of binary black holes (BBHs). Instead of prespecifying the
formalisms of mass distribution, we introduce a more flexible and nonparametric
model with which the distribution can be mainly determined by the observed
data. We first test our method with simulated data, and find that it can well
recover the injected distribution. Then we apply this method to analyze the
data of BBHs' observations from LIGO-Virgo Gravitational-Wave Transient Catalog
2. By reconstructing the chirp mass distribution, we find that there is a peak
or a platform locating at $20-30\,M_{\odot}$ rather than a
single-power-law-like decrease from low mass to high mass. Moreover, one or two
peaks in the chirp mass range of $\mathcal{M}<20\,M_{\odot}$ may be favored by
the data. Assuming a mass-independent mass ratio distribution of $p(q)\propto
q^{1.4}$, we further obtain a distribution of primary mass, and find that there
is a feature locating in the range of $(30, 40)\,M_{\odot}$, which can be
related to \textsc{Broken Power Law} and \textsc{Power Law + Peak}
distributions described in The LIGO Scientific Collaboration et al. (2020).
Besides, the merger rate of BBHs is estimated to $\mathcal{R} =
25.53^{+13.87}_{-8.76}~{\rm Gpc^{-3}~yr^{-1}}$ supposing there is no redshift
evolution.
