Let $E$ be a supersingular elliptic curve in characteristic $p$. Then
$A = \mathrm{End}(E)$ is a maximal order in the quaternion algebra
$D/\mathbf{Q}$ ramified exactly at $p$ and $\infty$. If $\alpha \in A \smallsetminus \mathbf{Z}$, then Deuring's theorem guarantees that the pair $(E,\alpha)$ can be lifted to a characteristic zero, and the lift will necessarily be a CM elliptic curve. If $K$ is an imaginary quadratic field then $\mathcal{O}_K = \mathbf{Z}[\alpha]$ for some $\alpha$, so the question as to whether $E$ can be lifted to an elliptic curve with CM by the full ring of integers $\mathcal{O}_K$ is equivalent to asking that there is an inclusion
$$\mathcal{O}_K \rightarrow A.$$
A necessary condition is that $K \hookrightarrow D$. This is equivalent to the condition that $K$ splits $D$, which is equivalent to the condition that $p$ does not split in $K$. So from now on restrict to such $p$.
Consider for convenience of exposition the case when $K$ has even discriminant so one can choose $\alpha = \sqrt{-\Delta_K}$. Then $\alpha \in A$ if and only if there exists a trace zero element in $A$ with norm $\Delta_K$. The trace zero elements form a rank three module (over $\mathbf{Z}$), and so the condition that $\Delta_K$ be a norm is equivalent to asking that $\Delta_K$ is represented by the corresponding ternary quadratic form $P_A(x,y,z)$. Understanding what numbers are represented by ternary quadratic forms is a tricky problem in general. In this setting, the ternary quadratic forms $P_A$ depend on $A$. However, these forms are all locally the same for all finite places. So the easy thing to compute is how many ways an integer is represented by $P_A(x,y,z)$ for some $A$. But this is just the number of CM elliptic curves with endomorphisms by
$\mathcal{O}_K$ which is the class number $h_K$. This is why things are much simpler when there is only one supersingular point which occurs for only finitely many $p$. When $p=2$, for example, then $D$ is the Hamilton quaternions and $P_A(x,y,z) = x^2 + y^2 + z^2$, and the number of ways of representing a prime in this form is directly related to class numbers, as was understood by Gauss. For larger $p$, there will be a similar formula for the number of ways of representing numbers in terms of some $P_A$, but understanding the individual factors is more complicated.
In more general settings of ternary definite quadratic forms,
The best one can hope for is that an integer $\Delta_K$ is represented by such a form as long as it is locally representable and $D$ is sufficiently large. (Small $D$ are always going to be difficult.) Even this is too much to ask for general ternary quadratic forms (a good introductory reference is here https://personal.math.ubc.ca/~cass/siegel/hanke-ternary.pdf).
Fortunately it turns out for this particular case that we are in the best situation, and the conclusion is indeed that, fixing $p$, lifts exist for any $K$ where $p$ does not split as long as $\Delta_K \gg_p 1$. But even more is true. We know there are $h_K \gg D^{1/2 - \varepsilon}$ CM elliptic curves by the ring of integers of $K$. It turns out that there reductions are distributed uniformly among the supersingular points as $\Delta_K \rightarrow \infty$. This follows from Theorem 3 of Michel's paper (https://annals.math.princeton.edu/wp-content/uploads/annals-v160-n1-p05.pdf)
Some Variations There are a number of other variations that could be asked. The question (somewhat artificially) considered supersingular elliptic curves over $\mathbf{F}_p$. One could insist on finding CM lifts which are themselves defined over $\mathbf{Z}_p$. To be precise, since $p$ is inert in $K$, these lifts do not actually have CM over the base field, but they
are lifts over $\mathbf{Z}_p$ which have CM over some extension. For convenience let us also restrict to $K/\mathbf{Q}$ of prime discriminant. Then one can ask:
Problem: For a fixed $p$, let $K/\mathbf{Q}$ range over imaginary quadratic fields of prime discriminant. What is the distribution of which supersingular elliptic curves over $\mathbf{F}_p$ lift to a (potentially) CM elliptic curve over $\mathbf{Z}_p$?
To be more specific, we know that when $K$ has prime discriminant then $h_K$ is odd. Let $H_K$ be the class field. The decomposition group at $p$ of $\mathrm{Gal}(H_K/\mathbf{Q})$
surjects onto $\mathrm{Gal}(K/\mathbf{Q})$ because $p$ is inert in this field; since $h_K$ is odd, this implies that the decomposition group in $\mathrm{Gal}(H_K/\mathbf{Q})$ is (any of the) reflections. Since the $j$-invariant of any CM elliptic curve is fixed by exactly one reflection, it follows that there will be a unique elliptic curve over $\mathbf{Z}_p$ which has (potential) CM by $\mathcal{O}_K$. The mod-$p$ reduction of this curve then gives a unique supersingular elliptic curve over $\mathbf{F}_p$ determined by $\Delta_K$. How are these reductions distributed as $\Delta_K$ ranges over all admissible primes?
To make this even more explicit, let $p=11$, so there are two supersingular elliptic curves with $j=0$ and $j = 1728$.
Let $p \equiv 3 \bmod 4$ be a prime such that $11$ is inert in $\mathbf{Q}(\sqrt{-p})$, so $p \not\equiv 2,6,7,8,10 \bmod 11$. There is a class polynomial $P(x)$ of degree $h_K$, and this will have a unique root defined over $\mathbf{Z}_{11}$. For what primes $p$ is this root $0 \bmod 11$ and for what primes is it $1728 \bmod 11$? (This question came up naturally from discussions with David Speyer in the comments). For example, with $\Delta_K = -3, -23, -31, -47$, the reductions are $j = 0,0,0,1728$. The first case is obvious, the modular equations in the next two cases are
$$t^3 + 3491750t^2 - 5151296875t + 12771880859375,$$
$$t^3 + 39491307t^2 - 58682638134t + 1566028350940383,$$
and the last case is computed in the comments). The default guess in this case is presumably that they each occur roughly half the time.
The question of understanding supersingular isogenies is something that seems to be intensively studied for cryptography,
so these questions may well be known.
