I've never dealt with fixing that particular problem, but you could use some kind of slope-scaled depth threshold.
e.g. threshold /= abs(normal_vs.z)
If you don't have per-pixel normals in this pass, you can guesstimate them with ddx/ddy on the depth values (or just guesstimate the slope directly rather than a full normal).
Regarding linearization, I've usually done this in a pass beforehand so that it occurs once per pixel, not once per blur sample.
I do have access to per pixel normals in this case but I was hoping to avoid the extra texture sampling. Depending on how I implement the slope method It may be cheaper, I am only going along one axis anyway.
EDIT: Just realized that the slope method wouldn't work very well. To calculate slope I would essentially do what I have already calculating:
slope = 1.0 / abs(centerdepth - depth)
step(abs(centerdepth - depth), maxdiff * abs(centerdepth - depth));
which would always result in 0 because abs(centerdepth - depth) > 0.01 * abs(centerdepth - depth)