数学建模-球谐函数

本文主题：已知二维空间的不均匀散点，基于球谐函数拟合曲面！

1.球谐函数

球谐函数模型的表达式见式 ( 1 ) − ( 4 ) (1)-(4) (1)−(4)。式中， β \\beta β为纬度， s s s为经度， N N N为球谐函数的最大展开阶数， C ~ n m \\widetilde{C}_{nm} C nm​和 S ~ n m \\widetilde{S}_{nm} S nm​为未知的球谐函数系数 (这是我们要计算的变量) ， P ~ n m ( sin ⁡ β ) \\widetilde{P}_{nm}(\\sin\\beta) P nm​(sinβ)为n度m阶勒让德函数， M C ( n , m ) MC(n,m) MC(n,m)为归化函数， P n m ( sin ⁡ β ) P_{nm}(\\sin\\beta) Pnm​(sinβ)为经典勒让德函数，， δ 0 m \\delta_{0m} δ0m​是kronecker型 δ \\delta δ函数。
V T E C ( β , s ) = ∑ n = 0 N ∑ m = 0 n P ~ n m ( sin ⁡ β ) ( C ~ n m cos ⁡ ( m s ) + S ~ n m sin ⁡ ( m s ) ) VTEC(\\beta,s)=\\sum_{n=0}^N\\sum_{m=0}^n\\widetilde{P}_{nm}(\\sin\\beta)(\\widetilde{C}_{nm}\\cos(ms)+\\widetilde{S}_{nm}\\sin(ms)) VTEC(β,s)=n=0∑N​m=0∑n​P nm​(sinβ)(C nm​cos(ms)+S nm​sin(ms))

P ~ n m ( sin ⁡ β ) = M C ( n , m ) ⋅ P n m ( sin ⁡ β ) \\widetilde{P}_{nm}(\\sin\\beta)=MC(n,m)\\cdot P_{nm}(\\sin\\beta) P nm​(sinβ)=MC(n,m)⋅Pnm​(sinβ)

M C ( n , m ) = ( n − m ) ! ( 2 n + 1 ) ( 2 − δ o m ) / ( n + m ) ! MC(n,m)=\\sqrt{(n-m)!(2n+1)(2-\\delta_{om})/(n+m)!} MC(n,m)=(n−m)!(2n+1)(2−δom​)/(n+m)! ​

δ 0 m = { 0 , m = 0 1 , m ≠ 0 \\delta_{0m}=\\left\\{ \\begin{aligned} 0,m=0 \\\\ 1,m\\ne0 \\end{aligned}\\right. δ0m​={0,m=01,m​=0​

2.方程求解

现在，我们建立计算 V T E C VTEC VTEC的等式，见式 ( 5 ) (5) (5)，其中， L ^ \\hat{L} L^为观测值， E E E为系统偏差。
V T E C = L ^ + E VTEC=\\hat{L}+E VTEC=L^+E

现在，建立式(5)的误差方程,见式 ( 6 ) (6) (6)。
V = V T E C − E − L V=VTEC-E-L V=VTEC−E−L

本文以二阶二次多项式为例，利用球谐函数展开 V T E C VTEC VTEC,见式 ( 7 ) (7) (7)。

V T E C = P ~ 00 ( s i n β ) ⋅ C ~ 00 + P ~ 10 ( s i n β ) ⋅ C ~ 10 + P ~ 11 ( s i n β ) ⋅ C ~ 11 cos ⁡ ( s ) + P ~ 11 ( s i n β ) ⋅ S ~ 11 sin ⁡ ( s ) + P ~ 20 ( s i n β ) ⋅ C ~ 20 + P ~ 21 ( s i n β ) ⋅ C ~ 21 cos ⁡ ( 2 s ) + P ~ 21 ( s i n β ) ⋅ S ~ 21 sin ⁡ ( s ) + P ~ 22 ( s i n β ) ⋅ C ~ 22 cos ⁡ ( 2 s ) + P ~ 22 ( s i n β ) ⋅ S ~ 22 sin ⁡ ( 2 s ) \\begin{aligned} VTEC&=\\widetilde{P}_{00}(sin\\beta)\\cdot\\widetilde{C}_{00}\\\\ &+\\widetilde{P}_{10}(sin\\beta)\\cdot\\widetilde{C}_{10} +\\widetilde{P}_{11}(sin\\beta)\\cdot\\widetilde{C}_{11}\\cos(s) +\\widetilde{P}_{11}(sin\\beta)\\cdot\\widetilde{S}_{11}\\sin(s)\\\\ &+\\widetilde{P}_{20}(sin\\beta)\\cdot\\widetilde{C}_{20} +\\widetilde{P}_{21}(sin\\beta)\\cdot\\widetilde{C}_{21}\\cos(2s) +\\widetilde{P}_{21}(sin\\beta)\\cdot\\widetilde{S}_{21}\\sin(s)\\\\ &+\\widetilde{P}_{22}(sin\\beta)\\cdot\\widetilde{C}_{22}\\cos(2s) +\\widetilde{P}_{22}(sin\\beta)\\cdot\\widetilde{S}_{22}\\sin(2s) \\end{aligned} VTEC​=P 00​(sinβ)⋅C 00​+P 10​(sinβ)⋅C 10​+P 11​(sinβ)⋅C 11​cos(s)+P 11​(sinβ)⋅S 11​sin(s)+P 20​(sinβ)⋅C 20​+P 21​(sinβ)⋅C 21​cos(2s)+P 21​(sinβ)⋅S 21​sin(s)+P 22​(sinβ)⋅C 22​cos(2s)+P 22​(sinβ)⋅S 22​sin(2s)​

将展开的 V T E C VTEC VTEC代入到式 ( 6 ) (6) (6)，然后将式 ( 6 ) (6) (6)简记为 V = A X − L V=AX-L V=AX−L，其中， A A A和 X X X见式 ( 8 ) (8) (8)。

A = [ P ~ 00 ( s i n β ) , P ~ 10 ( s i n β ) , P ~ 11 ( s i n β ) ⋅ cos ⁡ ( s ) , P ~ 11 ( s i n β ) ⋅ sin ⁡ ( s ) , P ~ 20 ( s i n β ) , P ~ 21 ( s i n β ) ⋅ cos ⁡ ( 2 s ) , P ~ 21 ( s i n β ) ⋅ sin ⁡ ( s ) , P ~ 22 ( s i n β ) ⋅ cos ⁡ ( 2 s ) , P ~ 22 ( s i n β ) ⋅ sin ⁡ ( 2 s ) , − 1 ] . X = [ C ~ 00 , C ~ 10 , C ~ 11 , S ~ 11 , C ~ 20 , C ~ 21 , S ~ 21 , C ~ 22 , S ~ 22 , E ] . \\begin{aligned} A=[&\\widetilde{P}_{00}(sin\\beta),\\widetilde{P}_{10}(sin\\beta),\\widetilde{P}_{11}(sin\\beta)\\cdot\\cos(s),\\widetilde{P}_{11}(sin\\beta)\\cdot\\sin(s),\\widetilde{P}_{20}(sin\\beta),\\\\ &\\widetilde{P}_{21}(sin\\beta)\\cdot\\cos(2s),\\widetilde{P}_{21}(sin\\beta)\\cdot\\sin(s),\\widetilde{P}_{22}(sin\\beta)\\cdot\\cos(2s),\\widetilde{P}_{22}(sin\\beta)\\cdot\\\\ &\\sin(2s),-1].\\\\ X=[&\\widetilde{C}_{00},\\widetilde{C}_{10},\\widetilde{C}_{11},\\widetilde{S}_{11},\\widetilde{C}_{20},\\widetilde{C}_{21},\\widetilde{S}_{21},\\widetilde{C}_{22},\\widetilde{S}_{22},E]. \\end{aligned} A=[X=[​P 00​(sinβ),P 10​(sinβ),P 11​(sinβ)⋅cos(s),P 11​(sinβ)⋅sin(s),P 20​(sinβ),P 21​(sinβ)⋅cos(2s),P 21​(sinβ)⋅sin(s),P 22​(sinβ)⋅cos(2s),P 22​(sinβ)⋅sin(2s),−1].C 00​,C 10​,C 11​,S 11​,C 20​,C 21​,S 21​,C 22​,S 22​,E].​

现在，采用间接平差计算球谐系数，见式 ( 9 ) (9) (9)，其中， P P P为权阵，若数据质量相同，可以设置为单位阵。
X = ( A T P A ) − 1 ( A T P L ) X=(A^TPA)^{-1}(A^TPL) X=(ATPA)−1(ATPL)
现在我们得到了球谐系数，因此，可以利用式 ( 7 ) (7) (7)计算 V T E C VTEC VTEC。

参考文献：

[1] Schaer S . Mapping and Predicting the Earth\'s Ionosphere Using the Global Positioning System[J]. Geod Geophys.arb.schweiz, 1999.
[2]李新星. 基于球谐函数构建VTEC模式与精度分析[D]. 中国地震局地震研究所, 2017.