矩諧分析表達區域地球重力場的基本思想是在局部直角坐標系下求解拉普拉斯方程， 將引力位展開成矩諧級數， 以矩諧係數的集合來表徵引力位。矩諧分析是在局部直角坐標系中進行的， 因此需選取數據區域的中心點作為局部直角坐標系的原點，將輸入數據點的大地坐標或球坐標轉換到局部直角坐標系中。
基本介紹
- 中文名：矩諧分析
- 外文名：Moment of harmonic analysis
矩諧分析數學模型,數值計算與分析,結論,
經濟高效的航空重力測量能夠很好地填補衛星重力和地面重力測量所得重力場信息的信號盲區， 可覆蓋地面重力測量難以達到的地區， 數據精度在全波長為5～10km 的空間解析度尺度上能夠達到1～3mGal （ 1Gal=1cm / s 2）。航空重力測量獲取的是距地面達數千米飛行高度處的重力觀測值， 大地測量和地球物理等領域需要的是地面或大地水準面上的重力值， 因此需將飛行高度處的重力觀測值向下調和延拓至地面或大地水準面上。
在航空重力觀測值向下延拓過程中噪聲會得到放大， 可能造成重力場信號的嚴重失真， 若不採用合適的向下延拓方法， 將無法得到穩定可靠的重力延拓解。現有的向下延拓方法主要包括逆泊松積分法、 最 小 二 乘 配 置（ LSC）、 快 速 傅 里葉變換（ FFT）、 直接代表法和直接法，其中比較常用的是逆泊松積分法， 求解泊松積分方程時一般需引入正則化處理。一種基於矩諧分析的航空重力向下延拓新方法， 矩諧分析最開始被文獻用於區域地磁場建模， 在地磁學研究中套用廣泛。文獻將矩諧分析用於局部重力場逼近， 數值結果驗證了其可行性。該方法的基本原理是在局部直角坐標系中將擾動位或其泛函展開成矩諧級數， 利用飛行高度處的重力異常或重力擾動進行矩諧分析以求解表征該區域內重力場的擾動位係數， 再進行矩諧綜合計算得到地面或大地水準面上的重力異常和重力擾動， 從而完成重力觀測值從飛行高度處到地面或大地水準面的向下調和延拓。
在航空重力觀測值向下延拓過程中噪聲會得到放大， 可能造成重力場信號的嚴重失真， 若不採用合適的向下延拓方法， 將無法得到穩定可靠的重力延拓解。現有的向下延拓方法主要包括逆泊松積分法、 最 小 二 乘 配 置（ LSC）、 快 速 傅 里葉變換（ FFT）、 直接代表法和直接法，其中比較常用的是逆泊松積分法， 求解泊松積分方程時一般需引入正則化處理。一種基於矩諧分析的航空重力向下延拓新方法， 矩諧分析最開始被文獻用於區域地磁場建模， 在地磁學研究中套用廣泛。文獻將矩諧分析用於局部重力場逼近， 數值結果驗證了其可行性。該方法的基本原理是在局部直角坐標系中將擾動位或其泛函展開成矩諧級數， 利用飛行高度處的重力異常或重力擾動進行矩諧分析以求解表征該區域內重力場的擾動位係數， 再進行矩諧綜合計算得到地面或大地水準面上的重力異常和重力擾動， 從而完成重力觀測值從飛行高度處到地面或大地水準面的向下調和延拓。
矩諧分析數學模型
諧展開模型採用周期性函式傅立葉級數表達重力場信號， 其成立的前提是假定待求信號在計算區域內也是周期性的。顯然有限區域的重力場信號不具有完全的周期性， 用周期函式表示非周期性信號， 在區域的邊界處， 重力位近似等於兩邊界處重力位的平均值， 會產生周期延拓邊界效應。在矩諧分析中， 擴展計算區域的範圍， 使計算區域內的待求重力場信號滿足周期性條件，可降低周期延拓邊界效應的影響。
航空重力數據在向下延拓過程中觀測噪聲會被放大。隨著數據格網間隔的減小和延拓高度的增加， 設計矩陣A的復共線性增強， 有可能呈病態甚至奇異， 若採用經典最小二乘方法求解，即使較小的觀測誤差也會導致解的嚴重失真甚至錯誤。此外， 在矩諧分析中引入擴展參數Δ 也會使法方程的條件數增加， 加重法方程的病態性， 同樣會影響待求擾動 位 系 數 的 求 解 精 度。因此， 必要時需引入正則化方法求解病態法方程， 以獲取待求位係數的最優解， 可採用 Tikhonov正則化方法， 則正則化解為＾ x=（ A TPA+α I） －1A TP y
式中， I為單位陣； α>0為正則 化 參 數， 至 於 最 優正則化參數的確定可採用廣義交叉檢驗（ GCV）。利用飛行高度處的重力擾動觀測值進行矩諧分析求得表示該區域內重力場的擾動位係數後，再按照式（ 1） 進行矩諧綜合計算即可得到大地水準面上的重力擾動， 至此便完成了航空重力觀測值的向下延拓。
式中， I為單位陣； α>0為正則 化 參 數， 至 於 最 優正則化參數的確定可採用廣義交叉檢驗（ GCV）。利用飛行高度處的重力擾動觀測值進行矩諧分析求得表示該區域內重力場的擾動位係數後，再按照式（ 1） 進行矩諧綜合計算即可得到大地水準面上的重力擾動， 至此便完成了航空重力觀測值的向下延拓。
數值計算與分析
為評價基於矩諧分析的向下延拓方法的可靠性、 精度和穩定性，利用 EGM2008重力位模型設計模擬數值試驗，將矩諧分析與直接法和基於廣義嶺估計的逆泊松積分法這兩種下延方法進行數值比較與分析。首先利用 EGM2008重力位模型的2～2190階位係數計算得到飛行高度處的重力擾動， 在模擬觀測值中加入數學期望為0、 標準差為1．5mGal的 高 斯 白 噪 聲， 移 去 由 GOCE衛星重力位模型2～120 階位係數計算的參考重力擾動得到殘餘重力擾動， 將殘餘重力擾動向下延拓至大地水準面得到延拓值， 再恢復參考模型值得到大地水準面上的重力擾動延拓值。將由EGM2008模型2～2190階位係數計算的大地水準面上的重力擾動作為真實值， 比較重力擾動延拓值與真實值， 可以檢驗向下延拓方法的可靠性與精度並比較其優劣。為模擬真實航空重力測量條件， 選定4個飛行高度 H（ 正高）， 分別是2km、3km、 4km、 5km， 數據點範圍為緯線[ 30°， 33°] 和經線[ 107°， 110°] 圍成的3°×3°區域， 格網間隔取2．5′×2．5′， 共 包 含5184 個 數 據 點。分 別 采 用 矩 諧 分 析、直接法和基於廣義嶺估計的逆泊松積分法， 將飛行高度處的重力擾動觀測值向下延拓至大地水準面。直接法和逆泊松積分法採用 EGM96重力位模 型 計 算 遠 區 流 動 點 的 貢 獻， 最 大 階 數 截 斷 至360階， 積分球 冠 區 半 徑 選 為1°。經 多 次 對 比 測試，確定了矩諧分析的最優展開階數和擴展參數，展開階數 取 為 N =M =20， 擴 展 參 數 取 為 Δ=50km， 選擇 GCV 方法確定正則化參數。
由於計算區域外沒有數據點分布， 直接法、 逆泊松積分法和矩諧分析的向下延拓結果均會受到邊界效應影響， 矩諧分析中還存在周期延拓邊界效應， 首 先 對 比 這 3 種 下 延 方 法 的 邊 界 效 應。直接法、 基於廣 義 嶺 估 計 的 逆 泊 松 積 分 法和矩諧分析的向下延拓誤差分布圖， 延拓高度為4km， 計算點範圍與數據點範圍保持一致。可以看出， 在這3種方法中， 矩諧分析的下延結果受邊界效應的影響最小， 誤差極值點的數量和量級是最小的， 主要分布在距四周邊界15′以內， 中心區域（ 圖2中虛線圍成的2°×2° 矩形） 計算點的結果不受邊 界 效 應 影 響。直 接 法 的 邊 界 效 應 最 為 嚴重， 誤差極值點的數量和量級均大大超過其他兩種方法， 即使在中心區域誤差也非常明顯。基於廣義嶺估計的逆泊松積分法的邊界效應介於矩諧分析和直接法之間， 中心區域的計算點基本不受邊界效應影響。
為消除邊界效應 的 影 響， 應 以2°×2°中 心 區域的計算結果來評定3 種向下延拓方法的精度，共包含2304個計算點。 取 不 同延拓高度時直接法、 逆泊松積分法和矩諧分析的向下延拓誤差統計信息， 向下延拓誤差定義為延拓值與 由 EGM2008 模 型 計 算 的 真 實 值 之 差， 其中誤差比定義為ρ=RMSe（ 11）
式中， RMS為延拓值相對於真實值的均方根誤差；e為白噪聲的標準差。可以看出， 在不同延拓高度處， 矩諧分析所得延拓解的精度都是最高的， RMS隨延拓高度的增大幅度也是最小的， 最大誤差比僅為2．12， 噪聲放大程度遠低於直接法， 略低於基於廣義嶺估計的逆泊松積分法。當延拓高度為5km時， 矩諧分析解的RMS為3．183mGal ， 占直接法所得延拓值 RMS（ 5．665mGal ） 的56%， 略低於逆泊松積分法所得延拓解的 RMS（ 3．245mGal ）。2°×2° 中心區域內3種方法的向下延拓誤差分布圖， 延拓高度為4km。從中可以看出， 矩諧分析的向下延拓誤差均在±10mGal以內， 大部分計算點處的誤差不超過±3mGal ， 其空間分布與逆泊松積分法相似， 但誤差量級略低。相比之下， 直接法的向下延拓效果最差， 中心區域受到邊界效應的影， 出現了較多誤差極值點。
為消除邊界效應 的 影 響， 應 以2°×2°中 心 區域的計算結果來評定3 種向下延拓方法的精度，共包含2304個計算點。 取 不 同延拓高度時直接法、 逆泊松積分法和矩諧分析的向下延拓誤差統計信息， 向下延拓誤差定義為延拓值與 由 EGM2008 模 型 計 算 的 真 實 值 之 差， 其中誤差比定義為ρ=RMSe（ 11）
式中， RMS為延拓值相對於真實值的均方根誤差；e為白噪聲的標準差。可以看出， 在不同延拓高度處， 矩諧分析所得延拓解的精度都是最高的， RMS隨延拓高度的增大幅度也是最小的， 最大誤差比僅為2．12， 噪聲放大程度遠低於直接法， 略低於基於廣義嶺估計的逆泊松積分法。當延拓高度為5km時， 矩諧分析解的RMS為3．183mGal ， 占直接法所得延拓值 RMS（ 5．665mGal ） 的56%， 略低於逆泊松積分法所得延拓解的 RMS（ 3．245mGal ）。2°×2° 中心區域內3種方法的向下延拓誤差分布圖， 延拓高度為4km。從中可以看出， 矩諧分析的向下延拓誤差均在±10mGal以內， 大部分計算點處的誤差不超過±3mGal ， 其空間分布與逆泊松積分法相似， 但誤差量級略低。相比之下， 直接法的向下延拓效果最差， 中心區域受到邊界效應的影， 出現了較多誤差極值點。
結論
矩諧分析通過在局部直角坐標系下求解關於引力位的拉普拉斯方程， 建立位係數與重力觀測值之間的函式關係， 其推導過程嚴格滿足地球重力場的位理論， 具有明確的物理意義。實際計算時矩諧函式展開至一定的截斷階次， 求得表征區域重力場結構的有限階次矩諧係數， 在譜域內對重力場信號具有一定的平滑作用， 在求解關於擾動位係數的線性方程組時， 又引入了正則化處理， 故矩諧分析能夠有效抑制向下延拓過程中觀測噪聲的放大。
直接法和逆泊松積分法屬於空域內的向下延拓方法， 由於其函式模型的固有局限性， 下延結果的邊界效應比較明顯， 特別是直接法的結果受邊界效應影響非常嚴重。矩諧分析方法同樣存在邊界效應， 但擴展計算區域的範圍， 求解法方程時進行正則化處理， 可有效控制邊界效應的量級和影響範圍， 使得中心區域不致受到邊界效應的影響。基於 EGM2008重力位模型設計的空中重力數據向下延拓數值試驗表明： 在下延精度、 穩定性和邊界效應等方面， 矩諧分析都要優於直接法和基於廣義嶺估計的逆泊松積分法，能夠實現航空重力數據的穩定下延， 是一種可靠的航空重力向下延拓方法。
直接法和逆泊松積分法屬於空域內的向下延拓方法， 由於其函式模型的固有局限性， 下延結果的邊界效應比較明顯， 特別是直接法的結果受邊界效應影響非常嚴重。矩諧分析方法同樣存在邊界效應， 但擴展計算區域的範圍， 求解法方程時進行正則化處理， 可有效控制邊界效應的量級和影響範圍， 使得中心區域不致受到邊界效應的影響。基於 EGM2008重力位模型設計的空中重力數據向下延拓數值試驗表明： 在下延精度、 穩定性和邊界效應等方面， 矩諧分析都要優於直接法和基於廣義嶺估計的逆泊松積分法，能夠實現航空重力數據的穩定下延， 是一種可靠的航空重力向下延拓方法。
