廖彬彬,徐建橋，陳曉東*,孫和平,周江存

1 中國科學院精密測量科學與技術創(chuàng)新研究院 大地測量與地球動力學國家重點實驗室，武漢 430077 2 中國科學院大學，北京 100049

0 引言

月球是離地球最近的天體，是人類邁向太空的第一站.研究月球內(nèi)部結構對人類進行空間探索、了解地月系統(tǒng)的演化、開采并利用月球資源等方面具有重要意義(Smith et al.,2017;Zhao et al.,2018；Sun et al.,2019).20世紀阿波羅計劃在月球表面開展了一系列月震觀測實驗，這些月震觀測數(shù)據(jù)使人類對月球內(nèi)部結構有了初步的認識.早期月震的分析結果表明月球內(nèi)部主要存在以下幾個構造圈層：一個厚度約為50～70 km的月殼、一個厚度較大的月幔且在月幔深處存在一個低速層(Nakamura et al.,1974;Goins et al.,1981).后來月球極移和月球磁場的觀測結果證明了月球存在一個液態(tài)的外核(Williams et al.,2001；Shimizu et al.,2013).在阿波羅計劃之后，高分辨率遙感衛(wèi)星技術的應用以及近年來對早期月震數(shù)據(jù)的重新處理使得對月球結構的認識有了新的突破(Hao et al.,2018).Sugano 和 Heki(2004)利用LP(Lunar Prospector，月球勘探者)衛(wèi)星的重力觀測數(shù)據(jù)推斷月殼厚度約為20～60 km.Khan等(2000)利用非線性貝葉斯方法對阿波羅月震數(shù)據(jù)進行了重新處理，結果表明在月殼內(nèi)深度為20 km以及45±5 km處存在速度間斷面，并且月幔內(nèi)深度為500 km處P波速度與S波速度顯著增加.Lognonné等(2003)利用體波走時數(shù)據(jù)通過隨機搜索方法發(fā)現(xiàn)月幔內(nèi)深度500 km以及750 km處均存在速度間斷面.對于月球更深部分的分層結構，Konopliv等(1998)基于LP衛(wèi)星得到的月球質(zhì)量、平均轉動慣量以及重力數(shù)據(jù)推測：如果月核完全由固態(tài)鐵組成，則半徑為220 km到370 km之間；而如果月核為液態(tài)的Fe-FeS共晶體，則半徑為330～590 km之間.Garcia等(2011)基于大地測量觀測數(shù)據(jù)(包括月球質(zhì)量、轉動慣量以及勒夫數(shù))以及月震數(shù)據(jù)得到了月球結構模型VPREMOON，其中的月核半徑為380±40 km，月核平均密度為5200±1000 kg·m-3.Yan等(2015)通過Monte Carlo方法得到月核半徑為370 km.而Weber等(2011)通過反演得到月核由固態(tài)內(nèi)核以及液態(tài)外核組成，其中固態(tài)內(nèi)核半徑為240 km，液態(tài)外核半徑范圍為240～330 km．Williams等 (2014)以Weber等(2011)月球內(nèi)部結構模型為初始模型，利用GRAIL(Gravity Recovery and Interior Laboratory,重力重建與內(nèi)部結構實驗室)重力場數(shù)據(jù)計算了月球2階固體潮勒夫數(shù)、平均密度以及平均轉動慣量等物理量，通過調(diào)整低速區(qū)彈性參數(shù)、低速區(qū)厚度以及內(nèi)外核厚度，得到了滿足GRAIL觀測數(shù)據(jù)的5個月球內(nèi)部結構模型 (GPM 1—5).然而Matsuyama等(2016)認為Williams等(2014)采用的勒夫數(shù)k2=0.02422±0.000012是粘彈性月球的引潮位響應，經(jīng)過改正后給出了彈性月球的勒夫數(shù)為k2=0.0220±0.0017.并且Matsuyama等(2016)利用Markov Chain-Monte Carlo方法進行反演，結果顯示月幔底部低速區(qū)存在的概率約為50%.綜上所述，從上面目前國內(nèi)外對月球結構研究結果可以看出，現(xiàn)階段對月球結構模型的認識，特別是月球深部結構的認識，仍然存在許多的不確定性，這也是至今沒有一個月球內(nèi)部結構模型被廣泛認可和采用的根本原因.

目前可以用來反演月球內(nèi)部結構的數(shù)學方法有很多，各種方法或多或少都存在一定的局限性.傳統(tǒng)線性反演方法要求模型參數(shù)與觀測數(shù)據(jù)之間的映射是(擬)線性的，然而實際中觀測數(shù)據(jù)與模型的關系往往是非線性的，甚至可能都不是一一對應的映射關系(Tarantola,2005).已有的研究結果表明利用非線性Bayes反演方法可以很好解決傳統(tǒng)線性反演方法的不足(De Wit et al.,2013).Bayes反演方法基于Bayes理論，利用觀測數(shù)據(jù)d、模型的先驗分布σ(m)和似然函數(shù)L(d|m)計算模型參數(shù)后驗概率密度分布σ(m|d)，從而獲得模型的最優(yōu)解，即最大后驗概率對應的模型(De Wit et al.,2013).由于采用了更加復雜的概率模型與參數(shù)估計方法，Bayes反演方法完全擺脫最小二乘方法要求觀測數(shù)據(jù)和模型參數(shù)之間存在線性關系的限制.Bayes反演方法通常利用Monte Carlo方法來實現(xiàn)(K?ufl et al.,2013)，然而Monte Carlo方法受限于“維度詛咒”(Bellman,1961)，往往計算量巨大，需要較長的計算時間.另一種可行的Bayes反演方法是混合密度神經(jīng)網(wǎng)絡方法(Mixture Density Neural Network，MDN)，其基本原理是建立一個以觀測數(shù)據(jù)為輸入，以模型參數(shù)的后驗概率密度分布為輸出的神經(jīng)網(wǎng)絡(Bishop,1995).與Monte Carlo方法相比，MDN方法可以避免“維度詛咒”，從而極大地提高計算效率(De Wit et al.，2013).此外MDN方法還可以得到觀測數(shù)據(jù)與模型參數(shù)之間的函數(shù)關系，為進一步分析月球內(nèi)部結構模型的精度提供了可能.

本研究采用Williams等(2014)基于GRAIL觀測任務給出的平均密度ρ和平均轉動慣量MOI以及Matsuyama等(2016)給出的彈性勒夫數(shù)k2作為觀測約束，通過MDN方法獲得了月球內(nèi)部結構模型的低速層、外內(nèi)核的半徑、P波和S波速度等8個模型參數(shù)的后驗概率密度分布，計算了平均月球內(nèi)部結構模型Mean、最大后驗概率密度對應的MAP模型以及滿足1-σ準則的月球內(nèi)部結構模型的參數(shù)范圍.最后利用觀測值與月球內(nèi)部結構模型的參數(shù)之間的函數(shù)關系，采用不同觀測值觀測誤差范圍內(nèi)月球內(nèi)部結構模型參數(shù)的變化幅度作為度量，定量地分析了不同觀測值的精度對月球內(nèi)部結構模型的影響.

1 數(shù)據(jù)和方法

對月球內(nèi)部結構的反演工作主要由兩部分組成：(1)以Weber等(2011)的月球內(nèi)部結構模型為初始模型，對內(nèi)外核及低速區(qū)的厚度、P波及S波速度進行均勻采樣，得到所有可能的月球內(nèi)部結構模型，然后正演計算得到對應的固體潮勒夫數(shù)、平均密度及平均轉動慣量；(2)利用混合密度神經(jīng)網(wǎng)絡方法和上面計算的三個觀測值，得到對應月球內(nèi)部結構模型的后驗概率密度.

1.1 觀測數(shù)據(jù)的正演計算

假設月球是一個非自轉、完全彈性、各向同性以及橫向均勻的分層球體，并且各個圈層的密度以及彈性參數(shù)均為常數(shù).Matsumoto等(2015)利用Bayes方法反演了月球內(nèi)部結構，他們的結果表明：月幔底部低速區(qū)以及月核的半徑、彈性參數(shù)和密度都存在很大的不確定性，而低速區(qū)以上的部分模型參數(shù)的不確定性相對較小.因此本研究假設在低速區(qū)以上部分月球內(nèi)部結構模型參數(shù)為定值，均采用Weber等(2011)給出月球內(nèi)部結構模型中的數(shù)值，而低速區(qū)以下的模型參數(shù)則為反演的自由度.此外由于觀測值個數(shù)較少，出于簡化反演問題的考慮，低速區(qū)與內(nèi)外核的密度同樣也采用Weber等(2011)給出月球內(nèi)部結構模型中的對應數(shù)值.并且為了與嫦娥1號、Clementine等其他月球任務給出的月球半徑觀測數(shù)據(jù)吻合(Ping et al.,2009;Zuber et al.,2013;Smith et al.,2017)，Williams等(2014)建議將Weber等(2011)的月球內(nèi)部結構模型中的半徑調(diào)整為R=1737.151 km.最終在反演問題中需要考慮的模型參數(shù)總共8個，包括內(nèi)核半徑ri、內(nèi)核P波速度VPi、內(nèi)核S波速度VSi、外核半徑ro、外核P波速度VPo、低速區(qū)半徑rl、低速區(qū)P波速度VPl以及低速區(qū)S波速度VSl，所有參數(shù)先驗的變化范圍如表2中“先驗模型范圍”一列所示.假設8個模型參數(shù)的先驗概率分布滿足均勻分布，在表2給出的先驗模型范圍內(nèi)對它們進行隨機采樣，總共生成100000個月球內(nèi)部結構模型樣本.

研究中觀測值固體潮勒夫數(shù)k2、平均密度ρ以及平均轉動慣量MOI的具體數(shù)值以及誤差范圍如表1所示.本文中固體潮勒夫數(shù)k2采用的是Matsuyama等(2016)給出的彈性勒夫數(shù)k2=0.0220±0.0017.為了獲得上述100000個月球內(nèi)部結構模型樣本的模擬觀測值，首先通過正演數(shù)值計算方法獲得每個月球內(nèi)部結構模型對應固體潮勒夫數(shù)、平均密度以及平均轉動慣量的模擬觀測值.其中在求解二階潮汐勒夫數(shù)時，需要求解一個由動量平衡方程與泊松方程組成的二階偏微分方程組邊值問題，本研究采用廖彬彬等(2019)給出的譜元法進行數(shù)值計算.對于平均密度和平均轉動慣量，則用Gauss-Legendre數(shù)值積分方法進行計算(Quarteroni，2009).此外由于實際觀測值存在誤差，因此需要對上面正演計算得到的模擬觀測值加上觀測誤差才能得到實際觀測值.研究中假設三個觀測值在誤差范圍內(nèi)服從均勻分布，對其進行隨機采樣得到相應的觀測誤差，然后加上對應的觀測值便可得到之后訓練神經(jīng)網(wǎng)絡所需要的模擬觀測值.

表1 反演計算采用的觀測值及其精度Table 1 Observations and their accuracies used in the inversion calculation

1.2 用混合密度神經(jīng)網(wǎng)絡方法反演后驗概率密度

神經(jīng)網(wǎng)絡反演方法主要特點在于它可以解決反演過程中的“非線性”和“不唯一性”等問題.該方法與一般采樣反演方法相比，計算效率更高，關鍵是還可得到模型參數(shù)后驗概率密度分布與觀測值間的映射關系.這種映射關系可用來進一步研究觀測值精度對模型精度的影響，也是本文研究模擬觀測值精度對月球內(nèi)部結構模型精度影響采用的主要數(shù)學方法.神經(jīng)網(wǎng)絡的數(shù)學本質(zhì)是實現(xiàn)從輸入層(Input layer)到輸出層(Output layer)之間任意映射關系的近似表達，其映射關系如圖1所示，可用下面公式(1)和(2)來描述(De Wit et al.,2013)：

圖1 神經(jīng)網(wǎng)絡映射關系示意圖Fig.1 Schematic diagram of the neural network mapping relationship

(1)

(2)

混合密度神經(jīng)網(wǎng)絡是以觀測數(shù)據(jù)為輸入層，以其反演的結構模型后驗概率密度函數(shù)為輸出層的兩層神經(jīng)網(wǎng)絡(Bishop,2006).MDN方法輸出層的后驗概率密度函數(shù)為高斯混合分布，其定義如下：

(3)

式中w表示公式(1)和(2)中的權重和偏置，αi是Softmax函數(shù)，它代表第i個高斯核函數(shù)φi的權重.高斯核函數(shù)φi定義為

φi(m|d;w)=

(4)

式中c為模型參數(shù)m的維度，μi和σi分別為第i個高斯核函數(shù)的均值及標準差，μi、σi和αi可以分別控制第i個高斯核函數(shù)的中心位置、寬窄、高低(權重).不同于傳統(tǒng)線性反演方法中以單個高斯分布模型作為概率模型，本研究月球內(nèi)部結構模型參數(shù)所用的概率模型(見公式(3))是M個高斯分布的線性組合.理論上隨著高斯核函數(shù)個數(shù)M趨近于無窮大，公式(3)給出的概率模型可以近似地表示任意概率分布函數(shù)(Bishops,1995).利用公式(3)給出的概率模型，可實現(xiàn)從一個觀測值到多個可能模型值(平均值μi)之間的映射關系，從而解決了反演的“不唯一性”問題.此外通過混合密度神經(jīng)網(wǎng)絡方法，還實現(xiàn)了由觀測數(shù)據(jù)到概率模型參數(shù)μi、σi和αi的函數(shù)映射，而不僅限于線性映射，因此該方法解決了傳統(tǒng)最小二乘方法只能處理觀測值與模型參數(shù)為線性映射的問題，這也是本研究選用該方法的主要原因.

在使用MDN反演方法的過程中，構建了一個以固體潮勒夫數(shù)、平均密度及平均轉動慣量的觀測值為輸入層，以低速層半徑、內(nèi)外核的厚度、P波和S波速度的聯(lián)合后驗概率分布(即公式(3))的參數(shù)μi、σi和αi為輸出層的兩層神經(jīng)網(wǎng)絡.該神經(jīng)網(wǎng)絡輸入層節(jié)點個數(shù)為I=3，其激活函數(shù)為“Rule”函數(shù)；隱藏層節(jié)點個數(shù)為50，其激活函數(shù)為g(a)=a；而對于輸出層，設高斯核個數(shù)M為10，節(jié)點個數(shù)為K=170.各項超參數(shù)(訓練步數(shù)、批訓練大小、優(yōu)化器類型、激發(fā)函數(shù)類型以及節(jié)點個數(shù)等)的具體選擇方法將在下面進行敘述.

神經(jīng)網(wǎng)絡的訓練過程是一個不斷更新節(jié)點權值和每一層的偏置(即公式(1)和(2)中的w與b)，從而使得損失函數(shù)接近最小的過程.MDN方法的損失函數(shù)loss是后驗分布函數(shù)(公式(3))的似然函數(shù)(Bishop，1995)，它可以用來衡量所輸出的概率模型與訓練數(shù)據(jù)之間的相似程度，其定義如下：

(5)

其中N為訓練樣本總數(shù).訓練所采用的最優(yōu)化方法是隨機梯度下降法，通過這種方法可將訓練集數(shù)據(jù)分批逐次地代入損失函數(shù)當中，而不必一次性計算所有訓練集數(shù)據(jù)損失函數(shù)的和(Goodfellow et al.,2016).本研究采用的隨機梯度下降算法是“Adam”算法.為了避免最優(yōu)化解陷入局部最優(yōu)解，需要對神經(jīng)網(wǎng)絡進行多次訓練，從中選擇損失函數(shù)最小的解，并且在每次訓練之前需要對權值進行不同的初始化(Bishop,1995).此外，為了避免輸出的月球內(nèi)部結構模型中各層半徑出現(xiàn)負值，本文采用log(ri)、log(ro-ri)和log(rl-ro)代替內(nèi)核半徑ri、外核半徑ro以及低速區(qū)半徑rl作為神經(jīng)網(wǎng)絡的輸出.最后通過z-socre標準化方法對觀測值與模型參數(shù)進行預處理可以有效避免訓練過程中出現(xiàn)“梯度飽和”(De Wit et al.,2013).

在神經(jīng)網(wǎng)絡訓練完后，還要對其泛化能力進行評估，神經(jīng)網(wǎng)絡的泛化能力是指其能正確預測未參與訓練數(shù)據(jù)的能力.為了度量神經(jīng)網(wǎng)絡的泛化能力，本研究以8∶1∶1的比例將100000組觀測數(shù)據(jù)與對應模型參數(shù)成對組成的數(shù)據(jù)集分成訓練集、測試集以及驗證集，三者中只有訓練集參與神經(jīng)網(wǎng)絡的訓練，后面兩者作為檢驗標準分別在訓練過程中超參數(shù)的選擇階段和訓練結束后結果評價階段對神經(jīng)網(wǎng)絡的泛化能力進行評估(Goodfellow et al.,2016).與傳統(tǒng)最優(yōu)化問題不同，為了提高訓練得到神經(jīng)網(wǎng)絡的泛化能力，避免“過擬合”的發(fā)生，神經(jīng)網(wǎng)絡最優(yōu)化的目標函數(shù)并非參與訓練的訓練集損失函數(shù)，而是未參與訓練的驗證集與測試集損失函數(shù)，訓練集的損失函數(shù)僅僅為最優(yōu)化過程提供了方向(Goodfellow et al.,2016).在訓練過程中，利用驗證集損失函數(shù)的大小來對訓練過程中的各種超參數(shù)進行抉擇.以訓練步數(shù)的選擇為例，圖2中展示了1000個步長學習過程中訓練集與驗證集損失函數(shù)隨訓練步數(shù)的變化曲線.在上文給出的各種超參數(shù)的條件下，經(jīng)過1000個步長時訓練集和驗證集的損失函數(shù)均趨近于收斂，且無明顯過擬合現(xiàn)象.而經(jīng)過1000個步長后，驗證集的損失函數(shù)不再繼續(xù)減小，因此我們將訓練步長選為1000是合理的.在訓練結束階段，用測試集損失函數(shù)的大小對訓練好的神經(jīng)網(wǎng)絡的泛化能力進行評估，最終得到的神經(jīng)網(wǎng)絡訓練集的平均損失函數(shù)值為-1.88，驗證集的平均損失函數(shù)值為-1.84，測試集的平均損失函數(shù)值為-1.82，三個值相差很小，說明未參與訓練數(shù)據(jù)的損失函數(shù)與參與訓練數(shù)據(jù)的損失函數(shù)相當，足以證明訓練得到的神經(jīng)網(wǎng)絡模型具有足夠強的泛化能力.

圖2 學習曲線(訓練集與驗證集損失函數(shù)隨訓練步數(shù)的變化曲線)Fig.2 Learning curves (loss function′s variation curves of training set and verification set with the number of training steps)

2 計算結果與分析

為了驗證混合密度神經(jīng)網(wǎng)絡(MDN)方法的有效性，需要將MDN方法的結果其與其他獨立方法的結果進行比較，本研究將其與隨機采樣方法的結果進行比較.并且根據(jù)MDN方法得到模型參數(shù)后驗概率密度分布對月球內(nèi)部結構模型進行采樣，并正演計算月球內(nèi)部結構模型參數(shù)樣本對應的模擬觀測值，從而分析了MDN方法的結果與實際觀測值誤差的符合程度．然后基于模型參數(shù)的后驗概率密度分布，計算平均月球內(nèi)部結構模型(Mean模型)、最大后驗概率密度對應的模型(MAP模型)以及基于1-σ準則給出的月球內(nèi)部結構模型(1-σ模型).最后依據(jù)得到的不同觀測值與模型參數(shù)間的函數(shù)關系，并利用模型參數(shù)在不同觀測值誤差范圍內(nèi)的變化幅度為度量，研究了不同觀測值的觀測精度對月球內(nèi)部結構模型的影響.

2.1 MDN方法計算結果有效性的獨立驗證

由于似然函數(shù)不是一個絕對的評價指標，其對應的損失函數(shù)值為-1.82，并不能說明MDN方法的反演結果能有效地反映觀測數(shù)據(jù)中所包含月球內(nèi)部結構的全部信息.因此為了驗證MDN方法結果的有效性，需要將MDN方法計算的結果與其他獨立方法的結果進行比較.本研究用來對比的獨立方法是隨機采樣方法，其主要步驟如下：(1)利用隨機采樣方法獲得100000個月球結構模型樣本，它們參數(shù)的變化范圍如表2所示；(2)利用1.1節(jié)所述正演方法計算每個樣本對應的觀測值；(3)假設這些觀測值在各自誤差范圍內(nèi)滿足均勻分布，對于步驟(2)中觀測值落在誤差范圍內(nèi)的樣本加以保留，對于落入誤差范圍外的樣本則拒絕，從而獲得符合觀測誤差的月球結構模型樣本.所采用的觀測值及其誤差范圍如表1所示.

圖3給出了MDN方法計算結果與隨機采樣方法獲得的月球內(nèi)部結構模型樣本的比較.圖3中藍色為MDN方法的結果，其右側和上側是單個月球內(nèi)部結構模型參數(shù)的邊緣概率密度分布.由于隨機采樣方法受限于“維度詛咒”，最終100000個樣本中符合觀測誤差的樣本只有36個，在圖中用紅色“×”表示，并且這36個樣本大部分都落在MDN方法聯(lián)合概率密度分布的高密度區(qū)域，因此MDN方法的計算結果能夠有效地表達隨機采樣方法結果的主要特征.以內(nèi)核半徑與外核半徑聯(lián)合概率密度分布為例，MDN方法得到二者的皮爾遜相關系數(shù)(Pearsonr)為-0.95，p值(顯著性水平)為0，隨機采樣方法得到二者皮爾遜相關系數(shù)為-0.91，p值為5.4×10-24；由于二者的皮爾遜相關系數(shù)均接近-1，說明二者存在線性負相關關系，而p值均遠小于0.05說明這種線性相關是顯著的.MDN方法與隨機采樣方法的比較結果可證明MDN方法能有效解決月球內(nèi)結構模型的反演問題，得到的解是非平凡解.

圖3 MDN方法結果與隨機采樣方法結果的比較藍色為MDN方法的結果，紅色叉號為隨機采樣方法的結果．(a)內(nèi)核半徑與外核半徑聯(lián)合概率分布；(b)外核半徑與低速區(qū)半徑聯(lián)合概率分布；(c)低速區(qū)P波與S波速度聯(lián)合概率密度分布．Fig.3 Comparison of results obtained by the MDN method and the random sampling methodThe blue areas are results of the MDN method,and the red crosses are results of the random sampling method.(a)Joint probability distribution of the inner core radius and the outer core radius.(b)Joint probability distribution of radii of the outer core and the LVZ (Low-Velocity Zone).(c)Joint probability distribution of P-wave and S-wave velocities of LVZ.

為了評價MDN方法得到月球內(nèi)部結構模型參數(shù)的精度，需要研究這些月球內(nèi)部結構模型對應的模擬觀測值與實際觀測值的符合程度.對MDN方法得到月球內(nèi)部結構模型后驗概率分布進行1000次獨立采樣，同樣利用正演方法計算它們對應的模擬觀測值.圖4給出了這些模擬觀測值與實際觀測值的比較，其中橙色豎線是實際觀測值，橙色橫線是其誤差范圍.圖4a中平均密度模擬觀測值的平均值為3345.64 kg·m-3,標準差為0.31 kg·m-3，與實際觀測值3345.56±0.4 kg·m-3相比差別非常小，模擬觀測值的平均值偏差與標準差的和(0.39 kg·m-3)在實際觀測值的誤差范圍以內(nèi)，說明反演結果對應的模擬觀測值相對于與實際觀測值的偏離程度以及反演結果的離散程度在誤差允許范圍內(nèi).圖4b中平均轉動慣量MOI的平均模擬觀測值為0.393100，標準差為0.000029，其實際觀測值為0.393112±0.000012，因此反演結果對應的平均值相對于實際觀測值的偏差(0.000012)在誤差范圍內(nèi)，然而反演結果對應的標準差要大于1倍中誤差.圖4c中勒夫數(shù)k2模擬觀測值的平均值為0.0228，標準差為0.0008，其實際觀測值為0.0220±0.0017；與實際觀測值相比，模擬觀測值平均值有0.0008的偏離，但是這個偏離加上標準差(0.0016)小于實際觀測值的一倍中誤差.造成三個模擬觀測值精度與實際精度不一致的主要原因是MDN方法的目標函數(shù)是月球內(nèi)部結構模型參數(shù)的似然函數(shù)，如公式(5)所示，并未利用任何觀測值的精度信息.

圖4 MDN方法獲得的月球內(nèi)部結構模型對應模擬觀測值的統(tǒng)計直方圖橙色實線表示實際觀測值以及其誤差范圍．(a)平均密度;(b)平均轉動慣量;(c)勒夫數(shù)k2.Fig.4 The statistical histogram of the simulated observation values corresponding to the lunar model obtained by the MDN methodThe orange solid lines are the real observations and their error range.(a)Mean density;(b)Mean solid MOI;(c)Love number k2.

MDN方法的反演結果是月球內(nèi)部結構模型的后驗概率密度分布，對于所有滿足該分布的月球內(nèi)部結構模型的模擬觀測值而言，它們與觀測值中的平均密度以及勒夫數(shù)k2的符合程度較好；而對于觀測值中的平均轉動慣量而言，它們的平均值的符合程度較好，而標準差則稍大于實際誤差范圍.但是本文給出的平均月球內(nèi)部結構模型(Mean模型)和最大后驗概率密度模型(MAP模型)對三個實際觀測值觀測誤差都符合，進一步說明本文給出的最優(yōu)化月球內(nèi)部結構模型滿足了所有觀測值的觀測誤差.

2.2 月球內(nèi)部結構模型的反演結果

表2給出了月球內(nèi)部結構模型參數(shù)先驗變化區(qū)間以及反演結果，其中Mean模型是服從后驗概率密度分布月球內(nèi)部結構模型的期望，它是所有反演得到月球內(nèi)部結構模型的平均值；MAP模型是最大后驗概率密度對應的月球結構模型，它是所有反演得到月球內(nèi)部結構模型當中可能性最大的模型.兩個模型相差很小，原因是反演得到的后驗概率密度分布是一個簡單的單峰分布，因此可以將MAP模型當作反演得到的最優(yōu)化月球結構模型.“1-σ”一列是在1-σ準則下月球內(nèi)部結構模型參數(shù)的變化范圍，它對應著在Mean模型參數(shù)附近概率密度積分為68.26%的區(qū)域內(nèi)所有月球內(nèi)部結構模型.

表2 月球內(nèi)部結構模型參數(shù)先驗變化區(qū)間和反演結果Table 2 The priori variation intervals and inversion results of the lunar internal structure model parameters

圖5是月球內(nèi)部結構模型中各層半徑、P波和S波速度的垂直剖面圖，其中灰色線條表示W(wǎng)eber等(2011)的月球內(nèi)部結構模型，藍色虛線表示表2中的MAP模型.在圖5a中的紅色誤差棒表示表2中1-σ模型中的低速區(qū)、外核以及內(nèi)核半徑的變化范圍，它們值分別為212～251 km、305～336 km、334～467 km.Matsuyama等(2016)給出1-σ準則下的半徑變化范圍分別為6～275 km，194～362 km，345～505 km.本文結果與他們的結果相比，給出的變化區(qū)間更窄，這是由于本文假設月球結構模型中的密度以及低速區(qū)LVZ以外圈層的其他模型參數(shù)均為定值，而在Matsuyama等(2016)的工作中這些參數(shù)都是變量，因此模型的不確定性更大.圖5(b，c)表示P波速度VP以及S波速度VS的垂直剖面，其中紅色虛線的區(qū)間表示表2中1-σ準則下P波速度與S波速度的變化范圍.Matsumoto等(2015)給出低速區(qū)內(nèi)的P波及S波速度分別為7.1±1.0 km·s-1、2.9±0.5 km·s-1；本文1-σ準則下給出的6.49～7.70 km·s-1、2.61～3.21 km·s-1與Matsumoto等(2015)的結果相比區(qū)間更窄，原因與上面半徑類似.MDN方法的反演結果與前人的反演結果基本吻合，但是由于簡化了反演參數(shù)個數(shù)，因此本文反演的模型不確定度小于前人的結果.此外Matsuyama等(2016)的結果顯示低速區(qū)LVZ存在的可能性約為50%，本文結果給出的低速區(qū)S波速度為2.63 ～ 3.21 km·s-1，明顯低于Weber等(2011)的月幔S波速度4.50 km·s-1，因此本文結果支持月幔底部存在低速區(qū)的觀點.

圖5 月球內(nèi)部結構模型各層半徑(a)、P波(b)及S波(c)速度的垂直剖面灰色實線為Weber等(2011)的月球內(nèi)部結構模型，左圖中紅色的誤差棒和右邊兩圖的紅色虛線為本研究中1-σ準則下月球內(nèi)部結構模型變化區(qū)間，藍色虛線為MAP模型.Fig.5 Vertical profiles of the radii (a),P-wave (b)and S-wave (c)velocities of lunar modeled layersThe gray solid lines are the lunar model given by Weber et al.(2011).The red error bars in the left figure and the red dotted lines in the two right figures are the ranges of the lunar model under the 1-σ criterion in this study.

利用模型參數(shù)間的聯(lián)合概率密度分布可以對模型參數(shù)之間的相關性進行研究.圖6給出了利用MDN方法得到的內(nèi)外核半徑以及外核半徑與低速區(qū)半徑的聯(lián)合概率密度分布，其中黃色菱形表示表2中MAP模型,紅色五角星則對應于Williams等(2014)的月球內(nèi)部結構模型GPM 1—5，其數(shù)值見表2.如2.1節(jié)所述，內(nèi)外核半徑之間存在明顯的線性負相關關系，這是平均密度、平均轉動慣量以及勒夫數(shù)k2的觀測值與誤差范圍共同作用的結果.從2.3節(jié)的研究結果來看，在三個觀測值現(xiàn)有的觀測誤差范圍內(nèi)，k2允許模型參數(shù)在較大的范圍內(nèi)變化，因此它對模型參數(shù)的約束較弱，因此月球內(nèi)部結構模型的變化范圍主要受到平均密度和平均轉動慣量約束.在本研究問題的假設背景下，低速區(qū)密度(3.4 kg·m-3)與月幔相同，因此低速區(qū)半徑的改變對平均密度和平均轉動慣量沒有任何影響，而月核中內(nèi)核密度(8.0 kg·m-3)大于外核密度(5.1 kg·m-3).因此為了保持平均密度不變，二者的半徑必然存在線性負相關關系；此外對于平均轉動慣量，由于采用的是歸一化的數(shù)值，因此其變化趨勢與平均密度一致.內(nèi)外核半徑聯(lián)合概率密度分布的結果與GPM 2—5模型符合程度較好，并且Mean模型和MAP模型幾乎與GPM 3模型重合.而月球內(nèi)部結構模型中的低速區(qū)半徑低于GPM 1—5模型，造成這種差異的主要原因是所采用的固體潮勒夫數(shù)k2為Matsuyama等(2016)給出的0.0220±0.0017，而非Williams等(2014)采用的0.02422±0.00022,該結果進一步說明了固體潮勒夫數(shù)k2的觀測值及其精度對月球內(nèi)部結構模型影響較為顯著.

圖6 月球內(nèi)部結構模型參數(shù)聯(lián)合概率密度分布(a)內(nèi)核半徑與外核半徑；(b)外核半徑與低速區(qū)半徑.黃色菱形為MAP模型,紅色五角星分別對應了月球內(nèi)部結構模型GPM 1—5．Fig.6 Joint probability density distribution (JPDD)of lunar model parameters(a)The JPDD of the inner core radius and the outer core radius；(b)The JPDD of the outer core radius and the LVZ radius.The yellow diamond refers to the MAP model and the red stars are lunar models GPM 1—5,respectively.

2.3 觀測值精度對反演的月球內(nèi)部結構模型的影響

與其他反演方法相比，MDN方法除了可以得到月球內(nèi)部結構模型的后驗概率密度分布之外，還可以得到觀測值與月球內(nèi)部結構模型參數(shù)的函數(shù)關系，這種函數(shù)關系為我們評估觀測值精度對所反演月球內(nèi)部結構模型的影響提供了可能.

圖7給出了在平均密度與平均轉動慣量為定值(取值如表1所示)條件下，勒夫數(shù)k2與月球內(nèi)部結構模型參數(shù)之間的關系.橙色實線為MDN方法得到的Mean模型與勒夫數(shù)k2的函數(shù)關系，藍色實線是利用正演數(shù)值計算方法得到二者函數(shù)關系.MDN方法與正演數(shù)值計算方法結果的差異反映了MDN方法作為一種函數(shù)近似方法的近似誤差.在一倍觀測誤差(0.0017)以內(nèi)，函數(shù)的梯度有著明顯的變化，模型參數(shù)與勒夫數(shù)k2之間函數(shù)關系呈現(xiàn)出明顯的非線性．在Williams等(2014)給出的固體潮勒夫數(shù)k2=0.02422±0.00022附近，內(nèi)核半徑與外核半徑的斜率趨近于0，而低速區(qū)半徑的斜率則是一個較大的正數(shù)，因此勒夫數(shù)k2在誤差范圍內(nèi)的變化對低速區(qū)半徑的精度影響較大，而對內(nèi)外核半徑的精度影響較小.而對于Matsuyama等(2016)給出的固體潮勒夫數(shù)0.0220±0.0017而言，勒夫數(shù)k2的觀測誤差對上述三個模型參數(shù)精度的影響程度相當.此外由于內(nèi)外核半徑以及低速區(qū)半徑存在ri

圖7 月球內(nèi)部結構模型參數(shù)與勒夫數(shù)k2的函數(shù)關系(a)內(nèi)核半徑；(b)外核半徑；(c)低速區(qū)半徑.Fig.7 Function relationship between lunar model parameters and the Love number k2(a)Radius of the inner core；(b)Radius of the outer core；(c)Radius of the LVZ.

圖8給出了觀測誤差為0時，利用MDN方法以及表1中的觀測數(shù)據(jù)得到月球內(nèi)部結構模型的后驗概率密度分布，其中紅色五角星為GPM1～5模型.圖8中內(nèi)外核半徑以及低速區(qū)半徑均在275 km附近，而2.2節(jié)中的平均月球內(nèi)部結構模型Mean內(nèi)核半徑為231 km，外核半徑為321 km，低速區(qū)半徑為414 km.造成二者差異的原因主要是因為觀測值k2=0.0220趨近于其取值范圍的下邊界，因此對于Matsuyama等(2016)給出的固體潮勒夫數(shù)k2=0.0220±0.0017,其有效變化范圍為k2=0.0220～0.0237，如果k2為0.0220時反演結果表明內(nèi)核半徑為275 km而外核與低速區(qū)厚度為0，而如果k2為0.0237時反演得到內(nèi)核半徑為218 km，外核半徑為318 km,低速區(qū)半徑為480 km；在假設k2服從均勻分布的條件，其平均值為0.02285，則反演結果便如圖6中所給出.這也是圖4反演結果中勒夫數(shù)k2模擬觀測值平均值相對于實際觀測值有明顯偏離的主要原因.

表3 觀測值在其誤差范圍內(nèi)變化時導致的月球內(nèi)部結構模型參數(shù)變化幅度
Table 3 Variation range of lunar model parameters caused by variations of observations within their errors

觀測值精度對月球內(nèi)部結構模型的影響主要由兩個方面決定：第一個方面是問題的物理背景，月球內(nèi)部結構模型參數(shù)通過正演過程得到對應的模擬觀測值，因此模擬觀測值與模型參數(shù)之間通過控制方程建立起如圖7所示的函數(shù)關系.在不同觀測值附近函數(shù)斜率不同，因此月球內(nèi)部結構模型對模擬觀測值的敏感程度也不同.第二個方面是實際觀測值的概率密度分布.如果實際觀測值觀測誤差為0，那么反演得到的模型參數(shù)會是一個準確值；而當用概率模型描述觀測值的精度時，實際觀測值并不是常數(shù)而是一個隨機變量.實際觀測值作為隨機變量在誤差范圍內(nèi)變化時，對應的模型參數(shù)會在函數(shù)的某一段區(qū)間內(nèi)上下浮動，其變化幅度的大小取決于誤差范圍的大小以及函數(shù)在各點上的斜率.

從圖7可以看出勒夫數(shù)k2與模型參數(shù)的函數(shù)在其誤差范圍內(nèi)具有明顯的非線性，因此利用函數(shù)在觀測值處的斜率分析模型參數(shù)可能的變化范圍是不合理的.此外考慮到其他模擬觀測值與模型參數(shù)的函數(shù)可能是非單調(diào)性的，因此用在觀測值誤差范圍內(nèi)模型參數(shù)最大值與最小值的差來度量模型參數(shù)的變化范圍是一種更為合理的方法.為了比較不同觀測值對模型精度的影響，在其他兩個觀測值為確定值的條件下，計算單個觀測值在其誤差范圍內(nèi)變化時對應月球內(nèi)部結構模型參數(shù)的變化幅度.表3給出了三個觀測值在各自一倍中誤差范圍內(nèi)對應月球內(nèi)部結構模型參數(shù)的變化幅度.表3結果顯示在觀測值的誤差范圍內(nèi)，勒夫數(shù)k2對應的月球內(nèi)部結構模型參數(shù)變化幅度遠遠大于其他兩個觀測值，例如對于低速區(qū)S波速度而言，勒夫數(shù)k2的影響幅度為202.79 m·s-1，約是MOI影響幅度1.91 m·s-1的200倍；對于其他月球內(nèi)部結構模型參數(shù)而言，k2的影響同樣遠大于MOI和平均密度的影響.

內(nèi)核外核低速區(qū)半徑ri(km)P波速度VPi(m·s-1)S波速度VSi(m·s-1)半徑ro(km)P波速度VPo(m·s-1)半徑rl(km)P波速度VPl(m·s-1)S波速度VSl(m·s-1)勒夫數(shù)k260.5952.9735.3648.4319.46198.50105.42202.79平均轉動慣量MOI0.357.681.430.281.030.354.691.91平均密度ρ(kg·m-3)2.5927.66.072.183.62.3821.03.21

3 結論

利用混合密度神經(jīng)網(wǎng)絡(MDN)方法以Williams等(2014)基于GRAIL任務給出的月球平均密度和平均轉動慣量、Matsuyama等(2016)給出的月球彈性勒夫數(shù)k2作為觀測約束，對月球內(nèi)外核以及低速區(qū)的半徑、P波和S波速度共計8個參數(shù)進行了反演，得到了月球內(nèi)部結構模型的后驗概率密度分布，獲得了最優(yōu)化月球內(nèi)部結構模型，即最大后驗概率密度對應的MAP模型.除此之外還計算了后驗概率密度分布的平均月球內(nèi)部結構模型(Mean模型)以及基于1-σ準則給出的月球內(nèi)部結構模型范圍.Mean模型與MAP模型差別較小且后驗概率密度分布是單峰分布說明以MAP模型作為最優(yōu)化月球內(nèi)部結構模型是合理的.反演得到的月球內(nèi)部結構模型中，內(nèi)核半徑與外核半徑呈現(xiàn)出明顯的線性負相關關系，這是因為在三個觀測值觀測精度內(nèi)勒夫數(shù)k2允許模型參數(shù)在較大的區(qū)間內(nèi)變化，而與之相比平均密度和平均轉動慣量在誤差允許范圍內(nèi)幾乎可以認為是不變的.反演結果中低速區(qū)S波速度(2.63～3.2 km·s-1)明顯低于月幔S波速度(4.50 km·s-1)，因此本文結果支持月幔底部存在低速區(qū).

相比于其他反演方法，MDN方法的優(yōu)勢主要體現(xiàn)在三個方面：(1)與傳統(tǒng)的最小二乘法反演方法相比，MDN方法可以解決反演的不唯一性與非線性性；(2)與反演中經(jīng)常采用的蒙特卡洛方法相比，MDN方法避免了“維度詛咒”，因此計算效率更高；(3)MDN方法不僅可以獲得一個實際觀測值對應的月球內(nèi)部結構模型參數(shù)，還可以獲得實際觀測值與月球內(nèi)部結構模型參數(shù)之間的函數(shù)關系，為進一步分析不同觀測值觀測誤差對月球內(nèi)部結構模型的影響提供了可能.

通過研究勒夫數(shù)k2與模型參數(shù)間的函數(shù)關系，數(shù)值結果表明勒夫數(shù)k2存在一個約為0.022的下限值.利用在觀測誤差范圍內(nèi)模型參數(shù)的變化分布作為度量，研究了不同觀測值精度對模型的影響，結果顯示現(xiàn)有觀測精度條件下勒夫數(shù)k2的精度對月球內(nèi)部結構模型的影響明顯要大于其他兩個觀測值.在將來的月球觀測任務中，期望能首先考慮提高固體潮勒夫數(shù)的觀測精度，這樣就能對月球內(nèi)部結構進行更加準確的約束.