I have a random field $h(\vec{r})$ that depends on $\vec{r}=(x,y)$, such that
\begin{equation} \langle h(\vec{r})h(\vec{r}+\vec{r}') \rangle \sim \exp(-||\vec{r}-\vec{r}'||/a^2) \end{equation}
where the average is over all $(x,y)$. $h(\vec{r})$ is also isotropic in $(x,y)$, so you can also choose a one-dimensional profile, and it has mean $0$.
If I want to determine
\begin{equation} \langle \frac{\partial h(\vec{r})}{\partial x}\frac{\partial h(\vec{r})}{\partial y} \rangle, \langle \frac{\partial h(\vec{r})}{\partial y}\frac{\partial h(\vec{r})}{\partial y} \rangle, \langle \frac{\partial h(\vec{r})}{\partial x}\frac{\partial h(\vec{r})}{\partial x} \rangle, \end{equation}
what is the proper way to approach it? For the terms without mixed derivatives, it makes sense to assume $\partial y = \partial r$ and compute the autocorrelation for the mean squared derivative, but how would one tackle the problem given the mixed derivative?
Like I said in the comment I think you actually mean random field.
Also I'm not sure how the right hand side of your expectation can still depend on $\vec{r}$ if an average has been taken over $(x,y)$.
It would be helpful to actually write out what the averages mean. I think perhaps you mean something like this
$$C(\vec{r}') = \langle h(\vec{r})h(\vec{r}+\vec{r}')\rangle = \int_{-\infty}^\infty \int_{-\infty}^\infty h((x,y)) h((x+x',y+y'))dx dy \\ = \exp(-\|\vec{r}'\|/a^2)$$
Then you are interested in quantities like $$C_{(1,1)}(\vec{r}') = \langle h_x(\vec{r})h_y(\vec{r}+\vec{r}')\rangle = \int_{-\infty}^\infty \int_{-\infty}^\infty h_x((x,y)) h_y((x+x',y+y'))dx dy = \int_{-\infty}^\infty\left(h(x,y) h_y(x+x',y+y')|_{x = -\infty}^{x=\infty} - \int_{-\infty}^\infty h(x,y)h_{x',y'}(x+x',y+y')dx\right)dy$$
You can probably make some argument that the boundary type terms $h(x,y) h_y(x+x',y+y')|_{x = -\infty}^{x=\infty}$ go away.
$$C_{(1,1)}(\vec{r}') = -\int_{-\infty}^\infty\int_{-\infty}^\infty h(x,y)h_{x',y'}(x+x',y+y')dx dy $$
So then of course you would like to commute the integrals and derivatives.
$$ -\int_{-\infty}^\infty \int_{-\infty}^\infty h(x,y)h_{x',y'}(x+x',y+y')dx dy = \\ -\frac{\partial}{\partial x'}\frac{\partial}{\partial y'} \int_{-\infty}^\infty \int_{-\infty}^\infty h((x,y)) h((x+x',y+y'))dx dy$$
In order to justify the above you need to know that the LHS and RHS integrals converge uniformly and absolutely. Presuming that you can prove that then you would have:
$$ C_{(1,1)}(\vec{r}') = -\frac{\partial}{\partial x'}\frac{\partial}{\partial y'}\exp(-\sqrt{x'^2 + y'^2}/a^2)$$