Given a set $\mathcal{C}=\{C_i\}_{i=1}^m$ of square matrices, the matrix blind joint block diagonalization problem (BJBDP) is to find a full column rank matrix $A$ such that $C_i=A\Sigma_iA^\text{T}$ for all $i$, where $\Sigma_i$'s are all block diagonal matrices with as many diagonal blocks as possible. The BJBDP plays an important role in independent subspace analysis (ISA). This paper considers the identification problem for BJBDP, that is, under what conditions and by what means, we can identify the diagonalizer $A$ and the block diagonal structure of $\Sigma_i$, especially when there is noise in $C_i$'s. In this paper, we propose a ``bi-block diagonalization'' method to solve BJBDP, and establish sufficient conditions under which the method is able to accomplish the task. Numerical simulations validate our theoretical results. To the best of the authors' knowledge, existing numerical methods for BJBDP have no theoretical guarantees for the identification of the exact solution, whereas our method does.