Now the answer is almost complete: modulo some extra work on the strictness of relevant inequalities, we do have the uniqueness. The additional ideas used to come to this conclusion are these: (i) to kind of reduce approximating a convex set by an ellipse by, vice versa, approximating an ellipse (and even a round disk) by shear-shifted versions of the convex set; (ii) shear-symmetrization; and
(iii)
minimizing the symmetric difference given the same areas is equivalent to maximizing the area of the intersection (the latter observation borrowed from Matt F.).
How to show that the best possible ellipse is unique?
Suppose that $E_1$ and $E_2$ are distinct best possible ellipses. By an appropriate rescaling of everything in the directions of principal axes of $E_1$, without loss of generality (wlog) $E_1$ is a round disk, say $D$, of radius $1/\sqrt\pi$.
Then the width of $E_2$ in some direction is the same as that of the round disk $E_1$ (that is, $2/\sqrt\pi$), where the width of a set in a given direction is defined as the width of the narrowest infinite band that contains the set and goes in
the direction perpendicular to the given one. This follows because the product of the widths of ellipse $E_2$ in the directions of its principal axes is $(2/\sqrt\pi)^2$.
So, there exist the following: (i) a real $t$; (ii) a vector $b\in\mathbb{R}^2$; and (iii) an orthonormal basis $(e_1,e_2)$ of $\mathbb{R}^2$ such that,
for the (shear-shift) affine operator $A=A_t=A_{b,t}$ on $\mathbb{R}^2$ defined by the conditions $A\mathbf0=b$, $Ae_1=b+e_1$, and $Ae_2=b+e_2+te_1$, we have $AE_1=AD=E_2$.
Let then $E_0:=\dfrac{I+A}2\,D$, where $I$ is the identity operator. Then $E_0$ is an ellipse of area $1$.
For real $y$, let $[u,v]=[u(y),v(y)]=K_y:=K(y):=\{x\in\mathbb R\colon xe_1+ye_2\in K\}$ be the $y$-"cross-section" of $K$. Similarly define the $y$-"cross-sections" $E_1(y)=[r_1,s_1]$ and $E_2(y)=[r_2,s_2]$ of $E_1$ and $E_2$.
Then the $y$-"cross-section" $E_0(y)$ of $E_0$ is $[r_0,s_0]=[\frac{r_1+r_2}2,\frac{s_1+s_2}2]$. Let $|\cdot|$ denote the Lebesgue measure on $\mathbb{R}$ or $\mathbb{R}^2$.
Then $|E_j\cap K|=\int_{\mathbb{R}}\delta_j(y)\,dy$ for $j=0,1,2$, where $\delta_j(y):=|E_j(y)\cap K(y)|$.
So, if we could show that
\begin{equation}
\delta_0(y)\ge\tfrac12\,\delta_1(y)+\tfrac12\,\delta_2(y) \tag{*}
\end{equation}
and this inequality is strict for some $y$ given that $E_1$ and $E_2$ are distinct, then we would obtain the desired contradiction: $|E_0\cap K|>\tfrac12\,|E_1\cap K|+\tfrac12\,|E_2\cap K|=|E_1\cap K|=|E_2\cap K|$.
Here we have used the mentioned additional idea (iii).
Unfortunately, inequality $(*)$ does not hold in general, without any assumptions on the convex sets $E_j$. Here we need the mentioned additional ideas (i) and (ii).
Consider first the round disk $E_1=D$, which is an optimal approximation of $K$. I then claim that for all the $y$-"cross-sections" of $E_1$ and $K$ we have $\frac{r_1+s_1}2=\frac{u+v}2$. Indeed, since $D$ is an optimal approximation of $K$, it is also an optimal approximation of $K^-$, where $K^-$ is obtained from $K$ by reflection about the diameter of $D$ parallel to $e_2$. Wlog $D$ is centered at the origin. So, for all $y$ we have $r_1=-s_1$, and the $y$-"cross-section" of $K^-$ is $K^-_y=[-v,-u]=-K_y$. Then for $K^0:=\frac12\,(K+K^-)$ we can check that
$|D\cap K^0_y|\ge\tfrac12\,|D\cap K_y| +\tfrac12\,|D\cap K^-_y|$ for all $y$. So,
the disk $D$ approximates $K^0$ no worse than it does $K$ (or $K^-$). Extra work is needed here to show that $D$ approximates $K^0$ (strictly) better than it does $K$ (or $K^-$) unless $K^-=K$.
But $K^0$ is the image of $K$ under an (area-preserving) shear-shift transformation (say $B$). So, modulo the mentioned extra work, we will conclude that $B^{-1}E_1=B^{-1}D$ is a better approximation of $K$ -- unless $K^-=K$. So, wlog $K^-=K$, which verifies the claim $\frac{r_1+s_1}2=\frac{u+v}2$. Similarly, by shear-shifting, $\frac{r_2+s_2}2=\frac{u+v}2$.
Given these two conditions, one can verify that $(*)$ holds (in that verification, wlog $\frac{r_1+s_1}2=\frac{r_2+s_2}2=\frac{u+v}2=0$).
It will remain to show that inequality $(*)$ will be strict for some $y$ unless $E_1=E_2$.
For readers' convenience(?), I am retaining below, under the horizontal line, the previous version of the answer.
This is not a complete answer, but I think it has a chance to lead to one.
