In complex stochastic high-dimensional reliability studies, polynomial chaos expansion (PCE) has been widely used to build surrogate models in lieu of prohibitively expensive Monte Carlo simulation (MCS). PCE relies on parametric distributions for associated variables and appropriate basis functions. However, incomplete or imperfect information on the stochastic variables can limit its use; accepted parametric forms for variable distributions, for instance, may not be justified when variables display multimodal character or mixed discrete-continuous support. Also, the dependency structure among the random variables may be complex, which can make probabilistic mapping or transformation to independent variables needed for PCE cumbersome. Nonlinearities in such transformations can affect the accuracy of PCE surrogate models and lead to slower convergence relative to truth system computations of desired QoIs (quantities of interest). To address these challenges, a distribution-free PCE approach is proposed. We compute joint raw moments of underlying random input variables for Gram-Schmidt orthogonalization in developing surrogate models. Using illustrative examples, we demonstrate the proposed approach as an efficient and accurate surrogate model-building alternative to traditional PCE.