There is a lot of detail in the question and I am not sure of all the requirements and desired results that would affect the processing after the 16 KHz decimated samples are produced. However I can provide details and comments up to that point that may be helpful.
To answer first, there is nothing fundamentally incorrect in the approach the OP is taking: ensuring any aliasing bands are filtered out prior to down-sampling and this filtering is done in the case using biquad IIR filters. The use of biquads could be an advantage if delay is of upmost concern, otherwise I would change this to a linear phase multiband filter designed with the least squares algorithm. Both can be made to work and the important points detailed below are to understand where the rejection is needed and how much.
Decimation is just a form of sampling, and all considerations and requirements for aliasing apply, so in the interest of understanding decimation, I start with an understanding of aliasing due to the sampling process.
The following post demonstrates how without a band-selection filter, the analog signal when sampled at 10 Hz, could have been an 11 Hz sine wave or a 1 Hz sinewave. If both were present (and neither filtered out prior to sampling) the resulting sampled waveform would be the interference result of the two.
Sampling and the creation of multiple images at integer multiples of the sampling frequency
Further insights with an analogy to frequency mixing (for those familiar with that analog process): Higher order harmonics during sampling
Below shows the spectrums for a typical sampling process of a low pass continuous-time waveform. The top graph is the continuous-time waveform with spectral content limited to less than half the sampling rate. The middle graph is the spectrum of the sampling process: when "sampling", we multiply the time domain waveform with an impulse train spaced at the sampling period. The Fourier Transform of an impulse train in time is an impulse train in frequency, spaced by the sampling rate. Multiplication in time is convolution in frequency, and thus we get the bottom graph as the resulting discrete time waveform. Here, the frequency axis is extended to plus or minus infinity, but due the periodicity, we only need to show the discrete-time waveform's spectrum from $-f_s/2$ to $+f_s/2$ since it repeats periodically everywhere else. For purposes of understanding resampling operations and to cohesively work with combined analog and digital signal processing, I find it convenient to consider the "unrolled digital spectrum" extending to $\pm \infty$ as I have done here.
To note, if the continuous time waveform had spectral content limited to $f_s \pm f_s/2$ as I have depicted below (and yes if the waveform is complex, then it can have spectrum in the positive frequency axis alone- however getting into complex signal processing is not important to this answer), the exact same discrete time waveform will result!
The important item to see here is from where aliasing occurs specifically in the frequency domain, and how to optimize filtering to prevent it. The plot below shows all the frequency zones in red that would result in irreversible aliasing (without other information about the signal): any spectral content within these red zones would alias into the same fundamental blue zone in our discrete-time spectrum representing our waveform of interest (in this case the low pass waveform in the continuous-time domain that is in blue). Signal anywhere else is of less of a concern since it would alias to be outside of our frequency band of interest, and therefore we have opportunity to filter it later, downstream with very efficient and effective digital filtering. Anything that lands in band however is problematic.
This suggests an optimum filter would concentrate it's rejection at these alias bands specifically. In the continuous-time domain, it is often simpler to implement a low pass filter but additional nulls could be added at these frequency locations to improve rejection if needed. However, in the discrete-time domain, such multiband filters are easily synthesized and commonly used for resampling applications.
What is to be noted from the plots above, that the decimation approximation and the aliasing mechanism is equivalent! Below shows the same considerations when going from the higher sampling rate for 48 KHz and resampling to 16 KHz. The purple lines on the second graph show the passband and rejection bands for the optimized anti-alias filter (the filter like the other spectrums is periodic so it typically defined in the range of $\pm fs/2$ where we see there is a defined passband and define rejection band. The rest of the regions are "don't care"
Such multiband filters for resampling applications are easily designed as linear phase FIR filters using the least-squares algorithm, provided optimum rejection and passband performance (in a least-squares error sense) for a given filter complexity. Such filters are designed using Matlab, Octave and Python scipy.signal using the firls function. The amount of rejection required is based on distortion requirements and amount of signal present (including quantization noise) in each of the alias bands.
