To save energy and alleviate interference, connected dominating set (CDS) was proposed to serve as a virtual backbone of wireless sensor networks (WSNs). Because sensor nodes may fail due to accidental damages or energy depletion, it is desirable to construct a fault tolerant virtual backbone with high redundancy in both coverage and connectivity. This can be modeled as a k-connected m-fold dominating set (abbreviated as (k, m)-CDS) problem. A node set C ⊆ V (G) is a (k, m)-CDS of graph G if every node in V(G)\C is adjacent with at least m nodes in C and the subgraph of G induced by C is k-connected. Constant approximation algorithm is known for (3, m)-CDS in unit disk graph, which models homogeneous WSNs. In this paper, we present the first performance guaranteed approximation algorithm for (3, m)-CDS in a heterogeneous WSN. In fact, our performance ratio is valid for any topology. The performance ratio is at most γ, where γ = α + 8 + 2 ln(2α - 6) for α ≥ 4 and γ = 3α +2 ln 2 for α <; 4, and α is the performance ratio for the minimum (2, m)-CDS problem. Using currently best known value of α, the performance ratio is ln δ +o(ln δ), where δ is the maximum degree of the graph, which is asymptotically best possible in view of the non-approximability of the problem. Applying our algorithm on a unit disk graph, the performance ratio is less than 27, improving previous ratio 62.3 by a large amount for the (3, m)-CDS problem on a unit disk graph.