One poor-man's approach to deblending multiple bands and visits is
- Deblend in one band (or a combination, e.g. chi^2 band)
- Fit models in that band
- In all bands separately fit those models simultaneously to all object in the blend, allowing only a minimum of coefficients to float (ideally only amplitudes, but given the realities of real galaxies the bulge and disk will need to both be fit, or the Sloan Swindle components, or the Lensfit Swindle, or ...). The correct PSF for each image would be used.
It will probably be necessary to include undetected objects in this fit (cf. HSC-1023), using point-source models.
It is not clear how much of a poor-man's solution this actually is, or whether it will work quite well.