Write $\beta = \lambda_n$, the top row of your GT patterns.  It's a theorem of [Baryshnikov] that if we choose a uniformly random point in the polytope GT${}_\lambda$, it's equivalent to choosing a Haar-random Hermitian matrix with spectrum $\beta$ and then taking its "principal minors".  (I've also seen this fact credited to Weyl, and others.)  More precisely, let $B = \mathrm{diag}(\beta)$, and form a matrix $X = U B U^\dagger$, where $U$ is a random unitary.  Then let $\lambda_{11}$ be the top-left entry of $X$, let $\lambda_{21}, \lambda_{22}$ be the eigenvalues of the top-left $2 \times 2$ submatrix of $X$, ..., and let $\lambda_{n1}, \dots, \lambda_{nn}$ be the eigenvalues of the top-left $n \times n $ submatrix of $X$ (namely, $\beta$).  Then $\lambda$ is uniformly random in the polytope GT${}_\lambda$.

This probability distribution on $\lambda$ is basically your integral, but we have to divide by the volume of the polytope, which is $V(\lambda)/[(n-1)! (n-2)! \cdots 2! 1!]$.  I guess this is standard?  If not, it's also in Baryshnikov.  

Having done so, your identity is the Harish-Chandra--Itzykson-Zuber identity applied to the matrices $A = \mathrm{diag}(\alpha)$ and $B$.  This follows by inferring the diagonal entries of $X$ from the Gelfand-Tsetlin pattern $\lambda$, which you can do because the Gelfand-Tsetlin pattern gives you the traces of all the top-left submatrices.

(By the way, I think the [Faraut] paper referenced below has a good exposition of some related things.)

<cite authors="Baryshnikov, Yu.">_Baryshnikov, Yu._, [**GUEs and queues**](http://dx.doi.org/10.1007/s004400000101), Probab. Theory Relat. Fields 119, No.2, 256-274 (2001). [ZBL0980.60042](https://zbmath.org/?q=an:0980.60042).</cite>

<cite authors="Faraut, Jacques">_Faraut, Jacques_, [**Rayleigh theorem, projection of orbital measures and spline functions**](http://dx.doi.org/10.1515/apam-2015-5012), Adv. Pure Appl. Math. 6, No. 4, 261-283 (2015). [ZBL1326.15058](https://zbmath.org/?q=an:1326.15058).</cite>