How to explain these results of integration of DiracDelta?

Here a way to get the correct symbolic result in Mathematica

Integrate[DiracDelta[1 - x^2 - y^2], {x, -1 - delta,1 + delta}, {y, -1 - delta, 1 + delta}
, Assumptions -> delta > 0]
(*Pi*)

numerical verification

eps = 0.005;
delta = .1
NIntegrate[1/Sqrt[2 Pi eps] Exp[-((1 - x^2 - y^2)^2/(2 eps))], {x, -1 - delta,1 + delta}, {y, -1 - delta, 1 + delta} ]
(*3.14091*)

A very simple limiting representation of the Dirac $delta$-function is

f[ε_, x_] = Piecewise[{{1/(2*ε), -ε < x < ε}}];

with $lim_{epsilonto0^+} f_{epsilon}(x)=delta(x)$ in the sense of a distribution.

The integral in question is then

Assuming[0 < ε < 1/2, 
  Integrate[f[ε, 1-x^2-y^2], {x, -1, 1}, {y, -1, 1}] // FullSimplify]

(*    (2 Sqrt[ε] + (1+ε)*ArcCot[Sqrt[ε]] + (ε-1)*ArcCsc[Sqrt[1/ε-1]] +
      (ε-1)*ArcTan[Sqrt[1/ε-2]] - ArcTan[Sqrt[ε]] - ε*ArcTan[Sqrt[ε]])/ε    *)

The limit is $pi$:

Limit[%, ε -> 0, Direction -> "FromAbove"]
(*    π    *)

The Mathematica expression

Integrate[
  DiracDelta[(x^2 + y^2) - 1], {x, -1, 1}, {y, -1, 1}]

evaluates (wrongly) to Pi/2 since your integration limits coincide with (x, y) points such that the argument of DiracDelta is equal zero (correct value: Pi, since we integrate over unit circle).

Thus, increasing the integration limit solves the problem:

Integrate[
  DiracDelta[(x^2 + y^2) - 1], {x, -2, 2}, {y, -2, 2}]

evaluates to Pi.

On the other hand, also integration over the upper half-plane yields the correct result:

Integrate[
  DiracDelta[(x^2 + y^2) - 1], {x, -1, 1}, {y, 0, 1}]

evaluates correctly to Pi/2.

Furthermore, by changing from cartesian to polar coordinates, Mathematica gives you a hint on why it's acting strange:

Integrate[r*DiracDelta[r^2 - 1], r, {[CurlyPhi], 0, 2 [Pi]}]

evaluates to Pi HeavisideTheta[r^2 - 1]. Note that HeavisideTheta != UnitStep, i.e. undefined behavior for r=1.

As a conclusion, your integration limit may not be element of your region of integration, which I think is a wrong behavior of mathematica.

I will focus on only one issue raised in the OP, namely, why MA produces the $pi/2$ result. I avoid on purpose any discussion on the mathematical justification of the integrals as I think this

I have never seen double integrals of distributions over bounded sets in math literature and don't know any definitions of such integrals. I would be very thankful for accessible and serious references.

goes beyond the scope of this forum. Coming back to this integral:

Integrate[DiracDelta[1 - x^2 - y^2], {x, -1, 1}, {y, -1, 1}]

1. Mathematica treats it as an iterated integral. Thus, it is sufficient to consider only the inner integral

Integrate[DiracDelta[1 - x^2 - y^2], {x, -1, 1}]

producing a wrong answer

(* ConditionalExpression[Boole[-1 < -Sqrt[1 - y^2] < 1]/(2 Sqrt[Abs[1 - y^2]]), -1 < Sqrt[1 - y^2] < 1] *)

2. The correct answer can obtained by the integration in the limits $[-a,a]$

 z=Integrate[DiracDelta[1 - x^2 - y^2], {x, -a, a}]
(*(Boole[-a < -Sqrt[1 - y^2] < a || a < -Sqrt[1 - y^2] < -a] + Boole[-a < Sqrt[1 - y^2] < a || a < Sqrt[1 - y^2] < -a])/(2 Sqrt[Abs[1 - y^2]])*)

and setting $a=1$ (no limit is necessary)

z1=z/.{a->1};
Integrate[z1, {y, -1, 1}]
(*    π    *)

3. Alternatively we can FullSimplify the intermediate result

za=FullSimplify[z, Assumptions -> a > 0 && -1 <= y <= 1]
(* Boole[Sqrt[1 - y^2] < a]/Sqrt[1 - y^2] *)

Integration over y can be performed in full generality

Integrate[Boole[Sqrt[1 - y^2] < a]/Sqrt[1 - y^2], {y, -1, 1}]

$$begin{array}{cc} Bigg{ & begin{array}{cc} pi & a>1 2 sin ^{-1}(a) & 0<aleq 1 end{array} end{array} tag{1}$$

Conclusion The reason for the $pi/2$ answer seems to be that mathematica forgot one root of the

$$1-x^2-y^2=0 tag{2}$$

equation. It certainly knows how to perform integrals of the $delta$-function of a function

$$ intdelta(f(x)),dx=int sum_kfrac{delta(x-alpha_k)}{left|f'(alpha_k)right|},dxtag{3} $$

where $alpha_k$ are the real roots of $f(x)=0$, however, in the OP example $a=1$ it misses one root and does not miss it otherwise.

It has been suggested by other answers that because the circle touches the boundaries of the integration domain this should somehow influence the result. It should not because these are points of measure zero: integration domain is $2D$, $delta$-function cuts out a $1D$ manifold out of it, whereas ill-defined integrand is on the $0D$ manifold. This seems to be just a coincidence that for $a=1$ mathematica misses one root of (3).