I refer to [this paper][1]. There is a straightforward approach in this case. You can separate the refraction index in the following way
$$
   n^2=n_0^2+\xi
$$
where $\xi$ is the random part. Then, the equation is
$$
  \nabla^2U+k^2n_0^2U=-\xi k^2U.
$$
Now, turning to the paper I just cited, take the Green function to be ($d=3$ for the sake of simplicity)
$$
   G(|{\bf r}|) =\frac{e^{ikn_0|{\bf r}|}}{4\pi |{\bf r}|}
$$
and the differential equation turns into an integral one
$$
   U({\bf r}) = U_0({\bf r}) - k^2\int d^3x'G(|{\bf r}-{\bf r'}|)\xi({\bf r'})U({\bf r'}).
$$
One can iterate this equation and obtain the following solution series
$$
   U({\bf r}) = U_0({\bf r}) - k^2\int d^3x'G(|{\bf r}-{\bf r'}|)\xi({\bf r'})U_0({\bf r'})+k^4\int d^3x'G(|{\bf r}-{\bf r'}|)\xi({\bf r'})U_0({\bf r'})\times
\int d^3x''G(|{\bf r'}-{\bf r''}|)\xi({\bf r''})U_0({\bf r''})+\ldots.
$$
By averaging on the distribution for $\xi$, a normal one by hypothesis, one gets a solution series in terms of the correlation function of the random variable. A simple case is this
$$
   \langle\xi({\bf r})\rangle =0 \qquad \langle\xi({\bf r})\xi({\bf r'})\rangle=\xi_0^2\delta^3({\bf r}-{\bf r'})
$$
yielding
$$
   U({\bf r}) \stackrel{?}{=} U_0({\bf r})+k^4\xi_0^2G(0)\int d^3x'G(|{\bf r}-{\bf r'}|)U_0^2({\bf r'})+\ldots.
$$
that is not well defined mathematically. You can consider a case with a finite volume and recover from this situation.



  [1]: https://www.researchgate.net/publication/228448963_Stochastic_Partial_Differential_Equations_as_priors_in_ensemble_methods_for_solving_inverse_problems