Why the imaginary part of green function is the spectral function? In the zero temperature Green's function, 
$$G_{\alpha\beta}(xt,x't')=-i\langle \Psi_H|\hat{T}[\psi_{H\alpha}(xt)\psi^{\dagger}_{H\beta}(x't')]|\Psi_H\rangle .$$
In Lehmann Representation, 
$$G_{\alpha\beta}(k,\omega)=\sum_n \left\{\frac{\langle \Psi_H|C_{S\alpha}(k)|n\rangle \langle n|C^{\dagger}_{S\beta}(k)]|\Psi_H\rangle }{\omega-(E^{N+1}_n-E^N_G)+i0^+}+\frac{\langle \Psi_H|C_{S\beta}^{\dagger}(k)|n\rangle \langle n|C^{\dagger}_{S\alpha}(k)]|\Psi_H\rangle }{\omega+(E^{N-1}_n-E^N_G)-i0^+}\right\}$$
So far so good.
Then I was told that the imaginary part of the Green's function is the spectral function and can be used to change Green function as:
$$G_{\alpha\beta}(k,\omega)=\int_0^{\infty}\left(\frac{A_{\alpha\beta}(k\omega')}{\omega-\mu-\omega'+i0^+}+\frac{B_{\alpha\beta}(k\omega')}{\omega-\mu+\omega'-i0^+}\right)\mathop{}\!\mathrm d\omega'$$
where 
\begin{align}
A_{\alpha\beta}(k\omega')&=\sum_n \langle\Psi_H|C_{S\alpha}(k)|n\rangle\langle n|C^{\dagger}_{S\beta}(k)]|\Psi_H\rangle\delta(\omega'+\mu-E_n^{N+1}+E_G^N)\\
B_{\alpha\beta}(k\omega')&=\sum_n \langle\Psi_H|C_{S\beta}^{\dagger}(k)|n\rangle\langle n|C_{S\alpha}(k)]|\Psi_H\rangle\delta(\omega'-\mu-E_n^{N+1}+E_G^N)
\end{align}
How is the third Green's function derived? I ask for detailed derivation process. Why is the imaginary part of the Green function the spectral function? The spectral function is for what? What is the relationship between the above thing and the Kramers-Kronig relationship?
 A: To take things in order

How is the third Green's function derived? I ask for detailed derivation process.

To reduce clutter and focus on ideas, I will drop a lot of subscripts and such. Note: the clutter I'm removing (e.g. $k,\alpha,\beta$) is still there, I'm just not writing it. I will also drop $\mu$ since I usually don't see it appearing the way you have it (see: https://en.wikipedia.org/wiki/Green%27s_function_(many-body_theory)#Zero-temperature_limit) and it can always be handled by absorbing it into the Fermi function or in shifting the energies. Begin with the second equation
$$G(\omega)=\sum_n \left\{\frac{\langle \psi|c|n\rangle \langle n|c^{\dagger}|\psi\rangle }{\omega-(E_n-E_m)+i0^+}+\frac{\langle \psi|c^\dagger|n\rangle \langle n|c|\psi\rangle }{\omega+(E_n-E_m)-i0^+}\right\}$$
Insert $1=\int d\omega'\delta(\omega'-E_n+E_m)$ to get
$$
G(\omega)=\int d\omega'\sum_n \delta(\omega'-E_n+E_m)\left\{\frac{\langle \psi|c|n\rangle \langle n|c^{\dagger}|\psi\rangle }{\omega-(E_n-E_m)+i0^+}+\frac{\langle \psi|c^\dagger|n\rangle \langle n|c|\psi\rangle }{\omega+(E_n-E_m)-i0^+}\right\}
$$
The delta function enforces $\omega'=E_n-E_m$ so we can substitute this into the expression. Note also that now the numerators in the integral are $A(\omega')$ and $B(\omega')$ respectively. This gives
$$
G(\omega)=\sum_n \left\{\int d\omega'\frac{A(\omega')}{\omega-\omega'+i0^+}+\int d\omega'\frac{B(\omega')}{\omega+\omega'-i0^+}\right\}
$$
Which is precisely your third expression.

Why is the imaginary part of the Green function the spectral function?

Well let's take the imaginary part of $G$.
$$\Im G(\omega) = \Im \int d\omega'\sum_n \left\{\frac{A(\omega')}{\omega-\omega'+i0^+}+\frac{B(\omega')}{\omega+\omega'-i0^+}\right\}
$$
For compactness let's drop the sum on $n$ and the $B(\omega')$ term - they all work the same way. To make the next part clearer, let's also replace $0^+$ by $\lim\limits_{\varepsilon\to0^+}\varepsilon$. So we are trying to analyze
$$
\Im \lim\limits_{\varepsilon\to0^+}\int d\omega'\frac{A(\omega')}{\omega-\omega'+i\varepsilon}
$$
Multiply the numerator and denominator the by $\omega-\omega'-i\varepsilon$ to get
$$
\Im \lim\limits_{\varepsilon\to0^+}\int d\omega'\frac{A(\omega')(\omega-\omega'-i\varepsilon)}{(\omega-\omega')^2+\varepsilon^2}=
$$
$$\Im \lim\limits_{\varepsilon\to0^+}\int d\omega'A(\omega')\left[\frac{-i\varepsilon}{(\omega-\omega')^2+\varepsilon^2}+\frac{\omega-\omega'}{(\omega-\omega')^2+\varepsilon^2}\right]
$$
(Note: taking the limit here gives a version of a Kramers–Kronig relation you ask about)
$A(\omega)$ is real and non-negative for Fermions (exercise left for reader) so the second term has no imaginary part and can be dropped. The first term is $-i\pi$ times a nascent delta function (specifically the Poisson kernel). This allows us to rewrite the above as
$$
\Im \lim\limits_{\varepsilon\to0^+}\int d\omega'A(\omega')\frac{-i\varepsilon}{(\omega-\omega')^2+\varepsilon^2}=\Im -i\pi\int d\omega'A(\omega')\delta(\omega-\omega')=-\pi A(\omega)
$$ 
Which gives us the known relation $-\frac{1}{\pi}\Im G(\omega,\vec{k})=A(\omega,\vec{k})$ (other conventions may give $-2$ instead of $-1/\pi$). I've brought back the $\vec{k}$ to emphasize that this doesn't hold for $G(\omega)=\sum_{\vec{k}} G(\omega,\vec{k})$ and $A(\omega)=\sum_{\vec{k}} A(\omega,\vec{k})$

The spectral function is for what?

The spectral function gives information about the excitations (particles) in your system. As you can see from the form of the spectral function $A(\omega,k)=\sum_n \delta(\omega-E_n+E_m)\langle \psi|c_k|n\rangle \langle n|c_k^{\dagger}|\psi\rangle $ the spectral function will be peaked at the energies and wavevectors associated with the species in your system. As we showed above, the spectral function allows us to get the Green's function. It can be used to get the filling of the system and information about the density of states. (Note that this applies to noninteracting systems which is what you've presented here. The interacting case is more subtle but the idea holds that the spectral function tells you about the excitations in your system).

What is the relationship between the above thing and the Kramers-Kronig relationship?

There are a few ways to understand the Kramers-Kronig relations. One is that given a function that's analytic in UHP and decays at least as fast as $1/|\omega|$ for large $|\omega|$) then the Kramers-Kronig relations allow you to reconstruct the whole function from knowing just the real part or just the imaginary part. This is essentially how we can use the spectral function (the imaginary part of $G$) to get the full $G$! This is why the third equation makes sense and how we show that the spectral function is (up to a constant) the imaginary part of the Green's function. These relations are a special case of the Sokhotski–Plemelj (https://en.wikipedia.org/wiki/Sokhotski%E2%80%93Plemelj_theorem) and is a useful thing to know for anyone working in condensed matter physics.
Hope that clears things up, cheers!
