There is an extensive literature on wave equations with variable coefficients, which are equivalent to wave equations on curved spacetime - as in the monograph by Friendlander. Depending on what you are willing to accept as 'explicit', the answer is that yes, there is an explicit generalisation of the Kirchhoff integral. However, it is very difficult to calculate the coefficients of the data $f,g$ in this integral, which are determined by the retarded Green function $G_R$ of the wave operator: $$\square = \partial_t^2 - c(x)\Delta.$$ This Green function is intimately related to the geometry of the spacetime with line element $ds^2=dt^2-c(x)d\vec{x}\cdot d\vec{x}$. In order to calculate $G_R$, you essentially need to fully solve the geodesic equations obtained by extremising the action $\int ds$.
The retarded Green's function satisfies
$$ \square G_R(t,x;t',x')=-4\pi\delta_4(t,x;t',x'),$$ where $\delta_4$ is the 4-dimensional Dirac distribution. Then the generalized Kirchhoff formula allows us to explicitly write down the value of a solution $\Psi$ of the wave equation at a point $(t,x)$ to the future of an initial data hypersurface $\Sigma^\prime=\{(t',x'):x'\in\mathbb{R}^3\}$:
$$ \Psi(t,x)=-\frac{1}{4\pi}\int_{\Sigma^\prime}(G_R(t,x;t',x')g(x')-f(x')\partial_{t'}G_R(t,x;t',x'))d\Sigma^\prime,$$
where $d\Sigma^\prime=c^{3/2}(x')d_3x'$ is the volume element on $\Sigma^\prime$.
These things (Green's function, waves in curved spacetime) are of great interest in General Relativity, and the GR literature is a good place to look for more details. In particular, I recommend starting with Poisson's Living Review article which covers the geometrical background in a very readable way. You'll find details here (and in the citations) on how to calculate $G_R$, which is required for practical applications of the Kirchhoff formula.
You are given a very special initial data. Let $f(x,y,z) = 3x - y + z$, note that this is a linear function on $\mathbb{R}^3$. Your initial data is
$$ \phi = f e^f $$
Why is linearity important? that means we can change coordinates with isometries.
Let $v_1 = \frac{1}{\sqrt{11}}\begin{pmatrix} 3 \\ -1 \\ 1\end{pmatrix}$ and complete this to an orthonormal basis $\{v_1, v_2, v_3\}$. Let $O$ be the orthogonal matrix whose column vectors are $v_1, v_2, v_3$. Now define the coordinates $x',y',z'$ by
$$ \begin{pmatrix} x' \\ y' \\ z' \end{pmatrix} = O \begin{pmatrix} x \\ y \\ z\end{pmatrix} $$
we have that in this coordinate, the function $f$ can be written as $f(x',y',z') = \sqrt{11} x'$.
Since this change of coordinates is orthogonal in $\mathbb{R}^3$, it preserves the Laplace operator. So in this new coordinate system your function $u$ still solves
$$ \frac{\partial^2}{\partial t^2} u = 4 \left( \frac{\partial^2}{\partial x'^2} u + \frac{\partial^2}{\partial y'^2} u + \frac{\partial^2}{\partial z'^2} u \right) $$
But now you are given initial data $\phi$ which is independent of $y'$ and $z'$. This implies that the solution will also be independent of $y'$ and $z'$. This means that $u$ in fact solves
$$ \partial^2_t u = 4 \partial^2_{x'} u $$
a one dimensional wave equation with initial data $\phi(x',y',z') = \phi(x') = \sqrt{11}x' e^{\sqrt{11} x'} $. Now you can use the one-dimensional representation formula for the solution of the wave equation to get precisely that
$$ u(t,x') = \frac12\phi( x' + 2 t) + \frac12\phi(x' - 2t) $$
which after you change the variables back to $x,y,z$ is precisely what you are looking for.
Best Answer
Because now we are considering the $u$ as a funcion of the two new variables then
$$u(x,t)=u\bigg(\frac{\xi+\eta}{2}, \frac{\xi-\eta}{2c}\bigg)=f(\xi)+g(\eta)=u(\xi,\eta)$$
The confusion arises because the same symbol is used for $u$ but it is justified since, even if the expressions for the two function are different, they represent the same function.