How to simulate a rough gaussian field?

I want to simulate a random Gaussian field $h(x)$ such that:
$$
sqrt{langleleft(h(x+r)-h(x)right)^2rangle}simleft(frac{r}{gamma}right)^u,
$$
and
$$langle h(x)rangle =0,$$
where $u$ is a parameter and $r$ a generic distance between points.
How I can do that? My plan was to obtain a Gram Matrix and then plug it into the numpy function:

numpy.random.multivariate_normal

But I am a lot confused, do you have suggestions about getting the Gram matrix?