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

左廷英，陳仲兒，宋迎春

（中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙410083）

在測量平差中，通常以G－M模型為平差模型，同時(shí)顧及非線性模型線性化的誤差［1－2］，求得的簡化模型往往是不準(zhǔn)確的，因而參數(shù)估值也非其準(zhǔn)確值，即參數(shù)存在不確定性［3－4］。尤其是在一些特殊情況下（如高精度GPS變形監(jiān)測網(wǎng)），由于模型給的不準(zhǔn)、觀測誤差、參數(shù)的不確定性等不確定性因素的影響，從而使得參數(shù)的最小二乘估值可能是與實(shí)際相違背的，即最小二乘的平滑、去噪性掩蓋了形變的真實(shí)信息，其原因在于沒有考慮系統(tǒng)的物理特性或先驗(yàn)信息，運(yùn)用純數(shù)學(xué)方法來處理，導(dǎo)致估值“失真”。

在參數(shù)估計(jì)問題上，整體最小二乘法是對以上不確定性的一個(gè)改進(jìn)。其將系統(tǒng)中不確定因素?cái)U(kuò)展到了模型問題上，綜合考慮設(shè)計(jì)矩陣與觀測噪聲誤差，估計(jì)得到的參數(shù)值在某些場合下［5－7］優(yōu)于最小二乘解。

整體最小二乘參數(shù)估計(jì)方法存在不足，其最大限制因素是必須了解觀測信息的統(tǒng)計(jì)特性和隸屬函數(shù)，而由于不確定性因素的出現(xiàn)，往往較難獲得上述統(tǒng)計(jì)信息。最小二乘法、整體最小二乘法都是通過人為的手段對其確定，這樣顯然是不合適的。相反，較之統(tǒng)計(jì)特性，獲取參數(shù)的上下界、誤差的極限值等［8］相對較容易些。國內(nèi)外學(xué)者對此類問題也做探討和研究 并成功摸索出一條求解上述不確定問題的道路［9－11］。就目前研究情況而言，較為成熟的解決辦法是區(qū)間分析，其已在巖土力學(xué)、結(jié)構(gòu)力學(xué)等地質(zhì)、材料方面取得了一些成績［9－11］。但由于區(qū)間分析的缺點(diǎn)和不足，多數(shù)學(xué)者更多關(guān)注的是其理論研究［12－13］。在測繪領(lǐng)域方面，陶本藻、王新洲和楊元喜等都對測量不確定度理論進(jìn)行研究，拓展了測量平差數(shù)據(jù)處理的理論和方法。本文基于參數(shù)有界性，構(gòu)建了參數(shù)有界約束下的平差模型，利用最小二乘原理，將問題轉(zhuǎn)化為附有箱型約束的二次規(guī)劃問題，運(yùn)用最優(yōu)化理論知識(shí)［14－15］來處理，進(jìn)而得到參數(shù)的一個(gè)可行解，并將其應(yīng)用到沉降監(jiān)測網(wǎng)實(shí)例中，得出較為合理的結(jié)果。

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

1．1 建立平差模型

參數(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．2 轉(zhuǎn)化為附有箱型約束的二次規(guī)劃問題

式（1）的解不只一組，故需附加最小二乘約束來求得參數(shù)的最優(yōu)解。

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

當(dāng)網(wǎng)中各觀測量為等權(quán)觀測量時(shí)，上式轉(zhuǎn)化為

其中，BΤPB是一個(gè)對稱正定（半正定）矩陣，式（3）中可簡化為

上式是一個(gè)凸函數(shù)。其中，N＝BΤPB，W＝BΤPL。故式（3）可等價(jià)為一個(gè)凸二次規(guī)劃問題：

其中，F(xiàn)（X）＝f（X）－LΤPL。

1．3 由分解算法計(jì)算參數(shù)估值

求解上述凸二次規(guī)劃問題的算法有很多，如內(nèi)點(diǎn)算法、積極集法、對偶方法等［16－18］，需要具有深厚的最優(yōu)化理論及數(shù)學(xué)知識(shí)，且沒有一個(gè)通用、簡單的算法來實(shí)現(xiàn)。本文采用文獻(xiàn)［19］中提到的一種分解算法來處理，以下稱“分解算法”。具體內(nèi)容如下：

N∈Rt×t是對稱矩陣，且N滿足下述條件：

N＝N1＋N2；N1－N2是對稱正定矩陣，其中N1，N2∈Rt×t。

則稱（N1，N2）為N的一正則分裂。

設(shè)＾X為問題（1）的最優(yōu)解，顯然有X＝＾X＋（X－＾X），因此

故而轉(zhuǎn)化為等價(jià)問題：

其中，Ω是式（3）的解集。

箱型約束二次規(guī)劃問題的分解算法：

2 算例與分析

圖1所示為南方某城市工業(yè)區(qū)部分沉降監(jiān)測網(wǎng)略圖，由于工業(yè)需要，大量開采地下水造成較為嚴(yán)重的地面下沉，為避免安全事故及不必要的人員和財(cái)產(chǎn)損失，在近幾年已響應(yīng)相關(guān)部門要求，禁止了地下水的過度開采 同時(shí)采用二等水準(zhǔn)測量方法獲得了前后兩期高差觀測值，其中BM1點(diǎn)為已知點(diǎn)（國家二等水準(zhǔn)點(diǎn)），高程為27．329 m，觀測周期為1年，觀測精度為0．86 mm。表1為兩期高差觀測成果。通過前后兩期水準(zhǔn)測量獲得的高差觀測值來分析地面沉降量。在具體處理過程中，以網(wǎng)中BM1點(diǎn)為固定點(diǎn)，其他點(diǎn)為待定點(diǎn)，由此建立參數(shù)有界約束下的平差模型：

相應(yīng)的附有箱型約束的凸二次規(guī)劃問題為

其中，N對稱正定，F(xiàn)（X）為嚴(yán)格凸函數(shù)。

圖1 沉降監(jiān)測網(wǎng)

表1 兩期高差觀測成果

解決上述問題關(guān)鍵在于如何確定參數(shù)（高程）的上下界，即上述模型中的c，d值。經(jīng)過對Ⅰ期水準(zhǔn)網(wǎng)做經(jīng)典自由網(wǎng)平差處理 可以得到參數(shù)的最或然值，在做Ⅱ期平差處理時(shí)，以Ⅰ期平差得到的參數(shù)值作為本次平差參數(shù)初值。根據(jù)已有資料顯示，在過去幾年，該地區(qū)年沉降量最大值不超過10 mm，由變形機(jī)理可知，由于禁止令的實(shí)施，其沉降量不應(yīng)超過往年最大值，故可近似認(rèn)為前后兩期沉降量不超過一個(gè)10 mm，如此就確定了參數(shù)的上下界。即

為了比較傳統(tǒng)最小二乘平差法與參數(shù)有界約束下的最小二乘平差算法的優(yōu)劣性，在具體處理過程中，本文還對前后兩期水準(zhǔn)網(wǎng)做了經(jīng)典自由網(wǎng)平差處理。同時(shí)為便于描述，分別將經(jīng)典自由網(wǎng)平差法和參數(shù)有界約束下的最小二乘平差算法記為LS、BLS。經(jīng)過MATLAB編程實(shí)現(xiàn)，兩種方法求得的參數(shù)估值及沉降量如表2、表3所示。

表2 兩種方法的參數(shù)估計(jì)結(jié)果 mm

表3 兩種方法得到的沉降量 mm

由上述結(jié)果可見，當(dāng)限定前后兩期網(wǎng)中待定點(diǎn)沉降量不超過10 mm時(shí)，由經(jīng)典自由網(wǎng)平差得到的參數(shù)估值已明顯不能滿足該要求，如P4點(diǎn)的沉降量已超出約束界限。而在同一個(gè)水準(zhǔn)網(wǎng)中，各待定點(diǎn)沉降量應(yīng)該相差不大，故參數(shù)有界約束下的最小二乘平差算法得到的結(jié)果較符合真實(shí)情況。

3 結(jié)束語

經(jīng)典最小二乘平差法基于傳統(tǒng)高斯—馬爾科夫模型 具有諸多優(yōu)良特性 可以得到參數(shù)的最優(yōu)估值，但由于不確定性問題的出現(xiàn)，依靠其單一平滑去噪的特點(diǎn)已然不能滿足高精度變形監(jiān)測數(shù)據(jù)處理要求，因此亟待構(gòu)建一個(gè)與真實(shí)系統(tǒng)更為接近的模型并提出相應(yīng)算法，較為有效的解決辦法是利用已知的先驗(yàn)信息或通過其他相關(guān)途徑，對模型附加約束條件，如參數(shù)有界，建立參數(shù)有界約束下的平差模型，并運(yùn)用相關(guān)知識(shí)來求解。實(shí)踐證明，經(jīng)過最優(yōu)化處理，得到的模型全局最優(yōu)點(diǎn)為參數(shù)的最佳估值，且滿足在參數(shù)有界下殘差平方和最小的原則，得到的估計(jì)結(jié)果在某些情況下可能比最小二乘更能反映系統(tǒng)真實(shí)狀態(tài)。對于小型變形監(jiān)測網(wǎng)數(shù)據(jù)處理較為適用。

本文以某一沉降監(jiān)測網(wǎng)為算例，分別采用經(jīng)典最小二乘法和基于參數(shù)有界約束下的最小二乘平差算法來處理，新算法的估值更符合實(shí)際。但在實(shí)際工作中，往往無法獲知系統(tǒng)的先驗(yàn)信息，因此如何更為有效確定參數(shù)初值及上下界還有待深入研究。事實(shí)上，多數(shù)情況下的變形監(jiān)測網(wǎng)為GPS變形監(jiān)測網(wǎng)，此時(shí)參數(shù)為三維且方程維數(shù)較大，對于該類問題，需進(jìn)一步延伸。同時(shí)，分解算法雖簡單實(shí)用，但在具體執(zhí)行時(shí)，效率不高，需迭代多次才接近收斂，故迫切需要提出一種新的最優(yōu)化算法來解決參數(shù)有界問題。

［1］ 王新洲．非線性模型能否線性化的實(shí)用判據(jù)［J］．武漢測繪科技大學(xué)學(xué)報(bào)，1999，24（2）：145－148．

［2］ 張朝玉，陶本藻．平差系統(tǒng)的模型誤差及其識(shí)別方法研究［J］．武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2005，30（10）：897－899．

［3］ CHANDRASEKARAN S，GOLUB C H，GU M，et al．Parameter Esti mation in the Presence of Bounded Data Uncertainties［J］．SIA M Jour nal on Matrix Analysis and Application，1998，19（1）：235－252．

［4］ CHANDRASEKARAN S，GOLUB C H，GU M，et al．Parameter Estimation in the Presence of Bounded Mod－eling Errors J．IEEE Signal Pr ocessing Letters 1997 4（7）：195－197．

［5］ 丁克良．整體最小二乘法及其在測量數(shù)據(jù)處理中的若干應(yīng)用研究［D］．武漢：中國科學(xué)院測量與地球物理研究所，2006．

［6］ 楚彬，范東明．基于比例整體最小二乘的GPS高程擬合［J］．測繪工程，2014，23（4）：37－39．

［7］ 袁豹，岳東杰，張亮，等．基于總體最小二乘的相似變換模型及其在地圖掃描數(shù)字化中的應(yīng)用［J］．測繪工程，2013，22（4）：45－47．

［8］ 劉世君，徐衛(wèi)亞，王紅春．不確定性巖石力學(xué)參數(shù)的區(qū)間反分析［J］．巖石力學(xué)與工程學(xué)報(bào)，2004，23（6）：885－888．

［9］ 羅麟．區(qū)間分析法在滑坡災(zāi)害風(fēng)險(xiǎn)評價(jià)中的應(yīng)用研究［D］．北京：中國地質(zhì)大學(xué)，2008．

［10］邱天．區(qū)間分析在邊坡工程中的應(yīng)用［D］．南京：河海大學(xué)，2006．

［11］于生飛．基于區(qū)間不確定方法的邊坡穩(wěn)定性分析及非概率可靠度評價(jià)研究［D］．南京：南京大學(xué)，2012．

［12］邱志平，顧元憲．有界不確定參數(shù)結(jié)構(gòu)位移范圍的區(qū)間攝動(dòng)法［J］．應(yīng)用力學(xué)學(xué)報(bào)，1999，16（1）：1－9．

［13］蘇靜波．工程結(jié)構(gòu)不確定性區(qū)間分析方法及其應(yīng)用研究［D］．南京：河海大學(xué)，2006．

［14］袁亞湘，孫文瑜．最優(yōu)化理論與方法［M］．北京：科學(xué)出版社，1997：422－451．

［15］寇述舜．凸分析與凸二次規(guī)劃［M］．天津：天津大學(xué)出版社，1994：148－182．

［16］張藝．框式約束凸二次規(guī)劃問題的內(nèi)點(diǎn)算法［J］．高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào)，2002（2）：163－168．

［17］COLEMAN T F，LAURIE A H．A Direct Active Set Algorith m f or Lar ge Sparse Quadratic Programs with Simple Bound［J］．Mathematical Programming，1989，45：373－406．

［18］馬圣容．楊正豪．一類凸二次規(guī)劃的對偶方法［J］．南京師大學(xué)報(bào)：自然科學(xué)版，2003，26（1）：39－44．

［19］盧戰(zhàn)杰，魏紫鑾．邊界約束二次規(guī)劃問題的分解方法［J］．計(jì)算數(shù)學(xué)，1999，21（4）：475－482．