问题来源
在火灾报警电话时，基站一般难以确定与呼叫手机的距离，但是可以确定呼叫手机到两个临近基站的距离差值（原理待考证）。问题在于：能否根据呼叫手机到三个临近基站的两两距离差值确定呼叫手机的位置。
简化描述为：根据目标点到三个锚点间两两的差值确定空间中目标点的位置（轨迹）。
$$\left\{\begin{aligned} PA-PB =& d_1\\ PB - PC=& d_2 \\ PC - PA=& -(d_1+d_2) \end{aligned}\right.$$
注意：(1) $d_1,d_2$均可以是负数；(2) 以上三个方程实际上等价为两个方程，因为任意一个方程均可以由另外两个导出。
设$P$坐标为$(x,y,z)$，$A,B,C$坐标分别为$(a_1,b_1,c_1)$，$(a_2,b_2,c_2)$，$(a_3,b_3,c_3)$，则以上问题转化为：
$$
\left\{
\begin{aligned}
& \sqrt{(x-a_1)^2+(y-b_1)^2+(z-c_1)^2} - \sqrt{(x-a_2)^2+(y-b_2)^2+(z-c_2)^2} = d_1 \\
& \sqrt{(x-a_2)^2+(y-b_2)^2+(z-c_2)^2} - \sqrt{(x-a_3)^2+(y-b_3)^2+(z-c_3)^2} = d_2
\end{aligned}
\right.
$$
注意：以上方程包含三个未知数，而只有两个方程，因此该方程并不只有唯一解。
简化描述为：根据目标点到三个锚点间两两的差值确定空间中目标点的位置（轨迹）。
数学描述
设目标点为$P$，三个锚点分别为$A,B,C$，目标点到三个锚点的距离差满足如下关系：$$\left\{\begin{aligned} PA-PB =& d_1\\ PB - PC=& d_2 \\ PC - PA=& -(d_1+d_2) \end{aligned}\right.$$
注意：(1) $d_1,d_2$均可以是负数；(2) 以上三个方程实际上等价为两个方程，因为任意一个方程均可以由另外两个导出。
设$P$坐标为$(x,y,z)$，$A,B,C$坐标分别为$(a_1,b_1,c_1)$，$(a_2,b_2,c_2)$，$(a_3,b_3,c_3)$，则以上问题转化为：
$$
\left\{
\begin{aligned}
& \sqrt{(x-a_1)^2+(y-b_1)^2+(z-c_1)^2} - \sqrt{(x-a_2)^2+(y-b_2)^2+(z-c_2)^2} = d_1 \\
& \sqrt{(x-a_2)^2+(y-b_2)^2+(z-c_2)^2} - \sqrt{(x-a_3)^2+(y-b_3)^2+(z-c_3)^2} = d_2
\end{aligned}
\right.
$$
注意：以上方程包含三个未知数，而只有两个方程，因此该方程并不只有唯一解。
需求分析
- 可视化描述目标点到锚点距离差值一定的轨迹方程
- 计算出目标点的轨迹（数值解或解析解）
解决思路一
Mathematica提供了
如图所示，三种颜色的曲面表示分别表示到两锚点距离差一定的轨迹，其中蓝色曲线即为三个曲面的交线。（橙色曲线为对称的交线）
Solve函数可以用于求解以上方程，但计算结果极为冗长，几乎无法拷贝至其他软件中实现。并且在Mathematica中改变参数进行计算时计算效率也很低，同时存在“除零”的情况。该方法可以用于Demo，但几乎不具备实用性。Demo效果图如图所示。
解决思路二
1. 设三锚点为$A, B, C$， 以$A, B, C$形成的平面为$XOY$平面，且以$\overrightarrow{AB}$为$x$轴正方向，$AB$的中点为原点，中垂线为$y$轴正方向；
2. 根据锚点间距离，确定参数$c$（两锚点间距离为$2c$）；
3. 根据动点到两锚点间距离差， 确定参数$a$（距离差为$2a$）；
4. 相应地确定参数$b = \sqrt{c^2 - a^2}$；
5. 确定两锚点方向相对于$x$轴正方向的旋转角度与平移量；
6. 根据$2 \sim 4$确定相应双叶双曲面的标准方程，然后根据$5$中确定的旋转角度与平移量确定旋转平移后的双叶双曲面方程；
7. 固定$z_0$，即选定平面$Z=z_0$，代入双曲面方程，确定双曲面与该平面的交线（为双曲线）方程；
8. 求两双曲线方程的交点坐标（将$x(y)$代入，转化为关于$y$的一元四次方程[3, 4]）；
9. 利用一元四次方程的求根公式计算出求出$y$的解，带回原方程验证，确定最终解即为两双曲面在平面$Z=z_0$上的交点坐标；
10. 遍历$z_0$， 重复步骤$7 \sim 9$。
2. 根据锚点间距离，确定参数$c$（两锚点间距离为$2c$）；
3. 根据动点到两锚点间距离差， 确定参数$a$（距离差为$2a$）；
4. 相应地确定参数$b = \sqrt{c^2 - a^2}$；
5. 确定两锚点方向相对于$x$轴正方向的旋转角度与平移量；
6. 根据$2 \sim 4$确定相应双叶双曲面的标准方程，然后根据$5$中确定的旋转角度与平移量确定旋转平移后的双叶双曲面方程；
7. 固定$z_0$，即选定平面$Z=z_0$，代入双曲面方程，确定双曲面与该平面的交线（为双曲线）方程；
8. 求两双曲线方程的交点坐标（将$x(y)$代入，转化为关于$y$的一元四次方程[3, 4]）；
9. 利用一元四次方程的求根公式计算出求出$y$的解，带回原方程验证，确定最终解即为两双曲面在平面$Z=z_0$上的交点坐标；
10. 遍历$z_0$， 重复步骤$7 \sim 9$。
思路二中的关键步骤
坐标转换顺序
按照以上解决思路中给出的顺序，其中涉及到的坐标转换顺序如图所示。
 |
| 坐标转换顺序 |
平面旋转变换
由前面的分析可以知道，三锚点位于$XOY$平面， 那么相应的形成的双叶双曲面的分割面均垂直于$XOY$平面。因此旋转变换仅限于平面$XOY$中，对$z$方向不涉及旋转变换。下面介绍平面旋转变换。
平面直角坐标系中左边平移与旋转变换
如图所示，给出了原始坐标系$XOY$到旋转平移后的坐标系$X'O'Y'$之间坐标的相互转换关系。
注意：从$XOY$到$X'O'Y'$需要先进行平移变换，再进行旋转变换，旋转的角度为顺时针$\alpha$，因此在图中标识为$-\alpha$；反之，从$X'O'Y'$到$XOY$坐标系需要先进行旋转变换，旋转角为逆时针$\alpha$，然后进行平移转换。
 |
| 平面直角坐标系中左边平移与旋转变换 |
注意：从$XOY$到$X'O'Y'$需要先进行平移变换，再进行旋转变换，旋转的角度为顺时针$\alpha$，因此在图中标识为$-\alpha$；反之，从$X'O'Y'$到$XOY$坐标系需要先进行旋转变换，旋转角为逆时针$\alpha$，然后进行平移转换。
子问题
平面（空间）中到两定点距离差固定的动点轨迹为双曲线（双曲面）
证明：假设两定点间的举例为$2c$，而动点到两定点的举例差为$2a$（根据三角形两边之差小于第三边，可知：当$2a>2c$时， 是不存在的；而当$2a=2c$时，动点只可能位于两定点中的任意一个；以下的讨论均假定$2c>2a$）
为了便于讨论，假设两个定点（焦点）的坐标分别位于$(-c, 0)$和$(c, 0)$。设动点坐标为$(x, y)$， 则已知： $$\left| \sqrt{(x-c)^2 + y^2} - \sqrt{(x+c)^2 + y^2}\right| = 2a$$确定$(x, y)$的轨迹方程。
因为$\sqrt{(x-c)^2 + y^2} - \sqrt{(x+c)^2 + y^2} = \pm 2a$，移项可得
$$\sqrt{(x-c)^2 + y^2} = \pm 2a + \sqrt{(x+c)^2 + y^2}$$
因为$\sqrt{(x-c)^2 + y^2} - \sqrt{(x+c)^2 + y^2} = \pm 2a$，移项可得
$$\sqrt{(x-c)^2 + y^2} = \pm 2a + \sqrt{(x+c)^2 + y^2}$$
两边平方可得
$$ (x-c)^2 + y^2 = 4a^2 + (x+c)^2 + y^2 \pm 4a \sqrt{(x+c)^2 + y^2}$$
化简可得
$$ -4cx - 4a^2= \pm 4a \sqrt{(x+c)^2 + y^2}$$
即
$$ -cx - a^2= \pm a \sqrt{(x+c)^2 + y^2}$$
两边平方
$$ (cx + a^2)^2= a^2 \left[(x+c)^2 + y^2\right]$$
化简可得
$$x^2(c^2 - a^2) - a^2 y^2 = a^2 (c^2 - a^2)$$
另$b^2 = c^2 - a^2$，可得
$$x^2 b^2 - a^2 y^2 = a^2 b^2$$
两边同时除以$a^2 b^2$，即可得到双曲线的标准方程：
$$\frac{x^2}{a^2} - \frac{y^2}{b^2} = 1$$
将以上的问题推广到空间中是类似的，假设焦点坐标位于$(-c, 0, 0)$与$(c, 0, 0)$，动点坐标为$(x, y, z)$，则动点满足如下关系：
$$\left| \sqrt{(x-c)^2 + y^2 + z^2} - \sqrt{(x+c)^2 + y^2 + z^2}\right| = 2a$$
经过与上述类似的变换（将其中的$y^2$替换为$y^2+z^2$即可）可以得到如下方程：
$$\frac{x^2}{a^2} - \frac{y^2+z^2}{b^2} = 1$$
该方程为双叶双曲线的标准方程（焦点位于$x$轴），以下给出一个示例。
该方程为双叶双曲线的标准方程（焦点位于$x$轴），以下给出一个示例。
 |
| 双叶双曲面：$x^2-y^2-z^2=1$ |
平面与双叶双曲面所形成的交线为双曲线？
如图，给出一个双叶双曲面的切割demo。
图中分别给出了一个双叶双曲面和三个平面，其中：红色和绿色的平面是双叶双曲面的“渐进平面”即在$XOY$平面上的投影对应为该双叶双曲线在$XOY$平面投影（双曲线）的渐近线。而蓝色平面为切割平面，为了便于观察，绘制蓝色平面与双曲面的交线如下图所示。
从图中可以看出，蓝色交线为双曲线，下面给出证明以及给出该交线的方程。
首先这一结论是有前提条件的，考虑的“切割”平面是垂直于$XOY$平面的（对于焦点坐标在$X$轴上时，即标准方程为$\frac{x^2}{a^2} - \frac{y^2+z^2}{b^2} = 1$）
证明：切割平面在本例中垂直于$XOY$平面，平面方程可以表述为：$y = kx+m$其中“斜率$k$”的范围在上述“渐进平面”对应的“渐近线”的斜率范围内，即有：$|k| < \frac{b}{a}$
将平面方程$y=kx+m$代入双曲面标准方程$\frac{x^2}{a^2} - \frac{y^2+z^2}{b^2} = 1$，可得：
$$\left(\frac{1}{a^2} - \frac{k^2}{b^2}\right)x^2 - \frac{2km}{b^2}x + \frac{m^2}{b^2} = 1 + \frac{z^2}{b^2}$$
进一步整理可得：
$$\left(\frac{1}{a^2} - \frac{k^2}{b^2}\right)\left(x - \frac{a^2 km}{b^2 - a^2 k^2}\right)^2 - \frac{z^2}{b^2} = 1 + \frac{m^2}{b^2 - a^2 k^2}$$
注意到：$|k| < \frac{b}{a}$，故：$b^2 - a^2k^2 > 0$，从而上式中$\left(\frac{1}{a^2} - \frac{k^2}{b^2}\right)>0$，等式右侧$1 + \frac{m^2}{b^2 - a^2 k^2} > 1$，故：上式可以写为：
$$\frac{(x-m')^2}{(a')^2} - \frac{z^2}{(b')^2} = 1$$
其中：$m' = \frac{a^2 km}{b^2 - a^2 k^2}$, $a' = \sqrt{\frac{a^2 b^2}{b^2-a^2 k^2}\left(1+\frac{m^2}{b^2-a^2k^2}\right)}$，$b' = \sqrt{b^2 \left(1 + \frac{m^2}{b^2 - a^2 k^2}\right)}$
同时还需要满足$y=kx+m$。
示例：以$x^2 - y^2 - z^2 = 1$被平面$x+2y=0$切割为例，可以绘制出解如下
图中绿色曲面即为$\frac{(x-m')^2}{(a')^2} - \frac{z^2}{(b')^2} = 1$表示的曲面。
 |
| Demo：双叶双曲面“切割” |
 |
| 双叶双曲面切割平面及其交线 |
首先这一结论是有前提条件的，考虑的“切割”平面是垂直于$XOY$平面的（对于焦点坐标在$X$轴上时，即标准方程为$\frac{x^2}{a^2} - \frac{y^2+z^2}{b^2} = 1$）
证明：切割平面在本例中垂直于$XOY$平面，平面方程可以表述为：$y = kx+m$其中“斜率$k$”的范围在上述“渐进平面”对应的“渐近线”的斜率范围内，即有：$|k| < \frac{b}{a}$
将平面方程$y=kx+m$代入双曲面标准方程$\frac{x^2}{a^2} - \frac{y^2+z^2}{b^2} = 1$，可得：
$$\left(\frac{1}{a^2} - \frac{k^2}{b^2}\right)x^2 - \frac{2km}{b^2}x + \frac{m^2}{b^2} = 1 + \frac{z^2}{b^2}$$
进一步整理可得：
$$\left(\frac{1}{a^2} - \frac{k^2}{b^2}\right)\left(x - \frac{a^2 km}{b^2 - a^2 k^2}\right)^2 - \frac{z^2}{b^2} = 1 + \frac{m^2}{b^2 - a^2 k^2}$$
注意到：$|k| < \frac{b}{a}$，故：$b^2 - a^2k^2 > 0$，从而上式中$\left(\frac{1}{a^2} - \frac{k^2}{b^2}\right)>0$，等式右侧$1 + \frac{m^2}{b^2 - a^2 k^2} > 1$，故：上式可以写为：
$$\frac{(x-m')^2}{(a')^2} - \frac{z^2}{(b')^2} = 1$$
其中：$m' = \frac{a^2 km}{b^2 - a^2 k^2}$, $a' = \sqrt{\frac{a^2 b^2}{b^2-a^2 k^2}\left(1+\frac{m^2}{b^2-a^2k^2}\right)}$，$b' = \sqrt{b^2 \left(1 + \frac{m^2}{b^2 - a^2 k^2}\right)}$
同时还需要满足$y=kx+m$。
示例：以$x^2 - y^2 - z^2 = 1$被平面$x+2y=0$切割为例，可以绘制出解如下
参考
4. 维基百科 - 四次方程
No comments:
Post a Comment