Extending the notion of second-order correlations, we define the cumulants of stationary non-Gaussian random fields, and demonstrate their potential for modeling and reconstruction of multidimensional signals and systems. Cumulants and their Fourier transforms called polyspectra preserve complete amplitude and phase information of a multidimensional linear process, even when it is corrupted by additive colored Gaussian noise of unknown covariance function. Relying on this property, phase reconstruction algorithms are developed using polyspectra, which can be computed via a 2-D FFT-based algorithm. Additionally, consistent ARMA parameter estimators are derived for identification of linear space-invariant multidimensional models which are driven by unobservable, i.i.d., non-Gaussian random fields. Contrary to autocorrelation based multidimensional modeling approaches, when cumulants are employed, the ARMA model is allowed to be non-minimum phase, asymmetric non-causal or non-separable.