人機(jī)純方位無(wú)源定位:從三組方位角反推坐標(biāo)的MATLAB實(shí)現(xiàn))
簡(jiǎn)介本資源為2022年全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽本科組B題的完整論文PDF聚焦無(wú)人機(jī)純方位無(wú)源定位問(wèn)題適合備戰(zhàn)數(shù)模競(jìng)賽的本科生及指導(dǎo)教師參考。論文圍繞圓形與錐形兩種編隊(duì)系統(tǒng)建立被動(dòng)接收信號(hào)無(wú)人機(jī)的定位模型與位置調(diào)整模型運(yùn)用二維直角坐標(biāo)系、余弦定理與“三圓共點(diǎn)”原理推導(dǎo)發(fā)射信號(hào)無(wú)人機(jī)數(shù)量需求并給出基于0.01米步長(zhǎng)的遍歷調(diào)整方案全部求解過(guò)程均以MATLAB編程實(shí)現(xiàn)。壓縮包內(nèi)僅含1個(gè)PDF文件約2.63MB完整呈現(xiàn)摘要、問(wèn)題重述、模型假設(shè)、符號(hào)說(shuō)明及建模求解全過(guò)程便于讀者直接研讀建模思路與代碼邏輯。該論文曾獲河北省一等獎(jiǎng)目前已有8039人學(xué)習(xí)下載可為參賽者提供從模型構(gòu)建到編程求解的完整參考范例。1. 從三組方位角反推坐標(biāo)這套無(wú)人機(jī)純方位無(wú)源定位方案到底能跑通什么編隊(duì)飛行里最怕的不是風(fēng)是電磁信號(hào)被干擾后整個(gè)隊(duì)形散架。純方位無(wú)源定位的思路很直接無(wú)人機(jī)不主動(dòng)發(fā)波只靠接收其他已知位置無(wú)人機(jī)發(fā)來(lái)的方向信息反推自己在哪、該往哪挪。這套 2022 年全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽本科組 B 題的完整解答核心就是把這個(gè)過(guò)程拆成可計(jì)算的模型——用二維直角坐標(biāo)系加余弦定理把三個(gè)方位角變成坐標(biāo)再用步長(zhǎng)遍歷把偏差位置拉回理想位置。它適合兩類(lèi)人正在備戰(zhàn)數(shù)模競(jìng)賽、需要一套能直接跑 MATLAB 的完整代碼框架的選手以及做無(wú)人機(jī)編隊(duì)控制、想理解純方位定位最小可解條件的工程師。整套資料包含問(wèn)題重述、模型建立、MATLAB 程序附錄和結(jié)果分析不是純理論推導(dǎo)附錄里四個(gè)完整程序覆蓋了從定位到調(diào)整的全流程。2. 定位模型怎么建從三個(gè)方位角到二維坐標(biāo)的完整推導(dǎo)鏈2.1 為什么選余弦定理而不是最小二乘純方位無(wú)源定位的常見(jiàn)做法有兩種一種是基于最小二乘的偽線性估計(jì)把非線性角度方程線性化后迭代求解另一種是直接利用幾何關(guān)系構(gòu)造解析解。這套方案選了后者原因很實(shí)際——圓形編隊(duì)里發(fā)射信號(hào)的無(wú)人機(jī)位置是精確已知的接收信號(hào)的無(wú)人機(jī)與它們之間的角度關(guān)系可以用余弦定理直接串起來(lái)不需要迭代也不存在初值敏感的問(wèn)題。具體來(lái)說(shuō)把第 i 架無(wú)人機(jī)的準(zhǔn)確位置固定在 x 軸上圓形編隊(duì)半徑記為 r相鄰無(wú)人機(jī)間隔 40 度。三架發(fā)射信號(hào)無(wú)人機(jī)分別記為 FY0i、FY0j、FY0k接收信號(hào)無(wú)人機(jī)為 FY0x。接收到的三個(gè)方位角 α1、α2、α3 分別對(duì)應(yīng)接收機(jī)與三架發(fā)射機(jī)連線之間的夾角。在三角形中反復(fù)使用余弦定理就能把三條邊的長(zhǎng)度解出來(lái)再反解出接收機(jī)的坐標(biāo)。這里有個(gè)容易翻車(chē)的地方當(dāng)接收機(jī)位于發(fā)射機(jī)連線的延長(zhǎng)線上時(shí)角度信息退化方程會(huì)出現(xiàn)多解甚至無(wú)解。這也是后面問(wèn)題二里第二輪定位時(shí) 4、7、9、11、13、14 號(hào)無(wú)人機(jī)無(wú)法定位的根本原因。2.2 被動(dòng)接收信號(hào)無(wú)人機(jī)定位模型的 MATLAB 實(shí)現(xiàn)附錄一的程序分兩段第一段求解當(dāng)前位置坐標(biāo)第二段求解準(zhǔn)確位置坐標(biāo)。先看第一段的核心邏輯。% 被動(dòng)接收信號(hào)無(wú)人機(jī)定位模型的建立 % 求解當(dāng)前位置的坐標(biāo) alpha1 input(第一個(gè)角度); % 0k i alpha2 input(第二個(gè)角度); % 0k j alpha3 input(第三個(gè)角度); % j k i FY0i input(發(fā)射信號(hào)無(wú)人機(jī) 1); FY0j input(發(fā)射信號(hào)無(wú)人機(jī) 2); r input(圓形編隊(duì)半徑); % 計(jì)算兩架發(fā)射信號(hào)無(wú)人機(jī)之間的距離 a1 sqrt(2*r^2 - cosd(abs(FY0i - FY0j)*40)*2*r^2); syms k0 a2 a3 % 三個(gè)余弦定理方程三個(gè)未知量 eq1 (k0^2 a3^2 - r^2)/(2*k0*a3) - cosd(alpha1); eq2 (k0^2 a2^2 - r^2)/(2*k0*a2) - cosd(alpha2); eq3 (a2^2 a3^2 - a1^2)/(2*a2*a3) - cosd(alpha3); [X, Y, Z] solve(eq1, eq2, eq3, k0, a2, a3); K0 double(X); % 求得 k0 長(zhǎng)度 % 計(jì)算接收機(jī)相對(duì)于 x 軸的偏轉(zhuǎn)角 k0i acosd(1 - 2*cosd(alpha1)^2); k01 k0i - 40*abs(i - 1); x real(K0*cosd(k01)) % 橫坐標(biāo) y real(K0*sind(k01)) % 縱坐標(biāo)這段代碼的關(guān)鍵參數(shù)有三個(gè)alpha1、alpha2、alpha3是接收到的三個(gè)方位角單位是度FY0i和FY0j是兩架發(fā)射信號(hào)無(wú)人機(jī)的編號(hào)用來(lái)計(jì)算它們之間的基準(zhǔn)距離r是圓形編隊(duì)半徑題目里是 100 米。syms聲明符號(hào)變量后solve一次性解出三個(gè)未知邊長(zhǎng)再用double轉(zhuǎn)成數(shù)值。最后兩行把極坐標(biāo)轉(zhuǎn)成直角坐標(biāo)real是防止符號(hào)計(jì)算產(chǎn)生微小虛部。第二段程序反過(guò)來(lái)已知接收機(jī)編號(hào)和發(fā)射機(jī)編號(hào)直接算出三個(gè)方位角的準(zhǔn)確值% 求解無(wú)人機(jī)接收到的角度信息 FY00 0; FY0i input(發(fā)射信號(hào)無(wú)人機(jī) 1); FY0j input(發(fā)射信號(hào)無(wú)人機(jī) 2); FY0k input(接收信號(hào)無(wú)人機(jī)); r input(圓形編隊(duì)半徑); % 計(jì)算三條邊的長(zhǎng)度 a1 sqrt(2*r^2 - cosd((FY0j - FY0i)*40)*2*r^2); a3 sqrt(2*r^2 - cosd((FY0k - FY0i)*40)*2*r^2); a2 sqrt(2*r^2 - cosd((FY0k - FY0j)*40)*2*r^2); % 用余弦定理求三個(gè)角 alpha1 acosd((r^2 a3^2 - r^2)/(2*r*a3)); alpha2 acosd((r^2 a2^2 - r^2)/(2*r*a2)); alpha3 acosd((a3^2 a2^2 - a1^2)/(2*a3*a2));注意alpha1和alpha2的公式里出現(xiàn)了r^2 a3^2 - r^2看起來(lái)像冗余實(shí)際上是因?yàn)閳A心無(wú)人機(jī) FY00 到圓周無(wú)人機(jī)的距離恒為 r代入余弦定理后 r 平方項(xiàng)抵消了一部分。這是圓形編隊(duì)的特殊性質(zhì)換成錐形編隊(duì)就不能這么簡(jiǎn)化。提示運(yùn)行第一段程序前需要先確定接收機(jī)在圓周上的大致方位否則k01的計(jì)算會(huì)引入 40 度整數(shù)倍的偏移誤差。我一般會(huì)先用第二段程序算出準(zhǔn)確角度再拿第一段反解驗(yàn)證。3. 額外需要幾架發(fā)射機(jī)枚舉法確定最小可解條件3.1 為什么一架額外發(fā)射機(jī)就夠了問(wèn)題一第二小問(wèn)的核心結(jié)論是除 0 號(hào)和 1 號(hào)外只需要額外增加 1 架發(fā)射信號(hào)無(wú)人機(jī)就能實(shí)現(xiàn)有效定位。這個(gè)結(jié)論的推導(dǎo)依賴(lài)一個(gè)幾何事實(shí)——一個(gè)動(dòng)點(diǎn)到兩個(gè)定點(diǎn)的夾角為定值時(shí)動(dòng)點(diǎn)軌跡是上下對(duì)稱(chēng)的兩段圓弧組成的封閉曲線。當(dāng)夾角為鈍角時(shí)軌跡像一只“眼睛”當(dāng)夾角為銳角時(shí)軌跡像一個(gè)“葫蘆”。三架發(fā)射機(jī)產(chǎn)生三條這樣的軌跡三條軌跡相交于一點(diǎn)該點(diǎn)就是接收機(jī)的實(shí)際位置。如果只有兩架發(fā)射機(jī)接收機(jī)只能得到一個(gè)方位角軌跡是一條曲線而非一個(gè)點(diǎn)無(wú)法唯一確定位置。所以額外增加一架是理論下限。但這里有個(gè)隱藏條件三條軌跡必須確實(shí)相交于同一點(diǎn)。當(dāng)接收機(jī)位于某兩架發(fā)射機(jī)連線的延長(zhǎng)線上時(shí)對(duì)應(yīng)的軌跡退化成直線交點(diǎn)不唯一。這就是為什么在問(wèn)題二的錐形編隊(duì)里第二輪定位時(shí)會(huì)有 6 架無(wú)人機(jī)無(wú)法定位。3.2 枚舉法求解額外發(fā)射機(jī)編號(hào)的 MATLAB 實(shí)現(xiàn)附錄二的程序用枚舉法遍歷所有可能的發(fā)射機(jī)組合找出使方位角偏差平方和最小的那一組。% 求出所有的可能性 FY00 0; FY0i 1; % 兩個(gè)已知的定點(diǎn) r input(圓形編隊(duì)半徑); % 遍歷所有可能的發(fā)射機(jī)組合預(yù)計(jì)算角度表 for FY0j 2:9 for FY0k 2:9 if FY0k FY0j continue; end a1 sqrt(2*r^2 - cosd((FY0j - FY0i)*40)*2*r^2); a3 sqrt(2*r^2 - cosd((FY0k - FY0i)*40)*2*r^2); a2 sqrt(2*r^2 - cosd((FY0k - FY0j)*40)*2*r^2); alpha1(FY0j, FY0k) acosd((r^2 a3^2 - r^2)/(2*r*a3)); alpha2(FY0j, FY0k) acosd((a3^2 a2^2 - a1^2)/(2*a3*a2)); alpha3(FY0j, FY0k) acosd((r^2 a2^2 - r^2)/(2*r*a2)); end end % 反解根據(jù)接收到的角度信息找最匹配的發(fā)射機(jī) beta1 input(第一個(gè)角度); beta2 input(第二個(gè)角度); beta3 input(第三個(gè)角度); FY0k input(接收信號(hào)的無(wú)人機(jī)編號(hào)); min_val 1000; mini 0; for i 2:9 if i FY0k continue; end m (beta1 - alpha1(i, FY0k))^2 ... (beta2 - alpha2(i, FY0k))^2 ... (beta3 - alpha3(i, FY0k))^2; if m min_val min_val m; mini i; end end mini % 輸出需要額外增加的發(fā)射信號(hào)無(wú)人機(jī)編號(hào)這段代碼的邏輯分兩步先建一張角度查找表把所有可能的發(fā)射機(jī)組合對(duì)應(yīng)的三個(gè)方位角算出來(lái)存進(jìn)alpha1、alpha2、alpha3三個(gè)矩陣然后對(duì)于給定的接收角度beta1、beta2、beta3遍歷所有候選發(fā)射機(jī)找平方和最小的那個(gè)。參數(shù)說(shuō)明FY0i固定為 1因?yàn)轭}目規(guī)定 0 號(hào)和 1 號(hào)已經(jīng)發(fā)射信號(hào)FY0j和FY0k從 2 到 9 遍歷跳過(guò)相同編號(hào)m是三個(gè)角度偏差的平方和越小說(shuō)明匹配度越高。最終輸出的mini就是需要額外增加的發(fā)射機(jī)編號(hào)。注意枚舉法在 8 架無(wú)人機(jī)里找 1 架額外發(fā)射機(jī)計(jì)算量是 C(7,1)7 種組合非常小。但如果問(wèn)題擴(kuò)展到幾十架無(wú)人機(jī)這種暴力枚舉就不合適了需要換成基于幾何約束的剪枝搜索。4. 位置調(diào)整怎么收斂步長(zhǎng)遍歷與極坐標(biāo)更新4.1 向北還是向南調(diào)整方向的分類(lèi)邏輯問(wèn)題一第三小問(wèn)要求把有偏差的無(wú)人機(jī)調(diào)整到理想位置。方案設(shè)定步長(zhǎng)為 0.01 米向四周全方位遍歷當(dāng)與準(zhǔn)確位置相差小于 0.01 米時(shí)認(rèn)為到達(dá)理想位置。調(diào)整方向分成向北和向南兩類(lèi)本質(zhì)上是極坐標(biāo)下角度增量的正負(fù)。具體做法無(wú)人機(jī)當(dāng)前位置記為 B向北移動(dòng)一個(gè)微小步長(zhǎng)后到達(dá) C。已知 B 點(diǎn)的方位信息可以算出 BC 的長(zhǎng)度和角度變化再用余弦定理求出 C 點(diǎn)相對(duì)于圓心的極坐標(biāo)。重復(fù)這個(gè)過(guò)程直到偏差小于閾值。這里有個(gè)工程上很實(shí)用的簡(jiǎn)化調(diào)整時(shí)間不計(jì)意味著無(wú)人機(jī)知道自己需要調(diào)整位置的瞬間就已經(jīng)到位。這在實(shí)際系統(tǒng)里當(dāng)然不成立但在建模階段可以大幅簡(jiǎn)化動(dòng)力學(xué)部分把注意力集中在定位精度上。4.2 圓形編隊(duì)位置調(diào)整的 MATLAB 實(shí)現(xiàn)附錄三的程序用 while 循環(huán)實(shí)現(xiàn)迭代調(diào)整。% 從無(wú)人機(jī)接收的方位信息去定位無(wú)人機(jī)找到無(wú)人機(jī)的位置 FY00 0; FY01 1; njl 0.01; % 最小距離閾值 r 100; % 圓形編隊(duì)半徑 bzrow 100; % 準(zhǔn)確位置 row bztheta [0 40 80 120 160 200 240 280 320]; % 準(zhǔn)確位置 theta FY0irow(2) input(FY02 的初始位置(row):); FY0itheta(2) input(FY02 的初始位置(theta):); jl 100; while(jl njl) mjl 10000; for j 0:359 % 遍歷所有可能的角度找距離最小的方向 % 此處省略具體角度計(jì)算核心是極坐標(biāo)更新 % ... end % 更新位置縮小 jl jl mjl; end這段代碼的框架是外層 while 控制迭代直到距離小于閾值內(nèi)層 for 遍歷 0 到 359 度所有方向找使距離最小的調(diào)整角度。njl是收斂閾值設(shè)成 0.01 米r是編隊(duì)半徑bztheta是 9 架圓周無(wú)人機(jī)的理想角度分布間隔 40 度。實(shí)際運(yùn)行時(shí)每次迭代需要根據(jù)當(dāng)前極坐標(biāo)和方位角信息用余弦定理算出移動(dòng)一個(gè)步長(zhǎng)后的新極坐標(biāo)再與理想位置比較。如果偏差沒(méi)有縮小說(shuō)明方向選反了需要反向遍歷。提示步長(zhǎng) 0.01 米在 100 米半徑的編隊(duì)里相當(dāng)于萬(wàn)分之一精度迭代次數(shù)可能上百次。如果只做仿真驗(yàn)證可以先用 0.1 米步長(zhǎng)快速收斂最后再用 0.01 米精調(diào)。5. 錐形編隊(duì)怎么套用從圓形到錐形的模型遷移與邊界5.1 錐形編隊(duì)的幾何差異與定位策略問(wèn)題二的錐形編隊(duì)和圓形編隊(duì)本質(zhì)區(qū)別在于圓周上相鄰無(wú)人機(jī)等距但錐形編隊(duì)里相鄰三架無(wú)人機(jī)構(gòu)成等邊三角形直線上的間距相等。這意味著圓形編隊(duì)里“相鄰間隔 40 度”的簡(jiǎn)化不再成立需要重新建立坐標(biāo)系。方案假設(shè) 1 號(hào)和 5 號(hào)無(wú)人機(jī)位置無(wú)偏差先定位 2 號(hào)。具體做法將 1 號(hào)、2 號(hào)作為發(fā)射機(jī)5 號(hào)接收得到一個(gè)方向信息再將 2 號(hào)、5 號(hào)作為發(fā)射機(jī)1 號(hào)接收得到另一個(gè)方向信息。在三角形 1-2-5 中兩個(gè)角已知第三個(gè)角隨之確定2 號(hào)的位置就解出來(lái)了。但這里有個(gè)硬邊界當(dāng)某架無(wú)人機(jī)位于 1、2、5 號(hào)構(gòu)成的三角形的邊的延長(zhǎng)線上時(shí)只能傳送一個(gè)方向信息無(wú)法準(zhǔn)確定位。所以第二輪定位時(shí)4、7、9、11、13、14 號(hào)無(wú)人機(jī)依然無(wú)法確定位置只能先定位 3、6、8、10、12、15 號(hào)。第三輪選取 2、6、8 為發(fā)射機(jī)剩下的未定位無(wú)人機(jī)才都能解算。5.2 錐形編隊(duì)位置調(diào)整的 MATLAB 實(shí)現(xiàn)要點(diǎn)附錄四的程序結(jié)構(gòu)與附錄三類(lèi)似但坐標(biāo)系和理想位置表不同。錐形編隊(duì)的理想位置表如下無(wú)人機(jī)編號(hào)X 坐標(biāo)Y 坐標(biāo)1100√350250√3503004-50√3505-100√31006-100√307-100√3-1008-50√3-509-50√35010-50√3-15011-100√320012-100√310013-100√3014-100√3-10015-100√3-200調(diào)整邏輯與圓形編隊(duì)一致設(shè)定步長(zhǎng)向四周遍歷直到偏差小于閾值。區(qū)別在于錐形編隊(duì)的理想位置不是均勻分布在圓周上而是按等邊三角形網(wǎng)格排列所以極坐標(biāo)更新時(shí)需要先轉(zhuǎn)成直角坐標(biāo)再計(jì)算距離。% 錐形編隊(duì)位置調(diào)整核心邏輯 target_x -100*sqrt(3); % 示例15 號(hào)機(jī)目標(biāo) X target_y -200; % 示例15 號(hào)機(jī)目標(biāo) Y step 0.01; current_x input(當(dāng)前 X 坐標(biāo)); current_y input(當(dāng)前 Y 坐標(biāo)); while sqrt((current_x - target_x)^2 (current_y - target_y)^2) step % 計(jì)算當(dāng)前距離 dist sqrt((current_x - target_x)^2 (current_y - target_y)^2); % 沿梯度方向移動(dòng)一個(gè)步長(zhǎng) dx (target_x - current_x) / dist * step; dy (target_y - current_y) / dist * step; current_x current_x dx; current_y current_y dy; end這段代碼比圓形編隊(duì)的極坐標(biāo)遍歷更直接已知目標(biāo)坐標(biāo)直接沿梯度方向走步長(zhǎng)。參數(shù)step控制精度target_x和target_y從理想位置表里查。實(shí)際運(yùn)行時(shí)每架無(wú)人機(jī)需要先通過(guò)純方位定位確定自己的當(dāng)前位置再查表找到目標(biāo)位置然后執(zhí)行梯度下降。注意錐形編隊(duì)里 4、7、9、11、13、14 號(hào)無(wú)人機(jī)在第二輪無(wú)法定位必須等第三輪 2、6、8 號(hào)調(diào)整到位后才能解算。如果強(qiáng)行在第二輪對(duì)這些無(wú)人機(jī)執(zhí)行調(diào)整會(huì)因?yàn)槲恢霉烙?jì)錯(cuò)誤導(dǎo)致越調(diào)越偏。6. 從跑通到跑對(duì)四個(gè)容易翻車(chē)的參數(shù)與驗(yàn)證習(xí)慣第一個(gè)坑是角度單位。MATLAB 的acosd返回的是度但cos和sin默認(rèn)是弧度。附錄一的程序里混用了cosd和cos如果輸入的角度沒(méi)有統(tǒng)一算出來(lái)的坐標(biāo)會(huì)差一個(gè) π/180 的系數(shù)。我一般會(huì)在程序開(kāi)頭加一行clear; clc;然后把所有三角函數(shù)統(tǒng)一成帶 d 的版本。第二個(gè)坑是符號(hào)計(jì)算的虛部。solve解方程組時(shí)可能產(chǎn)生復(fù)數(shù)解double轉(zhuǎn)換后虛部雖然很小但不為零。real函數(shù)能去掉虛部但如果虛部量級(jí)超過(guò) 1e-6說(shuō)明方程本身有問(wèn)題需要檢查角度輸入是否合理。第三個(gè)坑是枚舉法的邊界。附錄二里FY0j和FY0k從 2 到 9 遍歷但alpha1矩陣的索引是(FY0j, FY0k)當(dāng)FY0j或FY0k等于FY0k接收機(jī)編號(hào)時(shí)對(duì)應(yīng)的角度值沒(méi)有物理意義。程序里用continue跳過(guò)了相同編號(hào)但沒(méi)有跳過(guò)接收機(jī)編號(hào)實(shí)際使用時(shí)需要手動(dòng)排除。第四個(gè)坑是步長(zhǎng)與收斂速度的權(quán)衡。0.01 米步長(zhǎng)在 100 米半徑下需要約 10000 次迭代才能從 100 米外收斂到目標(biāo)每次迭代都要重新計(jì)算方位角計(jì)算量不小。我一般會(huì)先用 1 米步長(zhǎng)粗調(diào)再用 0.01 米精調(diào)總迭代次數(shù)能降到十分之一。驗(yàn)證方法很簡(jiǎn)單跑完附錄一算出坐標(biāo)后代回附錄二反解發(fā)射機(jī)編號(hào)如果得到的編號(hào)與輸入一致說(shuō)明定位模型自洽。再用附錄三調(diào)整位置檢查最終坐標(biāo)與理想位置表的偏差是否小于 0.01 米。三個(gè)程序串起來(lái)跑通這套純方位無(wú)源定位方案就算真正落地了。從那以后我每次拿到新的編隊(duì)模型都強(qiáng)制先跑一遍“正解-反解-調(diào)整”的閉環(huán)驗(yàn)證確認(rèn)三個(gè)環(huán)節(jié)的數(shù)值能對(duì)上再往下做。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取