  A significant fraction of observed galaxies in the Rubin Observatory Legacy
Survey of Space and Time (LSST) will overlap at least one other galaxy along
the same line of sight, in a so-called "blend." The current standard method of
assessing blend likelihood in LSST images relies on counting up the number of
intensity peaks in the smoothed image of a blend candidate, but the reliability
of this procedure has not yet been comprehensively studied. Here we construct a
realistic distribution of blended and unblended galaxies through high-fidelity
simulations of LSST-like images, and from this we examine the blend
classification accuracy of the standard peak-finding method. Furthermore, we
develop a novel Gaussian process blend classifier model, and show that this
classifier is competitive with both the peak-finding method as well as with a
convolutional neural network model. Finally, whereas the peak-finding method
does not naturally assign probabilities to its classification estimates, the
Gaussian process model does, and we show that the Gaussian process
classification probabilities are generally reliable.
