I'd like to give a partial answer (or a partial failure, depending on perspective:)) to my own question. The strategy used for the homotopy classification is derived from the classic paper by Avron, Seiler and Simon, discussing the homotopy perspective of IQHE, and the method is very succinctly described in the lecture notes by Moore(see discussions around equation (39)). Let me summarize it here, following the lines by Moore:

Mathematically, to classify n-band IQHE systems, we are homotopically classifying (gapped) maps from $T^2$ to $H_n$, where $H_n$ is the space of $n\times n$ Hermitian matrices(i.e. Hamiltonians). We can diagonalize the Hamiltonian $H$ using unitary matrices, i.e. $H=UDU^\dagger$. Since the spectrum of the Hamiltonian is gapped, we can homotopically flatten and shift the bands within $D$, keeping $U$ fixed, without band crossing, to a standard matrix, say diag$[1,2,\ldots, n]$. Hence, all the information about homotopy class should be encoded in the unitary matrices $U$. However, we can always right-multiply $U$ by any diagonal unitary matrix without changing $H$, hence such equivalence must be quotiented(not proper English, I know) away. In the end, we need to homotopically classify maps $T^2\to U(n)/\sim$. Such maps are precisely classified by n-1 Chern numbers.

Now come to my own problem, I slightly enlarge the space of maps to make standard techniques applicable, namely, instead of the specific maps in the form of (1) in the main post, I consider the maps of the form

$$U(k)=\begin{bmatrix} 0&U_1(k)\\ \\U_2(k)&0\end{bmatrix}\cdots\cdots (2),$$

with the only requirement that $U(k)$ is unitary. Now as a spoiler to the readers not patient enough for the would-be longwinded discussion, there exist nontrivial homotopy classification for such maps, and the classification is strictly finer than a Chern-number classification, however the classification is not good enough to distinguish the topological phase established numerically in my collegues' work.

Using the observation made in the main post that $\gamma_5 U=-U\gamma_5$, it is easily derived that $U$ can be diagonalized as

$$U=\begin{bmatrix}a_1&b_1&a_1&b_1\\a_2&b_2&a_2&b_2\\a_3&b_3&-a_3&-b_3\\a_4&b_4&-a_4&-b_4 \end{bmatrix}\begin{bmatrix}d_1& & & \\ &d_2& & \\ & & -d_1 & \\ & & & -d_2 \end{bmatrix}\begin{bmatrix}a_1&b_1&a_1&b_1\\a_2&b_2&a_2&b_2\\a_3&b_3&-a_3&-b_3\\a_4&b_4&-a_4&-b_4 \end{bmatrix}^{\dagger},$$

By the same reasoning as outlined in Moore's notes, the diagonal matrix diag$[d_1,d_2,-d_1,-d_2]$ can be deformed to a standard one, say diag$[1,i,-1,-i]$. Now we only need to worry about the similarity transformations. Note our similarity transformations have the form

$$\begin{bmatrix}A&A\\B&-B\end{bmatrix}\cdots\cdots(3),$$

where

$$A=\begin{bmatrix}a_1&b_1\\a_2&b_2\end{bmatrix}, \ B=\begin{bmatrix}a_3&b_3\\a_4&b_4\end{bmatrix}.$$

Hence A and B are independent matrices, satisfying $AA^\dagger=BB^\dagger=\frac{1}{2}I$ because of the unitarity of the similarity transformation. Also, the right-multiplications on (3) allowed are of the form

$$\begin{bmatrix}\Lambda& \\ &\Lambda\end{bmatrix},$$

to maintain the form (2), where $\Lambda$ is any 2 by 2 unitary diagonal matrix.

Therefore, we need to deal with maps $T^2\to M$, where $M$ is the quotient space of the pair of $2\times 2$ matrices $(A,B)$, with the equivalence relation $(A,B)\sim(A\Lambda,B\Lambda)$, that is, $M=(U(2)\times U(2))/\sim$. In fact, we can perform a diffeomorphism before doing quotient: $(A,B)\mapsto (A,BA^{-1})$, then the equivalence relation is simplified to $(A,BA^{-1})\sim (A\Lambda,BA^{-1})$, and so $M\approx (U(2)/\sim)\times U(2)$. (I learnt this trick from this mathoverflow post)

In the end, we have the homotopy classes $[T^2, (U(2)/\sim)\times U(2)]=[T^2, U(2)/\sim]\times [T^2, U(2)]$, the former is characterized by a Chern number as usual and the latter is characterized by two winding numbers. Therefore, in this setup, we have a $\mathbb{Z}^3$ classification scheme, and it is easy to prove such classification is finer than the Chern number classification implemented on (2). The honeycomb lattice has more or less the same mathematical structure, and it has a $\mathbb{Z}^5$(1 Chern number and two pairs of winding numbers) classification since the size of the matrix is $6\times 6$.

As I said, this scheme is not able to distinguish some topological phase in my college's work(the AFI phase mentioned in their paper). It is not very hard to understand the failure: In my collegues' work, the ring resonator system posses edge states which are robust against the missing of rings, whereas the constraint I imposed by $\gamma_5$ (which in position space represent a $\pi$ phase delay on 2 of the 4 amplitudes in a ring) is not, hence I'm probably barking at the wrong "symmetry" constraint. It seems the symmetry constraint of ring resonator model is less straightfoward to find than that of electronic condensed matter system, since it is from the beginning constructed on a very phenomenlogical level, while in electronic condensed matter system one can think of the fundamental symmetries that we are used to, like time reversal etc.