暑期數(shù)學建模培訓課件空洞探測

1、主講: 廖川榮(2000) 空洞探測問題空洞探測問題 山體、隧洞、壩體等的某些內(nèi)部結(jié)構山體、隧洞、壩體等的某些內(nèi)部結(jié)構 可用彈性波測量來確定。一個簡化問題可可用彈性波測量來確定。一個簡化問題可 描述為，一塊均勻介質(zhì)構成的矩形平板內(nèi)描述為，一塊均勻介質(zhì)構成的矩形平板內(nèi) 有一些充滿空氣的空洞，在平板的兩個鄰有一些充滿空氣的空洞，在平板的兩個鄰 邊分別等距地設置若干波源，在它們的對邊分別等距地設置若干波源，在它們的對 邊對等地安放同樣多的接收器，記錄彈性邊對等地安放同樣多的接收器，記錄彈性 波由每個波源到達對邊上每個接波由每個波源到達對邊上每個接 收器的時收器的時 間，間，根據(jù)彈性波在介質(zhì)中和在空氣

2、中不同根據(jù)彈性波在介質(zhì)中和在空氣中不同 的傳播速度，來確定板內(nèi)空洞的位置。的傳播速度，來確定板內(nèi)空洞的位置。 一一.問題重述問題重述 現(xiàn)考察如下的具體問題： 一塊240（m）240（m）的平板(如圖) 1.在 AB邊等距地設置7個波源 Pi (i =1, 7)， CD邊對等地安放7個接收器 Qj (j =1,7)， 記錄由Pi發(fā)出的彈性波到達Qj的時間tij (秒); 2.在AD邊等距地設置7個波源 Ri (i =1,7)， BC邊對等地安放7個接收器 Sj (j =1,7)， 記錄由Ri發(fā)出的彈性波到達Sj的時間 (秒)。 ijQjABCDPiRiSj240（m）240（m）ijijtQ1Q

3、2Q3Q4Q5Q6Q7P10.06110.08950.19960.20320.41810.49230.5646P20.09890.05920.44130.43180.47700.52420.3805P30.30520.41310.05980.41530.41560.35630.1919P40.32210.44530.40400.07380.17890.07400.2122P50.34900.45290.22630.19170.08390.17680.1810P60.38070.31770.23640.30640.22170.09390.1031P70.43110.33970.35660.195

4、40.07600.06880.1042 ijtS1S2S3S4S5S6S7R10.06450.06020.08130.35160.38670.43140.5721R20.07530.07000.28520.43410.34910.48000.4980R30.34560.32050.09740.40930.42400.45400.3112R40.36550.32890.42470.10070.32490.21340.1017R50.31650.24090.32140.32560.09040.18740.2130R60.27490.38910.58950.30160.20580.08410.070

5、6R70.44340.49190.39040.07860.07090.09140.0583 ij 已知彈性波在介質(zhì)和空氣中的傳播速度 分別為2880（米/秒）和320（米/秒），且彈 性波沿板邊緣的傳播速度與在介質(zhì)中的傳播 速度相同。 要求:1）確定該平板內(nèi)空洞的位置。 2）只根據(jù)由Pi發(fā)出的彈性波到達Qj的時 間tij (i,j=1,7)，能確定空洞的位置嗎？ 3.）討論在同樣能夠確定空洞位置的前提 下, 減少波源和接受器的方法。 二、符號說明二、符號說明1.).ijp 彈性波經(jīng)過介質(zhì)的長度 5.).ij 觀測誤差ijijijpq.2880320ijt注: ij2.).ijq 彈性波經(jīng)過空氣

6、的長度ij3.).ijt 或傳播時間觀測值4.).ijm 無空洞時傳播時間理論值1i0iix第 個方格為空氣第 個方格為介質(zhì)6.)ijd 第i條波線被第j個網(wǎng)格所截線段的長度.kk 1 242;7.)l =Pk7iji j 016;ijijQ = , ,.,第條波線., = ,.12369.)X,.,Tx xx全體方格的特征向量.10.).ix 第i方格的特征量k43,44,.84;8.)742.i, j0,1,.6;ijkijlS Rkij-為條波線 ,0iP i0,iRi6,jSj0ABCD,6jQj11121 36212223684 18428436.ddddddDddd10102061

7、01216606165.;.;.;.TBq qqq qqq qq3 63 63 611228 48 4111, .,.kkkkkkldldld2010206101216606165.;.;.;.TBq qqq qqq qqijijqq或彈性波經(jīng)過空氣的長度 網(wǎng)格單元的編號:31A32A33A36A7A8A9A12A1A2A3A6A 三三.問題分析問題分析 1)已知均勻平板的尺寸為240(m)240(m),彈性波在介質(zhì)和空氣中的傳播速度分別為2880(m/s)和320(m/s). 如果介質(zhì)中無空洞,則彈性波從 到 或從 到 (i=1,2,.7)的傳播時間 應為t=240(m) 2880(m/s)

8、=0.833(s). （沿網(wǎng)格線的傳播時間為最短時間）ipiQiRiS 通過對表中給出的測量時間的分析 發(fā)現(xiàn), 彈性波從 到 或從 到 (i=0,1,.6)的傳播時間 ipiQiRiS 0.0583 , 0.5895, iiiit 由于0.0583與理論計算值0.0833(s)不符. 可得最大偏差為: =0.0833(s)-0.0583(s)=0.025(s).t 顯然,實際測量值不應小于彈性波在 純介質(zhì)中傳播時的理論計算值. 故 可視為彈性波傳播時間測量值的最大誤差.t 以 乘以彈性波在空氣中傳播速度,得 320 = 8(m). tt 可見,由于測量誤差的限制,只能探測到沿傳播方向上洞徑大于

9、8m的空洞. 2)所給問題為一個反問題,不妨先看它的正問題: 已知空洞的位置及彈性波在介質(zhì)和空氣中的傳播速度, 求由波源 發(fā)出的彈性波到達接收器 的時間 .ipijtjQ 正問題容易求解: 算出波線經(jīng)歷的每個小方格的長度, 根據(jù)該方格是介質(zhì)還是空洞,分別除以波在介質(zhì)和空氣中的傳播速度,對經(jīng)歷的所有方格求和,即得 .ijPQijt 正問題等價于: 算出波線 經(jīng)歷的每個小方格的長度, 根據(jù)該方格是介質(zhì)還是空洞,分別乘以0和1, 對經(jīng)歷的所有方格求和,即可得到波線 經(jīng)歷的空洞的總長度.ijPQijPQ 于是,反問題反問題可以這樣求解:先由波源 發(fā)出的彈性波到達接收器 的時間 ,及彈性波在介質(zhì)和空氣中

10、的傳播速度, 算出波線 經(jīng)歷的空洞的總長度. ipijPQjQijt 再由 經(jīng)歷的每個方格的長度,即可解出每個方格是介質(zhì)還是空洞.(即是0還是1)ijPQ四四.模型假設模型假設1) 觀測數(shù)據(jù)有測量誤差. 觀測數(shù)據(jù)除測 量誤差外是可靠的.2) 波在傳播過程中沿直線單向傳播,且不考慮波的反射、折射以及干涉等現(xiàn)象. 3) 空氣密度和介質(zhì)密度都均勻. 4) 彈性波在傳播過程中沒有能量損失.其 波速僅與介質(zhì)有關,且在同一均勻介質(zhì) 中波速不變. 彈性波沿平板邊緣的傳播 速度與在介質(zhì)中的傳播速度相同. 5) 假設平板可劃分為網(wǎng)格, 空洞定位于每 個網(wǎng)格單元內(nèi), 空洞大小大致相同.6) 若某路徑上彈性波傳播時

11、間的實際測 量值與該路徑上無空洞時彈性波傳播 時間的理論計算值之差在測量誤差范圍 內(nèi),則認為該路徑上無空洞分布.五五. 模型的建立與求解模型的建立與求解1.圖形排除法模型 建立直角坐標系如圖,x,y軸分別等距地放置7個波源,將平板劃分為36個方格;方格的邊長為40(m). ,0iP i0,iRi6,jSj0ABCD,6jQj 由于 = /2880+ /320+ , ijtijpijqijij 所以,當 時,取 = 0, 即認為在波線 或 經(jīng)過的途中不存在空洞;否則認為存在空洞.ijqijPQiiRS 于是,通過計算得到不經(jīng)過空洞的波線:ijijtm-傳播時間觀測值.ijt-無空洞時傳播時間理論

12、值.ijm01104556646500112233445566010210365663646500112233445566,.P QPQP QP QP QP QP QPQP QP QP QP QP QR SR SR SR SR SR SR SR SR SR SR SR SR SR SR S 顯然,空洞出現(xiàn)在沒有這些波線穿過的方格. 由此得出如下判斷準則: 圖中凡有1條以上連線通過的方塊區(qū)域均無空洞存在. 8A9A14A15A17A21A22A26A據(jù)此可知: 僅有8個方格區(qū)域可能存在 空洞. 2.線性回歸模型 為第i條波線被第j個網(wǎng)格所截線段的長度.ijd 設 為第 k=7i+j 條波線及長度

13、, (k=1,42; i, j=0,1, . 6; i j) 為第k=7i+j+42 條波線, (k=43, ,84; i, j=0,1,. 6; i j)kijlPQkijlS Rijpijq 記 為彈性波 經(jīng)歷介質(zhì)的總長 度, 為彈性波 經(jīng)歷空氣的總長度, 則 和可由下式二元線性方程組解出:ijpijqijPQijPQ2880320ijijijpqt22406ijijpqij,0,1,.6i j ijPQ1.計算彈性波經(jīng)歷空氣的總長度ijq2.計算第i條波線被第j個網(wǎng)格所截線段的長度.（i=1,2,84;j=1,2,36)ijdij波 線 P Q 的 方 程 為 :600yxiji6yi6

14、x+ i-j,0,1,.6i j ijRS波線的方程為:06jiyix 66ij xyi,0,1,.6i j 計算波線 , 與網(wǎng)格線x=n , y=m 的交點,從而算出波線所經(jīng)歷的每個方格的長度 . ijPQ,0,1,6ij ijRS6yixk6x+ i-j或 6yiyk6x+ i-j 即解下列方程組: ijd,0,1,6i j 3 63 63 611228 48 4111, .,.kkkkkkldldld1010206101216606165.;.;.;.TBqqqq qqqqq11121 3621222 3684 184 284 36.ddddddDddd X= 全體方格的特征向量1236

15、,Tx xx2010206101216606165.;.;.;.TBqqqqqqqqq由于 DX=B11121 3610121222 3620284 184 284 363665.dddxqdddxqdddxq1010206101216606165,TBq qqq qqq qq126781 27 98 08 4.;.; .;.l lll lllll2010206101216606165;TBqqqqqqqqq 注:以上線性方程組稱為超定方程組.應 用最小二乘法求解此超定方程組,即 求 min . 此為一個 非線性0-1規(guī)劃模型. TDXBDXB 當Rank(D)=36 時可求出該方程最小 二乘

16、意義下的解. X=1TTD DD B 3)由于 Rank(D)=36,所以可用數(shù)學軟 件(如:Lingo,Matlab中非負最小二乘 法 X=lsqnonneg(D,B). )直接得到方 程組 DX=B的解如下:X1=-0.0365X2=0.0373X3=-0.0730X4=-0.0617X5=0.0383X6=-0.0016X7=0.0056X8=0.9847X9=1.0674X10=0.0008X11=-0.0013X12=-0.0128X13=0.0580X14=0.9512X15=0.9588X16=0.0707X17=0.9149X18=0.0313X19=0.1042X20=-0.

17、0099X21=1.0209X22=0.9968X23=0.0303X24=-0.0586X25=-0.0110X26=1.0184X27=0.0058X28=-0.0046X29=-0.0625X30=-0.0040X31=-0.0762X32=0.0184X33=-0.0046X34=-0.0577X35=0.0416X36=0.03030,取適的 令0111iiixXx 1ix 則 的小方格為介質(zhì), 的小方 格為空洞,其余的 (如果有的話) 無法確 定.0ix 若取 =0.11,可得下表: X10X20X3 0X14 0 X50 X60 X70 X81 X91X10 0 X110 X12

18、0 X130 X141X15 1 X160 X171X18 0 X190 X200 X211 X221 X230 X240 X250 X261 X270 X280 X290X30 0 X310X32 0 X330 X340X35 0 X360 注: 相應于X的排列,上下與平板的空洞 位置正相反. 對應于平板上的空洞位置如下圖. 3.對第2問的討論 只根據(jù)由 發(fā)出的彈性波到達 的 時間 ,可得長度距陣 . 若方程組 有解,則能確定空洞位置, 但 iPjQijt142 36D1142 36DXB14 23 6DRank( )=29,所以方程組: 無解,故不能確定空洞的位置. 114236DXB 在

19、同樣能夠確定空洞位置的前提 下,減少波源和接收器的方法. 設在BC邊仍安放7個接收器 ,而AD 邊設置波源 的數(shù)量一個一個地增加(由 下而上),發(fā)現(xiàn)設置2個波源時就有 Rank(D)=36.iSjR 雖 DX=B 有解: X= , 但由該解1TTD DD B1TTA AA B 無法選取 得到 .x 直到設置5個波源時, 才由X= 得到較滿意的解:1TTA AA BX1=-0.0324X2=0.0485X3=-0.0889X4=-0.0959X5=0.0231X6=0.0661X7=0.0181 X8=0.9766X9=1.0591X10=0.0345X11=0.0201X12=-0.0688X1