Magnetic resonance electrical impedance tomography (MREIT) is a non-invasive method to visualize cross-sectional electrical conductivity and/or current density by measuring a magnetic flux density signal when an electrical current is injected into a subject. In the MREIT system, it is crucial to reduce the scan duration while maintaining spatial resolution and contrast for practical in vivo implementation. The purpose of the study is to optimize the measured magnetic flux density using an interleaved multiple-echo pulse sequence (injected current nonlinear encoding) that acquires each spatial position multiple times, although these pixels vary between echoes in their signal-to-noise ratio due to (a) T*₂ decay and (b) the current density passing through the pixel. Using the acquired multiple measured magnetic flux densities, the noise level for the measured magnetic flux density B_z at each pixel is estimated using the relationship between the intensity of the magnitude and the width of the injected current. We determine an optimal combination of the multiply acquired magnetic flux densities, which optimally reduces the random noise artifacts. We develop a new denoising technique and apply it to a recovered conductivity distribution with a known noise level of the recovered magnetic flux density, which is designed to provide a stable internal conductivity distribution, while sustaining resolution. The proposed method uses three key steps: the first step is optimizing the magnetic flux density by using the multiple-echo magnetic flux densities at each pixel, the second step is estimating the noise level of this optimized magnetic flux density and the third step is applying a denoising technique using the pixel-specific estimated noise level. Numerical simulation experiments using a three-dimensional cylindrical phantom model validated the proposed method. Multiple-echo B_z data were generated, including in short T*₂ or low spin-density regions, as a function of T*₂ and the temporal extent of the injected current. In the simulation experiment, comparing between an equally averaged and the optimized B_z methods, relative L₂-mode errors were 0.053 and 0.024, respectively. In an actual imaging experiment of an agarose gel filled with objects of various conductivities and shapes, we acquired six echoes per repetition time. The optimal weighting factors minimized the effects of noise in B_z, and provided reconstructed conductivity maps that suppressed noise artifacts that normally accumulate in the low signal-to-noise-ratio defect regions.