This paper presents a robust and efficient algorithm for analyzing and restoring the observability of distribution grids. A linearized load flow model is used to characterize three-and multi-phase distribution grids under balanced and unbalanced conditions. Assuming that metering devices are installed to measure nodal magnitudes at some buses, the rank of the matrix associated with the non-metered magnitudes is determined by using a greedy strategy combined with a Gram-Schmidt orthonormalization procedure. Unlike other numerical-based techniques, the proposed algorithm avoids common errors caused by denominators close to zero. As an additional salient feature, a mathematically sound stopping criterion is provided to guarantee the robustness of the proposed algorithm. Furthermore, for unobservable grids, the algorithm identifies the linearly dependent rows in the aforementioned matrix, which correspond to the buses that should be metered to restore grid observability. The effective performance of the proposed algorithm is illustrated using several medium-and low-voltage benchmarks with up to 906 buses.