  The cosmic expansion rate can be directly measured with gravitational waves
(GWs) of the compact binary mergers, by jointly constraining the mass function
of the population and the cosmological model via the so called spectral sirens.
Such a method relies on the features in the mass functions, which may originate
from some individual sub-populations, and hence become blurred/indistinct due
to the superposition of different sub-populations. In this work we propose a
novel approach to constrain the cosmic expansion rate with sub-populations of
GW events, named multi-spectral sirens. We illustrate the advantage of the
multi-spectral sirens compared to the traditional spectral sirens by simulation
with mock data. The application of this approach to the GWTC-3 data yields
$H_0=73.25^{+29.87}_{-25.55}~{\rm Mpc}^{-1}~{\rm km}~{\rm s}^{-1}$ (median and
symmetric 68.3\% credible level). The incorporation of the bright standard
siren GW170817 with a uniform prior in $ [10,200] ~{\rm Mpc}^{-1}~{\rm km}~{\rm
s}^{-1}$ gives $H_0=72.38^{+15.03}_{-9.13}~{\rm Mpc}^{-1}~{\rm km}~{\rm
s}^{-1}$ (68.3\% confidence level), corresponding to an improvement of
$\sim28\%$ with respect to the measurement from sole GW170817.
