Are there alternatives to the bilinear transform?
Analog filters are stable if the poles are in the left half of the s-plane (figure on the left) and digital filters are stable if the poles are inside the unit circle (figure on the right). So mathematically all that is needed to convert from analog to digital is a mapping (conformal?) from the half-space to the unit disk and the $\jmath\Omega$ axis to the unit circle $\vert z\vert=1$. Any transformation that does this is a possible candidate for being an alternative to the bilateral transform.

Two of the well known methods are the impulse invariance method and the matched Z-transform method. Conceptually, both of these are similar to sampling a continuous waveform that we're familiar with. Denoting the inverse Laplace transform by $\mathcal{L}^{-1}$ and the Z transform as $\mathcal{Z}$, both these methods involve calculating the impulse response of the analog filter as
$$a(t)=\mathcal{L}^{-1}\{A(s)\}$$
and sampling $a(t)$ at a sampling interval $T$ that is high enough so as to avoid aliasing. The transfer function of the digital filter is then obtained from the sampled sequence $a[n]$ as
$$D_a(z)=\mathcal{Z}\{a[n]\}$$
However, there are key differences between the two.
Impulse invariance method:
In this method, you expand the analog transfer function as partial fractions (not in the matched Z transform as mentioned by Peter) as
$$A(s)=\sum_m \frac{C_m}{s-\alpha_m}$$
where $C_m$ is some constant and $\alpha_m$ are the poles. Mathematically, any transfer function with a numerator of lesser degree than the denominator can be expressed as a sum of partial fractions. Only low-pass filters satisfy this criterion (high-pass and bandpass/bandstop have at least the same degree), and hence impulse invariant method cannot be used to design other filters.
The reason why it fails is also quite clear. If you had a polynomial in the numerator of the same degree as in the denominator, you will have a free standing constant term, which upon inverse transforming, will give a delta function that cannot be sampled.
If you carry out the inverse Laplace and forward Z transforms, you'll see that the poles are transformed as $\alpha_m \to e^{\alpha_m T}$ which means that if your analog filter is stable, the digital filter will also be stable. Hence it preserves the stability of the filter.
Matched Z-transform
In this method, instead of splitting the impulse response as partial fractions, you do a simple transform of both the poles and the zeros in a similar manner (matched) as $\beta_m\to e^{\beta_m T}$ and $\alpha_m\to e^{\alpha_m T}$ (also stability preserving), giving
$$A(s)=\frac{\prod_m (s-\beta_m)}{\prod_n (s-\alpha_n)}\longrightarrow \frac{\prod_m \left(1-z^{-1}e^{\beta_m T}\right)}{\prod_n \left(1-z^{-1}e^{\alpha_n T}\right)}$$
You can easily see the limitation of both these methods. Impulse invariant is applicable only if your filter is low pass and matched z-transform method is applicable to bandstop and bandpass filters (and high pass up to the Nyquist frequency). They are also limited in practice by the sampling rate (after all, you can only go up to a certain point) and suffer from the effects of aliasing.
The bilinear transform is by far the most commonly used method in practice and the above two are rather more for academic interests. As for conversion back to analog, I'm sorry but I do not know and can't be of much help there as I hardly ever use analog filters.