  Line intensity mapping (LIM) has emerged as a promising tool for probing the
3D large-scale structure through the aggregate emission of spectral lines. The
presence of interloper lines poses a crucial challenge in extracting the signal
from the target line in LIM. In this work, we introduce a novel method for LIM
analysis that simultaneously extracts line signals from multiple spectral
lines, utilizing the covariance of native LIM data elements defined in the
spectral--angular space. We leverage correlated information from different
lines to perform joint inference on all lines simultaneously, employing a
Bayesian analysis framework. We present the formalism, demonstrate our
technique with a mock survey setup resembling the SPHEREx deep field
observation, and consider four spectral lines within the SPHEREx spectral
coverage in the near infrared: H$\alpha$, $[$\ion{O}{3}$]$, H$\beta$, and
$[$\ion{O}{2}$]$. We demonstrate that our method can extract the power spectrum
of all four lines at the $\gtrsim 10\sigma$ level at $z<2$. For the brightest
line, H$\alpha$, the $10\sigma$ sensitivity can be achieved out to $z\sim3$.
Our technique offers a flexible framework for LIM analysis, enabling
simultaneous inference of signals from multiple line emissions while
accommodating diverse modeling constraints and parameterizations.
