$\newcommand{\x}{\mathbf x}\newcommand{\tx}{\tilde{\mathbf x}}\newcommand{\0}{\mathbf0}\newcommand{\R}{\mathbb R}$You are citing Besag's paper (formulas (3.1) and (3.3) there) incorrectly. What the paper actually has is this:
\begin{equation*}
\begin{aligned}
Q(\x)&:=\ln\frac{P(\x)}{P(\0)}
\\
&= \sum_{1\le i\le n} x_iG_i(x_i)
+\sum_{1\le i<j\le n}x_ix_jG_{i,j}(x_i,x_j)\\
&+ \sum_{1\le i<j<k\le n} x_ix_jx_kG_{i,j,k}(x_i,x_j,x_k) +\cdots \\
&+x_1\cdots x_nG_{1,\dots,n}(x_1,\ldots,x_n).
\end{aligned}
\tag{1}\label{1}
\end{equation*}
This is indeed a version of the inclusion-exclusion formula. To show this, identify the vector $\x=(x_1,\dots,x_n)\in\R^n$ with the function $[n]\ni i\mapsto x_i$, which will also be denoted by $\x$; here, as usual, $[n]:=\{1,\dots,n\}$. For each $J\subseteq[n]$, let $\tx_J$ be the vector/function on $[n]$ such that $(\tx_J)(i):=x_i$ for $i\in J$ and $(\tx_J)(i):=0$ for $i\notin J$. Let $\x_J:=\x|_J=(x_i\colon i\in J)$, the restriction of the function $\x$ to $J$, so that $\tx_J$ is the extension by zeroes of the the function $\x_J$ from $J$ to $[n]$. Let $|J|$ denote the cardinality of $J$.
For any $K\subseteq[n]$, let
\begin{equation*}
H_K(\x_K):=\sum_{J\subseteq K}(-1)^{|K|-|J|}Q(\tx_J);
\end{equation*}
note that the last expression depends on $\x$ only through $\x_K$, so that $H_K(\x_K)$ is well defined.
Then
\begin{equation*}
\begin{aligned}
&\sum_{K\subseteq[n]}H_K(\x_K) \\
&=\sum_{K\subseteq[n]}\sum_{J\subseteq K}(-1)^{|K|-|J|}Q(\tx_J) \\
&=\sum_{J\subseteq[n]}Q(\tx_J)
\sum_{K: J\subseteq K\subseteq[n]}(-1)^{|K|-|J|} \\
&=\sum_{J\subseteq[n]}Q(\tx_J)
\sum_{k=|J|}^n(-1)^{k-|J|}\binom{n-|J|}{k-|J|} \\
&=\sum_{J\subseteq[n]}Q(\tx_J)
\sum_{m=0}^{n-|J|}(-1)^m\binom{n-|J|}m \\
&=\sum_{J\subseteq[n]}Q(\tx_J)
\,1(|J|=n) \\
&=Q(\tx_{[n]})=Q(\x).
\end{aligned}
\tag{2}\label{2}
\end{equation*}
Note also that $H_K(\x_K)=0$ if $x_j=0$ for some $j\in K$; this follows in view of the natural bijection between the (set of all subsets $J$ of the -- necessarily nonempty -- set $K$ such that $J\ni j$) and (the set of all subsets $J$ of $K$ such that $J\not\ni j$).
So, $H_K(\x_K)=0$ if $p_K(\x)=0$, where
\begin{equation*}
p_K(\x):=\prod_{j\in K}x_j.
\end{equation*}
So, letting $G_K(\x_K):=H_K(\x_K)/p_K(\x)$ if $p_K(\x)\ne0$ and $G_K(\x_K):=0$ if $p_K(\x)=0$, we get $H_K(\x_K)=p_K(\x)G_K(\x_K)$ for all $\x$ and all $K$. Thus, by \eqref{2},
\begin{equation*}
Q(\x)=\sum_{K\subseteq[n]}p_K(\x)G_K(\x_K), \tag{3}\label{3}
\end{equation*}
which is just another, more compact way of writing \eqref{1}.
Of course, one may refer to $G_K(\x_K)$ as a partial divided difference (of the function $Q$) of order $|K|$ with respect to $\x_K=(x_i\colon i\in K)$. In particular,
\begin{equation*}
G_{\{1\}}(x_1)=\frac{Q(x_1,0,\dots,0)-Q(0,0,\dots,0)}{x_1}
\Big(=\frac{Q(x_1,0,\dots,0)}{x_1}\Big)
\end{equation*}
if $x_1\ne0$ and
\begin{equation*}
G_{\{1,2\}}(x_1,x_2)=\frac{Q(x_1,x_2,0,\dots,0)-Q(x_1,0,0,\dots,0)
-Q(0,x_2,0,\dots,0)+Q(0,0,\dots,0)}{x_1x_2}
\end{equation*}
$x_1x_2\ne0$.
So, \eqref{1} and, equivalently, \eqref{3} may be considered difference analogues of the Maclaurin expansion for functions of $n$ variables.
Representations \eqref{1} and, equivalently, \eqref{3} are essentially unique: the values of $G_K(\x_K)$ are uniquely determined for all $\x$ and $K$ such that $p_K(\x)\ne0$. This can be verified by reasoning quite similar to \eqref{2}; this is left as an exercise for now.
