姚 娟,張 镠2,王坤朋,唐齊鴻2,劉忠仁

(1.西南科技大學(xué) 信息工程學(xué)院，四川 綿陽 621010; 2.重慶大學(xué) 自動化學(xué)院，重慶 400044)

0 引言

大型構(gòu)筑體表面形變監(jiān)測是保障其安全運行的主要手段之一。如久經(jīng)碾壓的橋梁在結(jié)構(gòu)上會發(fā)生一定的形狀變化；隧道內(nèi)壁會因為周圍土體擾動、建筑負(fù)荷變化或施工振動而產(chǎn)生形變[1]；高層建筑在地震、強風(fēng)等自然災(zāi)害后發(fā)生傾斜。物體的形變?nèi)粑幢患皶r檢測和采取應(yīng)對措施，將對人員和設(shè)備的安全構(gòu)成巨大的潛在威脅，因此對構(gòu)筑體表面沉降、曲率以及位移等形變參數(shù)的監(jiān)測具有重要意義和實用價值。

目前，三維激光掃描儀被廣泛應(yīng)用于形變監(jiān)測領(lǐng)域，涉及建筑、隧道、基坑、鐵路、橋梁、礦山、地形勘查[2-9]等，擁有快速、高效、高密度、非接觸[6-8]等優(yōu)點。但在諸多形變監(jiān)測應(yīng)用中，主要根據(jù)具體的監(jiān)測對象確定相應(yīng)的監(jiān)測參數(shù)和制定監(jiān)測方案，還未形成完善的監(jiān)測體系[10]，主要通過測量標(biāo)靶坐標(biāo)的改變達(dá)到監(jiān)測目的，存在監(jiān)測參數(shù)單一、靈活性差、表面形變測量精度低等缺點。近年來，激光測距儀雖然逐步應(yīng)用在形變監(jiān)測領(lǐng)域[4-5]，相對三維激光掃描儀而言，采用激光測距結(jié)合轉(zhuǎn)向云臺[5]對物體表面形變進(jìn)行測量的方法，具有成本低、易于實現(xiàn)和靈活性強等優(yōu)點，但因其只對少數(shù)特征點進(jìn)行沉降變化監(jiān)測，并受限于搭載平臺，對物體表面形變測量誤差較大。

本文針對于目前三維激光在物表形變監(jiān)測應(yīng)用中存在的不足，提出一種結(jié)合激光采樣和高精度建模對物體表面形變進(jìn)行監(jiān)測的方法。其中，岳天祥等提出的高精度曲面建模(high accuracy surface modeling, HASM)方法，已經(jīng)廣泛應(yīng)用于數(shù)字高程模型構(gòu)建、土壤污染物空間分布模擬、氣候要素(氣溫、降水、輻射等)時空變化模擬[11]等領(lǐng)域，其理論體系相對成熟，在建模精度上優(yōu)于各種傳統(tǒng)的經(jīng)典插值法[12-16]。利用激光測距儀對物體表面進(jìn)行柵格化采樣，通過高精度曲面建模得到表面沉降和曲率變化等形變參數(shù)，為實現(xiàn)大型構(gòu)筑體安全評價和預(yù)警提供數(shù)據(jù)支撐。

1 HASM4及形變參數(shù)求解

1.1 采樣點坐標(biāo)計算

根據(jù)激光非接觸、高效、精準(zhǔn)的特點，對檢測物體表面進(jìn)行柵格化采樣。不論是一維激光測距還是三維激光掃描，在采樣儀器坐標(biāo)系下，對于采樣點對應(yīng)的激光束，根據(jù)其偏航角ψ和俯仰角θ，以及距采樣點的距離S，可求得其坐標(biāo)系下空間坐標(biāo)(X,Y,Z)，如圖1所示，其中X=S·cosψ·cosθ，Y=S·sinψ·cosθ，Z=S·sinθ。

圖1 坐標(biāo)計算示意圖

1.2 坐標(biāo)轉(zhuǎn)換

將直接采樣的測量點映射到柵格坐標(biāo)系中，記為[x,y,z]T=β·A·[X,Y,Z]T+ [Δx,Δy,Δz]T，其中為柵格對應(yīng)坐標(biāo)系下的采樣坐標(biāo)值，Δx,Δy,Δz分別為3個方向上坐標(biāo)的偏移量，β為解算尺寸，A∈R3×3為旋轉(zhuǎn)矩陣[1]。

1.3 HASM4

HASM是一種空間插值法，HASM4為其第4代更新。若采樣基于I×J的柵格，設(shè)定插值密度為k，則插值柵格可選擇(k(I-1)+1)×(k(J-1)+1) 的大小，即N×M。先進(jìn)行簡單插值，其結(jié)果將作為HASM求解的初始值。

設(shè)物表曲面表示為函數(shù)f(x,y)。柵格頂點坐標(biāo)為(xi,j,yi,j)，其對應(yīng)的函數(shù)值為fi,j。以往的研究表明曲面模型可由Gauss方程模擬[16-18]。

(1)

(2)

對離散化的Gauss方程用矩陣形式表示，即HASM4為:

(3)

s.t.L×Fn=D

(4)

公式(4)可通過引入拉格朗日乘子λ，轉(zhuǎn)換為無條件約束的最小二乘問題，

(5)

PFn=Rn-1

(6)

1.4 Gauss-Seidel法求解

HASM的求解包含兩個循環(huán)嵌套迭代：更新Rn-1為外迭代；采用Gauss-Seidel法求解公式(6)中的線性方程組更新Fn為內(nèi)迭代?；诿恳淮蜶n-1更新后的線性方程表示為RF=R，正定的系數(shù)矩陣P可被拆分為P=D-L-U，其中:

(7)

其中:k為迭代次數(shù)，方程的矩陣表達(dá)形式為(D-L)Fk=UFk-1+R，即Fk=(D-L)-1UFk-1+(D-L)-1R，設(shè)T=(D-L)-1U，C=(D-L)-1R，則:

Fk=TFk-1+C(k=1,2,…)

(8)

1.5 HASM求解算法

求解算法偽代碼如表1所示，其中F1取值柵格插值結(jié)果，作為求解算法的輸入。設(shè)置最小優(yōu)化權(quán)重參數(shù)λ和內(nèi)外迭代結(jié)束閾值a、b，并計算出固定參數(shù)P、T，然后進(jìn)入內(nèi)外嵌套迭代。內(nèi)迭代為求解PFn=Rn-1過程，外迭代則是根據(jù)當(dāng)前內(nèi)迭代的解Fn更新Rn過程，最后輸出設(shè)定精度的模擬結(jié)果Fn。

表1 HASM求解算法

1.6 形變參數(shù)沉降、曲率求解

基于高精度曲面建模插值結(jié)果，通過局部多項式空間插值的方法實現(xiàn)對曲面上任意點的模擬值和第一、二基本量進(jìn)行近似求解，最終得到曲面形變參數(shù)：沉降、曲率或曲率半徑。

如圖2所示，在求解HASM得到的模擬值fi,j,i=1,…,N,j=1,…,M周圍選取n×n大小的柵格作為初始樣本，進(jìn)而通過多項式空間插值求解該點模擬值和曲率參數(shù)。根據(jù)曲面上樣本區(qū)域的彎曲程度合理確定其插值多項式階數(shù)m，并對樣本柵格的各行各列求解其插值多項式。

圖2 插值樣本

設(shè)物表曲面上監(jiān)測點O(x0,y0,z0)，若通過最小二乘法解得第i行插值多項式及其導(dǎo)數(shù)為:

(9)

第j列插值多項式及其導(dǎo)數(shù)為:

(10)

對于每一行，y=y0處模擬值:

fi,y0=fi(y0) (i=1,2,…,n)

(11)

得y=y0這一列的插值多項式:

(12)

同理，對于每一列，x=x0處模擬值:

fx0,j=fj(x0)(j=1,2,…,n)

(13)

得x=x0這一行的插值多項式:

(14)

則(x0,y0)處模擬值:

z0=f(x0,y0)≈fx0(y0)≈fy0(x0)

(15)

監(jiān)測點沉降則是物表曲面形變前后該處模擬值之差。其一階偏導(dǎo)：

(16)

二階偏導(dǎo)：

(17)

(18)

可求得混合偏導(dǎo):

(19)

(20)

可求得混合偏導(dǎo):

(21)

(22)

2 試驗測試

2.1 數(shù)值仿真實驗

為驗證該方法的可行性，選取已知表達(dá)式的兩組無量綱、非突變的曲面做數(shù)值仿真對照實驗。

f1(x,y)=3+2sin(2πx) sin(2πy)+

exp(-15(x-1)2-15(y-1)2)+exp(-10x2-15(y-1)2)

其中，曲面f1可模擬地表面貌，特征相對復(fù)雜；曲面f2為球面，模擬一般規(guī)則建筑表面，特征相對簡單；模擬區(qū)域為x,y∈[0,1]，在該區(qū)域內(nèi)隨機產(chǎn)生100個測試點，如圖3所示。通過該方法求取測試點處的模擬值和不同方向上的曲率半徑。分別在不同采樣密度/間距、插值密度、曲率方向、坡度條件下，以控制變量法對比、分析兩種曲面的模擬誤差。

圖3 曲面f1、f2和測試點分布

2.1.1 采樣密度對模擬誤差的影響

設(shè)采樣間距為d0，插值間距為d1，插值密度為K=d0/d1，模擬值均方差、曲半徑相對平均差分別為:

如圖4所示，隨采樣密度的倍數(shù)增加，兩個曲面模擬值誤差和曲率半徑誤差均迅速減小，表明兩個的曲面的模擬精度快速提高；曲面f1模擬誤差變化較大，而f2相對較小，表明更為復(fù)雜的物體表面曲面對采樣密度的敏感度較大；當(dāng)采樣間距為0.0125時，兩個曲面的模擬精度接近，表明足夠高的采樣密度能夠減小曲面的復(fù)雜程度對模擬精度的影響；曲率誤差比模擬值誤差更明顯，表明曲率精度對模擬值精度非常敏感。

2.1.2 插值密度對模擬誤差的影響

保持d0=0.1，插值密度K∈{2,5,10}，得如曲線圖5變化趨勢。

圖5 模擬誤差與插值密度關(guān)系

隨著插值密度的倍增，曲面模擬值誤差和曲率半徑誤差緩慢減小且速度變緩，表明增大插值密度在一定程度上能夠提高模擬精度，但插值密度對模擬精度的影響比采樣密度更小。

2.1.3 坡度對模擬誤差的影響

當(dāng)K=2，d0=0.025,分析測試點位置坡度對模擬精度的影響，獲取兩個曲面高距均為0.5的等值線圖，如圖6所示。

圖6 模擬誤差與坡度分布關(guān)系

由圖可知，曲面f1的整體坡度遠(yuǎn)大于曲面f2，且紅色標(biāo)記均主要分布于坡度相對較大部位。結(jié)果表明，大坡度是降低模擬精度的重要因素，也是導(dǎo)致在稀疏采樣條件下曲面f1模擬精度低于曲面f2的根本因素。

2.2 形變實測

為驗證該方法的有效性，設(shè)計半徑分別為0.4米和0.2米圓柱曲面模型模擬物體表面形變前、后狀態(tài)，通過針對該方法設(shè)計的實驗設(shè)備先后對兩個圓柱面柵格化采樣，在MATLAB環(huán)境下對采樣數(shù)據(jù)應(yīng)用該方法進(jìn)行建模計算，對10處測試點給出形變前后的高度和曲率半徑，并進(jìn)行實驗數(shù)據(jù)分析。

以0.025米為間距對模型進(jìn)行柵格化采樣，采樣區(qū)域為x∈[0,0.275]，y∈[0,0.125]測試模型和10個測試點A,B,…,I,J分布如圖7所示。

半徑的柱面模型模擬形變前物表曲面，測試結(jié)果如表2所示。

表中h0為直接測量高度，h為建模解算高度，Δh=|h0-h|為高度誤差；r為解算曲率半徑，Δr=|r0-r|/r0×100%為相對曲率半徑誤差，數(shù)據(jù)表明模擬高度平均誤差達(dá)0.000 39 m，模擬曲率半徑相對平均誤差達(dá)5.9%。

圖7 圓柱曲面模型模擬物體

數(shù)據(jù)表明模擬高度平均誤差達(dá)0.29毫米，模擬曲率半徑相對平均誤差達(dá)7.1%。

形變前后曲面狀態(tài)及其測試點形變參數(shù)及其誤差如圖8、表4所示。

表2 形變前測試點測量參數(shù)及誤差

表3 形變后測試點測量參數(shù)及誤差

表4 沉降、曲率半徑變化及誤差

圖8 模型形變前后狀態(tài)和測試點分布

3 結(jié)論

本文提出了一種基于激光采樣，并結(jié)合高精度曲面建模求解曲面模擬值和曲率參數(shù)，進(jìn)而達(dá)到從沉降和曲率角度監(jiān)測物體表面形變的方法。數(shù)值仿真實驗測試了采樣密度、插值密度和曲面坡度對該方法測量精度的影響；并用實物模型模擬物體表面形變，驗證了該方法的可行性和有效性。實驗表明，雖然坡度變化影響模擬精度，但提高采樣密度和插值密度能夠降低其影響，精確解算出不同復(fù)雜程度的物體表面曲面的模擬值和曲率參數(shù)，進(jìn)而精確監(jiān)測物表形變；在模擬形變監(jiān)測實驗中，沉降精度可達(dá)98.41%，曲率半徑變化精度可達(dá)81.82%，明顯地反映了物體表面形變趨勢，到達(dá)監(jiān)測目的，該方法能夠有效地適用于大型構(gòu)筑體表面形變的在線監(jiān)測。