Super-resolution ultrasound localization microscopy is an advanced microvessel imaging modality. Multiple studies have proved its significant preclinical and clinical application potential. As an important step of the super-resolution ultrasound localization microscopy, current microbubble localization methods need spatially isolated microbubble signal and overlapping microbubble signals are rejected or treated as single isolated microbubbles. Overlapping microbubbles would convey equally important information as isolated microbubbles. As such, conventional ultrasound localization microscopy techniques use reduced microbubble concentration to facilitate isolated microbubble signals. This would require multiple microbubble injections to construct robust super-resolution images and localization error for overlapping microbubbles, which elongates the overall data acquisition time for ultrasound localization microscopy. In this study, the multivariate Gaussian fitting (MGF) method is applied to achieve localization of spatially overlapping microbubble signals. The MGF-based localization technique provides a customizable size and shape for both isolated and overlapping microbubbles. This approach allows the use of high concentration microbubble injections which effectively reduces the data acquisition time and improves the temporal resolution of ultrasound localization microscopy. Simulation studies were designed and demonstrated improved microbubble localization accuracy with the MGF-based localization technique. In vivo chick embryo microvessel studies showed that the proposed method improved the microvessel image quality with more accurate microvessel depiction.