You need to distinguish two different regions. As long as you are not too close to either black hole their gravitational potentials will just add. This is the linear region. Close to either black hole this linearity breaks down and we enter the non-linear region.
In the linear region the time dilation (relative to an observer at infinity) is well approximated by the weak field expression:
$$ t_r = \sqrt{1 +\frac{2\Phi}{c^2}} \tag{1} $$
where $\Phi$ is the Newtonian gravitational potential. In fact for a single black hole the Newtonian potential is:
$$ \Phi = -\frac{GM}{r} $$
and substituting this in equation (1) gives us the relativistic expression that you cite:
$$ t_r = \sqrt{1 - \frac{2GM}{c^2r}} \tag{1} $$
However be cautious about assigning too much significance to this as it's just a coincidence. The meaning of the variable $r$ is different in the Newtonian and relativistic theories.
Anyhow, provided we are in the linear region to calculate the time dilation for two (or more) black holes we need only calculate the total gravitational potential, and this is just the sum of the two separate potentials:
$$ \Phi_t = -\frac{GM_1}{r_1} - \frac{GM_2}{r_2} $$
and substituting this in equation (1) is going to give the expression that you suggest in your question:
$$ t_r = \sqrt{1 +\frac{2\Phi_t}{c^2}} = \sqrt{1 -\frac{2G}{c^2}\left( \frac{M_1}{r_1} + \frac{M_2}{r_2} \right)} \tag{2} $$
This works for the linear region, but there is still the non-linear region to consider i.e. close to one or both of the black holes. The trouble is that there is no analytic solution for the metric that describes two orbiting black holes. All we can say is that equation (2) will become increasingly inaccurate the closer we approach. We would need to do a numeric calculation to get an accurate result.