My answer http://mathoverflow.net/questions/152912/wiener-measure-and-bochner-minlos/164658#164658 can help defining a free massless real scalar field as a random element of $S_0'(\mathbb{R}^2)$. Namely, take the Schwartz space of rapidly decaying test functions $S(\mathbb{R}^2)$ and consider the subspace $S_0(\mathbb{R}^2)=\{f| \widehat{f}(0)=0\}$ of "charge-neutral" test functions. The bilinear form
$$
B(f,g)=\int_{\mathbb{R}^2} \frac{d^2\xi}{(2\pi)^2}\ \frac{\overline{\widehat{f}(\xi)}\widehat{g}(\xi)}{|\xi|^2}
$$
is continuous and positive on $S_0(\mathbb{R}^2)$ and therefore by the Bochner-Minlos Theorem
there exists a unique centered Gaussian probability measure $\mu$ on the topological dual $S_0'(\mathbb{R}^2)$
for which
$$
\mathbb{E}(\phi(f)\phi(g))=B(f,g)
$$
where $\phi$ is the corresponding random element in $S_0'(\mathbb{R}^2)$.

I did not do the computation (which needs a lot of care) but I suspect that $B(f,g)$
should be a multiple of
$$
\int_{\mathbb{R}^2\times\mathbb{R}^2} d^2 x\ d^2y\ f(x)(-\log|x-y|)g(y)\ .
$$
The issue here is that you have both a UV and an IR problem to deal with.
The slow decay of the propagator for large $\xi$ makes it so that $X$ or $\phi$ is a random distribution rather than a random function and, in particular, punctual evaluations $X(x)$ do not make sense. Also, the divergence at $\xi=0$, or zero-mode, makes it so that there is an ambiguity of shifting the field by a constant so the field $X$ does not make sense but its "increments" do. Working with $\partial X$ is another way to circumvent this issue.


----------
Edit: I just did the computation and indeed
$$
B(f,g)=\frac{1}{2\pi}\int_{\mathbb{R}^2\times\mathbb{R}^2} d^2 x\ d^2y\ f(x)(-\log|x-y|)g(y)\ .
$$
This can be done by replacing $1/|\xi|^2$ by $1/|\xi|^{2-\epsilon}$ and taking the $\epsilon\rightarrow 0$ limit by dominated convergence. Then, use the representation
$$
\frac{1}{|\xi|^{2-\epsilon}}=\frac{1}{\Gamma\left(\frac{2-\epsilon}{2}\right)}
\int_{0}^{\infty}
d\alpha\ \alpha^{-\frac{\epsilon}{2}} e^{-\alpha |\xi|^2}\ ,
$$
write the Fourier transforms as integrals and integrate over $\xi$.
One then ends up with the computation of
$$
\lim_{\epsilon\rightarrow 0^+}
\int_{\mathbb{R}^2\times\mathbb{R}^2} d^2 x\ d^2y\ f(x)\left(\frac{1}{\epsilon}
|x-y|^{-\frac{\epsilon}{2}}\right)g(y)\ .
$$
Since $\int f(x) d^2 x=\int g(x) d^2 x=0$, the last expression is the same as
$$
\lim_{\epsilon\rightarrow 0^+}
\int_{\mathbb{R}^2\times\mathbb{R}^2} d^2 x\ d^2y\ f(x)\left\{\frac{1}{\epsilon}
\left(|x-y|^{-\frac{\epsilon}{2}}-1\right)\right\}g(y)\ .
$$
The resulting derivative at $\epsilon =0$ produces the wanted logarithm of $|x-y|$.


----------

Edit 2:
A substantial elaboration on the answer I gave above has appeared recently. See the review
<a href="http://arxiv.org/abs/1407.5605">"Log-correlated Gaussian fields: an overview"</a> by Duplantier, Rhodes, Sheffield and Vargas.
This is part of a wider program regarding the study of stationary, isotropic and self-similar Gaussian fields:

 1. <a href="http://arxiv.org/abs/1407.5598">"Fractional Gaussian fields: a survey"</a>
by Lodhia, Sheffield, Sun and Watson.
 2. Chapter 2 of <a href="http://libra.virginia.edu/catalog/libra-oa:7070">"Construction and Analysis of a Hierarchical Massless Quantum Field Theory"</a>, Ph.D. Thesis by Ajay Chandra.
 3. <a href="http://www.jstor.org/stable/2242835">"Gaussian and their Subordinated Self-similar Random Generalized Fields"</a> by Roland Dobrushin and its follow-up
<a href="http://link.springer.com/book/10.1007/BFb0094036">" Multiple Wiener-Itô Integrals"</a> by Péter Major.