顧勇為，歸慶明，張 璇，魏 萌

1.信息工程大學 理學院，河南 鄭州450001；2.中船重工集團 第七一三研究所，河南 鄭州450015

1 引 言

求解大規(guī)模超定線性方程組是控制網(wǎng)平差、GPS導航定位、重力場向下延拓等大地測量與地球物理應用中頻繁出現(xiàn)的問題，共軛梯度法（conjugate gradient method，CG）是目前求解此類方程組的有效方法之一，該方法在保證計算精度的前提下可顯著提高運算速度，且只需要幾百MB的內存空間，因此，應用十分廣泛［1－5］。然而，以上應用領域呈現(xiàn)的問題又往往具有病態(tài)性，在相應的線性方程組中則呈現(xiàn)出系數(shù)矩陣具有嚴重的復共線性特征，此時單純用CG方法解算將嚴重失真［1－2］。正則化方法是適合于消除或減弱復共線性危害的解算病態(tài)問題的一種主要方法，長期以來受到國內外學者的廣泛關注，其研究成果不斷涌現(xiàn)［4－12］。正則化理論的基本思想是將病態(tài)問題加上先驗條件，以恢復問題的穩(wěn)定性。通過正則化方法，可以將病態(tài)問題轉化為病態(tài)性相對較輕的問題或良性問題［4－5］。由于實際工作的巨大需求，正則化方法的理論與應用研究在測量界發(fā)展迅猛［10－18］。文獻［6］對測量領域中的不適定問題提出了正則化解的統(tǒng)一表達式。文獻［7］在GPS快速定位等實際工作中充分考慮參數(shù)的物理信息和先驗信息，并將其成功應用于正則化矩陣的構造和正則化參數(shù)的選取中。文獻［8］在航空重力測量數(shù)據(jù)向下延拓中比較了多種延拓方法，最后發(fā)現(xiàn)正則化方法效果較佳，并提出了以均方誤差極小值為目標確定正則化參數(shù)的方法。文獻［9］在GPS中長基線實時系統(tǒng)誤差估計中采用正則化方法，其模型中既有DD碼觀測又有載波相位（差）觀測，待求參數(shù)為位置坐標和對流延遲，正則化矩陣設計思想是僅對相對準確的位置參數(shù)實施控制，體現(xiàn)了對先驗信息已知的參數(shù)進行約束的選權擬合法思想。文獻［10—11］挖掘利用復共線性的關鍵信息構造正則化矩陣，提出了基于復共線性診斷的正則化方法。文獻［12］提出將廣義嶺估計用于求解航空重力向下延拓病態(tài)問題，通過確定多個正則化參數(shù)獲得更小的均方誤差。文獻［13］利用廣義正則化方法降低病態(tài)性對總體最小二乘求解的影響，以提高解算結果的穩(wěn)定性。文獻［14］利用GRACE數(shù)據(jù)解算重力場問題時采用了正則化方法，在構造正則化矩陣時利用了CSRRL04GRACE解誤差的信息，在選取正則化參數(shù)時采用了L曲線投影變量法，極大地減少了計算量。文獻［15］給出了4類可用于重力場解算的正則化矩陣（零次、一次、二次和Kaula），以及用于確定正則化參數(shù)的L曲線法和GCV方法的數(shù)學模型。文獻［16］為避免正則化參數(shù)對航空重力向下延拓過程中可靠成分的修正影響，提出Tikhonov雙參數(shù)正則化法。文獻［17］提出了基于均方誤差意義下的正則化參數(shù)的選取方法，具有較深入的理論意義。文獻［18］利用TSVD正則化方法有效地解決了水壩監(jiān)測模型中的復共線性問題。文獻［19］提出了修正偏差的正則化方法，并分析了其優(yōu)于普通正則化解的條件。

為了克服病態(tài)性對共軛梯度法的危害影響，文獻［5］根據(jù)正則化思想，針對法矩陣的最小特征值改進算法，使解的精度獲得提高。然而，實際工程技術中法矩陣往往具有多個小特征值對模型求解構成危害，僅考慮克服最小特征值對解算結果的影響是不夠的。因此本文在該方法的基礎上深入研究，從理論上提出了克服多個小特征值的解法，使解算方程的病態(tài)性降低到適當?shù)乃?，給出了具體的實施方法及步驟，最后通過算例驗證了新解法的有效性。

2 數(shù)學模型及共軛梯度法解法

大地測量與地球物理中大量的問題可歸結為對如下觀測方程的求解［1］

式中，L為n×1觀測向量為n×1觀測真值向量；Δ為n×1觀測誤差向量；A為n×t系數(shù)陣；rank（A）＝t，X＝［X1X2…Xt］T為t×1未知參數(shù)向量。式（1）具有如下統(tǒng)計特征

式中，σ20為未知單位權方差；單位陣I為觀測權陣（觀測權陣若不是單位陣，可將其單位化［7］）。由于式（1）為超定方程組，因此沒有常規(guī)解，只有最小二乘（least square，LS）解。在式（1）的兩邊同乘以AT，得到法方程

式（3）可簡化為

式中，N＝ATA；Y＝ATL。解之得未知參數(shù)X的LS估計

當線性方程組為大規(guī)模超定方程組時，共軛梯度法是一種常用的有效解法［1－3］。其步驟簡述如下。

步驟1：選取初值X0。

步驟2：計算殘差r0＝Y－NX0，令p0＝r0，α0＝，得到首次迭代結果X1＝X0＋α0p0。計算r1＝Y－NX1以備用。

步驟3：對于k＝1，2，…，n，重復以下迭代過程：計算pk＝rk＋βk－1pk－1，其中算得迭代結果Xk＋1＝Xk＋αkpk，計算rk＋1＝Y－NXk＋1。若rk＋12＜ε（式中，·2表示歐氏二范數(shù)，以下相同），停止迭代，轉入步驟4；否則繼續(xù)迭代。

步驟4：將最后的Xk＋1作為近似解，＝Xk＋1。

共軛梯度法在求解大規(guī)模超定的線性方程組時效果顯著，但是當線性方程組具有病態(tài)性時則出現(xiàn)失效、收斂速度太慢等問題。目前解算病態(tài)問題應用較廣的方法是Tikhonov正則化方法，它是針對不適定問題提出來的，對于線性模型式（1），其估計準則為

使Φ達到最小的參數(shù)X即為所要求的解。式（6）的第一項為殘差平方和，第二項Ω（X）稱為穩(wěn)定泛函，通常取Ω（X）＝X2H＝XTHX，式中，H為正則化矩陣。由此得到正則化解

正則化思想的核心是通過修正法矩陣，使解算方程的條件數(shù)大大降低，病態(tài)性大大減弱，因而其普適性強、應用廣泛。如果將正則化思想與共軛梯度法相結合，則改善的共軛梯度法就能夠用以求解病態(tài)性問題，下面將對此展開研究。

3 基于條件數(shù)控制的正則化迭代解法

3.1 干擾源向量與修正法矩陣

設法矩陣ATA的特征值為0＜λ1＜λ2＜…＜λt，相應的單位正交化特征向量為t×1向量e1、e2、…、et。法矩陣的病態(tài)性不僅表現(xiàn)在其最小特征值λ1相對于最大特征值λt很小，而且在多數(shù)情況下，法矩陣具有s個相對于最大特征值λt很小的特征值，即λ1＜λ2＜…＜λs＜＜λt。本文利用正則化思想修正或改造這些小特征值，獲得可以替換ATA的修正法矩陣，進而對法方程作改進，使修正法方程與原法方程相比，病態(tài)性減弱到便于穩(wěn)定解算的程度。為了對ATA的病態(tài)性進行干擾，構造以下n×1向量

式中，βi＞0，稱x0i為干擾源向量。利用干擾源向量x0i，設t×1向量

設

由e1，e2，…，et的正則性可得

式中，Bs是通過對法矩陣ATA改造獲得的，稱其為修正法矩陣。由式（11）可知，修正法矩陣Bs的t個特征值為λi＋βi2λi2（1≤i≤s）和λi（s＋1≤i≤t），特征向量仍為e1、e2、…、et。

由此可見，與ATA相比，Bs的條件數(shù)減少，由式（8）、式（9）可知，干擾源向量x0i的構造，或參數(shù)βi的選取是決定Bs的病態(tài)性的關鍵所在。干擾源向量的個數(shù)s的確定是另一個關鍵所在。

3.2 參數(shù)s及β1、β2、…、βs 的選取

干擾源向量的個數(shù)s以及β1、β2、…、βs的取法關系到法矩陣病態(tài)性被修正的程度，若s過小則修正不足，s過大則可能由于計算過程復雜引入過大的誤差，損害新舊方程的同解特性（下文詳述）。β1、β2、…、βs的選取則要考慮新方程的條件數(shù)被切實降低。為此對參數(shù)選取的總體策略是在滿足病態(tài)性降低到一定水平的條件下，干擾源向量的個數(shù)s盡量少。s以及β1、β2、…、βs的具體取法考慮如下。

（1）求出ATA的特征值為0＜λ1＜λ2＜…＜λt，計算cond（ATA），為修正法矩陣條件數(shù)的確定提供參考。

（2）確定修正法矩陣條件數(shù)K。根據(jù)實際問題，將法矩陣的條件數(shù)降低幾個數(shù)量級，即可達到明顯的效果，如本文針對航空測量數(shù)據(jù)向下延拓問題，選500較合適。

（3）計算條件指標［20］按其是否大于K分為兩類：第一類的條件指標過大，相應的特征值λ1、λ2、…、λs過小，應作修改，由此確定s。

（4）對于過小的特征值λ1、λ2、…、λs，利用其對應的特征向量按照式（8）構造相應的干擾源向量，其中參數(shù)β1、β2、…、βs的取法為

由式（12）可知λi＋βi2λi2＝λs＋1（1≤i≤s），又因修正法矩陣Bs的t個特征值為λi＋βi2λi2（1≤i≤s）和λi（s＋1≤i≤t），可知

3.3 修正法方程與基于條件數(shù)控制的正則化迭代解法

由式（8）、式（9）及式（1）可知

將式（3）與式（14）（1≤i≤s，共s個）相加，得到與式（3）同解的方程

簡記為

式中，Bs＝ATA＋y01yT01＋y02yT02＋…＋y0syT0s；Ls＝ATL＋y01xT01L＋y02xT02L＋…＋y0sxT0sL。式（16）是為改善法方程式（3）的病態(tài)性而對其進行的一種修正，稱其為修正法方程。至此，病態(tài)方程式（3）轉化為相對良態(tài)的方程式（16），用共軛梯度法求解即可生效。特別值得注意的是，式（3）和式（16）理論上是同解的，這是一般的正則化方法所不具有的良好性質。該結論是本文將文獻［5］的結論推演后得到的結論，文獻［5］給出了本算法當s取1時的一個特例，對本算法的形成提供了啟發(fā)。

在進行裝配過程中的信息追溯時，基于所建立的航空發(fā)動機裝配過程中的兩層數(shù)據(jù)世系，給出如下時間區(qū)間隊列和數(shù)據(jù)快照的定義

綜上所述，本文提出的算法主要分為3步：

步驟1：利用對法矩陣譜分解獲得的特征值與特征向量確定干擾源向量的個數(shù)、構造干擾源向量，進而構造修正法矩陣。

步驟2：利用干擾源向量及其特性，將法方程改造為修正法方程。

步驟3：用共軛梯度法解算修正法方程。稱該方法為基于條件數(shù)控制的正則化迭代解法（regularization based on controlling condition number，RBC）。

4 算例及分析

算例1：線性方程組L＝ˉL＋Δ＝AX，式中

未知參數(shù)的真值為X＝ ［1 1 1 1 1 1 1 1 1 1 ］T。為了比較幾種算法的解算效果，本文分別用LS估計、CG估計、Tikhonov估計和RBC估計對該問題進行解算，并計算其均方差作為精度，得到的結果見表1。

表1 4種算法的解算結果及精度Tab.1 The results and accuracies for the four solutions

從解算結果可以看出，LS估計嚴重失真，CG方法和兩種正則化方法均有效改善了LS估計求解的質量，而按新提出的RBC估計解算則達到了較高質量。

改正一個和多個特征值對估計結果的影響相差較大。從一個到多個一開始逐漸變好，但改正的個數(shù)達到一定程度后，繼續(xù)修改則效果不明顯（或效果開始變差）。通過試驗，按本文的修改方法效果較好，該問題將另文作進一步深入研究。

算例2：航空重力測量數(shù)據(jù)向下延拓模擬試驗。在2°×2°的某丘陵區(qū)域由EGM2008模型生成576個網(wǎng)格點的重力異常數(shù)據(jù)，將其作為真值，如圖1所示。將地面數(shù)據(jù)向上延拓至空中3400m［11］比地面區(qū)域稍大的區(qū)域內，并作正態(tài)隨機擾動，形成空中重力測量數(shù)據(jù)的模擬觀測值。本文按4種方式進行延拓解算：LS、CG、Tikhonov、RBC，并求出各種解值的均方差，作為精度指標考察求解質量［11］。解算結果見圖1—圖5及表2。

圖1 重力異常真值Fig.1 True value of gravity anomaly

圖2 重力異常的LS解Fig.2 LS solution of gravity anomaly

圖3 重力異常的CG解Fig.3 CG solution of gravity anomaly

圖4 重力異常Tikhonov正則化解Fig.4 Tikhonov solution of gravity anomaly

圖5 重力異常的RBC正則化解Fig.5 RBC solution of gravity anomaly

表2 4種解的統(tǒng)計分析Tab.2 Statistical analysis of the four solutions ms－2

可以看出，LS估計解算效果嚴重失真，CG、Tikhonov和RBC均可改善求解質量，而RBC改善求解質量效果最佳。

5 結 論

共軛梯度法是解算大型線性方程組的常用方法，具有解算精度高、數(shù)據(jù)存儲方便等優(yōu)點，但是大地測量與地球物理應用中經(jīng)常出現(xiàn)的問題具有強病態(tài)性，甚至是超強病態(tài)性，使得共軛梯度法在這種情形下失效。本文對此進行了深入的研究，提出了基于條件數(shù)控制的正則化迭代解法。通過對法矩陣作譜分析，構造干擾源向量，改造病態(tài)的法方程為修正法方程，使解算方程的病態(tài)性減弱到適當?shù)乃?，在此基礎上用共軛梯度法解算，得到了良好的效果。新算法較好地將正則化思想結合到共軛梯度法之中，在理論上具有一定的價值。通過算例表明了該算法的有效性及與相關算法比較的優(yōu)越性，特別在航空重力向下延拓的解算試驗中獲得成功，對于大地測量病態(tài)性問題的求解具有一定的應用價值。

［1］ ZHENG Wei，XU Houze，ZHONG Min，et al.Accurate and Rapid Determination of the GRACE Earth’s Gravitational Field Using the Improved Pre－conditioned Conjugate－gradient Approach and Three－dimensional Interpolation Method［J］.Progress in Geophysics，2011，26（3）：805－812.（鄭偉，許厚澤，鐘敏，等.基于改進的預處理共軛梯度法和三維插值法精確和快速解算GRACE地球重力場［J］.地球物理學進展，2011，26（3）：805－812.）

［2］ LIU C S，ATLURI S N.An Iterative Method Using an Optimal Descent Vector for Solving an Ill－conditioned System Bx＝b Better and Faster than the Conjugate Gradient Method［J］.Computer Modeling in Engineering and Sciences，2011，80（3）：275－298.

［3］ XU Shufang，GAO Li，ZHANG Pingwen.Numerical Linear Algebra［M］.Beijing：Beijing University Press，2000：139－160.（徐樹方，高立，張平文.數(shù)值線性代數(shù)［M］.北京：北京大學出版社，2000：139－160.）

［4］ XU Peiliang，YOICHI F，LIU Yumei.Multiple Parameter Regularization：Numerical Solutions and Applications to the Determination of Geopotential from Precise Satellite Orbits［J］.Journal of Geodesy，2006，80（1）：17－27.

［5］ LIU C S，HONG H K，ATLURI S N.Novel Algorithms Based on the Conjugate Gradient Method for Inverting Illconditioned Matrices，and a New Regularization Method to Solve Ill－posed Linear Systems［J］.Computer Modeling in Engineering and Sciences，2010，60（3）：279－308.

［6］ OU Jikun.Uniform Expression of Solutions of Ill－posed Problems in Surveying Adjustment and the Fitting Method by Selection of the Parameter Weights［J］.Acta Geodaetica et Cartographica Sinica，2004，33（4）：283－288.（歐吉坤.測量平差中不適定問題解的統(tǒng)一表達與選權擬合法［J］.測繪學報，2004，33（4）：283－288.）

［7］ WANG Zhenjie.Regularization of Ill－posed Problems in Surveying［M］.Beijing：Science Press，2006：66－85.（王振杰.測量中不適定問題的正則化解法［M］.北京：科學出版社，2006：66－85.）

［8］ WANG Xingtao，XIA Zheren，SHI Pan，et al.A Comparison of Different Downward Continuation Methods for Airborne Gravity Data［J］.Chinese Journal of Geophysics，2004，47（6）：1017－1022.（王興濤，夏哲仁，石磐，等.航空重力數(shù)據(jù)向下延拓方法比較［J］，地球物理學報，2004，47（6）：1017－1022.）

［9］ LUO Xiaowen，OU Jikun，YUAN Yunbin.An Improved Regularization Method for Estimating Near Real－time Systematic Errors Suitable for Medium－long GPS Baseline Solution［J］.Earth Planets Space，2008，60：793－800.

［10］ GU Yongwei，GUI Qingming，BIAN Shaofeng，et al.Comparison between Tikhonov Regularization and Truncated SVD in Geophysics［J］.Geomatics and Information Science of Wuhan University，2005，30（3）：238－241.（顧勇為，歸慶明，邊少鋒，等.地球物理反問題中兩種正則化方法的比較［J］.武漢大學學報：信息科學版，2005，30（3）：238－241.）

［11］ GU Yongwei，GUI Qingming.Study of Regularization Based on Signal－to－Noise Index in Airborne Gravity Downward to the Earth Surface［J］.Acta Geodaetica et Cartographica Sinica，2010，39（5）：458－464.（顧勇為，歸慶明.航空重力測量數(shù)據(jù)向下延拓基于信噪比的正則化方法的研究［J］.測繪學報，2010，39（5）：458－464.）

［12］ JIANG Tao，LI Jianchen，WANG Zhengtao，et al.Solution of Ill－posed Problem in Downward Continuation of Airborne Gravity［J］.Acta Geodaetica et Cartographica Sinica，2011，40（6）：684－689.（蔣濤，李建成，王正濤，等.航空重力向下延拓病態(tài)問題求解［J］.測繪學報，2011，40（6）：684－689.）

［13］ GE Xuming，WU Jicang.Generalized Regularization to Ill－posed Total Least Squares Problem ［J］.Acta Geodaetica et Cartographica Sinica，2012，41（3）：372－377.（葛旭明，伍吉倉.病態(tài)總體最小二乘問題的廣義正則化［J］.測繪學報，2012，41（3）：372－377.）

［14］ SAVE H，BETTADPUR S，TAPLY B.Reducing Errors in the GRACE Gravity Solutions Using Regularization［J］.Journal of Geodesy，2012，86（4）：695－711.

［15］ XU Xinyu，LI Jiancheng，WANG Zhengtao，et al.The Simulation Research on the Tikhonov Regularization Applied in Gravity Field Determination of GOCE Satellite Mission［J］.Acta Geodaetica et Cartographica Sinica，2010，39（5）：465－470.（徐新禹，李建成，王正濤，等.Tikhonov正則化方法在GOCE重力場求解中的模擬研究［J］.測繪學報，2010，39（5）：465－470.）

［16］ DENG Kailiang，HUANG Motao，BAO Jingyang，et al.Tikhonov Two－parameter Regulation Algorithm in Downward Continuation of Airborne Gravity Data［J］.Acta Geodaetica et Cartographica Sinica，2011，40（6）：690－696.（鄧凱亮，黃謨濤，暴景陽，等.向下延拓航空重力數(shù)據(jù)的Tikhonov雙參數(shù)正則化方法［J］.測繪學報，2011，40（6）：690－696.）

［17］ SCHAFFRIN B.Minimum Mean Squared Error Adjustment and the Optimal Tykhonov－Phillips Regularization Parameter via Reproducing Best Invariant Quadratic Uniformly Unbiased Estimates［J］.Journal of Geodesy，2008，82（2）：113－121.

［18］ XU Chang，DENG Chengfa.Solving Multicollinearity in Dam Regression Model Using TSVD ［J］.Geo－spatisl Information Science，2012，14（3）：230－234.

［19］ SHEN Yunzhong，XU Peiliang，LI Bofeng.Bias－corrected Regularized Solution to Inverse Ill－posed Models［J］.Journal of Geodesy，2012，86（8）：567－608.

［20］ GUI Qingming，YAO Shaowen，GU Yongwei，et al.A New Method to Diagnose Multicollinearity Based on Condition Index and Variance Decomposition Proportion［J］.Acta Geodaetica et Cartographica Sinica，2006，35（3）：210－214.（歸慶明，姚紹文，顧勇為，等.診斷復共線性的條件指標—方差分解比法［J］.測繪學報，2006，35（3）：210－214.）