$\newcommand\Si\Sigma\newcommand\si\sigma\newcommand\R{\mathbb R}\newcommand{\la}{\lambda}$
The calculations can be done as follows. 

Consider the spectral decomposition of $\Si$:
\begin{equation}
	\Si=QDQ^T,\tag{0}
\end{equation}
where $Q$ is an orthogonal matrix $n\times n$ matrix and $D$ is the diagonal matrix with diagonal entries $\si_1^2>0,\dots,\si_k^2>0,0,\dots,0$, so that $k$ is the rank of $\Si$. Let $Y:=Q^TX$, so that $X=QY$ and $Y=(Y_1,\dots,Y_k,0,\dots,0)\sim N(0,D)$. So, for any Borel-measurable function $g\colon\mathbb R^n\to\mathbb R$ such that $Eg(X)$ exists, we have 
\begin{equation}
\begin{aligned}
	Eg(X)&=Eg(Q[Y_1,\dots,Y_k,0,\dots,0]^T) \\ 
&=\int_{\mathbb R^k}g(Q[y_1,\dots,y_k,0,\dots,0]^T) \\
&\times f_{\si_1}(y_1)\cdots f_{\si_k}(y_k)dy_1\cdots dy_k,
\end{aligned}
\tag{1}
\end{equation}
where $f_\si$ is the pdf of $N(0,\si^2)$. 

Now you can compute the probabilities for $X$ as follows: 
$$P(X\in B)=E1(X\in B),$$
where $B$ is any Borel subset of $\mathbb R^n$.

--- 

For an illustration, here is a Mathematica notebook with the calculation of the probability $P(X_1>X_2+1)$ given $(X_1,X_2)\sim N\left((0,0),\left(
\begin{array}{cc}
 1 & 2 \\
 2 & 4 \\
\end{array}
\right)\right)$: 

[![enter image description here][1]][1]


  [1]: https://i.sstatic.net/OnnL0.png

---

As you insist on writing the expectation $Eg(X)$ as an integral with respect to the Lebesgue measure on the support of the distribution of $X$, here it is. 

Note that the support of the distribution of $X$ is the range(=column space) $\Si\R^n$ of $\Si$, which is the same as the range $q\R^k$ of $q$, where $q$ is the map 
\begin{equation*}
	\R^k\ni z=(z_1,\dots,z_k)\mapsto qz:=QJz\in q\R^k=\Si\R^n
\end{equation*}
and, in turn, $J$ is the map 
\begin{equation*}
	\R^k\ni z=(z_1,\dots,z_k)\mapsto Jz:=(z_1,\dots,z_k,0,\dots,0)\in\R^n. 
\end{equation*}
The equality $\Si\R^n=q\R^k$ follows because, in view of (0), the range $\Si\R^n$ of $\Si$ is the same as that of $QD^{1/2}$, which is the same as that of $q$ (since the range of $D^{1/2}$ is the same as that of $J$). Here, we identify matrices with the corresponding linear transformations. 

In view of (1), 
\begin{equation*}
\begin{aligned}
	Eg(X)&=\int_{\R^k}g(qz)\nu(dz)
	=\int_{\Si\R^n}g(x)\mu(dx),
\end{aligned}
\tag{2}
\end{equation*}
where $\nu(dz):=\pi(z)\la_k(dz)$, $\pi(z):=f_{\si_1}(z_1)\cdots f_{\si_k}(z_k)$, $\la_k(dz)$ is the Lebesgue measure on $\R^k$, and $\mu(dx):=\nu(q^{-1}(dx))=\pi(q^{-1}(x))\la_k(q^{-1}(dx))$. 

Since the map $q$ is an isometry, for the Lebesgue measure $\la_{\Si\R^n}$ over $\Si\R^n$ we have $\la_{\Si\R^n}(dx)=\la_k(q^{-1}(dx))$. So, by (2), 
\begin{equation*}
\begin{aligned}
	Eg(X)&=\int_{\Si\R^n}g(x)\pi(q^{-1}(x))\la_{\Si\R^n}(dx). 
\end{aligned}
\tag{3}
\end{equation*}
In particular, letting $g:=1_B$ for any Borel set $B\subseteq\Si\R^n$, we get 
\begin{equation*}
\begin{aligned}
	P(X\in B)&=\int_B \pi(q^{-1}(x))\la_{\Si\R^n}(dx). 
\end{aligned}
\tag{3a}
\end{equation*}

I think formulas (3) and, in particular, (3a) are completely useless, though -- because you will actually **compute** $Eg(X)$  (as in the above Mathematica notebook) by formula (1).