数学 / 三角函数 · 从角到波 / 三角函数模型与周期数据的拟合 待审核 7 / 13
简谐运动 · 最小二乘

三角函数模型与周期数据的拟合

潮汐的涨落、一年里日照时长的伸缩、气温的季节变化、交流电的电压，都是同一件事的不同外衣：同一个形态每隔一段时间重来一次。描述它们用的是同一个模型 y=Asin⁡(ωt+φ)+by = A\sin(\omega t + \varphi) + b。一般正弦型函数的图象变换 已把四个参数的几何效果说清，本页处理反方向的问题：给定一批带误差的观测值，如何把这四个数定出来。

1 · 简谐运动的参数含义

物理里最标准的周期运动是简谐运动 (simple harmonic motion)：质点受到的回复力与位移成正比、方向相反，解出的位移随时间的变化恰是 y=Asin⁡(ωt+φ)y = A\sin(\omega t + \varphi)。四个参数各有物理名字。

AA 是振幅，位移能达到的最大值，量纲与 yy 相同。ω\omega 是角频率，量纲是 rad/s，与周期、频率的换算是 T=2π/ωT = 2\pi/\omega、f=1/T=ω/2πf = 1/T = \omega/2\pi。φ\varphi 是初相，t=0t = 0 那一刻的相位；它不改变曲线形状，只决定从周期的哪一处开始计时。bb 是平衡位置，振动绕之往返的那条中线；理想的简谐运动里它是零，但一切实测数据都需要它，因为观测量的零点由测量方式而非物理定义。

ω\omega 与 φ\varphi 的地位并不对称。φ\varphi 换一个值，曲线只沿时间轴平移；ω\omega 换一个值，曲线的疏密就变了，任何平移都补不回来。这条不对称在 §3 会变成算法上的分水岭。

2 · 由图象读参数的顺序

手工建模有一条固定顺序，每一步只用前一步用不上的信息。

先读 bb 与 AA：b=ymax⁡+ymin⁡2b = \dfrac{y_{\max} + y_{\min}}{2}，A=ymax⁡−ymin⁡2A = \dfrac{y_{\max} - y_{\min}}{2}。两者只用到纵向极值，与时间轴无关。

再读 ω\omega：量出相邻两个同类点的时间间隔即是周期 TT，取 ω=2π/T\omega = 2\pi/T。同类点指两个峰、两个谷，或两个同方向穿过中线的点。峰与谷之间、相邻两次穿中线之间都只有半个周期，量成一个周期就差一倍。

最后读 φ\varphi：取一个已知点代回相位。习惯上取上行穿过中线的时刻 t0t_0，该处相位为 2kπ2k\pi，故 φ=−ωt0\varphi = -\omega t_0。φ\varphi 只在模 2π2\pi 的意义下确定，通常取 ∣φ∣<π|\varphi| < \pi 的那一个。

这套顺序对干净的图象够用，对带噪声的实测数据不好用：ymax⁡y_{\max} 与 ymin⁡y_{\min} 各只由一个点决定，一个离群值就能把 AA 整个抬起来。数据一多就该换成最小二乘。

3 · 频率已知时的线性最小二乘

最小二乘 (least squares) 要最小化 ∑i[yi−Asin⁡(ωti+φ)−b]2\displaystyle\sum_i \bigl[y_i - A\sin(\omega t_i + \varphi) - b\bigr]^2。这个式子对 AA 与 bb 是线性的，对 φ\varphi 与 ω\omega 都不是，看上去只能上非线性优化。和角公式把它拆开：

Asin⁡(ωt+φ)=(Acos⁡φ)sin⁡ωt+(Asin⁡φ)cos⁡ωtA\sin(\omega t + \varphi) = (A\cos\varphi)\sin\omega t + (A\sin\varphi)\cos\omega t

只要 ω\omega 已知，sin⁡ωti\sin\omega t_i 与 cos⁡ωti\cos\omega t_i 都是可以事先算出的常数，模型对三个新参数 p=Acos⁡φp = A\cos\varphi、q=Asin⁡φq = A\sin\varphi、bb 就是线性的。剩下的是一个普通的三元线性最小二乘：设计矩阵的三列分别是 sin⁡ωti\sin\omega t_i、cos⁡ωti\cos\omega t_i 与全 11 列，一组 3×33 \times 3 的正规方程一次解出，不需要初值、不需要迭代、没有局部极小。解完再合回去，A=p2+q2A = \sqrt{p^2 + q^2}，φ=atan2⁡(q, p)\varphi = \operatorname{atan2}(q,\ p)。

这正是 辅助角公式与三角式的最值 的实用兑现。那一页把 asin⁡x+bcos⁡xa\sin x + b\cos x 合成 Rsin⁡(x+φ)R\sin(x + \varphi)，本节是同一个变换反着用：先拆成两个能线性求解的系数，解完再合成振幅与初相。atan2⁡\operatorname{atan2} 而非 arctan⁡(q/p)\arctan(q/p) 的理由也照搬那一页，两个系数同时变号给出的相位相差 180°180°。

这条捷径值多少是量出来的。取北纬 39.9° 全年每 7 天一个日照读数共 53 个点、噪声上界 0.25 小时，在真频率处线性法一次解出的残差平方和是 0.989797。换成在 (A, φ)(A,\ \varphi) 平面上铺网格搜索、每格再由残差均值定 bb：40×4040 \times 40 的网格给出 1.100117，200×200200 \times 200 给出 1.018549，1000×10001000 \times 1000 共一百万次求值仍是 0.989947，依然高于线性解。网格加密只是越来越逼近同一个解，追不平。

警示 · 「线性」说的是模型对参数线性，不是对自变量线性。sin⁡ωti\sin\omega t_i 本身当然是 tit_i 的非线性函数，但它在正规方程里只是一列已知的数。判断一个拟合问题能否一次解出，看的永远是参数怎么进式子，不是自变量怎么进。

4 · 频率未知时的残差谱

ω\omega 未知时上一节的办法搬不动，它藏在 sin⁡ωti\sin\omega t_i 里面，当参数就回到非线性问题。可行的做法是把它单独拎出来做一维搜索：每个候选 ω\omega 做一次线性拟合，记下残差平方和，连成一条曲线取最低点。这条曲线称为残差谱，与周期图 (periodogram) 是同一件事的两种记法。

一维搜索比三维搜索便宜得多，而且这条曲线本身就是诊断图。真周期处是一道深谷，谷宽大致是观测跨度的倒数，观测跨得越长，频率定得越准。上面那 53 个日照读数的总变差是 206.5878，而残差谱的最高点是 206.566：试探周期离真值足够远时，拟合出的正弦几乎一点变化都吸收不掉，残差退回「只拟合一条水平线」的水平。谷底 0.989 只有它的 0.48%。

谱的两侧并不对称。同一批数据上，试探周期缩到 180 天时残差平方和是 205.87，放到 620 天却只剩 11.10。周期远长于观测跨度的正弦在这段数据上退化成一条缓慢的趋势线，仍能吸收掉不少变化，长周期那一侧的残差因此涨得慢。谱峰与谷底相差两个半数量级，画这条曲线时纵轴须取对数，线性纵轴会把整个谷压成贴底的一条直线。

图 4-1 · 上半幅是带噪声的观测点与当前试探周期下线性拟合出的曲线，下半幅是残差平方和随周期变化的谱，纵轴取对数。可调试探周期、噪声强度与采样点数，观察谷的深浅与两侧的不对称。

有一处容易误读：残差谱的最低点不一定落在真周期上。数据若不是纯正弦，最小二乘会拿频率去换一点残差。用几何地平口径的日照时长逐日数据做实验，它完全无噪声、真周期恰是 365 天，而残差谱的最低点落在 366.57 天，残差平方和从真频率处的 0.61897 降到 0.60268。挪掉 1.57 天换来 2.6% 的残差，这不是数值误差，是模型失配的必然结果。周期有独立来源时，宁可把 ω\omega 钉死再拟合另外三个参数。

5 · 日照时长的拟合

日照时长有解析近似，可以据此造出真实形状的周期数据。太阳赤纬 δ\delta 随年内第几天 NN 变化，Cooper 近似式给出 δ=23.45°sin⁡360°(284+N)365\delta = 23.45°\sin\dfrac{360°(284+N)}{365}[4]；日出方程给出半昼弧 HH 满足

cos⁡H=sin⁡h0−sin⁡ϕsin⁡δcos⁡ϕcos⁡δ\cos H = \frac{\sin h_0 - \sin\phi\sin\delta}{\cos\phi\cos\delta}

其中 ϕ\phi 是纬度，日照时长是 2H/152H/15 小时（HH 以度计）。h0h_0 是把日出定在太阳中心的哪个高度上：几何地平取 00，通行口径取 −0.833°-0.833°，那是大气折射约 34′34' 加上日面半径约 16′16'。

北纬 39.9° 全年逐日的这条曲线，用一年期正弦拟合（ω=2π/365\omega = 2\pi/365）得 A=2.777507A = 2.777507 小时、b=12.156520b = 12.156520 小时、φ=−79.89°\varphi = -79.89°。中线比 12 小时多出 9.39 分钟，来自那 −0.833°-0.833°：把日出定义放低，每天两头各多算一点。

图 5-1 · 上半幅是日照时长的逐日曲线、抽样得到的观测点与拟合出的正弦，下半幅是放大后的残差。可调纬度、采样间隔与噪声，并切换日出高度的口径，读数给出由观测点拟合出的振幅与中线，以及逐日残差里各次谐波的振幅。

残差的均方根是 2.47 分钟，最大残差 4.29 分钟出现在夏至那天。原先以为残差里的次要成分是二倍频，实测不是：把残差再分别对 2ω2\omega、3ω3\omega、4ω4\omega 各做一次线性拟合，三个振幅是 0.7356 分、3.4143 分、0.0327 分，三倍频比二倍频高 4.6 倍。

追下去能找到原因。几何地平口径下 cos⁡H=−tan⁡ϕtan⁡δ\cos H = -\tan\phi\tan\delta，日照时长减去 12 小时是 δ\delta 的奇函数，而 δ\delta 本身是时间的正弦，奇函数复合正弦只生奇次谐波，偶次项恒为零。实测证实了这一点：h0=0h_0 = 0 时 2ω2\omega 与 4ω4\omega 的振幅都是 0.0000 分，3ω3\omega 是 3.4101 分，中线 bb 也恰好是 12 小时。那 0.7356 分的二倍频整个来自折射修正，sin⁡h0\sin h_0 是 cos⁡H\cos H 分子里一个与 δ\delta 无关的常数项，它把奇对称破掉了。

一条正弦够不够用，取决于纬度。同样的拟合在赤道上残差均方根是 0.21 分，20°N 是 0.83 分，39.9°N 是 2.47 分，55°N 涨到 7.37 分，北极圈 66.5°N 是 55.10 分。最后这个纬度上已经出现极昼与极夜，曲线顶部被削成一段水平线，正弦无论如何配不上，要么换模型，要么只拟合中纬度那一段。

6 · 参考文献

  1. Simple harmonic motion. Wikipedia. 回复力与位移成正比的运动，解为正弦型函数。https://en.wikipedia.org/wiki/Simple_harmonic_motion
  2. Linear least squares. Wikipedia. 正规方程、设计矩阵，以及「对参数线性」的准确含义。https://en.wikipedia.org/wiki/Linear_least_squares
  3. Least-squares spectral analysis. Wikipedia. 逐频率做线性拟合得到的谱，与周期图的关系。https://en.wikipedia.org/wiki/Least-squares_spectral_analysis
  4. Cooper, P. I. (1969). The absorption of radiation in solar stills. Solar Energy, 12(3), 333–346.
  5. Sunrise equation. Wikipedia. 半昼弧公式与日出高度 h0h_0 的取法。https://en.wikipedia.org/wiki/Sunrise_equation
  6. Atmospheric refraction. Wikipedia. 地平附近约 34′ 的折射量，以及它对日出时刻的影响。https://en.wikipedia.org/wiki/Atmospheric_refraction
  7. Daytime. Wikipedia. 日照时长随纬度与季节的变化，以及分点前后各地并非恰好 12 小时的原因。https://en.wikipedia.org/wiki/Daytime