This paper is devoted to the study of estimates of heat kernel of generalized Schrödinger operators and the related Hardy spaces. Let \(\mathcal{L}=-\Delta+\mu\) be the generalized Schrödinger operator on \(\mathbb{R}^n\) with \(n\geq 3\), where \(\mu\not\equiv0\) is a nonnegative Radon measure satisfying the scale-invariant Kato condition: for all \(x\in\mathbb{R}^n\) and \(0<r<R<\infty\), \[ \mu(B(x,r))\leq C_0(\frac rR)^{n-2+\delta}\mu(B(x,R)), \] and also the doubling condition: for all \(x\in\mathbb{R}^n\) and \(0<r< \infty\), \[ \mu(B(x,2r))\leq C_1\{\mu(B(x,r))+r^{n-2}\}, \] where \(C_0,C_1,\delta\) are constants independent of \(x, r\) and \(R\). In this paper, the author prove that the heat kernel \(\mathcal{K}_t\) associated with \(\mathcal{L}\) satisfies the following upper estimate: \[ 0\leq \mathcal{K}_t(x,y)\leq Ch_t(x-y)e^{-\varepsilon d_\mu(x,y,t)}, \] for some \(\varepsilon>0\), where \(h_t(x)=(4\pi t)^{-n/2}e^{-|x|^2/(2t)}\) and \(d_\mu(x,y,t)\) is some parabolic type distance function with respect to \(\mu\). As applications, a Hardy space with respect to \(\mathcal{L}\) is introduced via maximal functions of heat semigroup \(e^{-t\mathcal{L}}\) and characterized via atoms and Riesz transforms. The dual space of this Hardy space is also obtained.