You can use the Fourier Transform theory to seek for a solution. I think you are wrong in making your inverse transforms.
The first task is to transpose your problem into Fourier space. We shall use vector notation letting first $\vec x = (x,y)$, $\vec k=(k_x,k_y)$ the wave vector conjugated in 2D Fourier plane with spatial position $\vec x$ in direct plane. Let's write $\delta(\vec x-\vec x_0)$ the 2D delta distribution in direct plane. We define the direct and inverse Fourier transform of a distribution of
density $f(x)$ by,
$$\hat f(\vec k)=\int f(\vec x) \exp(-i \vec k.\vec x) d\vec x$$
and its inverse (normalization of the FT is inessential in what follows),
$$f(x)=\int \hat{f}(\vec k) \exp(i \vec k.\vec x) d\vec k$$
You seek for solving the Poisson equation,
$$\nabla^2 v(\vec x) = -\frac{q}{\epsilon} \delta(\vec x-\vec x_0)$$
It is very easy to show that the action of the Laplacian operator $\nabla^2 v(\vec x)$ translates on Fourier plane into a multiplication of $\hat v(\vec k)$ with the opposite of the squared modulus of the wave vector $-k^2$. Transforming the RHS of the equation is also straightforward when using the basic properties of the Dirac function. Hence your equation translates on Fourier plane into,
$$\hat v(\vec k)=\frac{q}{\epsilon} \frac{1}{k^2} \exp(-i \vec k.\vec x_0)$$
Take the inverse Fourier transform of the preceding relation side by side,
$$ v(\vec x, \vec x_0) = \frac{q}{\epsilon} \int \frac{1}{k^2} \exp[-i \vec k.(\vec x-\vec x_0)] d\vec k$$
Make use of a system of polar coordinates on the Fourier plane writing the wave vector $\vec k = k \hat{e}_k(\theta)$ with $\hat e_k(\theta)$ the local radial unit vector of cartesian coordinates $(\cos\theta,\sin\theta)$ on the Fourier plane, where $\theta$ is the polar angle. Using this new system of coordinates, the inverse Fourier transforms into,
$$ v(\vec x, \vec x_0) = \frac{q}{\epsilon} \int_{0}^{\infty} dk \int_{0}^{2\pi} \frac{1}{k} \exp(-i k \hat e_k(\theta).(\vec x-\vec x_0) d\theta$$
Similarly, make use of a system of polar coordinates on the spatial plane to write down the displacement $(\vec x-\vec x_0)$ as $(\vec x-\vec x_0)=r \hat e_r(\phi)$ with unit radial vector $\hat e_r(\phi)=(\cos \phi,\sin\phi)$ where $\phi$ is the polar angle on the spatial plane. This yields,
$$v(\vec x,\vec x_0) \equiv v(r=|\vec x-\vec x_0|,\phi) = \frac{q}{\epsilon} \int_{0}^{\infty} dk \int_{0}^{2\pi} \frac{1}{k} \exp[-i kr (\hat e_k(\theta).\hat e_r(\phi))] \ d\theta$$
The scalar product $\hat e_k(\theta).\hat e_r(\phi)$ is simply $\cos(\theta-\phi)$. It is always possible to choose the referent cartesian frame of the polar system such that $\phi=0$. Hence the transform is isotropic (does not depend upon $\phi$) reducing to,
$$ v(r) = \frac{q}{\epsilon} \int_{0}^{\infty} dk \int_{0}^{2\pi} \frac{1}{k} \exp[-i kr \cos\theta] \ d\theta$$
Integration over polar angle $\theta$ is straightforward it gives the Bessel function of the first kind $2\pi J_0(kr)$. As a conclusion you obtain an integral representation of your solution.
$$v(r)=\frac{2\pi q}{\epsilon} \int_{0}^{\infty} \frac{ J_0(kr)}{kr} d(kr)$$
Note that this integral representation is divergent. Introducing a strictly positive lower bound $a>0$ in the low part of the wavenumber spectrum allows to recover a definite solution. Further integrating over new variable $s=kr$,
$$v(r,a)=\frac{2\pi q}{\epsilon} \int_{ar}^{\infty} \frac{ J_0(s)}{s} ds$$
For instance mathematica/wolfram gives a function which is equivalent at the vicinity of $a=0$ to,
$$v(r,a) \sim -\frac{2\pi q}{\epsilon} [\log(\frac{1}{2} ar) +\gamma]$$
where $\gamma$ is the Euler-Mascheroni constant.
Best Answer
It suffices to show \begin{align} \int_{\mathbb{R}^2} G(\textbf{r}, \textbf{r}') \nabla^2f(\textbf{r}')\ d^2\textbf{r}' = f(\textbf{r}). \end{align}
Observe we have \begin{align} \int \log|\textbf{r}-\textbf{r}'|\nabla^2f(\textbf{r}')\ d^2\textbf{r}'&=\ \int_{|\textbf{r}-\textbf{r}'|\leq\ \epsilon} \log|\textbf{r}-\textbf{r}'|\nabla^2f(\textbf{r}')\ d^2\textbf{r}'+\int_{|\textbf{r}-\textbf{r}'|>\epsilon} \log|\textbf{r}-\textbf{r}'|\nabla^2f(\textbf{r}')\ d^2\textbf{r}'\\ &=:\ I_1+I_2. \end{align}
Let us first focus on $I_1$. Note that we have the following estimate \begin{align} \left| I_1\right| \leq&\ \|f\|_{C^2(B(0, 1))} \int_{|\textbf{r}-\textbf{r}'|\leq\ \epsilon} \left|\log|\textbf{r}-\textbf{r}'|\right|\ d^2\textbf{r}'\\ \leq &\ -C\int^\epsilon_0 r\log r\ dr \leq C\epsilon \end{align} which goes to $0$ as $\epsilon \rightarrow 0$. For the $I_2$ term, observe we have that \begin{align} I_2 &=\ \int_{|\textbf{r}-\textbf{r}'|=\epsilon} \log|\textbf{r}-\textbf{r}'|\frac{\partial f}{\partial n}(\textbf{r}')\ dS(\textbf{r}')-\int_{|\textbf{r}-\textbf{r}'|>\epsilon} \nabla\log|\textbf{r}-\textbf{r}'|\cdot \nabla f(\textbf{r}')\ d^2\textbf{r}'\\ &=\ \int_{|\textbf{r}-\textbf{r}'|=\epsilon} \log|\textbf{r}-\textbf{r}'|\frac{\partial f}{\partial n}(\textbf{r}')- \frac{\partial \log|\textbf{r}-\textbf{r}'|}{\partial n}f(\textbf{r}')\ dS(\textbf{r}')\\ &=:\ J_1+J_2. \end{align} For the $J_1$ term we have the estimate \begin{align} |J_1| \leq |\log \epsilon| \int_{|\textbf{r}-\textbf{r}'| = \epsilon} \left|\frac{\partial f}{\partial n}(\textbf{r}')\right|\ dS(\textbf{r}') \leq C|\epsilon \log \epsilon| \end{align} which also goes to $0$ as $\epsilon \rightarrow 0$.
Lastly, observe \begin{align} -\int_{|\textbf{r}'|=\epsilon} \frac{\partial \log|\textbf{r}'|}{\partial n} f(\textbf{r}-\textbf{r}')\ dS(\textbf{r}')= - \int_{|\textbf{r}'| = \epsilon} \frac{f(\textbf{r}-\textbf{r}')}{|\textbf{r}'|} dS(\textbf{r}')\rightarrow -2\pi f(\textbf{r}'). \end{align} as $\epsilon \rightarrow 0$.
Note: If we use Fourier transform for the Poisson equation $\Delta u = \delta$, we get \begin{align} -4\pi^2|\xi|^2 \hat u= 1 \end{align} which means \begin{align} \hat u = \frac{-1}{4\pi^2|\xi|^2}. \end{align} Taking the Fourier inverse of $\hat u$ yields the desired result.