In a strongly nonlinear system the particle distribution in the phase space may develop long tails which contribution to the covariance (sigma) matrix should be suppressed for a correct estimate of the beam emittance. A method is offered based on Gaussian approximation of the original particle distribution in the phase space (Klimontovich distribution) which leads to an equation for the sigma matrix which provides efficient suppression of the tails and cannot be obtained by introducing weights. This equation is easily solved by iterations in the multi-dimensional case. It is also shown how the eigen-emittances and coupled optics functions can be retrieved from the sigma matrix in a strongly coupled system. Finally, the developed algorithm is applied to 6D ionization cooling of muons in HFOFO channel.