艾科纳尔方程

艾科纳尔方程（eikonal equation，亦称程函方程）是一类规定未知函数梯度大小的一阶非线性偏微分方程。在传播问题中，未知函数可以表示波前第一次到达各位置的时间：局部速度越慢，单位距离累积的时间越多，时间函数在空间中变化得也越快。

声音穿过不同介质、光线通过折射界面，或者机器人绕过障碍到达目的地，都可以产生这一类问题。方程描述的是在既定速度场与允许路径下的最早到达时间；它本身不描述波的振幅，也不保留同一点处较晚到达的其他传播分支。

到达时间与局部速度

设位置为 $x\in \Omega \subset {ℝ}^n$，局部各向同性速度为 $F(x)>0$，源集合为 $\Gamma$。在源之外的 $\Omega \setminus \Gamma$，到达时间 $T(x)$ 满足 $$F(x)|\mathrm{\nabla }T(x)|=1,T{|}_{\Gamma }=0.$$ 这里 $\mathrm{\nabla }T$ 是各空间偏导组成的向量，$|\mathrm{\nabla }T|$ 是它的欧氏长度。若时间用秒、距离用米，则 $\mathrm{\nabla }T$ 的单位为 $\mathrm{s}/\mathrm{m}$，恰是速度倒数的单位。也常把 $s(x)=1/F(x)$ 称为慢度，写成 $|\mathrm{\nabla }T|=s$；不同文献中的字母约定并不一致。Sethian 的概述使用速度形式，并说明其最早到达时间含义。

在 $T$ 光滑且梯度非零处，等时面 $T(x)=t$ 的单位法向为 $n=\mathrm{\nabla }T/|\mathrm{\nabla }T|$。沿法向移动一小段 $\mathrm{d}\ell$，到达时间增加约 $|\mathrm{\nabla }T|\mathrm{d}\ell$。它又等于“距离除以速度”，即 $\mathrm{d}\ell /F$，因此得到上述方程。这一推导也表明方程假设方向不影响局部传播速度；各向异性情形通常需要更一般的汉密尔顿—雅可比方程。

路径费用表示

对从源到 $x$ 的一条可行曲线 $\gamma$，传播时间为 $$𝒯(\gamma )={\int }_{\gamma }\frac{\mathrm{d}s}{F}.$$ 其中 $\mathrm{d}s$ 为弧长。在速度有正的上下界、相关位置可到达等适当条件下，最早到达时间由 $$T(x)=\underset{\gamma :\Gamma \to x}{inf}𝒯(\gamma )$$ 给出。障碍物通过限制可行路径进入问题；穿越禁区的直线即使更短，也不能参加这个最小化。将到达路径最后一小段的费用与此前最低费用相接，就得到方程背后的动态规划关系。

若速度处处为常数 $c>0$，且无障碍、源为一个点 $a$，任意路径的长度至少为 $|x-a|$，直线达到下界，所以 $$T(x)=\frac{|x-a|}{c}.$$ 在 $x\ne a$ 处，$\mathrm{\nabla }T=(x-a)/(c|x-a|)$，故 $c|\mathrm{\nabla }T|=1$。在源点梯度不存在，应结合源条件理解；不能因为源点没有经典导数就舍弃这一解。

二维常速 c=1 时，t=1、2、3 的等时线分别是半径 1、2、3 的圆。相邻波前的法向距离等于时间间隔。

两个源与不可微的交汇线

设两个源位于 $(-1,0)$ 和 $(1,0)$，速度为一，二者同时开始传播。对任一点比较两条直线路径即可： $$T(x,y)=min\{\sqrt{(x+1{)}^2+y^2},\sqrt{(x-1{)}^2+y^2}\}.$$ 在左半平面，左源先到；在右半平面，右源先到。直线 $x=0$ 上两种到达时间相同，但左右梯度的水平分量相反，因此这里通常不可微。沿连接两源的线段，时间简化为 $T(x,0)=1-|x|$，$-1\le x\le 1$；中点是到达时间的尖峰。

等时线由两组圆弧拼接而成。x=0 是两个源同时到达的位置；它是时间函数的不可微集合，不是传播无法到达的障碍。

这类尖点是几何传播本身的结果。最早到达时间应在黏性解意义下解释：用从上、下接触的光滑测试函数检查偏微分不等式，而不强迫所有点都有经典梯度。对于不连续的速度界面，还要结合路径费用或相应的界面条件选择解，不能任意拼接两侧满足方程的函数。

两种速度下的折射算例

考虑分界线 $x=0$，左侧速度为 $F_1=1$，右侧为 $F_2=2$，单位已统一。源点 $A=(-2,0)$，目标 $B=(2,3)$。每个半平面内速度恒定，给定穿越点后，区域内的最优路径必为直线。若过界点为 $C=(0,y)$，总时间为 $$\tau (y)=\sqrt{4+y^2}+\frac{1}{2}\sqrt{4+(3-y{)}^2}.$$ 第二段速度较大，因此同样长度只支付一半时间。求导得 $${\tau }^{\prime }(y)=\frac{y}{\sqrt{4+y^2}}-\frac{3-y}{2\sqrt{4+(3-y{)}^2}}.$$ 而 $${\tau }^{″}(y)=\frac{4}{(4+y^2{)}^{3/2}}+\frac{2}{(4+(3-y{)}^2{)}^{3/2}}>0.$$ 所以时间函数严格凸，驻点若存在就是唯一最小点。${\tau }^{\prime }(0)<0$、${\tau }^{\prime }(3)>0$，故最优点位于零和三之间。对这个区间二分求根，可得 $$y_{\ast }\approx 0.797018742,\tau (y_{\ast })\approx \mathrm{3.640670680.}$$ 连接 $A,B$ 的直线在 $y=1.5$ 处过界，两段长度均为 $2.5$，时间为 $2.5+2.5/2=3.75$。折线路径虽然总长度更长，却把更多路程放在快速区域，使总时间下降约 $2.92\mathrm{\%}$。

实线经 C=(0,0.7970) 到达 B，用时约 3.6407；虚线为直线路径，用时 3.75。界面法线沿水平方向，夹角应从这条法线量起。

设两段路径与界面法线的夹角分别为 ${\theta }_1,{\theta }_2$，驻点方程可写成 $$\frac{\sin {\theta }_1}{F_1}=\frac{\sin {\theta }_2}{F_2}.$$ 这就是速度形式的 Snell 折射关系。本例两角约为 $21.73^{\circ }$ 和 $47.76^{\circ }$。右侧传播更快，最优路径在右侧偏离法线更多。若改用折射率，折射率与速度成反比，关系也可写为两侧“折射率乘正弦”相等。

正交网格上的上风更新

在二维等间距网格上，设步长为 $h>0$，当前点速度为 $F>0$。从水平方向两邻点取较小的已知到达值，从竖直方向也取较小值，记为 $a,b$，并交换次序使 $a\le b$。一阶上风离散近似为 $$(max\{T-a,0\}{)}^2+(max\{T-b,0\}{)}^2={\delta }^2,\delta =\frac{h}{F}.$$ 正部保证只使用早于当前点的方向信息。若 $a+\delta \le b$，竖直或水平中的一个方向尚不能参与更新，故 $T=a+\delta$。否则两个方向都参与，把两个平方展开，取满足 $T\ge b$ 的根，得到 $$T=\{\begin{array}{cc}a+\delta ,& b-a\ge \delta ,\\ \frac{a+b+\sqrt{2{\delta }^2-(a-b{)}^2}}{2},& b-a<\delta .\end{array}$$ 第二分支的条件使根号有定义，并使所得值大于两个所用邻值；第一分支则给出单方向传播。更新值不早于提供信息的邻点，这一因果性使快速行进法可以从小到大确定到达时间。用堆维护候选点时，标准实现的复杂度为 $O(N\log N)$，其中 $N$ 是网格点数。

以下程序只实现两个有限、非负邻值的局部更新，不是完整快速行进求解器，不修改输入。尚未到达、以无穷表示的邻点应由外层求解器另作处理。

Python 3

import math


def eikonal_update(a, b, h=1.0, speed=1.0):
    """二维等间距网格；返回两个有限邻值的上风更新时间。"""
    values = (a, b, h, speed)
    if not all(math.isfinite(v) for v in values):
        raise ValueError("inputs must be finite")
    if min(a, b) < 0 or h <= 0 or speed <= 0:
        raise ValueError("nonnegative times and positive h, speed required")
    a, b = sorted((a, b))
    delta = h / speed
    if not math.isfinite(delta) or delta <= 0:
        raise OverflowError("h / speed is outside floating-point range")
    gap = b - a
    if gap >= delta:
        result = a + delta
    else:
        # 先将根与差各除以二，避免可表示结果的中间计算溢出。
        result = a + gap / 2 + delta * math.sqrt(
            0.5 - 0.25 * (gap / delta) ** 2
        )
    if not math.isfinite(result):
        raise OverflowError("updated time is outside floating-point range")
    return result


print(eikonal_update(2.0, 2.0))  # 2.7071067811865475
print(eikonal_update(2.0, 4.0))  # 3.0

当 $a=b=2,h=F=1$ 时，$T=2+1/\sqrt{2}$；两个方向共同提供信息。当 $a=2,b=4$ 时，$T=3<4$，较晚的邻值不能参与。单次更新只需常数时间和常数额外空间。该离散式有网格误差；它给出的不是任意网格上都精确的欧氏距离。

几何传播的边界

几何光学使用艾科纳尔方程时，通常假设波长相对介质变化尺度足够小。衍射、干涉、幅度变化及多次反射后的较晚到达，需要波动方程或额外的振幅、相位与路径信息。速度为零的障碍不能直接代入 $1/F$；应先从可行域中排除，或者明确采用怎样的退化方程与解概念。