Allow SpectralRegion to work on spectral axes that are arbitrarily ordered

#911 updated `SpectralRegion` to work for cases where the spectral axis is either increasing or decreasing. This covers most "in the wild" spectra, I believe. However, there are some edge cases where spectral axes are not monotonic. The most practical examples I can think of are where multiple spectra are mashed together into a spectral energy distribution, or e.g. echelle spectra where some of the pixels from adjacent orders overlap in spectral axis space.

Of course one can always manually make a spectrum monotonic be re-sorting the pixels. But it would be mildly more convenient if `SpectralRegion` did that instead of forcing the user to do that. https://github.com/astropy/specutils/pull/911#discussion_r789923462 has a rough sketch how that might be done.

Start by reading the SpectralRegion implementation and the approach sketched in the discussion on PR #911. Reproduce cases with non-monotonic spectral axes, including concatenated spectra and overlapping echelle orders, then verify that SpectralRegion handles them while retaining support for increasing and decreasing axes.

Written by the indexing model from the issue text.