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

陳仲兒 左廷英 宋迎春

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙市麓山南路932號(hào)，410083

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

若N正定，即模型（1）為滿(mǎn)秩問(wèn)題，式（10）等價(jià)于最小二乘解：

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

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

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

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

具體算法如下：

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

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

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

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

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

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

其中，B′＝B（t1＋1），P′＝P（t1＋1），D′＝D（t1＋1）分別表示第一部分觀(guān)測(cè)值的權(quán)陣和方差陣。因?yàn)槿〕?shù)，所以。

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

3 實(shí)例應(yīng)用

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

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

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

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

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

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

轉(zhuǎn)步驟3。

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

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

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

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

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

4 結(jié) 語(yǔ)

本文提出參數(shù)有界約束下的平差模型新算法，不僅克服了不等式單邊約束的缺陷，且新算法具有簡(jiǎn)單、可行、收斂速度快的優(yōu)勢(shì)。計(jì)算分析證明，附加約束條件后，估值優(yōu)于最小二乘解，在一定程度上減少了部分?jǐn)?shù)值的部分不確定性，使其結(jié)果更加真實(shí)。本文采用模擬測(cè)量數(shù)據(jù)來(lái)驗(yàn)證算法的可行性，但在實(shí)際應(yīng)用中，由于系統(tǒng)先驗(yàn)信息往往不是確切已知的，但可大致確定其模糊范圍，故可先確定參數(shù)模糊上下界，再用本文算法來(lái)計(jì)算。得到初步解后，按照解所在的范圍逐漸縮小模糊度范圍，從而可更精確地得到參數(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］.測(cè)繪學(xué)院學(xué)報(bào)，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］丁克良.整體最小二乘法及其在測(cè)量數(shù)據(jù)處理中的若干應(yīng)用研究［D］.武漢：中國(guó)科學(xué)院測(cè)量與地球物理研究所，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］彭軍還，張亞利，章紅平，等.不等式約束最小二乘問(wèn)題的解及其統(tǒng)計(jì)性質(zhì)［J］.測(cè)繪學(xué)報(bào)，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］朱建軍，謝建.附不等式約束平差的一種簡(jiǎn)單迭代算法［J］.測(cè)繪學(xué)報(bào)，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)亞，王紅春.不確定性巖石力學(xué)參數(shù)的區(qū)間反分析［J］.巖石力學(xué)與工程學(xué)報(bào)，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ū)間分析法在滑坡災(zāi)害風(fēng)險(xiǎn)評(píng)價(jià)中的應(yīng)用研究［D］.北京：中國(guó)地質(zhì)大學(xué)，2008（Luo Lin.Interval Analysis Applied and Studied in Landslide Hazard and Risk Assessment［D］.Beijing：China University of Geoscience，2008）

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

［9］邱志平，顧元憲.有界不確定參數(shù)結(jié)構(gòu)位移范圍的區(qū)間攝動(dòng)法［J］.應(yīng)用力學(xué)學(xué)報(bào)，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