(-\Delta)^\alpha  + h(u)}
  u=0& {\rm in}\quad{\mathbb{R}}^N\setminus\Omega,
\end{array}$$ where $h:{\mathbb{R}}_+\to{\mathbb{R}}_+$ is a continuous nondecreasing function and $(-\Delta)^\alpha $ denotes the fractional Laplacian of exponent $\alpha$ defined by $$(-\Delta)^\alpha  u(x)=\lim_{\epsilon\to0^+} (-\Delta)_\epsilon^\alpha u(x),$$ where for $\epsilon>0$, $$(-\Delta)_\epsilon^\alpha  u(x)=-\int_{{\mathbb{R}}^N}\frac{ u(z)-
u(x)}{|z-x|^{N+2\alpha}}\chi_\epsilon(|x-z|) dz$$ and $$\chi_\epsilon(t)=\left\{ \arraycolsep=1pt
\begin{array}{lll}
0,\quad & {\rm if}\quad t\in[0,\epsilon],\\[2mm]
1,\quad & {\rm if}\quad t>\epsilon.
\end{array}
\right.$$

In the pioneering work [@B12] (also see [@BB]), Brezis studied the existence of weak solutions to second order elliptic problem $$\label{eq003}
 \arraycolsep=1pt
\begin{array}{lll}
 -\Delta  u+h(u)=\nu \quad & {\rm in}\quad\Omega,\\[2mm]
 \phantom{ -\Delta  +g(u)}
u=0  \quad & {\rm on}\quad \partial\Omega,
\end{array}$$ where $\Omega$ is a bounded $C^2$ domain in ${\mathbb{R}}^N(N\ge3)$, $\nu$ is a bounded Radon measure in $\Omega$, and the function $h:{\mathbb{R}}\to{\mathbb{R}}$ is nondecreasing, positive on $(0,+\infty)$ and satisfies the subcritical assumption: $$\int_1^{+\infty}(h(s)-h(-s))s^{-2\frac{N-1}{N-2}}ds<+\infty.$$ In particular case that $0\in\Omega$, $h(s)=s^q$ and $\nu=k\delta_0$ with $k>0$, Brezis et al showed that (\[eq003\]) admits a unique weak solution $v_k$ for $1< q< N/(N-2)$, while no solution exists if $q\geq  N/(N-2)$. Later on, Véron in [@V0] proved that if $1< q< N/(N-2)$, the limit of $v_k$ is a strong singular solution of $$\label{local}
 \arraycolsep=1pt
\begin{array}{lll}
 -\Delta  u+u^q=0 \quad & {\rm in}\quad\Omega\setminus\{0\},\\[2mm]
 \phantom{ -\Delta  +u^q }
u=0  \quad & {\rm on}\quad \partial\Omega.
\end{array}$$ After 