Discrete implementations of integral operators are essentially matrix-vector multiplications. In this paper, we study the structure of Kirchhoff-type matrices where each row corresponds to a summation surface and each column corresponds to an impulse response. Due to irregularities in seismic coverage, the columns and rows are generally badly scaled. To balance the coefficients of the matrix, we propose two formulations for normalization: a data-space formulation based on row scaling of pull (sum) operators, and a model-space normalization based on column scaling of push (spray) operators. In both approaches, the final image is normalized by a reference model that is the operator's response to an input vector with all components equal to one. We apply this normalization technique to approximate a data covariance matrix based on the definition of a data-space pseudo inverse. This data covariance is an AMO matrix. It represents an equalization filter that corrects the imaging operator for the interdependencies among data parameters. We investigate the use of the normalization operators as preconditioners for the iterative solution of multichannel inversion. The diagonal transformation ensures common magnitude to all the variables and accelerates the convergence of the linear solver. Beyond the goal of fold normalization, the advantage of the iterative solution is to interpolate for missing data and reduce the artifacts related to data aliasing.