Suppose that $N = 3$ and $d = 2$ (so that $\Sigma$ is a surface in $\mathbb R^3$), and suppose that $$f: R^3 \to R$$ is $C^\infty$, and that $\operatorname{grad} f $ does not vanish on $\Sigma =f^{-1}(0)$. Then the normal vector $$\vec n(X) = \frac{\operatorname{grad} f(X)}{|\operatorname{grad} f(X)|}$$ can be defined everywhere in a neighbourhood of $\Sigma$, and one can take $$P(X) = I-\vec n(X)\cdot \vec n(X)^T,$$ where $P(X)$ is the orthogonal projection on to the tangent space $T_p \Sigma$.
How to prove this result? I know that the projection of the vector onto the subspace is given by $P=A(A^TA)^{-1}A^T$.
If I understand correctly, you want to prove that P is the orthogonal projection onto $T_p \Sigma$. Since both of these are linear maps, you only need to check that they agree on a basis. So it suffices to check it on two vectors spanning $T_p \Sigma$ and on $n(p)$. The first two should be preserved and the latter should be mapped to zero. Your formula does exactly that.