Your problem is you have a multiple input multiple output system which couples the elements of your sweep, thus we indeed expect the edge effects of your zero-phase digital filter to work their way inwards. You can consider that at each sweep you get a vector of inputs, rather than a timeseries of data, with the transfer function of your system being:
$$\mathbf{y}=\mathbf{x}-\left(\frac{1}{1-z^{-1}}\mathbf{I}\right)\mathbf{H}\left(z^{-1}\mathbf{I}\right)\mathbf{y}=\left(\mathbf{I}+\left(\frac{z^{-1}}{1-z^{-1}}\mathbf{I}\right)\mathbf{H}\right)^{-1}\mathbf{x}$$
Where lower case boldface are vectors, and upper case boldface are matrices ($\mathbf{I}$ being the identify matrix). That matrix $\mathbf{H}$ (your FIR filter) is coupling everything, which is why it's expected that the edge effects will work their way inwards. To assess stability, you'd have to solve:
$$\mathbf{det}\left(\mathbf{I}+\left(\frac{z^{-1}}{1-z^{-1}}\mathbf{I}\right)\mathbf{H}\right)=0$$
This is much more complex than things need to be.

As a solution, redesign your control loop as per my figure below, which will theoretically guarantee stability and prevent the edge effects creeping inwards (I say theoretically because this analysis assumes there is no inaccuracy in your system, which would be fine if everything was done digitally, but you're converting the correction to an analog signal to subtract it there, albeit the analysis could be extended to include a disturbance term relatively trivially):
$$\mathbf{e}=\mathbf{H}\left(L[z]z^{-1}\mathbf{I}\right)(\mathbf{e}+(\mathbf{x}-\mathbf{e}))=\mathbf{H}\left(L[z]z^{-1}\mathbf{I}\right)\mathbf{x}$$

$$\mathbf{y}=\mathbf{x}-\mathbf{e}=\left(\mathbf{I}-\mathbf{H}\left(L[z]z^{-1}\mathbf{I}\right)\right)\mathbf{x}$$
Where $L$ is the transfer function of a low pass filter to provide some averaging of uncorrected sweeps (probably a moving average filter will do, but I leave that to you). This system clearly has no poles, except in the low pass filter, $L$, so design a stable low pass filter, and you're guaranteed a stable system. Hopefully that final transfer function is quite intuitive in that it takes $\mathbf{x}$ and subtracts the component you don't want to give you $\mathbf{y}$. Note that the role of the low pass filter is just to give some averaging between sweeps to attenuate their noise, which I guess you are aware of that you need insofar that one sweep of data will not provide an accurate correction, which is why I you used that integrator. On the first iteration, I would initialise the low pass filter so that its output immediately is equal to the first sweep signal, like what your integrator did.

[![system redesign][1]][1]


  [1]: https://i.sstatic.net/xVU9ZGHi.png