One option would be to compute the weights as the reciprocal of the variance of the voxel time series -- those voxels that vary a lot (from whatever cause) will get relatively little weight. For example:

3dAutomask -prefix q_mask -dilate 5 epi+orig
3dTstat -stdev -prefix q_sigma epi+orig
3dcalc -prefix q_weight -datum float -a q_mask+orig -b q_sigma+orig -expr 'a/(b*b)'
/bin/rm -f q_mask+orig.* q_sigma+orig.*

Then use q_weight+orig'[0]' as the weight volume. You will want to make sure that you do NOT use the initial transient part of the epi+orig time series (if it has such a transient) in the 3dTstat calculation.

I tried this on one dataset only, just now. The results don't differ markedly from the standard weight, but then this dataset doesn't have the problem to which you allude. At least the concept isn't grotesquely wrong.