  In this paper, we consider the following problem $$ -\Delta u -\zeta
\frac{u}{|x|^{2}} = \sum_{i=1}^{k} \left( \int_{\mathbb{R}^{N}}
\frac{|u|^{2^{*}_{\alpha_{i}}}}{|x-y|^{\alpha_{i}}} \mathrm{d}y \right)
|u|^{2^{*}_{\alpha_{i}}-2}u + |u|^{2^{*}-2}u , \mathrm{~in~} \mathbb{R}^{N}, $$
where $N\geqslant3$, $\zeta\in(0,\frac{(N-2)^{2}}{4})$, $2^{*}=\frac{2N}{N-2}$
is the critical Sobolev exponent, and
$2^{*}_{\alpha_{i}}=\frac{2N-\alpha_{i}}{N-2}$ ($i=1,\ldots,k$) are the
critical Hardy--Littlewood--Sobolev upper exponents. The parameters
$\alpha_{i}$ ($i=1,\ldots,k$) satisfy some suitable assumptions. By using
Coulomb--Sobolev space, endpoint refined Sobolev inequality and variational
methods, we establish the existence of nontrivial solutions. Our result
generalizes the result obtained by Yang and Wu [Adv. Nonlinear Stud. (2017)].
