Support for averaging spectra using masks; proposed API

I'm attaching one notebook examining a few related topics, so please pardon me if this should be multiple tickets.

I worked through the documentation example for [Splicing/Combining Multiple Spectra](https://specutils.readthedocs.io/en/stable/manipulation.html#splicing-combining-multiple-spectra), and found I didn't initially appreciate that the example applies only to non-overlapping spectra. An important set of use cases (previously in the domain of IRAF's `stsdas.hst_calib.ctools.splice` function) is the appropriate handling of overlapping spectra during combination. Here I examine overlapping echelle orders, but this is more generally applicable.

In the attached notebook, I attempted to manually implement an error-weighted averaging to produce an output `Spectrum1D` object with correct flux and uncertainty. **I scope out a method interface to `SpectrumCollection` and `SpectrumList` to do this within the machinery of specutils.**

The rationale for natively supporting averaging is:
- Users won't have to recreate these algorithms for each application (and sometimes get them wrong).
- A workflow of nested averaging allows for a wide range of applications with a minimal core of code. A user could create a collection and collapse it down to a `Spectrum1D` object by a method call with different algorithm options. These resulting `Spectrum1D` objects can themselves be built up into further collections with possibly different options to reduce a hierarchy of data. E.g.:
- combining mostly-overlapping and/or dithered data taken in the same mode
- combining partially-overlapping & disjoint echelle orders
- combining neighboring cenwaves (tilts of a grating)
- combining neighboring gratings
- combining spectra from multiple instruments (with some support for mask definition dissimilarities)

As you work down this hierarchy, the differences in spectral binning become larger. Being able to specify parameters related to the output grid is a must!

Some major issues encountered:
- Resampling (with `FluxConservingResampler` here) does not preserve mask values. Applicable masks are needed for a later averaging step.
- `Spectrum1D` addition has some important distinctions from averaging:
- Adding non-overlapping regions destroys any useful information about the uncertainty there.
- Underlying data should typically only be averaged where the mask is `0`, masks should otherwise be bitwise-or combinations of data actually used. (More details are spelled out in the notebook.)
- Support for `SpectrumList` to be user-facing is needed for combining spectra of different dimensions.

Some ancillary issues encountered:
- `FluxConservingResampler(extrapolation_treatment='zero_fill')` in v1.1 produces some nans at spectrum edges.
There are fewer than in v1.0. I haven't tested with other Resampler objects.
- Resampled objects should not always give uncertainties as `InverseVariance`.
- The representation of [`SpectrumCollection.spectral_axis`, `SpectrumCollection.flux`, ...] appear to be transposed from expected shape.
- _Documentation:_ Does specifying the `Spectrum1D` wavelength data **and** the WCS give precedence to the WCS?
(My instrument's WCS is less accurate, and it would be good to know if this is the case.)

**Attached notebook, HTML rendering, and data:**
[specutils_averaging_notebook.zip](https://github.com/astropy/specutils/files/5394200/specutils_averaging_notebook.zip)

Thank you all for your hard work on this project! I can see this being very useful for HST/STIS users shortly.

Start with the attached averaging notebook, then inspect the SpectrumCollection and SpectrumList entry points and the FluxConservingResampler behavior it discusses. Establish the scope of mask-aware, error-weighted averaging and output-grid options before implementation; done should include a user-facing combination workflow that preserves appropriate flux, uncertainty, and mask information.

Written by the indexing model from the issue text.