Is $D(x)$ here considered as a probability distribution?
Yes. In fact, it is the probability of label or class given an input feature vector or image, denoted as $p(c|\mathbf{x})$. Here are a few resources to better understand this:
What is the intuition behind the expected value in it?
Unfortunately, in my opinion, the notation is a bit misleading. In general, the cross entropy between two probability distributions $p$ and $q$ with the same support $\mathcal{X}$ is:
$$
H(p,q) = -\mathbb{E}_{p(x)}[\log{q(x)}] = -\sum_{x\in\mathcal{X}} p(x) \log{q(x)} \tag{1} \label{eq1}
$$
This means that the expression:
$$
-\frac{1}{N} \sum_{i=1}^{N}{ \; y_i \; log(D(x_i))} \tag{2} \label{eq2}
$$
Is just the cross entropy multiplied by the scalar $\frac{1}{N}$. I think this notation is misused because in the context of classification, both distributions $p$ and $q$ are discrete, and so the cross entropy is just a summation. Instead, I think the multiplication by $\frac{1}{N}$ is used to show that expression in \ref{eq2} is an approximation of the expectation in equation \ref{eq1} as the sample mean. However, as I said, all distributions are discrete, so the multiplication by $\frac{1}{N}$ is unnecessary. Also, since we want to minimize the cross entropy with respect to $q$ (see here and the links above for why), then it does not make a difference whether we multiply the cross entropy by $\frac{1}{N}$ or not.
UPDATE
It turns out that although the scalar $\frac{1}{N}$ does not serve a theoretical purpose, it serves a practical one. The scalar $\frac{1}{N}$ is used to normalize the negative log-likelihood (cross entropy) function to make its gradients smaller, and therefore allows larger step sizes to be chosen if an iterative optimization procedure, such as gradient descent, is used to minimize the cross entropy function. These larger step sizes are desirable because if backtracking line search is used to find the optimal step size, then it will be easier/quicker to find this optimal step size if it is larger.
For example, consider the gradient descent update rule:
$$
\theta_{k+1} \leftarrow \theta_{k} - \alpha \cdot \nabla_{\theta} f(x;\theta)
$$
Where $\alpha$ is the step size and $f(x;\theta)$ is some function parameterized by $\theta$. Suppose that $f$ was replaced with the expression in $\ref{eq2}$. Also, notice that in $\ref{eq2}$, if the factor $\frac{1}{N}$ was not there, then as $N$ increases, the value of the summation will also increase. The gradient of this summation would also be a function of $N$, such that as $N$ increases, the gradient increases in value too. This would mean that the steo size $\alpha$ would need to be smaller in order to achieve the same overall step in $\theta$ according to the gradient descent update rule.
However, if we multiply this summation by the factor $\frac{1}{N}$, then the gradient of the summation will also be a function of $\frac{1}{N}$, such that if $N$ increases, then the gradient may not increase as much as before. This will help keep gradients small such that the step size $\alpha$ does not need to be as small to achieve the same overall step according to the gradient descent rule.
As further motivation for the use of $\frac{1}{N}$, suppose that we start with $\theta_0=1$ and we know that the optimal $\theta$ is $\theta^*=0$. Also, suppose that $N$ is larger (e.g. $N=1,000$). Without $\frac{1}{N}$, the term $\nabla_{\theta} f(x;\theta)$ in the gradient descent update rule would be very large, and if $\alpha\approx 0.05$, then the new value of $\theta$ after only 1 iteration would be extremely small (even negative), clearly overshooting our target of $\theta^*=0$. If we re-introduce $\frac{1}{N}$ while keeping $\alpha$ the same, then the value of $\theta$ after 1 iteration will still be smaller, but not as much as before.
