『Japan City Plans』古地図の投影系パラメータ推定にご協力お願いします

※ 本件、別途記事上げますが解決しました。ご協力ありがとうございました。

数日前から@Say_noさんと、米軍の70年前の日本都市地域図、Japan City Plans 1:12,500 U.S. Army Map Service, 1945-1946をタイル地図化するための謀略を行っております。

人手と時間さえかければ、GCP *1を与えまくってのむりくり投影変換=>タイル化、も可能なのですが、人も金もなかりけりなので、投影系の情報を読み解いて一気に変換したいと思っています。

とりあえず、地図上の凡例から、ONE THOUSAND YARD WORLD POLYCONIC GRID BAND IIIn, Zone A(東日本)又はB(西日本)というのが読み取れます。
それをヒントに検索しまくったところ、こちらの記事をみつけまして、

当時の米軍は米国内地図で多円錐図法を使っており、基準経線は8度刻みでZone AからGまであったらしい
基準緯線は北緯40.5度だったらしい
その座標系の単位はヤードで、基準経線/緯線での座標値はx:1000000、y:2000000だったらしい
測地系はNAD27（準拠楕円体はclrk66）だったらしい

といった情報がありました。

だったら日本の地図作成時も同じ基準流用してるだろ、と思って、全地図見直してみると、

X座標の1000000は、ZoneAでは東経143度、ZoneBでは東経135度付近にあるっぽい。両者の度間隔は8度
Y座標の2000000は、北緯40.5度付近にあるっぽい

というのが判って、これは間違いなく流用してるなと。

というわけで、この座標系のproj.4定義は、おそらく

ZoneA:
+proj=poly +lat_0=40.5 +lon_0=143 +x_0=914398.5307444400 +y_0=1828797.0614888800 +ellps=clrk66 +to_meter=0.9143985307444408 +no_defs
ZoneB:
+proj=poly +lat_0=40.5 +lon_0=135 +x_0=914398.5307444400 +y_0=1828797.0614888800 +ellps=clrk66 +to_meter=0.9143985307444408 +no_defs

だろうというのが判りました*2。

が、これでタイル化してやると、領域はほぼ正しく変換されるのですが、微妙に合わない。
場所によって違いますが、だいたい南東方向に数百mずれる。
これはきっと、米軍はあらたに測量し直したわけじゃなくて日本の旧版地形図の測量成果を流用してるでしょうから、旧日本測地系と世界測地系間の楕円体中心ズレがそのまま入ってきてるんじゃね？と思って、試しに旧日本測地系での楕円体ズレパラメータを入れてみました*3。

ZoneA:
+proj=poly +lat_0=40.5 +lon_0=143 +x_0=914398.5307444400 +y_0=1828797.0614888800 +ellps=clrk66 +to_meter=0.9143985307444408 +towgs84=-146.336,506.832,680.254,0,0,0,0 +no_defs
ZoneB:
+proj=poly +lat_0=40.5 +lon_0=135 +x_0=914398.5307444400 +y_0=1828797.0614888800 +ellps=clrk66 +to_meter=0.9143985307444408 +towgs84=-146.336,506.832,680.254,0,0,0,0 +no_defs

すると、数百mもズレるというような事はなくなり、場所によって異なるものの東西はほぼ一致、南北だけ数十mずれる感じにほぼおさまりました。

が、残りの数十mズレをどう抑え込めばいいのか判りません。
旧日本測地系の成果を流用してるのなら、楕円体をbesselにすればいいんじゃね？*4とか、楕円体ズレパラメータもヤードに直すのか？とかいろいろやってみましたが、いっこうに合う様子がない。
そして最終的に、これはもう楕円体が違う時点で、楕円体中心ズレパラメータもだいたい似た領域にいるのは間違いないけど、微妙に修正しないといけないんじゃね？という結論に達しました。

が、その算出方法がわからない。
正確にはパラメータ振って最小二乗法で検定して…とか判るんだけど、使いこなした経験がありません*5。
なので、もしそれが得意な方おられるならば、どなたかにサクッと出してもらえればなー、と思って記事共有します。

やっていただきたい事は下記の通りです。

評価する値を出力するためのproj.4スクリプトは、

cs2cs +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +to +proj=poly +lat_0=40.5 +lon_0=135 +x_0=914398.5307444400 +y_0=1828797.0614888800 +ellps=clrk66 +to_meter=0.9143985307444408 +towgs84=-146.336+δx,506.832+δy,680.254+δz,0,0,0,0 +no_defs -I -f '%10.5f'

です。
入力が米軍地図座標値（ZoneB）、出力がGoogleメルカトル座標値になります*6。
この中のδx,δy,δzの値をいろいろに振ってもらって、下記の米軍地図座標値を入力し、出力として出てくるメルカトル座標値を下記の比較値と比較し、誤差が最小になるδx,δy,δzの値を見つけてください。
パラメータは記事書き始めた時点で6点しかないですが、もうちょっと集めて追って追加します。

GISエキスパートの方、数学の得意な方の、ご協力をいただけるとありがたいです。

データ：
直江津起源データ

入力米軍地図座標値X	入力米軍地図座標値Y	比較メルカトルX	比較メルカトルY
1316000.25716	1601177.52710	15390049.84597	4463284.59138
1313461.53062	1600634.89841	15387001.91832	4462651.74175
1317747.73014	1603350.73837	15392113.70933	4465785.62674

下松起源データ

684380.42133	1215480.00090	14680025.19063	4028738.91679
682682.42843	1217635.20624	14678142.77804	4031064.79392
680005.65779	1218443.14815	14675188.35875	4031957.93863