參數(shù)有界約束下平差模型的一種新算法

陳仲兒 左廷英 宋迎春

1 中南大學地球科學與信息物理學院，長沙市麓山南路932號，410083

在測量數(shù)據(jù)處理中，通常以G－M 模型為平差模型，同時顧及非線性模型線性化的誤差，建立的簡化模型是不準確的，參數(shù)估值也非其準確值，即參數(shù)存在不確定性［1－2］。一個較為有效的解決辦法是整體最小二乘平差法［3］，其優(yōu)點在于同時考慮量測信息和設(shè)計矩陣的誤差，但是只考慮上述誤差是遠遠不夠的。另一種處理方式是，利用系統(tǒng)已知先驗信息，構(gòu)造不等式約束，建立附不等式約束的平差模型，采用相關(guān)方法來處理，其理論和應用在測量數(shù)據(jù)處理領(lǐng)域已有一定研究［4－5］，但是不等式約束也只是不確定因素中極小的一部分。就目前研究情況來看，綜合考慮上述所有不確定性問題幾乎是不可能的。而對于一個真實系統(tǒng)，有時可以得到量測誤差的極限值、參數(shù)的上下界等［6］，其綜合反映為相應的數(shù)值部分，可以通過區(qū)間分析［7－8］來解決。但區(qū)間運算也存在一些缺點和不足，可能導致在解算時出現(xiàn)法方程系數(shù)陣秩虧、解范圍擴大的不良結(jié)果，而且近幾年研究工作尚處于理論階段［6－10］。本文基于參數(shù)有界性，提出參數(shù)有界約束下的平差模型，該模型可轉(zhuǎn)化為相應的附不等式約束的平差模型，但由于不等式是單邊的，該轉(zhuǎn)化過程會造成不等式個數(shù)增加增加計算量。為此，利用最小二乘原理，將問題轉(zhuǎn)化為附有箱型約束的二次規(guī)劃問題，提出一種求解參數(shù)最優(yōu)估值的新算法，克服了不等式模型單邊約束的缺陷。

1 參數(shù)有界約束下的最小二乘平差算法

參數(shù)有界約束下的平差模型為：

式中，B是一個n×t的設(shè)計矩陣；X＝［X1，…，Xt］Τ是t×1 的參數(shù)向量；Δ是隨機誤差向量，Δ～N（0，σ2I），I是單位矩陣；c＝［c1，…，ct］Τ、d＝［d1，…，dt］Τ分別為參數(shù)的上、下界。將式（1）轉(zhuǎn)化為：

其中，G＝［－II］Τ，WC＝［－cd］Τ。顯然，式（2）為一個附不等式約束的平差模型。

式（2）所示的平差模型若無參數(shù)的邊界限制，可利用最小二乘求得參數(shù)最優(yōu)解。而當參數(shù)具有有界約束時，最小二乘解就不一定落入該范圍內(nèi)。于是，需要尋求附加約束條件的最優(yōu)解。利用最小二乘原理，可將式（1）表述為以下問題：

從而轉(zhuǎn)化為附有箱型約束的二次規(guī)劃問題：

其中，A＝BΤPB，b＝BΤPL，A是一個對稱正定（半正定）矩陣，因而，式（4）中第1式是一個凸函數(shù)，式（3）等價為一個凸二次規(guī)劃問題。

由于式（1）可轉(zhuǎn)化為式（2），因此也可采用附不等式約束的平差算法來求解，但是由于不等式為單邊的，該轉(zhuǎn)化過程會造成不等式個數(shù)增加，不便于計算。為此給出一個新算法，克服現(xiàn)有算法的不足。

式（3）中第1式是殘差平方和最小表達式，按照最小二乘準則，可得到相應的法方程：BΤ－BΤPL＝0，其中，為參數(shù)無約束下的最小二乘解。若對X作中心化處理，可將式（3）中第2式轉(zhuǎn)化為：

其中，Z＝X－z，z＝（c＋d）／2為各參數(shù)邊界范圍內(nèi)中心點，r＝（d－c）／2為中心點到邊界的長度。由于Z具有對稱的取值范圍，‖Z‖≤‖r‖，也即ZΤZ≤rΤr，對Z中每一分量Zi，均有ri。在求解過程中，若Z不滿足上述取值范圍，將使解偏離于問題（3）的解。令ZΤZ盡可能小，以使解盡量限制在其邊界范圍內(nèi)。為便于計算，可限制其達到最小，即對于 可行解，滿足min，也就是范數(shù)最小條件。因而，問題（3）的可行解滿足如下條件：

其中，K＝［k1k2…kt］Τ為聯(lián)系數(shù)向量。

將Φ對X＊求一階導數(shù)，并令其為零，得：

將上式代入式（6a），從中解出向量K的表達式：

同樣，將式（9）代入式（6a），反解出X＊：

若N正定，即模型（1）為滿秩問題，式（10）等價于最小二乘解：

為獲得，還需檢驗X＊是否落入X的取值區(qū)間。求得的結(jié)果可能存在以下兩種情況：

1）若對于X＊中每一分量均有，則根據(jù)式（10）求得的解X＊即為問題（3）的最優(yōu)解；

2）若X＊中至少存在一個分量不滿足（3）中所給參數(shù)的取值條件，即，則記下這些分量的下標，組成新集合，構(gòu)造一個空集合S，并作如下處理：

定理1［10］若N半正定，則可行解是問題（3）的最優(yōu)解的充要條件是，對i＝1，2，…，t，均有：

具體算法如下：

步驟1 對于一個具體平差問題，建立形如式（1）的模型，利用式（10）或式（11）求得無有界約束下的參數(shù)值X＊。令其為參數(shù)初值X（1），同時生成一空的解的集合X，X＝｛X（2），X（3），…，X（k），…｝，其中X（k）為第k次計算后的參數(shù)解。若X（1）滿足情況1），則轉(zhuǎn)步驟4，反之，轉(zhuǎn)步驟2。

步驟2 若X（k）出現(xiàn)情況2），則計算出指標s，判斷的大小。由于問題（3）的全局最優(yōu)解極有可能在參數(shù)無約束解的附近，而最優(yōu)解往往是未知的，因此，若比較接近其下界，可令，反之，取ds。判斷設(shè)計矩陣B（k）的s列中是否存在不為零的元素，記下該行下標i，將其寫入一個空集合同時，對L（k）中滿足下標為i∈I的每一行，減去常數(shù)，并刪去B（k）的第s列元素和P（k）的第s行與第s列所有元素，，重新按式（10）或式（11）計算X（k＋1）中其他未知分量，轉(zhuǎn)步驟3。

步驟3 若步驟2求得的X（k＋1）滿足定理1，則轉(zhuǎn)步驟4；反之，排除，對剩下的作情況2）處理，其中，下標“＼s”表示除去下標為s的所有分量，轉(zhuǎn)步驟2。

步驟4 計算終止，輸出當前的參數(shù)可行解。

定理2 設(shè)X（k）為式（10）的解，若X（k）滿足情況2），則稱相應的分量為自由變量。若令或1，2，…，t，則稱算法具有有限步收斂性。

2 新算法解PB的統(tǒng)計性質(zhì)

其中，B′＝B（t1＋1），P′＝P（t1＋1），D′＝D（t1＋1）分別表示第一部分觀測值的權(quán)陣和方差陣。因為取常數(shù)，所以。

在參數(shù)的方差－協(xié)方差陣中，對于i1，j1∈I1，，均有而對于i2，，均有且對于i1、i2，均有

3 實例應用

設(shè)參數(shù)有界約束下的平差模型為：

其中，B、L、c、d分別為平差模型的設(shè)計矩陣、觀測向量、參數(shù)的下界與上界，且

求參數(shù)X的最優(yōu)估值。

解：1）由式（10）計算最小范數(shù)解X＊。設(shè)P（1）＝I，

上述解為無約束下的最小二乘解。令其為參數(shù)初值X（1），可以看出小于其下界，即不滿足給出的邊界條件，于是轉(zhuǎn)步驟2。

2）取指標s為2，由于比較接近其下界，令，判斷B（1）中第2列中不為零的元素。對L（1）中第2行，減去常數(shù)，刪去B（1）中第2列元素和P（1）中第2行與第2列所有元素，重新按式（10）計算X（2）中其他未知分量

轉(zhuǎn)步驟3。

3）由計算出的X（2）＝［6／17 －2 －6／17－3／17］Τ，對于所有分量，均滿足模型所給參數(shù)的邊界條件，且對于Ni、（WB）i（i＝1，3，4），均有Ni X（2）－（WB）i＝0，對于N2、（WB）2，有N2X（2）－（WB）2＝2／17＞0，即對于所有分量，均滿足定理1，所以X的最優(yōu)估值即為X（2）＝［6／17－2－6／17－3／17］Τ。可見，在參數(shù)有界約束下，最優(yōu)解要明顯優(yōu)于最小二乘解，也在一定程度上減少了不確定性。

4）單位權(quán)方差估值計算：

的協(xié)因數(shù)陣：中誤差為：

如果將模型轉(zhuǎn)化為附不等式約束的平差模型：

其中，G＝［－II］Τ，WC＝［0 2 2 10 1 1 2 1］Τ。不等式個數(shù)由4個增加為8個。采用文獻［5］附不等式約束的平差模型迭代算法，同樣可求得，與本文算法的結(jié)果一致。通過比較，兩種算法都是收斂的，本文算法只需3步，而文獻［5］算法需要5步。

4 結(jié) 語

本文提出參數(shù)有界約束下的平差模型新算法，不僅克服了不等式單邊約束的缺陷，且新算法具有簡單、可行、收斂速度快的優(yōu)勢。計算分析證明，附加約束條件后，估值優(yōu)于最小二乘解，在一定程度上減少了部分數(shù)值的部分不確定性，使其結(jié)果更加真實。本文采用模擬測量數(shù)據(jù)來驗證算法的可行性，但在實際應用中，由于系統(tǒng)先驗信息往往不是確切已知的，但可大致確定其模糊范圍，故可先確定參數(shù)模糊上下界，再用本文算法來計算。得到初步解后，按照解所在的范圍逐漸縮小模糊度范圍，從而可更精確地得到參數(shù)的上下界。

［1］Chandrasekaran S，Golub G H，Gu M，et al.Parameter Estimation in the Presence of Bounded Data Uncertainties［J］.SIAM Journal on Matrix Analysis and Applications，1998，19（1）：235－252

［2］陶本藻.GIS質(zhì)量控制中不確定度理論［J］.測繪學院學報，2000，17（4）：235－238（Tao Benzao.Basic Theory of Uncertainty of Quality Control in GIS［J］.Journal of Institute of Surveying and Mapping，2000，17（4）：235－238）

［3］丁克良.整體最小二乘法及其在測量數(shù)據(jù)處理中的若干應用研究［D］.武漢：中國科學院測量與地球物理研究所，2006（Ding Keliang.Research on the Total Least Squares and Its Application in Surveying Data Processing［D］.Wuhan：Institute of Measurement and Geophysics，CAS，2006）

［4］彭軍還，張亞利，章紅平，等.不等式約束最小二乘問題的解及其統(tǒng)計性質(zhì)［J］.測繪學報，2007，36（1）：50－55（Peng Junhuan，Zhang Yali，Zhang Hongping，et al.The Solution of Inequality－constrained Least Squares Problem and Its Statistical Properties［J］.Acta Geodaetica et Cartographica Sinica，2007，36（1）：50－55）

［5］朱建軍，謝建.附不等式約束平差的一種簡單迭代算法［J］.測繪學報，2011，40（2）：209－212（Zhu Jianjun，Xie Jian.A Simple Iterative Algorithm for Inequality Constrained Adjustment［J］.Acta Geodaetica et Cartographica Sinica，2011，40（2）：209－212）

［6］劉世君，徐衛(wèi)亞，王紅春.不確定性巖石力學參數(shù)的區(qū)間反分析［J］.巖石力學與工程學報，2004，23（6）：885－888（Liu Shijun，Xu Weiya，Wang Hongchun.Interval Back Analysis on Uncertain Parameters in Rock Mechanics［J］.Chinese Journal of Rock Mechanics and Engineering，2004，23（6）：885－888）

［7］羅麟.區(qū)間分析法在滑坡災害風險評價中的應用研究［D］.北京：中國地質(zhì)大學，2008（Luo Lin.Interval Analysis Applied and Studied in Landslide Hazard and Risk Assessment［D］.Beijing：China University of Geoscience，2008）

［8］邱天.區(qū)間分析在邊坡工程中的應用［D］.南京：河海大學，2006（Qiu Tian.The Application of Interval Analysis in Slope［D］.Nanjing：Hohai University，2006）

［9］邱志平，顧元憲.有界不確定參數(shù)結(jié)構(gòu)位移范圍的區(qū)間攝動法［J］.應用力學學報，1999，16（1）：1－9（Qiu Zhiping，Gu Yuanxian.Perturbation Methods for Evaluating the Bounds on Displacements of Structures with Uncertain but Bounded Parameters［J］.Chinese Journal of Applied Mechanics，1999，16（1）：1－9）

［10］Liu Xiaojun，Jiao Yangchang，F(xiàn)ujishige S.An Algorithm for Strictly Convex Quadratic Programming with Box Constraints［J］.OR Transactions，1998，2（1）：8－22