There is a trigonometric identity:
$$\cos^2A+\cos^2B+\cos^2C+2\cos A\cos B\cos C\equiv 1\text{ when }A+B+C=\pi$$
It is easy to prove it in an algebraic way, just like that:
$\quad\cos^2A+\cos^2B+\cos^2C+2\cos A\cos B\cos C\\=\cos^2A+\cos^2B+\cos^2\left(\pi-A-B\right)+2\cos A\cos B\cos \left(\pi-A-B\right)\\=\cos^2A+\cos^2B+\cos^2\left(A+B\right)-2\cos A\cos B\cos \left(A+B\right)\\=\cos^2A+\cos^2B+\left(\cos A\cos B-\sin A\sin B\right)^2-2\cos A\cos B\left(\cos A\cos B-\sin A\sin B\right)\\=\cos^2A+\cos^2B+\cos^2A\cos^2B+\sin^2A\sin^2B-2\sin A\cos A\sin B\cos B-2\cos^2A\cos^2B+2\sin A\cos A\sin B\cos B\\=\cos^2A+\cos^2B-\cos^2A\cos^2B+\left(1-\cos^2A\right)\left(1-\cos^2B\right)\\=\cos^2A+\cos^2B-\cos^2A\cos^2B+1-\cos^2A-\cos^2B+\cos^2A\cos^2B\\=1$
Then, I want to find a geometric way to prove this identity, as $A+B+C=\pi$ and it makes me think of the angle sum of triangle. However, it is quite hard to prove it in a geometric way. Therefore, I hope there is someone who can help. Thank you!

Since the accent in the OP is put on a purely geometric solution, i can not even consider the chance to write $\cos^2 =1-\sin^2$, and rephrase the wanted equality, thus having a trigonometric function which is better suited to geometrical interpretations.
So this answer has two steps, first we reformulate the given identity in a mot-a-mot geometric manner, the geometric framework is introduced, some strictly geometrically transposed equivalent relations are listed, then we give a proof:
Proof: We have: $$ \sin \hat B =\sin \widehat{C'HA} =\frac{C'A}{AH} =\frac{AC\;\cos A}{AH} =\frac{2R\sin B\; \cos A}{AH} =\frac{\sin B\; \cos A}{AH} \ , $$ which implies $AH=\cos A$, and the similar relations. Then we express twice the area of $\Delta HBC$ as $$ HA'\cdot BC =2[HBC]=HB\cdot HC\cdot \sin\widehat{BHC}\ ,$$ thus getting $HA'=\cos B\cos C$.
We are in position to give a geometric mask to the given equality:
Proof: The relations above are equivalent:
$AH^2=4A^*H^2=4OM_A^2$, and $2 HA\cdot HA'=4 HA^*\cdot HA'$ is the power of $H$ in the circle $(N)$, so it can be rewritten using its radius $NA^*=\frac 12 R$ and the distance to its center, $NH=\frac 12 OH$ as $2 HA\cdot HA'=4 HA^*\cdot HA'=R^2-OH^2$.
From the triangle $OBM_A$, $4OM_A^2+BC^2 =4(OM_A^2+BM_A^2)=4OB^2=4R^2$.
Note that $G$ cuts the median $AM_A$ in the proportion $AG:GM_A=2:1$, so it projects on $BC$ in the same proportion. This also holds for the colinear points $H,G,O$, so $HG:GO=2:1$, so $HO=3GO$.
The last relation, $OG^2 = R^2-\frac 13(a^2+b^2+c^2)$, is a standard formula. We have in general the formula for an arbitrary point $P$: $$PA^2+PB^2+PC^2=GA^2+GB^2+GC^2+3GP^2\ .$$ We apply it for $P=O$, getting $3R^2=3OG^3+\sum AG^2=3OG^3+\frac 49\sum AM_A^2=3OG^3+\frac 49\sum \left(\frac 12b^2+\frac 12 c^2-\frac 14 a^2\right)=3OG^3+\frac 49\sum \frac 34a^2=3OG^3+\frac 13\sum a^2\ .$
$\square$