在模拟的3D空间中实现简单的三边测量算法

4
背景: 我正在为OpenComputers这个Minecraft mod添加的移动计算机实现导航系统。对于不熟悉这个mod的人来说,它基本上添加了各种可编程的Lua计算机,包括移动设备-即机器人、无人机和平板电脑。在尝试编写自主任务的机器人和无人机时经常出现的许多挑战之一是确保它们始终知道它们的坐标。

最简单的解决方案是使用Navigation升级,它确切地提供了计算机相对于制作它的地图中心的精确坐标。然而,它有两个主要缺点——它占用了Tier II升级插槽,这不是小事,而且限于地图范围。后者更或多或少可以接受,但仍使得该导航方法在某些使用情况下不可用。

另一个解决方案是让计算机记忆它们的坐标,然后跟踪它们的移动,但这也有许多潜在的注意事项-你必须通过自定义子例程控制所有移动或使用黑客拦截组件调用,你不能移动计算机而无需每次手动输入坐标,对于无人机存在一些精度误差,而这对平板电脑根本行不通。

第三种方法——我正在研究的方法——类似于现实生活中的GPS。它基于计算机可以升级为无线网络卡,以便能够在相当大的距离(400个块)内彼此发送消息,并且随着消息本身,它们接收到发件人和接收方之间的精确距离(浮点数,在块中)。如果我们将某些固定的计算机指定为“卫星”,它们不断广播它们的位置,我们可以让移动计算机能够使用4个或更多卫星的信息来三角定位其准确位置。

这种方法是可扩展的(您可以继续向网络中添加更多卫星以扩展其覆盖范围),不占用额外的导航升级插槽(因为许多移动计算机已经升级为无线网络卡),而且精确,这使得它比其他两种方法具有明显的优势。但是,它需要一些令人惊讶的复杂计算,这就是我卡住的地方。

问题: 我需要找到一个三角定位算法(最好附带代码示例),使任何移动计算机都能够在知道指定“卫星”的坐标和距离的情况下计算其位置(误差范围为 ~0.25个块)。我们假设所有计算机和卫星都配备了Tier II无线电卡(即它们可以在总范围为400个块内相互发送消息,并且以float32数字允许的精度知道发件人和自身之间的距离)。该解决方案将使用纯Lua编写,不访问任何第三方服务,因此像Mathematica这样的包是行不通的。目前,我打赌是一种适配方法,虽然我不知道如何实现它,以及它能否很好地适

在最基本的水平上,我们可以假设有4个卫星不断地并且正确地广播它们的位置,它们相互之间距离适当,并且不在同一个2D平面上。还有一些可选条件,理想情况下算法应该能够适应 - 请参见下面的部分。
额外加分项:
- 让算法足够小,以适应无人机的2KB内存(假设采用UTF8编码)。它应该比这个更小,以便主程序也能适应。越小越好。 - 制作一个算法,允许卫星非常接近彼此,并具有非整数坐标(以允许用一个不断移动的机器人或无人机代替多个固定卫星,或使移动计算机本身在从单个卫星获取测量时移动)。 - 制作一个算法,允许少于4颗卫星存在,假设已经可以确定位置 - 例如,如果所涉及的移动计算机是机器人,并且除了一种可能的位置在方块允许高度范围以下或以上(y<0或y>255)之外。如果有三颗卫星位于y=255的高度,这种设置是可能的。 - 制作一个算法,对某些卫星广播的略微错误的位置具有抗干扰性(设置中的小错误)。鉴于足够正确的测量存在,算法应该推导出正确的位置或直接报错。理想情况下,它还可以记录“关闭”卫星的位置。 - 制作一个算法,对同时存在于不同坐标系统中正确广播其位置的两个或多个卫星组具有抗干扰性(设置中的重大错误)。每个网络都有一个(可能独特的)识别符,允许区分由不同玩家(或者只是一个人)独立设置的不同网络。但是,如果他们没有费心正确设置标识符,则不同的信号可能会混合在一起,使移动计算机感到困惑。因此,抗干扰算法应该能够检测到这种情况,并且要么直接报错,要么区分不同的网络(然后可以进行微调以适应特定应用的目的 - 拒绝加载、选择最近的网络、选择最大的网络、提示用户或控制服务器等)。
我尝试过除了自己解决问题之外,还尝试在互联网上寻找合适的解决方案。然而,我找到的所有解决方案都不适合这个任务。
  • 通过谷歌搜索“三边测量算法”,我找到的大部分内容都与现实中的GPS系统有关,即仅使用2个坐标,强烈考虑误差并且通常不提供足够的精度。
  • 相反的是,一些纯数学的方法建议构建一系列方程来查找球体的交点。可惜的是,就我薄弱的数学背景而言,这种方法无法考虑浮点数的精度误差-圆圈并不完全相交,点并不完全在相同的位置,因此方程没有解。
  • 有些似乎解释了解决方案,但涉及了很多我无法理解的复杂数学,并没有包括确切的算法或至少一个代码示例。
  • 至少有一个使用像Mathematica这样的外部程序包,但在这种情况下不可用。

如果我留下了一些重要的点不清楚,请留言以便我改进问题。先谢谢了!


当有两组卫星(每组都有自己的坐标系)时,预期的输出是什么?程序是否应该仅输出您自己组的坐标?还是应该输出类似于这样的内容:“A组:(xa,ya,za);B组:(xb,yb,zb)”?不同坐标系的轴是否平行?A组中的x轴是否可以与B组中的x轴平行,但方向相反? - Egor Skriptunoff
感谢您的提问,@EgorSkriptunoff! 1)如果导航系统设置正确,则移动计算机应忽略所有它们不知道如何使用的网络。因此,是的,程序应仅在其自己的组中输出坐标。如果它没有正确设置 - 即检测到许多组并且不知道哪个是“自己的” - 那么理想情况下,它应该检测哪些卫星属于哪个组,并将其余部分留给实现(并默认抛出错误)。 2)是的,不同坐标系统的轴都以相同的方式对齐并且彼此平行。 - DeFazer
2个回答

3

这样的三边定位系统已经针对不同的模组(ComputerCraft)开发出来了。由于它可能与您的特定问题不兼容,因此您将需要修改和适应其逻辑,但算法本身应该是可行的。

以下是源代码

CHANNEL_GPS = 65534

local function trilaterate( A, B, C )
    local a2b = B.vPosition - A.vPosition
    local a2c = C.vPosition - A.vPosition
        
    if math.abs( a2b:normalize():dot( a2c:normalize() ) ) > 0.999 then
        return nil
    end
    
    local d = a2b:length()
    local ex = a2b:normalize( )
    local i = ex:dot( a2c )
    local ey = (a2c - (ex * i)):normalize()
    local j = ey:dot( a2c )
    local ez = ex:cross( ey )

    local r1 = A.nDistance
    local r2 = B.nDistance
    local r3 = C.nDistance
        
    local x = (r1*r1 - r2*r2 + d*d) / (2*d)
    local y = (r1*r1 - r3*r3 - x*x + (x-i)*(x-i) + j*j) / (2*j)
        
    local result = A.vPosition + (ex * x) + (ey * y)

    local zSquared = r1*r1 - x*x - y*y
    if zSquared > 0 then
        local z = math.sqrt( zSquared )
        local result1 = result + (ez * z)
        local result2 = result - (ez * z)
        
        local rounded1, rounded2 = result1:round( 0.01 ), result2:round( 0.01 )
        if rounded1.x ~= rounded2.x or rounded1.y ~= rounded2.y or rounded1.z ~= rounded2.z then
            return rounded1, rounded2
        else
            return rounded1
        end
    end
    return result:round( 0.01 )
end

local function narrow( p1, p2, fix )
    local dist1 = math.abs( (p1 - fix.vPosition):length() - fix.nDistance )
    local dist2 = math.abs( (p2 - fix.vPosition):length() - fix.nDistance )
    
    if math.abs(dist1 - dist2) < 0.01 then
        return p1, p2
    elseif dist1 < dist2 then
        return p1:round( 0.01 )
    else
        return p2:round( 0.01 )
    end
end

function locate( _nTimeout, _bDebug )
    -- Let command computers use their magic fourth-wall-breaking special abilities
    if commands then
        return commands.getBlockPosition()
    end

    -- Find a modem
    local sModemSide = nil
    for n,sSide in ipairs( rs.getSides() ) do
        if peripheral.getType( sSide ) == "modem" and peripheral.call( sSide, "isWireless" ) then   
            sModemSide = sSide
            break
        end
    end

    if sModemSide == nil then
        if _bDebug then
            print( "No wireless modem attached" )
        end
        return nil
    end
    
    if _bDebug then
        print( "Finding position..." )
    end
    
    -- Open a channel
    local modem = peripheral.wrap( sModemSide )
    local bCloseChannel = false
    if not modem.isOpen( os.getComputerID() ) then
        modem.open( os.getComputerID() )
        bCloseChannel = true
    end
    
    -- Send a ping to listening GPS hosts
    modem.transmit( CHANNEL_GPS, os.getComputerID(), "PING" )
        
    -- Wait for the responses
    local tFixes = {}
    local pos1, pos2 = nil, nil
    local timeout = os.startTimer( _nTimeout or 2 )
    while true do
        local e, p1, p2, p3, p4, p5 = os.pullEvent()
        if e == "modem_message" then
            -- We received a reply from a modem
            local sSide, sChannel, sReplyChannel, tMessage, nDistance = p1, p2, p3, p4, p5
            if sSide == sModemSide and sChannel == os.getComputerID() and sReplyChannel == CHANNEL_GPS and nDistance then
                -- Received the correct message from the correct modem: use it to determine position
                if type(tMessage) == "table" and #tMessage == 3 then
                    local tFix = { vPosition = vector.new( tMessage[1], tMessage[2], tMessage[3] ), nDistance = nDistance }
                    if _bDebug then
                        print( tFix.nDistance.." metres from "..tostring( tFix.vPosition ) )
                    end
                    if tFix.nDistance == 0 then
                        pos1, pos2 = tFix.vPosition, nil
                    else
                        table.insert( tFixes, tFix )
                        if #tFixes >= 3 then
                            if not pos1 then
                                pos1, pos2 = trilaterate( tFixes[1], tFixes[2], tFixes[#tFixes] )
                            else
                                pos1, pos2 = narrow( pos1, pos2, tFixes[#tFixes] )
                            end
                        end
                    end
                    if pos1 and not pos2 then
                        break
                    end
                end
            end
            
        elseif e == "timer" then
            -- We received a timeout
            local timer = p1
            if timer == timeout then
                break
            end
        
        end 
    end
    
    -- Close the channel, if we opened one
    if bCloseChannel then
        modem.close( os.getComputerID() )
    end
    
    -- Return the response
    if pos1 and pos2 then
        if _bDebug then
            print( "Ambiguous position" )
            print( "Could be "..pos1.x..","..pos1.y..","..pos1.z.." or "..pos2.x..","..pos2.y..","..pos2.z )
        end
        return nil
    elseif pos1 then
        if _bDebug then
            print( "Position is "..pos1.x..","..pos1.y..","..pos1.z )
        end
        return pos1.x, pos1.y, pos1.z
    else
        if _bDebug then
            print( "Could not determine position" )
        end
        return nil
    end
end

来自https://github.com/dan200/ComputerCraft/blob/master/src/main/resources/assets/computercraft/lua/rom/apis/gps.lua

如果您对源代码有任何具体问题,请提出。


1
谢谢!这段代码看起来是真实可行的,但我可能会选择@EgorSkriptunoff的答案,因为它更详细、先进,提供更多功能,并已经适用于OC。尽管如此,你提供的代码随附在官方的ComputerCraft发行版中,而且被广泛使用并良好运作,应该说明它相当不错和高效。虽然你可能应该提供源链接。 :) - DeFazer

1

函数trilateration需要卫星列表(它们的坐标和距离移动计算机的距离)和移动计算机的先前坐标。
仅收集您自己组的卫星,排除所有其他组的卫星。
您的一些卫星可能会发送错误数据,这没问题。

如果没有足够的卫星可访问,则该函数将返回nil,因为无法确定当前位置。
否则,该函数将返回移动计算机的当前坐标和被指责为不正确的卫星索引列表。
在存在歧义的情况下,新位置被选择为最接近移动计算机先前位置的位置。
输出坐标为整数,Y坐标限制在0..255范围内

以下条件应满足才能正确三边测量:

  • (number_of_correct_satellites)必须≥3
  • 如果存在至少一个不正确的卫星,则(number_of_correct_satellites)必须≥4
  • (number_of_correct_satellites)必须>(number_of_incorrect_satellites)

识别不正确的卫星是昂贵的CPU操作。
一旦将卫星识别为不正确,请将其存储在某个黑名单中,并将其从所有未来计算中排除。

do
   local floor, exp, max, min, abs, table_insert = math.floor, math.exp, math.max, math.min, math.abs, table.insert

   local function try_this_subset_of_sat(satellites, is_sat_incorrect, X, Y, Z)
      local last_max_err, max_err = math.huge
      for k = 1, math.huge do
         local oldX, oldY, oldZ = X, Y, Z
         local DX, DY, DZ = 0, 0, 0
         max_err = 0
         for j = 1, #satellites do
            if not is_sat_incorrect[j] then
               local sat = satellites[j]
               local dx, dy, dz = X - sat.x, Y - sat.y, Z - sat.z
               local d = (dx*dx + dy*dy + dz*dz)^0.5
               local err = sat.distance - d
               local e = exp(err+err)
               e = (e-1)/(e+1)/(d+1)
               DX = DX + dx*e
               DY = DY + dy*e
               DZ = DZ + dz*e
               max_err = max(max_err, abs(err))
            end
         end
         if k % 16 == 0 then
            if max_err >= last_max_err then
               break
            end
            last_max_err = max_err
         end
         local e = 1/(1+(DX*DX+DY*DY+DZ*DZ)^0.5/max_err)
         X = X + DX*e
         Y = max(0, min(255, Y + DY*e))
         Z = Z + DZ*e
         if abs(oldX - X) + abs(oldY - Y) + abs(oldZ - Z) <= 1e-4 then
            break
         end
      end
      return max_err, floor(X + 0.5), floor(Y + 0.5), floor(Z + 0.5)
   end

   local function init_set(is_sat_incorrect, len, ctr)
      for j = 1, len do
         is_sat_incorrect[j] = (j <= ctr)
      end
   end

   local function last_combination(is_sat_incorrect)
      local first = 1
      while not is_sat_incorrect[first] do
         first = first + 1
      end
      local last = first + 1
      while is_sat_incorrect[last] do
         last = last + 1
      end
      if is_sat_incorrect[last] == nil then
         return true
      end
      is_sat_incorrect[last] = true
      init_set(is_sat_incorrect, last - 1, last - first - 1)
   end

   function trilateration(list_of_satellites, previous_X, previous_Y, previous_Z)
      local N = #list_of_satellites
      if N >= 3 then
         local is_sat_incorrect = {}
         init_set(is_sat_incorrect, N, 0)
         local err, X, Y, Z = try_this_subset_of_sat(list_of_satellites, is_sat_incorrect, previous_X, previous_Y, previous_Z)
         local incorrect_sat_indices = {}
         if err < 0.1 then
            return X, Y, Z, incorrect_sat_indices
         end
         for incorrect_ctr = 1, min(floor((N - 1) / 2), N - 4) do
            init_set(is_sat_incorrect, N, incorrect_ctr)
            repeat
               err, X, Y, Z = try_this_subset_of_sat(list_of_satellites, is_sat_incorrect, previous_X, previous_Y, previous_Z)
               if err < 0.1 then
                  for j = 1, N do
                     if is_sat_incorrect[j] then
                        table_insert(incorrect_sat_indices, j)
                     end
                  end
                  return X, Y, Z, incorrect_sat_indices
               end
            until last_combination(is_sat_incorrect)
         end
      end
   end
end

使用示例:

-- assuming your mobile computer previous coordinates were 99 120 100
local previous_X, previous_Y, previous_Z = 99, 120, 100
-- assuming your mobile computer current coordinates are 111 112 113
local list_of_satellites = {
   {x=22, y=55, z=77, distance=((111-22)^2+(112-55)^2+(113-77)^2)^0.5},  -- correct satellite
   {x=35, y=99, z=42, distance=((111-35)^2+(112-99)^2+(113-42)^2)^0.5},  -- correct satellite
   {x=44, y=44, z=44, distance=((111-94)^2+(112-94)^2+(113-94)^2)^0.5},  -- incorrect satellite
   {x=10, y=88, z=70, distance=((111-10)^2+(112-88)^2+(113-70)^2)^0.5},  -- correct satellite
   {x=54, y=54, z=54, distance=((111-64)^2+(112-64)^2+(113-64)^2)^0.5},  -- incorrect satellite
   {x=91, y=33, z=15, distance=((111-91)^2+(112-33)^2+(113-15)^2)^0.5},  -- correct satellite
}

local X, Y, Z, list_of_incorrect_sat_indices = trilateration(list_of_satellites, previous_X, previous_Y, previous_Z)
if X then
   print(X, Y, Z)
   if #list_of_incorrect_sat_indices > 0 then
      print("Satellites at the following indices are incorrect: "..table.concat(list_of_incorrect_sat_indices, ","))
   end
else
   print"Not enough satellites"
end

输出:

111 112 113
Satellites at the following indices are incorrect: 3,5

谢谢!我目前无法在实际操作中检查此代码(Java在玩弄操作系统后会出问题),但它看起来是合法的。我会尽快尝试并报告结果。从我所看到的,它还可以检测和排除不正确的卫星,但是它是通过暴力破解来完成的,这意味着它的执行时间将呈指数增长,并且如果给定足够大的卫星列表,它将挂起整个计算机并显示“太长时间未进行让步”的错误。然而,对于特定应用程序,这个问题很容易解决,超出了本问题的范围。 :) - DeFazer
@DeFazer - 如果不正确的卫星一个接一个地出现,那么检测它们将非常快速。是的,因为经验值增长。 - Egor Skriptunoff
@DeFazer - 你在现实世界中测试了我的代码吗?我想听听一些批评 :-) - Egor Skriptunoff
@EgorSkriptunoff:抱歉,还没有。我目前遇到了一个问题,甚至无法启动Minecraft(似乎是这个问题http://www.minecraftforge.net/forum/topic/33291-fixed-minecraft-forge-freezes-on-quotreloading-model-manager-815quot/)。非常尴尬,因为我很快就收到了一个看起来不错的答案,但无法验证和接受它。我希望我能很快找到解决方案,然后能够尝试您的代码。^^ - DeFazer
请问这个算法使用的是什么方法?最好能提供一些包含算法内部数学知识的网址。能否解释一下这里正在发生的事情:local e = exp(err+err); e = (e-1)/(e+1)/(d+1); DX = DX + dxe; DY = DY + dye; DZ = DZ + dz*e,或者至少提供一个带有解释的网址。 - Prizoff
@Prizoff - 我使用了“逐步逼近法”:在每一步中,我们都会稍微移动我们的点(X,Y,Z)。每颗卫星都会告诉我们它是否希望我们离它更近或更远。只有当所有卫星都正确时,这个序列才会收敛。 (e-1)/(e+1)tanh(err),对于小误差是线性的,但对于大误差是有限制的,/(d+1) 使得靠近的卫星比远处的更重要。 - Egor Skriptunoff

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接