The linear inverse problem of estimating the input random field in a first-kind stochastic integral equation relating two random fields is considered. For a wide class of integral operators, which includes the positive rational functions of a self-adjoint elliptic differential operator on L2(ℝd), the ill-posed nature of the problem disappears when such operators are defined between appropriate fractional Sobolev spaces. In this paper, we exploit this fact to reconstruct the input random field from the orthogonal expansion (i.e. with uncorrelated coefficients) derived for the output random field in terms of wavelet bases, transformed by a linear operator factorizing the output covariance operator. More specifically, conditions under which the direct orthogonal expansion of the output random field coincides with the integral transformation of the orthogonal expansion derived for the input random field, in terms of an orthonormal wavelet basis, are studied.