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} = 26.29^{+14.21}_{-8.96}~{rm Gpc^{-3}~yr^{-1}}$ supposing there is no redshift evolution.
Download